Acoustic tweezing of microparticles in microchannels with sinusoidal cross sections

Acoustic tweezing of bioparticles has distinct advantages over other manipulation methods such as electrophoresis or magnetophoresis in biotechnological applications. This manipulation method guarantees the viability of the bio-particles during and after the process. In this paper, the effects of sinusoidal boundaries of a microchannel on acoustophoretic manipulation of microparticles are studied. Our results show that while top and bottom walls are vertically actuated at the horizontal half-wave resonance frequency, a large mono-vortex appears, which is never achievable in a rectangular geometry with flat walls and one-dimensional oscillations. The drag force caused by such a vortex in combination with the tilted acoustic radiation force leads to trapping and micromixing of microparticles with diameters larger and smaller than the critical size, respectively. Simulation results in this paper show that efficient particle trapping occurs at the intermediate sinusoidal boundary amplitudes. It is also indicated that in a square-sinusoidal geometry there are two strong vortices, instead of one vortex. Sub-micrometer particles tend to be trapped dramatically faster in such a geometry than in the rectangular-sinusoidal ones.

Two fundamental nonlinear acoustic phenomena are involved in all the above cases and compete with each other to move suspended particles inside a microchannel. First, the acoustic radiation force, caused by the scattering of acoustic waves from the surface of particles, tends to focus particles on the nodal or anti-nodal plane of the acoustic standing waves [33][34][35][36] . Then, the Stokes drag force of the acoustic streaming velocity field, caused by viscous stresses in acoustic boundary layers, appears as steady fluid flow vortices which tend to defocus and spread out the suspended particles [37][38][39][40] . A critical particle diameter is determined as a crossover from the radiation force dominated region to the acoustic streaming-induced drag force dominated region 41,42 .
The acoustic radiation force depends linearly on particle volume, acoustic energy and contrast factor (a parameter used to determine the contrast between density and compressibility of a particle to the surrounding medium). Particles with diameters larger than the critical diameter are enforced with the acoustic radiation force. On the other hand, the acoustic streaming drag-force depends linearly on the acoustic energy, viscosity of the fluid, radius of particles and a factor related to the channel geometry. Tiny particles smaller than the critical size are affected by the acoustic streaming flows. Any changes in the boundary conditions dramatically affect the shapes of this vortices and therefore on the acoustophoretic motion of tiny particles. This phenomenon can be troublesome in some cases but if treated properly, it can be helpful in patterning or mixing applications 43,44 .
The geometry plays a key role in the shapes of the acoustic streaming patterns. Some biotechnological applications have been reported by several groups, considering geometrical effects. Focusing of sub-micrometer particles and bacteria in the cross-section of a nearly-square microchannel has been reported by Antfolk et al. 10 . A single roll streaming flow has been observed experimentally using two-dimensional acoustic streaming phenomena. Nama et al. 32 have suggested a microchannel with the sharp-edged sidewalls as a micromixer. Huang et al. 29 have optimized the design of such sharp-edged microchannels and reported that mixing performance varies at www.nature.com/scientificreports/ different frequencies and tip angles . Shortly afterward, a programmable microfluidic pump has been established using oscillations of tilted sharp-edge structures 31 . Ahmed et al. 45 have found that acoustically driven sidewalltrapped microbubbles would act as a fast microfluidic mixer. As another application for cavitation microstreaming, Yazdi and Ardekani 12 have shown that microbubbles inside a horseshoe structure are able to trap bacterial aggregations. Attempts have also been made to theoretically investigate the resonance frequencies of liquid-gas interface trapped over horseshoe-shaped structures and compare with experimental observations 46 . A theoretical model has been proposed 47 to quantify the resonant frequencies and viscous dissipation factors for water-air bubbles trapped in the circular microcavities. Further investigation 48 has shown that microstreaming induced by air-water membranes on defended oscillating membrane equipped structures (DOMES) leads to an active micromixer on lab-on-a foil devices. As the micron scale surface profiles dramatically affect the inner streaming vortices, Lei et al. 49 have performed a numerical study on the effects of the roughness of the fluid channels rather than being perfectly flat. Sinusoidally shaped, semi-circular, triangular, square, and trapezoidal surfaces were included. The effects of sinusoidal boundaries of a microchannel on acoustic streaming patterns have studied by Jannesar and Hamzehpour 50 and a classification of some repetitive patterns have reported. There are many experimental methods for fabrication of microchannels 51,52 . Soft lithography techniques, especially rapid prototyping and replica molding are some appropriate methods to design arbitrary geometries. In rapid prototyping, a design is first established in a computer-aided design program (CAD), then cast on PDMS materials.
In our study, creation of a single strong dominant vortex in the cross section of a microchannel with sinusoidal boundaries in one-dimensional oscillatory condition is numerically studied using a very simple mechanism. Such a vortex can never be produced in rectangular geometry with flat boundaries using one-dimensional oscillations. Also, the trajectory of several particles inside the cross section is simulated to investigate the functionality of the system at different frequencies for special curvatures in the cross-section of the microchannel.
The rest of the paper is organized as follows. In "Theoretical background" section the theoretical background is presented. This is followed in "Numerical model and boundary conditions" section by a description on the numerical model. The results are presented and discussed in "Results and discussion". Some conclusions are stated in the last section.

Theoretical background
In the absence of external body forces and heat sources, there are three important governing equations in the microfluidic systems as 53,54 The continuity Eq. 1a expresses conservation of mass where the mass density is ρ , the Navier-Stokes equation 1b expresses conservation of momentum where momentum density is ρv , and the heat-transfer equation 1c expresses conservation of energy where energy density is ρ(ε + 1 2 v 2 ) . Noteworthy, ε is internal energy per unit mass, T is the temperature field, v is the velocity fluid in the medium and k th is the thermal conductivity. Also, σ is the stress tensor of the fluid (Cauchy stress tensor) as where p is the pressure field, I is the unit tensor and τ is viscous stress tensor, which is expressed in terms of dynamic shear viscosity η and bulk viscosity η B .
Thermal effects are neglected in this work, because the thermal boundary layer thickness (thermal diffusion length), δ t , in fluids is much smaller than viscous boundary layer thickness (viscous penetration depth), δ ν 56 that are defined as where D th is the thermal diffusion constant, ω is angular frequency of the acoustic field and ν = η ρ is the dynamic viscosity.
First-order perturbation approximation. The initial (zeroth-order) state of the fluid is quiescent, homogeneous, and isotropic. Considering the external acoustic field as a perturbation of the steady state of a fluid, all the fields in standard first-order approximation can be expanded in the form g = g 0 + g 1 where g 1 is much smaller than g 0 . The perturbation parameter is defined as the dimensionless acoustic Mach number 53 as www.nature.com/scientificreports/ where ρ 0 is the unperturbed density of the fluid and ρ 1 is the first order perturbation term of density. In this study, processes are considered to be adiabatic. As a result, entropy is conserved for any small fluid volume and the acoustic pressure is independent of the first-order fluid temperature. If an acoustic wave constitutes tiny perturbations, the (isentropic) first-order expansion of the equation of state is defined as: p(ρ) = p 0 + ∂p ∂ρ s ρ 1 . This means that the small change in the pressure is related to the small change in the density by p 1 = ∂p ∂ρ s ρ 1 . The isentropic compressibility in classical fluid mechanics is defined as where c 0 is the speed of sound in water at rest. As a result, in the adiabatic limit we have ρ 1 = ρ 0 κ s p 1 . Additionally, it can be assumed that all first-order fields have a harmonic time-dependent form resulted from an imposed ultrasound field. The first-order equations in the frequency domain become 54,55 where σ 1 is given by Subscript 0 denotes the equilibrium conditions and subscript 1 represents the first-order quantities.
The time averaged second-order perturbation approximation. Using the second order approximation of the perturbation theory, all the fields can be expressed as g = g 0 + g 1 + g 2 . Then, the time averaged second-order perturbation approximation of governing equations are 54,55 Note that the stress tensor of the fluid using the second-order perturbation theory is given by Acoustophoretic forces acting on a small particle in a viscous fluid. A suspending spherical particle with the diameter of a, in the long-wavelength limit, δ ν , a ≪ (i.e. particle size and viscous boundary layer are much smaller than the wavelength of the imposing acoustic wave, ) encounters with two acoustophoretic forces.
The the time-averaged radiation force, which appears as a result of scattering of the sound wave on the particles, is given by 56,57 where v in 1 and p in 1 are the first-order pressure and velocity fields of the incident acoustic wave evaluated at a particle position. Asterisks denote complex conjugations. The prefactors f 1 and f 2 are the so-called mono and dipole scattering coefficients. In viscous fluids these are calculated as where κ p and ρ p are the compressibility and density of the particles.
The acoustic streaming flow field v 2 is caused by the attenuation of the acoustic wave inside a microchannel due to the boundary effects. The time-averaged streaming-induced drag force on a spherical particle with the radius of a, moving with the velocity of u far from the channel walls, is given by www.nature.com/scientificreports/ The crossover from the dominance of each force is defined through a critical particle radius. For a spherical particle inside a rectangular microchannel, the crossover diameter is where is a factor related to the channel geometry. The acoustophoretic contrast factor is given by , that contains material parameters. The monopole scattering coefficient f 1 is real valued and depends only on the compressibility ratio κ . The viscosity-dependent dipole scattering coefficient f 2 is in general a complex-valued number, and its real value is abbreviated as f r 2 (ρ,δ ν ) = Re f 2 (ρ,δ ν ) . If the radius of the particle becomes less than the critical size, streaming flow effects will become more significant and the radiation force will be suppressed. Otherwise, for particles larger the crossover limit, radiation force will be dominant and make particles move towards the pressure nodes or anti-nodes due to their contrast factor.

Numerical model and boundary conditions
As shown in Fig. 1, a microchannel with a rectangular cross-section of the height, h = 100 µ m and width w = 150 µ m in the yz plane is considered, with the top and bottom walls having sinusoidal shapes given by where A is the amplitude of the sinusoidal boundaries.
The governing Eqs. (5,7) are solved using the finite element method. We have used the weak-form-PDE of mathematics for both first and second-order equations. The conducted steps are as follows: First, the flow equations are written as source-free flux formulation, ∇ · J + F = 0 ; Then they are converted to the weak-form; Finally, the weak-form equations are solved by mathematics-weak-form-PDE module.
Considering J as a vector and F as a scalar (e.g. Eqs. 5a and 7a), the weak form of a source-free flux equation is given by where ψ is a test function.
For J as a tensor and F a vector (e.g. Eqs. 5b and 7b), the weak form of a source-free flux equation is where m is the mth component of a test vector . In all cases, the zero-flux boundary condition, J · n = 0 , is supposed; where n is the normal vector to the boundary surface. Also, all boundaries are considered as hard walls. Top and bottom walls are actuated by an external acoustic field. Therefore, the boundary conditions of the first-order velocity field are where v bc = ωd , ω = 2πf , and d = 0.1nm. The displacement of the oscillating walls in the z direction, d, is small enough to justify the use of the perturbation theory.
A zero-mass-flux boundary condition is considered for the second-order velocity field as www.nature.com/scientificreports/ The maximum mesh size in boundaries until 10δ ν and bulk are considered 0.5µm and 5µm , respectively. In both cases, the mesh element growth rate is 1.3 . Figure 2 presents a sketch of spatial mesh of our computational domain. The fluid inside the microchannel is considered to be quiescent water. Additionally, the physical parameters for water at a temperature of T = 25 • C and pressure of p 0 = 0.1013 MPa are presented in Table 1.

Results and discussion
Extensive numerical calculations were carried out to study the effects of the geometry of the microchannel walls on two-dimensional acoustophoretic trapping of the microparticles. The results show that the geometry plays a key role in creating a dominant strong vortex of the streaming pattern. The drag force caused by such a vortex inside a microchannel in combination with the tilted acoustic radiation force results in trapping of microparticles with diameters larger than the critical particle diameter. In the following sub-sections we present our results and discuss their implications.

Creation of a mono-vortex pattern.
In this section our simulation results for a microchannel with sinusoidal cross-section are compared with a rectangular one. The results for the pressure field at the half-wave resonance frequency, i.e. f v = c 0 /2h = 7.48 MHz is shown in Fig. 3a. The actuated boundaries are top and bottom walls and the direction of the oscillation is vertical. At this frequency, eight Rayleigh-Schlichting acoustic streaming vortices are created inside the cross-section of a rectangular microchannel, as expected 57 . Four Rayleigh bulk streaming patterns are distinguishable in Fig. 3b, but four other Schlichting streaming flows inside the boundary layers cannot be recognized in this scale. If the top and bottom boundaries are actuated in the z direction with the frequency of f h = c 0 /2w = 5 MHz which is the horizontal half-wave resonance frequency of the microchannel, the creation of acoustic standing waves or acoustic streaming flows are not expected. Figure 3c, d show the first-order pressure field and the acoustic streaming patterns for this frequency in a rectangular crosssection.
As a validity test for our simulations, we have compared some of our results, such as shown in Fig. 3a, b with what was reported by Muller et al. 57 . The results are the same.
A different geometry of cross-section, i.e. sinusoidal walls with an amplitude of A = h/50 is considered as shown in Fig. 1. Meanwhile, Fig. 4a indicates that at the frequency of f v , an approximately vertical standing wave is established. Also, the acoustic streaming pattern presented in Fig. 4b is very similar to Fig. 3b with the same order of velocity magnitudes.
For the horizontal half-wave resonance frequency, f h , the results are different from the rectangular case. While the top and bottom walls are oscillating vertically with the frequency f h , a standing acoustic wave is almost created in the horizontal direction as shown in Fig. 4c. In addition, boundary layers appear in the vicinity of the oscillating sinusoidal walls instead of the no-slip fixed boundaries. As shown in Fig. 4d, the resulting acoustic  www.nature.com/scientificreports/ streaming pattern, in this case, is a large mono-vortex. Such a pattern is never achievable in a rectangular geometry with one-dimensional oscillation. A similar vortex has been reported by Antfolk et al. 10 , when applying two-dimensional acoustic actuation to a semi-square cross-section. In their study, all four side walls are involved in the oscillation, whereas in this study, only the top and bottom walls are oscillating in one direction.
To verify the correctness of the solution, we prove the mesh convergence in our study. As reported by Muller et al. 57 , the time-averaged second-order velocity field v 2 converges considerably slower than the first-order fields. So we consider the magnitude of the second-order velocity field for two cutlines, z = 0 and z = + 45 µm , of  www.nature.com/scientificreports/ Fig. 4d with various mesh sizes. The results in Fig. 5 show that v 2 values converge for mesh sizes smaller than 1 µm for boundary and 10 µm for bulk domains. In our simulations, we use 0.5 µm and 5 µm for boundary and bulk domains, respectively, which are smaller than the above mentioned values.
To study the effects of actuation frequencies on the streaming patterns, some frequencies around two resonance frequencies of the system, f h and f v , are investigated. The results for frequencies between 4 to 8 MHz are shown in Figs. 6a-i. The mono-vortex patterns appear at the frequencies around f h (see Fig. 6a-e). Clearly, the strongest one is created at the frequency of 5 MHz.
Particle tracing through the mono-vortex. In this section, 80 polystyrene microparticles are traced inside the sinusoidal microchannel at the actuation frequency of f h . The particles are spreading homogeneously inside the bulk of the domain. The physical properties are shown in Table 2. Figure 7 represents the acoustic radiation and streaming drag force fields and the corresponding particle trajectories after 100 s for various particle sizes. These trajectories are effectively controlled by the strength ratio of two acoustophoretic forces. According to Eqs. (9) and (11), F rad and F drag are functions of a 3 and a, respectively. The results in Fig. 7 confirm that increasing particle sizes by one order of the magnitude, causes F rad to increases by three order of the magnitude, while F drag increases linearly. Figure 7g-i show that the small particles move in closed paths similar to the streaming pattern and eventually get approximately aligned, while the larger ones tend to move in straight lines towards the pressure nodes (see Fig. 4c for the pressure profile). The reason is that tiny particles are more affected by the acoustic streaming drag force of the fluid flow while large particles are impacted by the acoustic radiation force. The combination of radiation and streaming drag forces at the same order of magnitude causes particles to move in spiral trajectories and finally be trapped in the center of the microchannel.
The particle size has a significant impact on the trapping time. In all geometries, the critical size is obtained when two acoustic forces are equal. In this study, the viscous boundary layer thickness is , δ ν = 0.24 µm , the geometry dependent factor, , is approximately equal to 0.375, which is valid for a planar wall, and the contrast factor is ≈ 0.165 for the polystyrene particles 57 . For the frequency, f h = 5 MHz the critical size, 2a c , approximately is calculated as 1.24 µm.
With a careful attention to the above results, these two distinct regimes can be extracted: (1) Particles much smaller than the critical size, which can be used as a micromixer of suspension fluids. (2) Particles larger than the www.nature.com/scientificreports/ critical size, which rapidly concentrate at the center of the microchannel. These phenomena can be considered in setups to trap or mix microparticles and biological cells.
Effects of the sinusoidal boundary amplitudes. In this section, effects of the sinusoidal boundary amplitude, A, on particle trajectories are studied. Figure 8 shows snap-shots with various A values for particles with d = 2 µm at 100 s. The results show that for small and large values of A trajectories of particles will be almost aligned, while for intermediate amplitude values, particles will be trapped at the center of the microchannel. In addition, the streaming velocities for intermediate amplitudes are faster. It is noteworthy that there is an optimal value for A about 2.5 µm , in which the particles trapping at the center of the microchannel occurs more efficiently. Figure 9 shows a typical particle distance from the center of the microchannel, R, after 100 s for various A.

Effects of the microchannel widths and heights. Another investigation has carried out focusing on
the different values of microchannel widths. Figure 10 illustrates some examples. Results show that at the fixed vertical oscillation frequency of f h = 5 MHz , spiral paths become less compressed and the concentration of microparticles slows down gradually as a microchannel becomes wider. However, decreasing the width divide the mono vortex into two strong vortices. The microchannel height is another parameter that worth to be investigated. Our results show that in the fixed vertical oscillation frequency of f h = 5 MHz , as the height of the microchannel is increased, the mono vortex turns into two strong vortices (see Fig. 11a-c). If we continue to increase the height of the microchannel, the trapping effect will be decreased as shown in Fig. 11d,e. Additionally, the movement of particles slows down.
Square-sinusoidal geometry. In previous results, we find that in a square-sinusoidal cross-section, two strong vortices emerge instead of one vortex. Figure 12 represents the trajectories for the particles with a radius of a = 250 nm inside rectangular-sinusoidal and square-sinusoidal cross-sections after definite time scales. The results show that the tiny particles tend to concentrate dramatically faster (about 10 times) in square-sinusoidal geometry than the rectangular ones. Therefore, it is suggested to use square-sinusoidal microchannels for effective trapping of sub-micrometer particles. The particles in a square-sinusoidal cross-section at the frequency of f h completely aggregate after 300 s. www.nature.com/scientificreports/

Conclusion
In this paper, two-dimensional trajectories of micro-particles inside an acoustically actuated microchannel with sinusoidal top and bottom boundaries have been studied. Our results have shown while the top and bottom walls are vertically oscillating with the horizontal half-wave resonance frequency, an approximately standing acoustic wave is created in the horizontal direction. The resulting acoustic streaming pattern is a large mono-vortex which could never be achieved in a rectangular geometry with flat walls and one-dimensional oscillations. The drag force caused by such vortex inside the microchannel in combination with the tilted acoustic radiation force, leads to trapping of micro-particles with diameters larger than the critical size, after moving on spiral paths. Particles smaller than the critical size, follow close trajectories for a long time. For larger particles, the acoustic radiation force is dominant, which makes particles to be trapped. Whereas, for smaller particles it became negligible, but the acoustic streaming drag force is dominant. In addition, the simulation results have shown that for small and large amplitudes of sinusoidal walls trajectories of particles would eventually be aligned, while for intermediate amplitudes the particles would be trapped at the center of the microchannel.
Moreover, results of the square sinusoidal geometry analysis have shown that there are two strong vortices in this case instead of one vortex. In this setup the tiny particles tend to concentrate dramatically faster than the rectangular sinusoidal geometry.  www.nature.com/scientificreports/