An unrecognized force in inertial microfluidics
Siddhansh Agarwal, Fan Kiat Chan, Mattia Gazzola, Sascha Hilgenfeldt
AAn unrecognized force in inertial microfluidics
Siddhansh Agarwal, ∗ Fan Kiat Chan, ∗ Bhargav Rallabandi, Mattia Gazzola,
1, 3, 4 and Sascha Hilgenfeldt † Mechanical Science and Engineering, University of Illinois, Urbana-Champaign, Illinois 61801, USA Mechanical Engineering, University of California, Riverside, USA National Center for Supercomputing Applications,University of Illinois, Urbana-Champaign, Illinois 61801, USA Carl R. Woese Institute for Genomic Biology, University of Illinois, Urbana-Champaign, Illinois 61801, USA (Dated: January 12, 2021)
Describing effects of small but finite inertia onsuspended particles is a fundamental fluid dynam-ical problem that has never been solved in fullgenerality [1–5]. Modern microfluidics has turnedthis academic problem into a practical challengethrough the use of high-frequency ( ω ∼ kHz–MHz) oscillatory flows, perhaps the most effi-cient way to take advantage of inertial effects atlow Reynolds numbers, to precisely manipulateparticles, cells and vesicles without the need forcharges or chemistry [6–8]. The theoretical un-derstanding of flow forces on particles has so farhinged on the pioneering work of Maxey and Ri-ley (MR in the following) [9], almost 40 years ago.We demonstrate here theoretically and computa-tionally that oscillatory flows exert previously un-explained, significant and persistent forces, thatthese emerge from a combination of particle iner-tia and spatial flow variation, and that they canbe quantitatively predicted through a generaliza-tion of MR. Oscillatory microfluidics is usually set up by or pasta localized object (e.g. a microbubble or a no-slip solid[7, 10]), resulting in spatially non-uniform flows charac-terized by strong variations on gradient L Γ and curvature L κ length scales. Such flows exert remarkably consistentand controllable forces on particles, and have been em-ployed with great success for guidance, separation, aggre-gation, and sorting [8, 11–16]. Nonetheless, it is preciselythis use of localized oscillations in modern microfluidicsthat is now pushing the envelope of the MR equation, ex-posing its limits through the observation of unexplained,significant and persistent forces. Here we provide a thor-ough revision of its theoretical foundations, but first, inlight of the importance of this work for applications, westate a major practical outcome: in any oscillatory back-ground flow field ¯ U associated with a localized object,density-matched ( ρ ) spherical particles (radius a p ) expe-rience an attractive force towards the object. The com-ponent of this force along the object-to-particle connector e takes the explicit form F Γ κ = m f (cid:10) a p ∇ ¯ U : ∇∇ ¯ U (cid:11) F ( λ ) · e , (1)where m f is the displaced fluid mass (= 4 πρa p /
3) andthe inner product represents the interaction of flow gra- dients and curvatures. Force (1) is steady, resulting froma time average h·i . The effect of oscillation frequencyis quantified by the universal, analytically derived func-tion F of the Stokes number λ . For harmonic oscillatoryflows, λ ≡ a p ω/ (3 ν ) and to excellent approximation F ( λ )reads F ( λ ) = 13 + 916 r λ , (2)valid over the entire range from the viscous λ (cid:28) λ (cid:29) r p from the localized object, so thatthe steady equation of motion becomes simply dr p dt = F Γ κ πa p νρ , (3)with ν the kinematic viscosity of the fluid. Generally, F Γ κ <
0, since the amplitude of ¯ U decays with distancefrom the oscillating object, indicating attraction. If anadditional steady flow component is present, (3) quanti-fies the deviation between particle and fluid motion.The above equations completely describe the particledynamics and stem from a rigorous, general formalismdeveloped here to respond to discrepancies observed ex-perimentally. As illustrated in Fig. 1ab, when neutrallybuoyant particles of moderate λ approach the surface ofoscillating bubbles (cf. [8, 15, 24, 25]), we find evidenceof significant radial attractive forces, even at a consider-able distance from the bubble. This observation is notexplained by existing theories [7, 12–14, 17–23, 26, 27]that either predict no attraction at all or a much tooweak effect (see caption of Fig. 1).Our goal here is to develop a unifying theory thatexplains observations, accounts for particle inertia, andseamlessly spans the full viscous-to-inviscid operationalflow spectrum. Accordingly, we revisit MR [9] and sys-tematically account for all leading-order terms in particleReynolds number Re p = a p U ∗ /ν , with U ∗ the velocityscale of the background flow. We then reveal their ef-fect through a specially constructed case: a bubble ofradius a b oscillating in pure volume (breathing) mode,with a spherical, neutrally buoyant particle placed at aninitial center-to-center distance r p (0). This scenario in-duces no rectified (streaming) flow in the absence of the a r X i v : . [ phy s i c s . f l u - dyn ] J a n FIG. 1. Particle attraction to oscillating bubbles. (a) A neutrally buoyant particle ( a p = 10 µm , λ ≈
4) is transported past anoscillating microbubble ( a b = 40 µ m, ω/ (2 π ) = 20kHz). (b) Close-up shows the particle trajectory (red) intersecting streamlines(blue), indicating a net attraction towards the bubble over fast time scales of a few ms, unexplained by existing theories: Inertialparticle migration due to shear gradients [4, 17, 18] is far slower; the secondary radiation force of acoustofluidics [12, 19–22] isproportional to the particle-fluid density contrast and thus vanishes here; an ad hoc theory for nearly inviscid flows ( λ (cid:29) (cid:15) . Top figure: instantaneous streamlines (colorbar is flow speed in units of U ∗ ); bottom figure: time-averaged streamlines (color bar is steady flow speed in units of (cid:15)U ∗ ). particle [28], and therefore allows for the precise eval-uation of the newly considered disturbance flow effectsintroduced by the particle itself. The analysis is com-plemented by direct numerical simulations (DNS) thatprovide first-principle solutions of flow field and parti-cle displacement. Figure 1c (upper half) shows that thecomputed oscillatory flow component closely resemblesthe background flow even in the presence of the particle,while time-averaging over an oscillation cycle (bottomhalf) reveals the much richer secondary steady distur-bance flow induced by the particle.Like MR, we wish to describe the hydrodynamic forceson a particle centered at r p using only information fromthe given undisturbed background flow ¯U . We fix a(moving) coordinate system at r p and non-dimensionalizelengths by a p , times by ω − , and velocities by U ∗ (us-ing lowercase letters for non-dimensional velocities). Aspherical particle exposed to a known (lab-frame) back-ground flow ¯ u and moving with velocity u p (neglectingeffects of rotation) then experiences the effects of theundisturbed flow w (0) = ¯u − u p and a disturbance flow w (1) . Following [9], the latter obeys ∇ w (1) − ∇ p (1) = 3 λ ∂ w (1) ∂t + Re p f , where (4) f = w (0) · ∇ w (1) + w (1) · ∇ w (0) + w (1) · ∇ w (1) with boundary conditions w (1) = u p − ¯u on r = 1, and w (1) = 0 as r → ∞ . This equation is exact and doesnot rely on small Re p . To obtain explicit results, we usetwo expansions: one, like MR, expands the backgroundflow around the particle position into spatial moments ofalternating symmetry: ¯u = ¯u | r p + r · E + rr : G + . . . , (5) where E = ( a p /L Γ ) ∇ ¯u | r p and G = ( a p /L κ ) ∇∇ ¯u | r p capture the background flow shear gradients and curva-tures, whose scales are, in practice, much larger than a p ,justifying (5).The other cornerstone of our theory is a regular pertur-bation expansion of all variables in (4), using subscriptsfor orders of Re p , e.g., w (1) = w (1)0 + Re p w (1)1 + O (Re p ).In contrast to MR, this retains a term Re p f in (4), where f = w (0) · ∇ w (1)0 + w (1)0 · ∇ w (0) + w (1)0 · ∇ w (1)0 is theleading-order nonlinear forcing of the disturbance flow.Note also that w (1)0 is purely oscillatory, while w (1)1 has anon-zero time-average, exemplified by the flow in Fig. 1c(bottom).Forces on the particle, as integrals of the fluid stresstensor over the particle surface S p , are also expandedin this fashion. Application of a reciprocal theorem [2]formally yields the inertial force components as volumeintegrals over the entire fluid domain without the needto explicitly compute the flow field at that order. Thereciprocal theorem employs a known test flow u = u ( t ) e in a chosen direction e . The component of the equationof particle motion in that direction, to O (Re p ), is then m p dU p dt = F (0)0 + F (1)0 +Re p ( F (0)1 + F (1)1 ) + O (Re p ) , (6a) F (0)0 = F S π Z V (3 λ∂ t ¯ u ) · e dV, (6b) F (1)0 = F S π L − (Z S p (cid:0) ˆ u p − ˆ¯ u (cid:1) ˆ u · ( ˆ σ · n ) dS ) , (6c) F (0)1 = F S π Z V (¯ u · ∇ ¯ u ) · e dV, (6d) F (1)1 = − F S π L − (cid:26) u Z V ˆ u · ˆ f dV (cid:27) , (6e) FIG. 2. Flow field simulation results. (a-e) Streamlines of the steady flow h w i = h w (1) i (Stokes streamfunction isolines) fordifferent λ ; color bar is velocity magnitude in units of (cid:15)U ∗ ; (f,h) The magnitude of Fourier-transformed quantities (indicated bytildes) evaluated at the driving frequency ω demonstrates that the flow field has no outer, inertia-dominated region. The ratiobetween oscillatory disturbance flow advective force ˜f ( ω ) and the Fourier component of the unsteady inertia ∂ w (1) /∂t remainssmall away from the bubble. (g,i) The Fourier component of vorticity at ω is confined to the oscillatory Stokes layer thickness δ S = p ν/ω (orange-dashed circle) around the particle. where σ is the stress tensor of the test flow, hats de-note Laplace transforms, and L − their inverse. Alldimensional forces have the common Stokes drag scale F S / π = νρa p U ∗ . Equations (6b) and (6d) are forces ex-erted by the background flow, while (6c) and (6e) stemfrom the disturbance flow. The original MR equationcontains F (0)0 and F (1)0 , but only part of F (0)1 , while F (1)1 is an entirely new term due to particle inertia. We shallshow that these unrecognized contributions are not smallcorrections, but are dominant in relevant applications,particularly the inertial disturbance force F (1)1 .This formalism is entirely general for arbitrary givenbackground flows and provides (see Methods) analyticalexpressions for the new forces F (0)1 and F (1)1 . The formerreads F (0)1 F S = 49 ( E : G ) · e F (0)1 , (7)where F (0)1 = 1 / F (1)1 simplifies con-siderably if the background flow around the particle ispotential (this is fulfilled in almost all cases, requiringonly that the distance between the particle and objectsurfaces is greater than the Stokes boundary layer thick-ness). If furthermore the particle is neutrally buoyant, we additionally obtain F (1)1 F S = 49 ( E : G ) · e F (1)1 ( λ ) , (8)where the function F (1)1 ( λ ) is determined analytically(see SI for details) and is universal, i.e., valid for arbitraryflow fields. While both (7) and (8) need non-zero gradi-ent and curvature terms of the background flow, F (1)1 ( λ )captures the nonlinear effect of inertia of the leading or-der unsteady disturbance flow w (1)0 on the particle. Formicron-size particles where λ ∼ F (1)1 is considerablylarger than F (0)1 , so that (8) is the dominant effect inpractical microfluidic applications. Adding both contri-butions (7) and (8), the resulting dimensional force is (1)before time-averaging.We now turn to the prototypical oscillatory flow exam-ple of Fig. 1c. This flow field’s unique scale is the bubbleradius ( L Γ = L κ = a b ). With an oscillation amplitudeof (cid:15)a b ( (cid:15) (cid:28) U ∗ = (cid:15)a b ω , and we anticipate that the relevant rectified(time-averaged) force will be proportional to (cid:15) (cf. [23]).It is advantageous to change the length scale to a b here,introducing α ≡ a p /a b , and to change the coordinate ori-gin to the bubble center, so that the background flow has FIG. 3. Comparison of theoretical (red) and simulated (blue) particle dynamics (radial displacements). (a) Full unsteadydynamics (solid lines) from DNS and theory Eq. (13) and time-averaged dynamics (dashed lines; theory uses Eq. (9) with (2)).The classical MR equation solutions (green) fail to even qualitatively capture the particle attraction to the bubble. (b-e) Steadydynamics from the uniformly valid asymptotic theory agrees with DNS for the entire range of λ values. Dashed lines show theinviscid-limit theory, demonstrating significant quantitative discrepancies even for the largest λ . only one component ¯ u = sin t/r in the direction e = e r .The oscillatory forces and the particle motion now followexplicitly (see Methods).Our ultimate goal is to predict the rectified trajectoryof the particle after time-averaging over the fast oscil-latory time scale, to provide practically useful guidancefor precision applications. Time scale separation usingthe slow time T = (cid:15) t analogous to [23] (see Methods)obtains the leading order equation for rectified particlemotion r p ( T ) dr p dT = − r p α λ F ( λ ) , (9)where F ( λ ) = F (1)1 ( λ )+ F (0)1 ; (9) is readily solved analyt-ically and is analogous to the result (3). Indeed, while theanalytical form of the universal function F (1)1 is compli-cated (see SI), one can Taylor expand in both the viscouslimit ( λ →
0) and the inviscid limit ( λ → ∞ ) to obtain F v = 916 r λ + O (1) , F i = 13 + O (1 / √ λ ) . (10)Simply adding leading terms yields the uniformly validexpression (2) for the total dimensionless force F ( λ ) onthe particle. Note that our derivation is based fundamen-tally on the presence of both viscous and inertial effects,so that even F v is a finite-inertia force. Its λ − / scalingfor small λ is analogous to Saffman’s lift force [30], but isobtained without decomposing the domain into viscousand inertial regions (see Methods). Remarkably, the op-posite limit F i exactly asymptotes to the result obtainedfrom the purely inviscid formalism of [23] as λ → ∞ . We now demonstrate that (2) is accurate over the en-tire range of Stokes numbers by comparing our theorywith independent, large-scale, 3D numerical simulations,previously validated in a range of streaming scenarios[31, 32]. Fig. 2a-e illustrate the rich time-averaged flow h w i at different λ , while Fig. 2g,i exemplify the expectedconfinement of vorticity around the particle. The sim-ulations also serve to justify our omission of an inertia-dominated outer region (Fig. 2f,h). In Fig. 3a, we com-pare analytical and simulated particle trajectories onboth the oscillatory and slow time scales. The classi-cal MR equation fails to capture any of the attractionobserved in DNS, while the present theory is in excellentagreement both for the instantaneous motion and the rec-tified drift of the particle. Moreover, it succeeds over theentire range of λ values, cf. Fig. 3b-e. We see here thatthe inviscid formalism of Ref. [23] (dashed lines) givesa much too weak attraction, particularly for practicallyrelevant λ ∼
1. This is an intuitive outcome of takingviscosity into account, as the Stokes boundary layer (cf.Fig. 2g,i) effectively increases particle size, so that forcesscaling with particle size (cf. (1)) become larger. Fig-ure 3 also illustrates the great benefit of the analyticaltheory (9), as individual DNS incur massive computa-tional cost up to ∼ ,
000 core-hours on the Stampede2supercomputer (see SI).Figure 4 summarizes the comparison between theoryand simulations: Time-averaged DNS trajectories (be-yond an initial transient – see SI for details) for differentvalues of λ were fitted to (9) to determine the dimension-less force F . Our analytical predictions are in quantita-tive agreement with DNS across the range of λ , exhibitingan average error of ≈ FIG. 4. Comparison of the overall inertial force magnitude F in theory (lines) and simulation (symbols), for various λ andinitial particle positions r p (0). The uniformly valid expres-sion (red) is extremely close to the full solution (orange) andin excellent agreement with all DNS data, while the inviscidtheory (black dashed) severely underestimates the forces. Finally, we emphasize that the novel forcing terms (7)and (8), investigated here in isolation, are not small cor-rections relative to the original MR terms in many com-monplace applications. Indeed, for customary microflu-idic settings, they are an order of magnitude strongerthan density contrast induced forces, acoustofluidic radi-ation forces, or Fax´en forces – see the Methods sectionfor a quantitative analysis.In summary, motivated by advancements in microflu-idics, the present work reveals previously unrecognizedforces acting on particles in viscous flows with finite in-ertial effects from oscillatory driving. These forces stemfrom flow gradients and curvatures, are attractive to-wards the oscillating object under mild assumptions, aremuch stronger than inviscid forces, and can be dominantover classically understood MR terms. They lead to sig-nificant displacements of cell-sized particles (1 − µm )over ms time scales, making them a promising tool forprecision manipulation strategies. Our analysis showsthat a surprisingly simple expression accurately predictsparticle motion, as quantitatively confirmed against first-principle, large-scale direct numerical simulations. Thetheory highlights the immense reduction in computa-tional effort between DNS and an explicit analytical the-ory, and as a generalization of the Maxey-Riley formalismis applicable to a wide variety of flow situations. Acknowledgments:
The authors thank KaitlynHood, Gabriel Juarez, and Howard A. Stone for fruit-ful discussions. The authors also acknowledge supportby the National Science Foundation under NSF CA-REER Grant No. CBET-1846752 (MG) and by the BlueWaters project (OCI-0725070, ACI-1238993), a joint ef-fort of the University of Illinois at Urbana-Champaignand its National Center for Supercomputing Applica-tions. This work also used the Extreme Science and En-gineering Discovery Environment (XSEDE) [33] Stam-pede2, supported by National Science Foundation grantno. ACI-1548562, at the Texas Advanced ComputingCenter (TACC) through allocation TG-MCB190004.
Methods
General solutions and the Reciprocal Theorem . Theleading-order oscillatory disturbance flow field w (1)0 is ob-tained by inserting (5) into the leading order of (4) and canbe formally expressed as a series solution [34, 35] w (1)0 = M D · u s − M Q · ( r · E ) − M O · ( rr : G ) + . . . , (11)where u s = u p − ¯u | r p is the slip velocity and M D,Q,O ( r , λ )are spatially dependent mobility tensors independent of theparticular background flow – see SI for explicit expressionsin the case of harmonic oscillatory flows. All informationabout the specific background flow is contained in the con-stant quantities u s , E , and G . The O (Re p ) flow field w (1)1 does not need to be computed explicitly; instead, we use areciprocal theorem. Denoting Laplace-transformed quantitiesby hats, application of the divergence theorem results in thefollowing symmetry relation: I S ( ˆ w (1)1 · ˆ σ − ˆ u · ˆ σ (1)1 ) · m dS = Z V h ∇ · ( ˆ w (1)1 · ˆ σ ) − ∇ · (ˆ u · ˆ σ (1)1 ) i dV. (12)As shown in the SI, the above expression yields the O (Re p )force on the particle captured by (6e). We note that thecomputation of the volume integral simplifies considerably:the integrand is proportional to f , in which only certainproducts are non-vanishing when the angular integrationaround the particle is performed. For instance, the firstterm in f is ( ¯ u − u p ) · ∇ w (1)0 = ( − u s + r · E + rr : G ) ·∇ ( M D · u s − M Q · ( r · E ) − M O · ( rr : G )). Due to the al-ternating symmetry of terms in the background flow and con-sequently w (1)0 , only products of adjacent terms survive inte-gration, while e.g. a term involving u s · ∇ ( M D · u s ) vanishesafter volume integration. Simplification for neutrally buoyant particles . Inthis work, we restrict our analysis to the case of neutrallybuoyant particles. Consequently, the slip velocity u s vanishes,so that only tensor products involving E and G contribute tothe volume integral (given in SI). Furthermore, in typical mi-crofluidic applications, the particle experiences an oscillatorybackground flow that is potential: to excellent approximation, this holds true if the particle is outside the Stokes boundarylayer of the oscillating object, that is if h p = r p − a p − a b & δ S ,which simplifies to the easily satisfied condition λ & ( a p /h p ) .Then, the only way to construct a force vector from a con-traction of the higher rank tensors E and G in a potentialflow is E : G (cf. [36, 37]), which can be pulled out of the in-tegral as a common factor. Thus, the volume integral can beevaluated generally for any background flow field and resultsin the single term (8) with the universal function F (1)1 ( λ ). Inner-Outer (wake) formalism . Often, the evaluationof forces on particles in a flow is complicated by the tran-sition between a viscous-dominated inner flow volume (nearthe particle) and an inertia-dominated outer volume, neces-sitating an asymptotic matching of the two limits (such asfor the Oseen [2] and Saffman [30] problems). The presentformalism, however, only employs an inner-solution expan-sion and still obtains highly accurate predictions (see also[5]). This behavior can be rationalized by invoking the anal-ysis of Lovalenti and Brady [2] who showed that an outerregion does not occur when the characteristic unsteady timescale ω − is shorter than the convective inertial time scale ν/ ( U ∗ w (0) ) , where w (0) is the dimensionless velocity scaleof the fluid as measured in the particle reference frame. Fordensity matched particles w (0) = O ( α ), so that this criterionreduces to (cid:15) λ (cid:28)
1, requiring the oscillation amplitude of theflow to be smaller than δ S α − , which is easily satisfied in mostexperimental situations. More directly, the Lovalenti-Bradycriterion relies on the magnitude of oscillatory inertia in thedisturbance flow ∂ w (1) /∂t being much larger than that of theadvective term f . DNS verifies that this relation holds for theentire range of λ treated here (see Fig. 2f,h). As a separate ef-fect, outer flow inertia due to the slow (steady) motion of theparticle will be present, but only results in O ( (cid:15) ) correctionsto the Stokes drag. Oscillatory equation of motion in radial flow . Forthe special case of the bubble executing pure breathing oscil-lations with the radial flow field ¯ u = sin t/r , it is straight-forward to compute E : G · e r = −
18 sin t/r p , where r p ( t ) isthe instantaneous particle position. Using (6), (7), (8), andnoting α Re p = 3 (cid:15)λ , the non-dimensional equation of motionfor r p ( t ) of a neutrally buoyant particle explicitly reads: λ d r p dt = (cid:15)λ (cid:18) cos tr p − (cid:15) sin tr p (cid:19) − λ (cid:15) α
18 sin tr p F (0) + (cid:20) sin tr p − dr p dt (cid:21) − (cid:20) λ (cid:15) α (18 sin t ) r p F (1)1 λ ) (cid:21) , (13)where the first line on the RHS represents contributions from F (0)0 and F (0)1 , while the first and second terms in squareparentheses represent F (1)0 and F (1)1 , respectively. Note that,for neutrally buoyant particles, the time-periodic character ofthe flow precludes memory terms that would otherwise emergefrom the inverse Laplace transforms [2, 9, 38]. Time scale separation and time averaging . Assuming (cid:15) (cid:28)
1, we introduce the slow time T = (cid:15) t , in addition to thefast time t . Using the following transformations r p ( t ) r p ( t, T ) , (14a) ddt ∂∂t + (cid:15) ∂∂T , (14b) d dt ∂ ∂t + 2 (cid:15) ∂ ∂t∂T + (cid:15) ∂ ∂T , (14c)we seek a perturbation solution in (cid:15) of the general form r p ( t, T ) = r p ( T ) + (cid:15) ˇ r p ( t, T ) + (cid:15) ˇˇ r p ( t, T ) + . . . , and separateorders in (13). The procedure is outlined in [23] and resultsin a leading-order equation for r p ( T ) given by (9), dependenton the slow time scale only (the scale t being averaged out). Simulation method and numerical implementation .Here we briefly describe the governing equations and numeri-cal technique used in our simulations. We consider two spheri-cal bodies (an oscillating microbubble and a neutrally buoyantparticle) immersed in an unbounded domain of incompress-ible viscous fluid. We denote the computational domain asΩ = Ω f ∪ Ω B , where Ω f is the fluid domain and Ω B = Ω b ∪ Ω p is the domain in which the bubble (Ω b ) and particle (Ω p ) re-side, and denote the interface between the fluid and the bod-ies as ∂ Ω B . The flow is then described by the incompressibleNavier–Stokes equation ∇· u = 0 , ∂ u ∂t + ( u ·∇ ) u = − ρ ∇ p + ν ∇ u x ∈ Ω \ Ω B (15)where ρ , p , u and ν are the fluid density, pressure, velocityand kinematic viscosity, respectively. We impose the no-slipboundary condition u = u B at ∂ Ω B , where u B is the bodyvelocity, and feedback from the fluid to the body is describedby Newton’s equation of motion. The system of equations issolved in velocity–vorticity form using the remeshed vortexmethod combined with Brinkmann penalization and a pro-jection approach [39]. This method has been extensively val-idated across a range of fluid–structure interaction problems,from flow past bluff bodies to biological swimming [39–43].Recently, the accuracy of this method has been demonstratedin rectified flow contexts as well, capturing steady streamingresponses from arbitrary shapes in 2D and 3D [31, 32]. Moredetails on method implementation and simulation techniquescan be found in the SI. Comparison with other hydrodynamic forces . Wehave discussed a special case of radial symmetry quantita-tively because it isolates the novel inertial forces reportedhere as the only effect, thus allowing us to assess the accu-racy of the theory. In more general flow situations, otherforces will compete with F Γ κ , and we assess their rela-tive magnitude here. If the particle density ρ p does notmatch ρ , a density contrast force [23] is induced, general-izing acoustofluidic secondary radiation forces. In order forthis force to exceed F Γ κ , the density contrast needs to fulfill ρ p /ρ − & a p /r p ) (1 + 2 / √ λ ). Appreciable forces only actwhen r p & a b and if λ is not very small; thus, ρ p /ρ − & . α . .
2. In most mi-crofluidic, and certainly in biomedical applications, the den-sity contrast is far less: even at 5% density difference (e.g.for polystyrene particles), F Γ κ is 5-30 times stronger than thedensity contrast force for 0 . < λ <
5. Other forces resultfrom steady flows: oscillation of an a b -sized object will gener-ically induce steady streaming flow at speed ∼ (cid:15) a b U ∗ , and itmay have transverse gradients of scale a b (in addition to ra-dial gradients). This situation induces a Saffman lift force L S [30] for particles with finite slip velocity V s (again because ofdensity mismatch) [2, 23]. L S and F Γ κ are of equal magnitudeif V s ∼ α (4 . √ λ ) U ∗ . In realistic settings, V s would need to exceed U ∗ , implying that the steady flow would overwhelmthe oscillatory motion, defeating the purpose of oscillatory-flow microfluidics. Lastly, flows with finite ∇ ¯U give rise toFax´en terms in added mass and drag. However, the oscilla-tory flows discussed here are (almost) potential flows as shownabove, so that the leading order effect of Fax´en terms comesfrom steady flow curvature and provides only an O ( α ) cor-rection to the steady-flow Stokes drag.We conclude that theinertial force described here is the dominant effect in manyrealistic oscillating microfluidics applications. ∗ S.A. and F.K.C. contributed equally to this work. † [email protected][1] B. Ho and L. Leal, Journal of fluid mechanics , 365(1974).[2] P. M. Lovalenti and J. F. Brady, Journal of Fluid Me-chanics , 561 (1993).[3] J. A. Schonberg and E. Hinch, Journal of Fluid Mechan-ics , 517 (1989).[4] D. Di Carlo, Lab on a Chip , 3038 (2009).[5] K. Hood, S. Lee, and M. Roper, Journal of Fluid Me-chanics , 452 (2015).[6] B. R. Lutz, J. Chen, and D. T. Schwartz, Analyticalchemistry , 5429 (2006).[7] P. Rogers and A. Neild, Lab on a Chip , 3710 (2011).[8] R. Thameem, B. Rallabandi, and S. Hilgenfeldt, PhysicalReview Fluids , 052001 (2017).[9] M. R. Maxey and J. J. Riley, The Physics of Fluids ,883 (1983).[10] B. R. Lutz, J. Chen, and D. T. Schwartz, Physics of Flu-ids , 023601 (2005).[11] C. Wang, S. V. Jalikop, and S. Hilgenfeldt, AppliedPhysics Letters , 034101 (2011).[12] L. Schmid, D. A. Weitz, and T. Franke, Lab on a Chip , 3710 (2014).[13] Y. Chen and S. Lee, Integrative and comparative biology , 959 (2014).[14] I. S. Park, J. H. Shin, Y. R. Lee, and S. K. Chung, Sensorsand Actuators A: Physical , 214 (2016).[15] R. Thameem, B. Rallabandi, and S. Hilgenfeldt, Biomi-crofluidics , 014124 (2016).[16] A. Volk, M. Rossi, B. Rallabandi, C. J. K¨ahler, S. Hilgen-feldt, and A. Marin, Physical Review Fluids , 114201(2020).[17] D. Di Carlo, D. Irimia, R. G. Tompkins, and M. Toner,PNAS , 18892 (2007).[18] M. E. Warkiani, A. K. P. Tay, B. L. Khoo, X. Xiaofeng,J. Han, and C. T. Lim, Lab on a Chip , 1101 (2015).[19] V. Bjerknes, Fields of force (General Books, 1906).[20] T. A. Hay, M. F. Hamilton, Y. A. Ilinskii, and E. A. Zabolotskaya, The Journal of the Acoustical Society ofAmerica , 1331 (2009).[21] W. T. Coakley and W. Nyborg, Ultrasound: Its applica-tions in medicine and biology , 77 (1978).[22] H. Bruus, Lab on a Chip , 1014 (2012).[23] S. Agarwal, B. Rallabandi, and S. Hilgenfeldt, PhysicalReview Fluids , 104201 (2018).[24] C. Wang, S. V. Jalikop, and S. Hilgenfeldt, Biomicroflu-idics , 012801 (2012).[25] C. Wang, B. Rallabandi, and S. Hilgenfeldt, Physics ofFluids , 022002 (2013).[26] A. Hashmi, G. Yu, M. Reilly-Collette, G. Heiman, andJ. Xu, Lab on a Chip , 4216 (2012).[27] Y. Chen, Z. Fang, B. Merritt, D. Strack, J. Xu, andS. Lee, Lab on a Chip , 3024 (2016).[28] M. S. Longuet-Higgins, Proceedings of the Royal Societyof London. Series A: Mathematical, Physical and Engi-neering Sciences , 725 (1998).[29] B. Rallabandi, (2020), preprint.[30] P. Saffman, Journal of fluid mechanics , 385 (1965).[31] T. Parthasarathy, F. K. Chan, and M. Gazzola, Journalof Fluid Mechanics , 647 (2019).[32] Y. Bhosale, T. Parthasarathy, and M. Gazzola, Journalof Fluid Mechanics , A13 (2020).[33] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither,A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D.Peterson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr,Computing in Science & Engineering , 62 (2014).[34] L. D. Landau and E. Lifshitz, Course of TheoreticalPhysics Vol. 6 Fluid Mechanies (Pergamon Press, 1959).[35] C. Pozrikidis et al. , Boundary integral and singularitymethods for linearized viscous flow (Cambridge Univer-sity Press, 1992).[36] S. Danilov and M. Mironov, The Journal of the Acousti-cal Society of America , 143 (2000).[37] A. Nadim and H. A. Stone, Studies in Applied Mathe-matics , 53 (1991).[38] A. B. Basset, A treatise on hydrodynamics: with nu-merous examples , Vol. 2 (Deighton, Bell and Company,1888).[39] M. Gazzola, P. Chatelain, W. M. Van Rees, andP. Koumoutsakos, Journal of Computational Physics , 7093 (2011).[40] M. Gazzola, C. Mimeau, A. A. Tchieu, and P. Koumout-sakos, Physics of fluids , 043103 (2012).[41] M. Gazzola, W. M. Van Rees, and P. Koumoutsakos,Journal of Fluid Mechanics , 5 (2012).[42] M. Gazzola, B. Hejazialhosseini, and P. Koumoutsakos,SIAM Journal on Scientific Computing , B622 (2014).[43] M. Gazzola, A. A. Tchieu, D. Alexeev, A. de Brauer, andP. Koumoutsakos, Journal of Fluid Mechanics , 726(2016). upplementary Information: An unrecognized force in inertial microfluidics Siddhansh Agarwal, Fan Kiat Chan, Bhargav Rallabandi, Mattia Gazzola,
1, 3, 4 and Sascha Hilgenfeldt ∗ Mechanical Science and Engineering, University of Illinois, Urbana-Champaign, Illinois 61801, USA Mechanical Engineering, University of California, Riverside, USA National Center for Supercomputing Applications,University of Illinois, Urbana-Champaign, Illinois 61801, USA Carl R. Woese Institute for Genomic Biology, University of Illinois, Urbana-Champaign, Illinois 61801, USA
I. THEORETICAL FORMALISM
In order to systematically account for the inertial forces on a sphere of radius a p centered at r p moving with withvelocity u p (neglecting effects of rotation) and exposed to a known (lab-frame) background undisturbed flow ¯ u , wesplit the Navier–Stokes equations that govern the flow field into an undisturbed flow w (0) = ¯ u − u p and a disturbanceflow w (1) (we adopt the same notation as [1]). Then, in a particle-centered (moving) coordinate system, we have ∇ w (0) − ∇ p (0) =3 λ ∂ w (0) ∂t + Re p (cid:16) w (0) · ∇ w (0) (cid:17) , (1a) ∇ w (1) − ∇ p (1) =3 λ ∂ w (1) ∂t + Re p (cid:20) ( ¯u − u p ) · ∇ w (1) + w (1) · ∇ ¯u + w (1) · ∇ w (1) (cid:21) , (1b) ∇ · w (0) =0 , ∇ · w (1) = 0 , (1c) w (1) = u p − ¯u on r = 1 and w (1) = 0 as r → ∞ , (1d)where Re p = U ∗ a p /ν is the particle Reynolds number. Quantities in these equations are non-dimensionalized byscaling velocities with U ∗ , lengths with a p , pressure with µU ∗ /a p , and time by ω − .The force contribution from the undisturbed flow is F (0) = ( F S / π ) H S n · σ (0) dS , like in the original Maxey–Riley(MR) formalism [1], where σ (0) = − p (0) I + ∇ w (0) + (cid:0) ∇ w (0) (cid:1) T is the stress tensor associated with the undisturbedflow field w (0) , and F S / π = νρa p U ∗ is the Stokes drag scale. The force contribution at the disturbance flow order isgiven by F (1) = ( F S / π ) H S n · σ (1) dS , where σ (1) = − p (1) I + ∇ w (1) + (cid:0) ∇ w (1) (cid:1) T is the stress tensor associated withthe disturbance flow field w (1) . The corresponding (dimensional) equation of motion for the particle then reads m p d U p dt = F (0) + F (1) . (2)Note that everything up to this point is exact and no assumptions have been made. MR [1] make the unsteady Stokesflow approximation in (1b) by setting Re p = 0, and compute F (1) without explicitly evaluating the disturbance flow,using a symmetry relation. While this assumption is plausible in many traditional microfluidic flow situations, fastoscillatory particle motion can give rise to large disturbance flow gradients so that the inertial terms on the RHS of(1b) are not necessarily negligible compared to the viscous diffusion term (typically Re p ∼ O (1), as in the experimentdescribed in Fig. 1 of the main text). A. Small Re p expansion In order to make analytical progress, following [2–4], we expand w (1) , p (1) , r p , u p and σ (1) (and consequently F (1) )in a regular asymptotic expansion for small Re p , w (1) = w (1)0 + Re p w (1)1 + . . . , (3a) p (1) = p (1)0 + Re p p (1)1 + . . . , (3b) ∗ [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] J a n r p = r p + Re p r p + . . . , (3c) u p = u p + Re p u p + . . . , (3d) σ (1) = σ (1)0 + Re p σ (1)1 + . . . , (3e) F (1) = F (1)0 + Re p F (1)1 + . . . . (3f)The leading-order equations for ( w (1)0 , p (1)0 ) are unsteady Stokes, ∇ w (1)0 − ∇ p (1)0 = 3 λ ∂ w (1)0 ∂t , (4a) ∇ · w (1)0 = 0 , (4b) w (1)0 = u p − ¯u on r = 1 , (4c) w (1)0 = 0 as r → ∞ . (4d)We note that in the original derivation of MR [1], a symmetry relation was used at this order to compute F (1)0 withoutexplicitly solving for w (1)0 . However, since we are interested in computing the force contribution at O (Re p ), we needan explicit solution for the leading-order disturbance flow w (1)0 . To obtain explicit results, as stated in the maintext, we expand the background flow field ¯u around the leading-order particle position r p into spatial moments ofalternating symmetry, ¯u = ¯u | r p + r · E + rr : G + . . . , (5)where E = ( a p /L Γ ) ∇ ¯u | r p and G = ( a p /L κ ) ∇∇ ¯u | r p with gradient L Γ and curvature L κ length scales. As aconsequence of (5), the boundary condition (4c) is also expanded around r p , so that in the particle fixed coordinatesystem w (1)0 = u p − ¯u = u p − ¯u | r p − r · E − rr : G + . . . on r = 1 , (6)where we have retained the first three terms in the background flow velocity expansion. Owing to the linearity of theleading order unsteady Stokes equation, the solution can generally be expressed as [5, 6] w (1)0 = M D · u s − M Q · ( r · E ) − M O · ( rr : G ) + . . . , (7)where M D,Q,O ( r, λ ) are spatially dependent mobility tensors. For oscillatory flows, they depend on the Stokes number λ . More explicit forms of these tensors will be given below.With the leading-order disturbance flow field known, the equations at O (Re p ) are as follows, ∇ w (1)1 − ∇ p (1)1 = ∇ · σ (1)1 = 3 λ ∂ w (1)1 ∂t + f , (8a) ∇ · w (1)1 = 0 , (8b) w (1)1 = u p on r = 1 , (8c) w (1)1 = 0 as r → ∞ , (8d)where f = w (0) · ∇ w (1)0 + w (1)0 · ∇ w (0) + w (1)0 · ∇ w (1)0 is the (explicitly known) leading-order nonlinear forcing of thedisturbance flow. In order to compute the force at this order, we employ a reciprocal relation in the Laplace domainsince the problem is time-dependent and, for oscillatory flows, the Laplace transform is explicitly obtained. B. Reciprocal theorem and test flow
A known test flow (denoted by primed quantities such as u ) is chosen around an oscillating sphere such that itsatisfies the following unsteady Stokes equation: ∇ u − ∇ p = ∇ · σ = 3 λ ∂ u ∂t , (9a) ∇ · u = 0 , (9b) u = u ( t ) e on r = 1 , (9c) u = 0 as r → ∞ , (9d)where the unit vector e is chosen to coincide with the direction in which the force on the particle is desired. Thesolution to this problem is of the same form as (7), but with only the first term, i.e., u = u ( t ) M D · e . (10)Denoting Laplace transformed quantities by hats (e.g., ˆ u ), one can write down the following symmetry relation usingthe divergence theorem (cf. [1, 4, 7]): I S ( ˆ w (1)1 · ˆ σ − ˆ u · ˆ σ (1)1 ) · m dS = Z V h ∇ · ( ˆ w (1)1 · ˆ σ ) − ∇ · (ˆ u · ˆ σ (1)1 ) i dV, (11)where m is the outward unit normal vector to the surface (pointing inward over the sphere surface), and ˆ σ = ∇ ˆ u + ( ∇ ˆ u ) T − ˆ p I . Substituting boundary conditions from (8) and (9), and setting the volume equal to the fluid-filleddomain, we obtainˆ u (1) p · Z S p ( ˆ σ · m ) dS − ˆ u e · Z S p ( ˆ σ (1)1 · m ) dS + Z S ∞ ( ˆ w (1)1 · ˆ σ ) · m dS − Z S ∞ (ˆ u · ˆ σ (1)1 ) · m dS = Z V h ˆ w (1)1 · ( ∇ · ˆ σ ) − ˆ u · ( ∇ · ˆ σ (1)1 ) + ∇ ˆ w (1)1 : ˆ σ − ∇ ˆ u : ˆ σ (1)1 i dV . (12)The third term on the LHS is 0 since the viscous test flow stress tensor decays to zero at infinity. Similarly, the integralin the fourth term vanishes in the far field if viscous stresses dominate inertial terms, and also in the case of inviscidirrotational flows (see [7, 8]). The third and fourth terms on the RHS also go to zero, owing to incompressibilty andsymmetry of the stress tensor: ∇ ˆ w (1)1 : ˆ σ − ∇ ˆ u : ˆ σ (1)1 = ∇ ˆ w (1)1 : ( ∇ ˆ u + ( ∇ ˆ u ) T ) − ˆ p ∇ · ˆ w (1)1 − ∇ ˆ u : ( ∇ ˆ w (1)1 + ( ∇ ˆ w (1)1 ) T ) − ˆ p (1) ∇ · ˆ u = 0 . (13)The divergence of the hatted stress tensors in the remaining two terms of the RHS can be obtained by taking theLaplace transforms of (8) and (9) and using the property [ f ( t ) = s d f ( t ) − f (0), so that ∇ · ˆ σ = ¯ λs ˆ u − u (0) , (14a) ∇ · ˆ σ (1)1 = ¯ λs ˆ w (1)1 − w (1)1 (0) + ˆ f . (14b)Now, the force on the sphere at this order is given by F (1)1 = R S p ( σ (1)1 · n ) dS = − R S p ( σ (1)1 · m ) dS , since m pointsinwards while n points outwards on the surface of the sphere. Assuming both flows start from rest, we have (cf. [7])ˆ u e · ˆ F (1)1 F S / (6 π ) = ˆ u p · Z S p ( ˆ σ · n ) dS − Z V ˆ u · ˆ f dV + O (Re p ) . (15)Adding the force contribution from the previous order, the net force on the particle due to its disturbance flow readsˆ u e · ˆ F (1) F S / (6 π ) = ˆ u e · (cid:16) ˆ F (1)0 + Re p ˆ F (1)1 (cid:17) + O (Re p )= Z S p (cid:0) ˆ u p − ˆ¯ u + Re p ˆ u p (cid:1) · ( ˆ σ · n ) dS − Re p Z V ˆ u · ˆ f dV + O (Re p ) (16a)= ⇒ e · F (1) = F S π L − (Z S p (cid:0) ˆ u p − ˆ¯ u (cid:1) ˆ u · ( ˆ σ · n ) dS − u Re p Z V ˆ u · ˆ f dV ) + O (Re p ) , (16b)where we have used u p = u p + Re p u p + O (Re p ), and L − denotes the inverse Laplace transform. The first term onthe RHS of (16b) is denoted as F (1)0 in the main text (and is the same as that obtained by MR), while the secondterm represents the O (Re p ) inertial force and is denoted as F (1)1 in the main text. II. EVALUATION OF THE O (Re p ) INERTIAL FORCE
In this section, we will explicitly evaluate the volume integral in (16b) representing the O (Re p ) inertial force. Thisrequires obtaining f from the leading-order oscillatory disturbance flow field w (1)0 . A. General solution to equation (4)
We already remarked that, given the background flow field expansion in uniform, linear, and quadratic parts aroundthe particle, w (1)0 is formally obtained as the linear combination (7). For harmonically oscillating, axisymmetricbackground flows (i.e., ¯u ( r ) = { ¯ u r , ¯ u θ , } in the spherical particle coordinate system, with all components ∝ e it ),general explicit expressions can be derived for the mobility tensors M D,Q,O , ensuring no-slip boundary conditionson the sphere order-by-order. A procedure obtaining M D is described in Landau–Lifshitz [5]; the other tensors aredetermined analogously. Using components in spherical coordinates, they read M D = a ( r ) r a ( r ) r
00 0 0 , M Q = b ( r ) r b ( r )3 r
00 0 0 , M O = − c ( r )3 r c ( r )3 r
00 0 0 , (17)where a ( r ) = 12 β r h β − iβ + 3 − e − iβ ( r − (1 + iβr ) i , (18a) b ( r ) = 1 β ( β − i ) r h β ( −
15 + β ( β − i )) + 15 i + 5 e − iβ ( r − ( βr (3 + iβr ) − i ) i , (18b) c ( r ) = − β ( β ( −
45 + β ( β − i )) + 105 i )) + 21 e − iβ ( r − (15 + βr ( − βr (6 + iβr ) + 15 i ))32 β ( − β ( β − i )) r , (18c)and β = q − ia p / ( ν/ω ) = √− iλ is the complex oscillatory boundary layer thickness. We emphasize that theseexpressions are the same for arbitrary axisymmetric oscillatory ¯u . Accordingly, only the expansion coefficients u s , E ,and G contain information about the particular flow.Similarly, the solution to the unsteady test flow is obtained directly as u = M D · cos θ − sin θ e it . (19)It is understood everywhere that physical quantities are obtained by taking real parts of these complex functions. B. Evaluation of F (1)1 In order to compute the volume integral in (16b), we first note that only certain products in f are non-vanishingwhen the angular integration over θ is performed. In particular, due to alternating symmetry of terms in the back-ground flow field expansion (5), and consequently in the leading order disturbance flow (7), only products of adjacentterms survive. This is because, in the Taylor expansion of the background flow field, the first and third terms aresymmetric ( u ( − r ) = u ( r )) while the second one is anti-symmetric ( u ( − r ) = − u ( r )). For example, the first term in f reads w (0) · ∇ w (1)0 = ( − u s + r · E + rr : G ) · ∇ ( M D · u s − M Q · ( r · E ) − M O · ( rr : G )) , (20)and the only terms that survive the angular integration are the symmetric ones (after a contraction with the symmetrictest flow u ), i.e., ( − u s + rr : G ) · ∇ ( − M Q · ( r · E )) + ( r · E ) · ∇ ( M D · u s − M O · ( rr : G )) . (21)Furthermore, in this paper we restrict ourselves to the case of neutrally buoyant particles and consequently the slipvelocity is u s = 0. In summary, only the following terms in f have non-trivial contributions to the volume integral: f = − ( rr : G ) · ∇ ( M Q · ( r · E )) − ( r · E ) · ∇ ( M O · ( rr : G )) − ( M Q · ( r · E )) · ∇ ( rr : G ) − ( M O · ( rr : G )) · ∇ ( r · E )+ ( M Q · ( r · E )) · ∇ ( M O · ( rr : G )) + ( M O · ( rr : G )) · ∇ ( M Q · ( r · E )) . (22)All information about the background flow field is contained in the constant tensors E and G , which are evaluatedat the particle position. If the particle is farther away from the surface of the oscillating object exciting the flowthan the Stokes layer thickness δ S , it is exposed to a pure potential flow; this will be the case in the overwhelmingmajority of realistic scenarios. For potential flows it can be shown that all non-zero terms of (22) are proportionalto E : G . Choosing a test flow in direction e , one obtains a surprisingly compact result for the e -component of the O (Re p ) inertial force: * F (1)1 F S + = − π (cid:28) L − (cid:26) u Z V ˆ u · ˆ f dV (cid:27) (cid:29) = 49 h E : G i · e F (1)1 ( λ ) . (23)We have here applied the required Laplace transforms as well as a time average to extract the steady part of theforce. Performing the volume integral leaves a universal dimensionless function F ( λ ), whose contributions stem from M D,Q,O . Explicitly, this function reads F (1)1 ( λ ) = (cid:20) π (cid:18) − λ / − λ / + 34005¯ λ / + 59790¯ λ / + 3312¯ λ / + 568¯ λ + 14078¯ λ + 97470¯ λ − λ − λ − λ − p ¯ λ − (cid:19) + e (1 − i ) √ ¯ λ π ¯ λ / (cid:18) (cid:18) π (4410 + 2033 i ) − e (1+ i ) √ ¯ λ Ei (cid:16) − p ¯ λ (cid:17) + e (2+2 i ) √ ¯ λ (5600 − i )Ei (cid:16) ( − − i ) p ¯ λ (cid:17) − (2033 + 4410 i ) e i √ ¯ λ Ei (cid:16) ( − − i ) p ¯ λ (cid:17) + e √ ¯ λ (5600 + 12600 i )Ei (cid:16) ( − i ) p ¯ λ (cid:17) − (2033 − i )Ei (cid:16) ( − i ) p ¯ λ (cid:17) + e (2+2 i ) √ ¯ λ (12600 + 5600 i ) π + e i √ ¯ λ (4410 − i ) π + e √ ¯ λ (12600 − i ) π (cid:19) ¯ λ / + 6 (cid:18) π (4195 + 3982 i ) − e (1+ i ) √ ¯ λ Ei (cid:16) − p ¯ λ (cid:17) − (3982 + 4195 i ) e i √ ¯ λ Ei (cid:16) ( − − i ) p ¯ λ (cid:17) − (3982 − i )Ei (cid:16) ( − i ) p ¯ λ (cid:17) + e i √ ¯ λ (4195 − i ) π (cid:19) ¯ λ / + 4 (cid:18) π (241 + 1714 i ) + 720 e (1+ i ) √ ¯ λ Ei (cid:16) − p ¯ λ (cid:17) − (1714 + 241 i ) e i √ ¯ λ Ei (cid:16) ( − − i ) p ¯ λ (cid:17) − (1714 − i )Ei (cid:16) ( − i ) p ¯ λ (cid:17) + e i √ ¯ λ (241 − i ) π (cid:19) ¯ λ / − (120 + 120 i ) (cid:18) π (cid:16) − i + e i √ ¯ λ (cid:17) − ie i √ ¯ λ Ei (cid:16) ( − − i ) p ¯ λ (cid:17) + Ei (cid:16) ( − i ) p ¯ λ (cid:17) (cid:19) ¯ λ / − (4 + 4 i ) (cid:18) π (cid:16) e i √ ¯ λ (248 + 127 i ) + ( − − i ) (cid:17) + e i √ ¯ λ (127 − i )Ei (cid:16) ( − − i ) p ¯ λ (cid:17) + (248 − i )Ei (cid:16) ( − i ) p ¯ λ (cid:17) (cid:19) ¯ λ − (6 + 6 i ) (cid:18) e i √ ¯ λ π (567 + 2134 i ) + e (1+ i ) √ ¯ λ (736 − i )Ei (cid:16) − p ¯ λ (cid:17) + e i √ ¯ λ (2134 − i )Ei (cid:16) ( − − i ) p ¯ λ (cid:17) + (567 − i )Ei (cid:16) ( − i ) p ¯ λ (cid:17) + ( − − i ) π (cid:19) ¯ λ + (cid:18) π (39033 + 25089 i ) − e (1+ i ) √ ¯ λ Ei (cid:16) − p ¯ λ (cid:17) − (25089 + 39033 i ) e i √ ¯ λ Ei (cid:16) ( − − i ) p ¯ λ (cid:17) − (25089 − i )Ei (cid:16) ( − i ) p ¯ λ (cid:17) + e i √ ¯ λ (39033 − i ) π (cid:19) ¯ λ + (315 + 315 i ) (cid:18) e (2+2 i ) √ ¯ λ π (420 + 60 i ) − (96 − i ) e (1+ i ) √ ¯ λ Ei (cid:16) − p ¯ λ (cid:17) + e (2+2 i ) √ ¯ λ (60 − i )Ei (cid:16) ( − − i ) p ¯ λ (cid:17) − (49 + 28 i ) e i √ ¯ λ Ei (cid:16) ( − − i ) p ¯ λ (cid:17) + e √ ¯ λ (420 − i )Ei (cid:16) ( − i ) p ¯ λ (cid:17) + (28 + 49 i )Ei (cid:16) ( − i ) p ¯ λ (cid:17) − (60 + 420 i ) e √ ¯ λ π + (49 − i ) π + e i √ ¯ λ (28 − i ) π (cid:19) ¯ λ + 15120 e √ ¯ λ (cid:18) − ie √ ¯ λ π + 5 ie (1+2 i ) √ ¯ λ π + 5 e (1+2 i ) √ ¯ λ Ei (cid:16) ( − − i ) p ¯ λ (cid:17) + 5 e √ ¯ λ Ei (cid:16) ( − i ) p ¯ λ (cid:17) − e i √ ¯ λ Ei (cid:16) − p ¯ λ (cid:17) (cid:19) + 945 (cid:18) π + 7 e i √ ¯ λ π − ie √ ¯ λ π + 160 ie (2+2 i ) √ ¯ λ π + 160 e (2+2 i ) √ ¯ λ Ei (cid:16) ( − − i ) p ¯ λ (cid:17) + 160 e √ ¯ λ Ei (cid:16) ( − i ) p ¯ λ (cid:17) − e (1+ i ) √ ¯ λ Ei (cid:16) − p ¯ λ (cid:17) − ie i √ ¯ λ Ei (cid:16) ( − − i ) p ¯ λ (cid:17) + 7 i Ei (cid:16) ( − i ) p ¯ λ (cid:17) (cid:19)p ¯ λ (cid:19)(cid:21)(cid:30) (cid:20) p ¯ λ (cid:16) λ / + 32¯ λ / + 8¯ λ + 64¯ λ + 72¯ λ + 36 p ¯ λ + 9 (cid:17) (cid:21) . (24)Here ¯ λ = (3 / λ and Ei is the exponential integral function. We show below that this lengthy expression is approxi-mated to great accuracy by two simple terms.We stress again here that the result is universal for any oscillatory potential flow; for the prototypical case of thevolumetrically oscillating bubble, h E : G i · e r = − /r p , as noted in the Methods section. C. Net inertial force
The time-averaged force contribution from the background flow at O (Re p ) is of the same form as (23), except that F ( λ ) is replaced by the simple constant F (0)1 = [9]. The two contributions F (1)1 and F (0)1 can thus be simply added.Transforming back to dimensional variables, we obtain the net time-averaged force on the particle as F Γ κ = m f a p (cid:10) ∇ ¯ U : ∇∇ ¯ U (cid:11) F ( λ ) , (25)where F = F (1)1 + F (0)1 and m f = 4 πa p /
3, as noted in the main text. This time-averaged inertial force on the particleis derived for a background flow that is symmetric about an axis e passing through the center of the particle.It was remarked above that the simple form of (25) is a consequence of the background flow being potential. Thiscan be backed up by symmetry arguments and dimensional analysis for an arbitrary oscillatory background flow thathas a harmonic scalar potential, ¯ U = ∇ ¯ ϕ with ∇ ¯ ϕ = 0. Such a flow is in fact generic since the background flowvorticity decays exponentially outside the Stokes boundary layer of the compact object driving the background flow.We are interested in a time-averaged force on the particle that is (i) quadratic in the oscillation amplitude and (ii)involves contractions of the flow gradient ∇ ¯ U = ∇∇ ¯ ϕ and the flow curvature tensor ∇∇ ¯ U = ∇∇∇ ¯ ϕ . The onlydimensionless parameter in the problem not already specified by ¯ U is the Stokes number λ . Collecting the abovestatements, the only way to construct the time-averaged force (a vector) from the higher rank tensors ∇∇ ¯ ϕ and ∇∇∇ ¯ ϕ is by their contraction ∇∇ ¯ ϕ : ∇∇∇ ¯ ϕ . All other combinations are either of insufficient tensor rank or areidentically zero (since ∇ ¯ ϕ = 0). See [10] for similar arguments for flows without curvature. Including the correctdimensions, the time averaged force for any oscillatory potential background flow thus has the form F Γ κ = m f a p h∇∇ ¯ ϕ : ∇∇∇ ¯ ϕ i F ( λ ) . (26)Note that although the background flow is irrotational, the disturbance flow has a finite vorticity within the Stokeslayer around the particle. Under this general setting there is no requirement of axisymmetry of the background flow,so (25) as well as (1) in the main text apply to the generic case of an oscillatory potential flow background, and withthe same universal function F ( λ ). III. ACCURACY OF THE UNIFORMLY VALID EXPRESSION FOR F As stated in the main text, while the explicit functional form (24) of F (1)1 ( λ ) is rather lengthy, we Taylor expandin both the viscously dominated limit ( λ →
0) and the inviscid limit ( λ → ∞ ) to obtain F v = 916 r λ + O (1) , F i = 13 + O (1 / √ λ ) . (27)We construct a uniformly valid solution by simply adding the two leading order results, yielding F ( λ ) = + q λ .In Fig. 1(a), we plot the uniformly valid F (red curve) and the full theory represented by Eq. (24) (orange), alongwith the viscous and inviscid limits denoted by dashed lines. Figure 1(b) shows that the relative error between thered and orange curves is small ( . λ , even those far smaller or larger than practically relevant values. IV. FITTING PROCEDURE TO OBTAIN F FROM DNS
The DNS outputs (unsteady) particle trajectories as a function of time, with an initial transient period when theparticle starts from rest before periodic motion is fully established. As shown in Fig. 3(a) of the main text, theseoscillatory trajectories were time-averaged to obtain the steady particle dynamics r p ( T ), which is a function of the FIG. 1. (a) Logarithmic plot of the overall inertial force magnitude F ( λ ): the uniformly valid expression (red) closely tracksthe full solution (orange) while the inviscid theory (gray dashed) severely underestimates the inertial force even for moderatelylarge λ . (b) The magnitude of the percentage error between the uniformly valid and full solutions is small throughout the entirerange of λ , with a maximum error of ∼
8% where the two limits blend, as one would expect. slow time T = (cid:15) t . We fit these trajectories to Eq. (9) in the main text with F as the fitting parameter in order toobtain the simulation points of Fig. 4 of the main text. This was done in two ways: i) Taking a derivative with respectto time of the time-averaged trajectories from simulations, one obtains F directly from Eq. (9) of the main text. ii)The explicit analytical solution to Eq. (9) of the main text, namely, r p ( T ) = (cid:0) r p (0) − α T F (cid:1) / , was fitted to thetime-averaged trajectories from DNS using the method of least squares over a time interval that excludes the periodof initial transient behavior in simulations. We found that both these strategies yielded virtually identical values for F , which are displayed in Fig.4 of the main text. V. SIMULATION METHODS AND DETAILS
In order to perform three dimensional flow–structure interaction simulations with deforming geometries, we usethe remeshed Vortex Method (rVM) described in [11]. Here, we list our simulation methodology and parameters forcompleteness and reproducibility, as well as convergence tests used to assess simulations accuracy.
A. Fluid–structure interaction
We briefly recap the governing equations and numerical method used for our simulation. We consider incompressibleviscous flows in an unbounded domain in which two density-matched spherical bodies (i.e. bubble and particle) areimmersed. We denote the computational domain as Ω = Ω f ∪ Ω B , where Ω f is the fluid domain and Ω B = Ω b ∪ Ω p is the domain in which the bubble (Ω b ) and particle (Ω p ) reside, and denote the interface between the fluid and thebodies as ∂ Ω B . Both the bubble and the particle are then represented by mollified characteristic functions χ b ( x ) and χ p ( x ), respectively, on a regular Cartesian grid mesh such that χ b ( x ) = 1 for x ∈ Ω b , χ p ( x ) = 1 for x ∈ Ω p , and χ b ( x ) = χ p ( x ) = 0 for x ∈ Ω f . In order to avoid discontinuities, for each of the bodies, we smoothly blend the χ values using the mollification function χ ( d ) = d < − (cid:15) m , (1 + d(cid:15) m + π sin( π d(cid:15) m )) | d | ≤ (cid:15) m , d > (cid:15) m , (28)where d is the signed-distance to the body–fluid interface, and (cid:15) m is a user-defined smoothing parameter.We then solve the incompressible Navier–Stokes equation eq. (29) in its velocity–vorticity form ∇ · u = 0 , D ω Dt = ( ω · ∇ ) u + ν ∇ ω + λ penal ∇ × ( χ ( u B − u )) x ∈ Ω (29)where ω is the vorticity field, u is the fluid velocity field, u B is the body velocity and ν is the kinematic viscosity.Here λ penal (cid:29) λ penal ∇ × ( χ ( u B − u )) is the Brinkmann penalization term usedto approximate the no-slip boundary condition [11]. We note that while a bubble has a no-stress boundary conditionat the interface, which in general is different from no-slip, in the case of a bubble oscillating in pure breathing mode,where tangential boundary conditions have no effect, both would result in the same flow response. Indeed, using thismethod, fluid velocity within a body is forced to approach the body velocity (i.e. u ( x ) = u B ( x ) for x ∈ Ω B ). Thebody velocity u B can be decomposed into its rigid components of motion and the body deformation velocity field as u B ( x , t ) = u T ( t ) + u R ( x , t ) + u def ( x , t ), where u T and u R are rigid translational and rotational velocities, and u def isthe (imposed) deformation velocity field. The body rigid velocity (as a result of action from the fluid) is obtained viaa projection approach [11] where u T and u R are computed through the conservation of momentum in the system. Theimposed deformation velocity field used to prescribe the bubble’s breathing mode is u def ( x , t ) = x a b (cid:15)a b ω sin( ωt ) for x ∈ Ω B . This methodology based on remeshed vortex methods, penalization and projection has been validated acrossa range of fluid–structure interaction problems involving both rigid and deformable bodies, from bluff body flows tobiological swimming [11–15]. Recently, it has also been demonstrated in resolving the spatio-temporal scales relatedto oscillatory flow problems, particularly in viscous streaming settings involving individual and multiple arbitrary-shaped bodies, both in two and three dimensions [16, 17]. For a more detailed description of the numerical method,we refer to [11]. B. Simulation details
We simulate both the bubble and the particle as spheres of radii a b = 0 .
01 and a p = 0 . a p /a b = 0 . (cid:15) m = √ x used in the characteristic functioneq. (28), where ∆ x is the simulation grid size. The computational domain is initialized with a physical size of[ − , . a b × [ − . , . a b × [ − . , . a b . We then discretized the domain with a mesh of N =560 × ×
272 nodes, resulting in a uniform grid size of ∆ x = 1 . × − = 1 . × − a b in each direction.The domain boundary conditions are set to free-space (unbounded) boundary conditions so that u → x → ∞ .We initialize the particle at [ r p (0) , ,
0] and the bubble at [0 , , u def ( x , t ) = x a b (cid:15)a b ω sin( ωt ) for x ∈ Ω B ,so that the interface of the bubble moves with radial velocity (cid:15)a b ω sin( ωt ). Throughout this paper, we set ω = 16 π .The viscosity ν is set based on λ = a p ω/ (3 ν ) and the simulation is allowed to run until the particle’s steady velocityis achieved (typically 40–200 oscillation cycles, depending on λ where larger λ require longer time for transient effectsto vanish). Finally, we note that the bubble, particle and fluid are density-matched and density is set to unity. C. Implementation and resources
The simulation algorithm is implemented in Fortan90 and relies on MPI for distributed memory parallelism. Thesoftware relies on Parallel Particle Mesh library (PPM) [18] which provides a convenient abstraction layer over MPIparticle–mesh operations, mapping on processors, processor communication and load-balancing. The software alsouses FFTW3 library for Poisson solves and HDF library for data output, visualization and post-processing. Thesimulations performed in this paper typically run for 48–96 hours on 16 nodes, each with 64 threads, on the Stampede2supercomputing facility.
D. Resolution convergence test
It is important that we capture the different length scales involved in order to properly resolve the physics atplay. We first identify the different physical (bubble oscillations (cid:15)a b , particle oscillations, Stokes boundary layerthickness δ S ) and numerical (mollification length (cid:15) m = √ x ) length scales in this problem. Taking these scales intoconsideration, we then need to ensure that (i) the oscillations are properly resolved (i.e. ∆ x < oscillation amplitudes)and (ii) δ S measured from the bubble interface is not embedded in or under-resolved relative to the mollificationregion (i.e. δ S > (cid:15) m ).We conduct a resolution convergence test where we run a series of separate simulations with increasing resolution(hence decreasing ∆ x ). We then track the particle’s trajectory (via center-to-center distance between bubble andparticle) and observe a convergence towards a fixed trajectory, beyond which decreasing the grid size further does notsignificantly affect the results while requiring considerably larger computational cost. We illustrate the convergencebehavior in Figure 2 for the case of r p (0) = 2 and λ = 20, deliberately chosen from the larger λ regime in the test FIG. 2. Resolution convergence: Trajectory of particle for simulations with different ∆ x for the case of r p (0) = 2 and λ = 20. cases explored in this paper so that δ S is thin (hence requiring finer ∆ x to resolve). Here we note that a grid sizeof ∆ x = 1 . × − a b provides a good compromise between computational cost and accuracy as it resolves thephysical and numerical length scales reasonably well, effectively ensuring (i) ∆ x is finer than oscillation amplitudesand (ii) δ S > (cid:15) m . Therefore, throughout this paper, we use ∆ x = 1 . × − a b for all our simulations. E. Domain size convergence test
FIG. 3. Domain convergence: Trajectory of particle for simulations with different domain size for the case of r p (0) = 2 and λ = 1. In order to perform the simulations within feasible computational costs while ensuring that all the length scalesinvolved are properly resolved (see section V D), we adjust our simulation domain to a reasonable size such that effectsfrom domain boundaries do not affect the computed results. We perform a simple test by fixing ∆ x = 2 . × − a b and explore the boundary effects for different domain sizes. The case of study here is r p (0) = 2 and λ = 1 (hencea thick δ S that might interact with domain boundaries). We note that while a lower resolution is used here forexploration (see section V D), the δ S is still resolved in the simulations (i.e. δ S > (cid:15) m ) and the test still serves todemonstrate the effects from domain boundaries. Figure 3 shows the time-averaged trajectories for different domainsizes 0.375 L, 0.5 L and L, where L = [ − , . a b × [ − . , . a b × [ − . , . a b . We observe that the0particle trajectories do not change when doubling the domain size from 0.5 L to L. Therefore, we use the domain sizeof [ − , . a b × [ − . , . a b × [ − . , . a b for simulations conducted throughout this work. [1] M. R. Maxey and J. J. Riley, Equation of motion for a small rigid sphere in a nonuniform flow, The Physics of Fluids ,883 (1983).[2] R. Cox and H. Brenner, The lateral migration of solid particles in poiseuille flow—i theory, Chemical Engineering Science , 147 (1968).[3] B. Ho and L. Leal, Inertial migration of rigid spheres in two-dimensional unidirectional flows, Journal of fluid mechanics , 365 (1974).[4] K. Hood, S. Lee, and M. Roper, Inertial migration of a rigid sphere in three-dimensional poiseuille flow, Journal of FluidMechanics , 452 (2015).[5] L. D. Landau and E. Lifshitz, Course of Theoretical Physics Vol. 6 Fluid Mechanies (Pergamon Press, 1959).[6] C. Pozrikidis et al. , Boundary integral and singularity methods for linearized viscous flow (Cambridge University Press,1992).[7] P. M. Lovalenti and J. F. Brady, The hydrodynamic force on a rigid particle undergoing arbitrary time-dependent motionat small reynolds number, Journal of Fluid Mechanics , 561 (1993).[8] H. Stone, J. Brady, and P. Lovalenti, Inertial effects on the rheology of suspensions and on the motion of individualparticles, preprint (2001).[9] B. Rallabandi, A note on the Maxey–Riley equation in nonuniform flows, (2020), preprint.[10] S. Danilov and M. Mironov, Mean force on a small sphere in a sound field in a viscous fluid, The Journal of the AcousticalSociety of America , 143 (2000).[11] M. Gazzola, P. Chatelain, W. M. Van Rees, and P. Koumoutsakos, Simulations of single and multiple swimmers withnon-divergence free deforming geometries, Journal of Computational Physics , 7093 (2011).[12] M. Gazzola, C. Mimeau, A. A. Tchieu, and P. Koumoutsakos, Flow mediated interactions between two cylinders at finitere numbers, Physics of fluids , 043103 (2012).[13] M. Gazzola, W. M. Van Rees, and P. Koumoutsakos, C-start: optimal start of larval fish, Journal of Fluid Mechanics ,5 (2012).[14] M. Gazzola, B. Hejazialhosseini, and P. Koumoutsakos, Reinforcement learning and wavelet adapted vortex methods forsimulations of self-propelled swimmers, SIAM Journal on Scientific Computing , B622 (2014).[15] M. Gazzola, A. A. Tchieu, D. Alexeev, A. de Brauer, and P. Koumoutsakos, Learning to school in the presence ofhydrodynamic interactions, Journal of Fluid Mechanics , 726 (2016).[16] T. Parthasarathy, F. K. Chan, and M. Gazzola, Streaming-enhanced flow-mediated transport, Journal of Fluid Mechanics , 647 (2019).[17] Y. Bhosale, T. Parthasarathy, and M. Gazzola, Shape curvature effects in viscous streaming, Journal of Fluid Mechanics , A13 (2020).[18] I. F. Sbalzarini, J. H. Walther, M. Bergdorf, S. E. Hieber, E. M. Kotsalis, and P. Koumoutsakos, Ppm–a highly efficientparallel particle–mesh library for the simulation of continuum systems, Journal of Computational Physics215