Sorting droplets beyond the Rayleigh limit under acoustic-gravity forces: a theoretical approach

1 Introduction

Acoustophoresis, the manipulation of matter using acoustic fields, has its foundations since the late 18th century (Chladni, 1787; Faraday, 1831; Kundt and Lehmann, 1874). Significant advancements in the theoretical understanding of the acoustic radiation force ( F a c ) on particles within the Rayleigh limit (d ≪ λ, where d is the droplet diameter and λ is the wavelength of the wave) were made by King (1934), Yosioka and Kawasima (1955), and Gor’kov (1962), with Gor’kov providing a simplified framework for calculating this force from the potential function. The acoustic radiation force acting on a particle ( d ≪ λ ) subjected to a standing wave along the y direction (pressure p = p a ⁡ sin ( k y ) is given by,where E a c = p a 2 / ( 4 ρ avg c avg 2 ) is acoustic energy density, r is the droplet radius, E a c is the acoustic energy density, k = 2 π / λ is the wavenumber, e y is the unit vector along the wave direction, and φ is the contrast factor, which is a function of the density ratio and the compressibility ratio of particle and fluid. The acoustic radiation force (Equation 1) described above has been used to explain the levitation/suspension of particles (Gould and Coakley, 1974; Weiser and Apfel, 1982), droplets (Crum, 1971; Foresti et al., 2013), bubbles (Eller, 1968), and red blood cells (Baker, 1972; Dyson et al., 1971; Coakley et al., 1989) against gravity. In microfluidic systems, this force has enabled the development of two extensively utilized particle (such as beads, biological cells and droplets) sorting methods using acoustic fields (Lenshof et al., 2012; Wu et al., 2019), namely, size-based sorting and contrast factor-based sorting. Contrast factor ( φ ) based sorting is achieved when the physical properties of the particles differ, allowing the contrast factor of one particle to be made positive and that of another negative by adjusting the properties of the suspending liquid medium. This results in the particles experiencing opposing forces, with the former being directed toward the pressure node and the latter toward the pressure antinode, enabling their separation. This technique has been employed for separating polyethene particles of different densities (Gupta et al., 1995), isolating lipids from blood (Petersson et al., 2004), sorting encapsulated cells (Nam et al., 2012), separated bacteria from blood cell (Ohlsson et al., 2018), separated cancer cells from WBC (Karthick et al., 2018; Cushing et al., 2018). Size-based sorting, on the other hand, involves exposing two particles of different sizes (but with contrast factors of the same sign) to acoustic fields. The larger particle experiences a stronger acoustic force than the smaller one while moving against the drag force, allowing it to reach the pressure node or antinode more quickly due to its acoustic mobility ( ∼ φ r 2 ). This difference in movement, combined with the timely bifurcation or trifurcation of the outlet fluid streams, allows for effective size-based sorting. This size-based approach has been utilized for separating polystyrene particle mixtures of different sizes (Petersson et al., 2007), isolating blood components such as erythrocytes, leukocytes, and platelets (Shi et al., 2009), and separating cancer cells from blood (Thévoz et al., 2010; Augustsson et al., 2012; Li et al., 2015).

F a c = 4 π r 3 E a c φ k s i n 2 k y e y , (1)

The acoustic radiation force ( F a c ) has been extended further to account for particles beyond the Rayleigh limit (d > ∼ λ, droplet size greater than or approaching the order of the wavelength) in subsequent studies (Hasegawa, 1979a; Hasegawa, 1979b; Marston, 2017; Baasch and Dual, 2018; Pazos Ospina et al., 2022; Rueckner et al., 2023). These studies demonstrated that the behavior of particles beyond the Rayleigh limit differs significantly from that of particles within the Rayleigh limit. In a recent experimental study, it was shown that the acoustic power or acoustic energy density needed to suspend droplets beyond the Rayleigh limit varies significantly with their size, in contrast to smaller particles, where it remains independent of size (Thirisangu et al., 2023). They have also demonstrated a droplet sorting technique by suspending smaller droplets while allowing larger ones to settle. Following this, Thirisangu et al. (2024b) developed a unified theoretical framework in Eulerian form to examine the behaviour of droplets of arbitrary size subjected to acoustic fields, predicting their experimental observations of droplet migration, deformation, and splitting within microchannels. Extending the theoretical framework above by incorporating gravity, Thirisangu et al. (2024a), theoretically demonstrated that the suspension dynamics of droplets beyond the Rayleigh limit become highly size-dependent and non-monotonous, and could predict the experimental outcomes of droplet suspension observed in Thirisangu et al. (2023).

Building on the theoretical understanding gained from our previous work (Thirisangu et al., 2024a), we present a novel acoustic sorting method for droplets beyond the Rayleigh limit, specifically focusing on selective sorting of droplets across a wide range of sizes, including intermediate sizes. By controlling the wavelength and acoustic energy density, we demonstrate the ability to suspend larger droplets while smaller ones settle, or vice versa. Additionally, we present the methodology for sorting intermediate-sized droplets from a broad size range. Using this method, we further show the ability to sort droplets of the same size but with differing interfacial tensions through droplet splitting. Unlike size-based sorting methods, our approach enables droplets of the selected size to remain suspended, while droplets of other sizes settle indefinitely until the critical acoustic force is applied.

2 Physics of the problem

When immiscible inhomogeneous fluid (water droplet in the mineral oil) is subjected to gravity and acoustic fields as shown in Figure 1, the dynamics is governed by the incompressible mass continuity, momentum conservation, and advection-diffusion equations (Landau and Lifshitz, 1987; Thirisangu et al., 2024a), which are given below, where ρ is the density, v is the velocity, η is the dynamic viscosity, p is the pressure, ϕ is the phase fraction (ϕ = 0 for fluid 1 and ϕ = 1 for fluid 2), and f g = − ρ g e y is the gravitational force. f σ is the interfacial force which is expressed as σ κ δ s n, where σ is the surface tension, κ is the curvature, δ s is a surface Dirac delta function that becomes non-zero only at the interface between the fluids, and n is the unit normal. The parameter ϵ defines the width of the transition region where ϕ varies gradually from 0 to 1. To avoid singularities and enhance numerical stability, the droplet interface is modeled as diffusive in the simulations. To improve numerical stability, the parameter γ regulates degree of reinitialization or stabilization of the level set function. For the standing acoustic wave applied along the y-direction (fast time scale pressure p f = p a ⁡ sin ( k y ) and fast time scale velocity v f = p a i ρ 0 c 0 cos ( k y ) e y), the acoustic force ( f a c ) responsible for the movement of the droplets (inhomogeneous liquids) is expressed as follows (Thirisangu et al., 2024a; Thirisangu et al., 2024b),

∇ ⋅ v = 0 , (2a)

ρ ∂ t v + v ⋅ ∇ v = − ∇ p + η ∇ 2 v + f g + f σ + f a c , (2b)

∂ t ϕ + v ⋅ ∇ ϕ = γ ∇ ⋅ ϵ ∇ ϕ − ϕ 1 − ϕ ∇ ϕ ∇ ϕ , (2c)

f a c = − 2 k E a c ⁡ sin 2 k y Z ̂ e y . (3)

Where E a c = p a 2 / ( 4 ρ avg c avg 2 ) is the acoustic energy density, k = 2 π / λ is the wavenumber, p a is the pressure amplitude, Z ̂ = Z / Z avg where Z denotes impedance which is the product of density ( ρ ) and speed of sound ( c ), Z avg = ( Z 1 + Z 2 ) / 2, ρ avg is the average density of the fluids, c avg is the average speed of sound of the fluids and e y is the unit normal vector along the y direction.

3 Numerical implementation

The numerical investigation of droplet behavior under acoustic fields is conducted by solving the governing equations simultaneously. These include the continuity equation (Equation 2a), the momentum equation (Equation 2b), and the advection-diffusion equation (Equation 2c), with the acoustic force equation (Equation 3) incorporated as a body force term. The simulations are performed using COMSOL Multiphysics 6.0, accounting for gravitational and interfacial forces within the computational framework. The interface of the droplet is modeled using phase field module, which are solved alongside the laminar flow equations. A time-dependent study is solved using the Implicit Backward Differentiation Formula. The flow field variables and the phase field variables are solved using the segregated Parallel Direct Solver. A two-dimensional channel with a width ( w ) of 12 mm and a height ( h ) of 6 mm is modeled, where a droplet is suspended in a continuous medium. An acoustic wave with wavelength λ is applied along the height (y-direction) of the channel, while a no-slip boundary condition is imposed on all walls. Due to the closed nature of the domain, a pressure constraint is applied at a single point to stabilize the pressure field. The computational domain is discretized using a total of 804,596 grid elements, with a maximum grid size of 15 microns. Further refinement beyond this resolution does not significantly alter the flow field, as demonstrated in Figure 2. The material properties of the water droplet (ρ water = 1000 k g / m 3, c water = 1480 m / s, and η water = 0.8 m P a . s) and mineral oil (ρ oil = 857.5 k g / m 3, c oil = 1440 m / s, and η oil = 26.5 m P a . s) are taken from Hemachandran et al. (2019). In Sections 4.1, 4.2, selective sorting of the droplets beyond the Rayleigh limit is studied, when the interfacial force is strong enough to retain the droplet shape, where the following three cases are considered. Case I: Suspending the bigger droplet while settling the smaller droplet. Case II: Suspending the smaller droplet while settling the bigger droplet. Case III: Suspending the intermediate droplet while settling both the smaller and larger droplets. In Section 4.3, the sorting of the droplets, based on the interfacial tension is studied.

4 Results and discussion

When an acoustic field is not present, a higher-density droplet (water) settles in a lower-density continuous medium (mineral oil) due to the balance between gravity and drag forces. In the presence of both gravity ( f g ) and a plane standing acoustic wave field in y direction ( f a c ), the droplet can either suspend ( f a c ≥ f g ) or settle ( f a c < f g ). In this section, different sorting techniques for droplets beyond the Rayleigh limit ( d > ∼ λ ) are presented, using the interplay between acoustic ( f a c ) and gravitational ( f g ) forces.

4.1 Acoustic sorting beyond the Rayleigh limit

As shown in Figure 3A, for the given two droplet sizes (d 1 = 1.3 mm and d 2 = 1.8 mm) beyond the Rayleigh limit ( d > ∼ λ ), at a frequency f = 720 KHz (λ = 2 mm) and E a c = 12 J/m 3, the larger droplet ( d 2 = 1.8 mm ) remains suspended as the net acoustic force acting on it dominates the gravitational force, whereas the smaller droplet ( d 1 = 1.3 mm ) settles as gravitational force dominates over net acoustic force. By changing the frequency to f = 480 KHz (λ = 3 mm) and the acoustic energy density to E a c = 20 J/m 3, the opposite of the above is observed where the bigger droplet (d 2 = 1.8 mm) settles, and the smaller droplet (d 1 = 1.3 mm) suspends as shown in Figure 3B. Remarkably, as illustrated in Figure 3C, droplets of intermediate size (d 2 = 1.8 mm) can be suspended while settling both the smaller (d 1 = 1.3 mm) and bigger (d 1 = 2.6 mm) droplets using an E a c = 12 J/m 3 and f = 720 KHz (λ = 2 mm). This behavior is attributed to the dynamics of droplets beyond the Rayleigh limit, where the minimum acoustic energy density required to suspend the droplet (i.e., critical acoustic energy density, E cr) is highly size-dependent and exhibits a non-monotonic pattern, as shown in Figure 3. In contrast, for droplets within the Rayleigh limit ( d ≪ λ ), such sorting is not feasible, as demonstrated by Crum (1971), that E cr is independent of droplet size, meaning that at a given E ac, droplets of all sizes either suspend or settle. The E c r line separates the settling and suspending regime as shown in Figure 4A. If the applied E a c is less than the E c r, the droplet settles, if E a c > E c r, the droplet suspends.

To explain the results obtained in Figure 3, it is essential to understand the dependence of E cr on droplet size, shown in Figure 4. Although this has previously been explained in detail (Thirisangu et al., 2024a), we provide a summary nonetheless. When the droplets of size d > ∼ λ are subjected to acoustic fields, both positive and negative acoustic force regions exist within the droplet due to the presence of multiple nodes and antinodes (Figure 1). The portion of the droplet placed in the positive acoustic force region (quarter-wavelength below the pressure node) experiences an upward force (f a c is positive), while the portion placed in the negative acoustic force region (quarter-wavelength above the pressure node) experiences a downward force (f a c is negative). The partial cancellation of these opposing forces within the droplet reduces the net acoustic force acting against gravity, making the suspension of the droplet highly size-dependent. Additionally, as the droplet size increases, the net volume contributing to the positive acoustic force required to overcome gravity decreases significantly. This causes the critical acoustic energy ( E c r ) needed for suspension to increase significantly, reaching a peak as shown in Figure 4A. The subsequent decrease in E c r after the peak is due to an increase in the available volume contributing to the positive net acoustic force acting on the droplet. Beyond this peak, E c r follows a pattern of alternating decreases and increases, with the magnitude of each successive peak or trough exceeding the previous one (Figure 4A). The above behavior of droplets reveals an intriguing phenomenon: larger droplets of some sizes require less E c r to suspend against gravity, whereas smaller droplets need higher E c r within the same wavefield ( f ) as shown in the Figure 3.

With this understanding, the results of Figure 3 can be explained as follows. In Figure 3A, for the applied E a c = 12 J/m 3 & λ = 2 mm, the bigger droplet (d 2 = 1.8 mm) suspends, as it is in the suspending regime (applied E a c = 12 J/m 3 > E c r ( d 2 ) = 8.8 J/m 3 as shown in Figure 4A), and the smaller droplet (d 1 = 1.3 mm) settles, as it is in the settling regime (applied E a c = 12 J/m 3 < E c r ( d 1 ) = 58 J/m 3 as shown in Figure 4A). For the same applied condition of E a c = 12 J/m 3 & λ = 2 mm as shown in Figure 3C, both the bigger droplet (d 3 = 2.6 mm), and smaller droplet (d 1 = 1.3 mm) settle, as both are in the settling regime (applied E a c = 12 J/m 3 < E c r ( d 3 ) = 16.5 J/m 3 < E c r ( d 1 ) = 58 J/m 3 as shown in Figure 4A), but the intermediate sized droplet (d 2 = 1.8 mm) suspends, as it is in the suspending regime (applied E a c = 12 J/m 3 > E c r ( d 2 ) = 8.8 J/m 3 as shown in Figure 4A). It is to be noted that the plot given in Figure 4A cannot be used directly since the applied frequency (wavelength) has changed in Figure 3B. Instead, a similar plot to Figure 4A can be obtained for the required frequency (wavelength) from the non-dimensional plot in Figure 4B. For E a c = 20 J/m 3 & λ = 3 mm in Figure 3B, the bigger droplet (d 2 = 1.8 mm) settles, as it is in the settling regime (applied E a c = 20 J/m 3 < E c r ( d 2 ) = 87 J/m 3), and the smaller droplet (d 1 = 1.3 mm) suspends, as it is in the suspending regime (applied E a c = 12 J/m 3 > E c r ( d 1 ) = 13.2 J/m 3).

4.2 Sorting of arbitrary intermediate droplet sizes from a broad range

The following steps can be utilized to selectively sort intermediate-sized droplets ( d int ) from a broad range d min to d max.

Step 1 Determine λ 1 such that the d int / λ 1 ratio falls before the first local maximum of f a c / f g (m a x 1 as shown in the Figure 4B) and matches the value of the first local minimum of f a c / f g (m i n 1 as shown in the Figure 4B). Set the frequency f 1 = c / λ 1 to generate λ 1 and apply the corresponding E c r for d int. As a result of this Step 1, droplets larger than d int will settle, while droplets of other sizes, including d int, will remain suspended.

Step 2 Determine a λ 2 such that d int / λ 2 corresponds to the first local maximum of f a c / f g (m a x 1 as shown in the Figure 4B). Set the frequency f 2 accordingly. Apply the E a c slightly less than the E c r corresponding to d int. As a result of this Step 2, smaller droplets remain suspended while the droplets of interest ( d int ) begin to settle.

Step 3 While d int is settling, increase the E a c slightly more than the E c r corresponding to d int. At the end of Step 3, d int is prevented from settling, allowing it to remain suspended within the fluid medium between the smaller droplets (stay at the top) and the larger droplets (settle at the bottom). This enables the sorting of arbitrary intermediate-sized droplets.

4.3 Sorting of the droplet based on the interfacial tension

Since the maximum Bond number B o = Δ ρ g r 2 / σ is 0.063 (for a radius r = 1.5 mm and σ = 50 mN/m) for the results discussed until the previous sections, which is much less than 1, the interfacial force is strong enough to prevent the droplet from deforming or splitting. When interfacial effects weaken or the acoustic force becomes stronger, the opposing acoustic forces within the droplet compete, potentially leading to the droplet splitting into daughter droplets under suitable conditions (Thirisangu et al., 2024a). A combination of droplet splitting and the acoustic sorting method beyond the Rayleigh limit, as presented in Section 4.1, facilitates the sorting of droplets of the same size but with differing interfacial tensions, which will be discussed in this section.

Two droplets of the same size (d A = d B = 2 mm), as shown in Figure 5, are subjected to an acoustic field with a wavelength of λ = 4 mm. Droplet A exhibits a strong interfacial force with an interfacial tension of σ A = 50 mN/m (Bond number B o = 0.028), and droplet B exhibits a weak interfacial force with an interfacial tension of σ B = 2 mN/m ( B o = 0.699 ), are sorted in 2 steps (Note: The interfacial effects in practical fluids can be reduced by introducing surfactants). Step 1: Upon applying E a c ( 200 J/m 3 ) such that the droplet B with the weak interfacial tension splits into two daughter droplets, while the droplet A with the strong interfacial tension remains intact. Both droplet A and the daughter droplets ( B ′ ) remain suspended in the medium as long as E a c ≫ E c r , A > E c r , B ′ is applied, as shown in the Figure 5B. Step 2: Upon reducing the E a c to 10 J/m 3 ( E a c , A > E a c > E c r , B ′ ), the droplet A settles, and the droplet B suspends in the fluid medium, enabling the sorting of the droplets of same size with differing interfacial tension as shown in the Figure 5B. It is important to note that for the given wavelength (λ = 4 mm), the E c r for the droplet A (d = 2 mm) is 14.5 J/m 3, and the E c r for the droplet B ′ (d ≈ 1 mm) is 4.8 J/m 3.

5 Conclusion

The theoretical findings presented in this study on the acoustic sorting of droplets beyond the Rayleigh limit hold promise for enabling the development of advanced droplet-sorting devices. However, due to the high acoustic energy density ( E a c ) or pressure amplitude ( p a ) involved, the experimental implementation of this technique requires pressurized systems to prevent cavitation. The study relies on a few assumptions: first, a predetermined E a c is employed without solving the first-order acoustic fields, and second, the acoustic body force f a c is derived using the Boussinesq approximation (Δ ρ / ρ avg ≪ 1 and Δ c / c avg ≪ 1) (Rajendran et al., 2022). In order to produce more accurate results on droplet dynamics, it is necessary to solve the bidirectionally coupled equations of fast-time acoustic fields and slow-time hydrodynamic motion, incorporating the divergence of the Reynolds stress tensor as the body force (Thirisangu et al., 2024b).

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.