## Abstract

Acoustic streaming, the generation of mean flow by dissipating acoustic waves, provides a promising method for flow pumping in microfluidic devices. In recent years, several groups have been experimenting with acoustic streaming induced by leaky surface waves: (Rayleigh) surface waves excited in a piezoelectric solid interact with a small volume of fluid where they generate acoustic waves and, as result of the viscous dissipation of these waves, a mean flow. We discuss the computation of the corresponding Lagrangian mean flow, which controls the trajectories of fluid particles and hence the mixing properties of the flows generated by this method. The problem is formulated using the averaged vorticity equation which extracts the dominant balance between wave dissipation and mean-flow dissipation. Particular attention is paid to the thin boundary layer that forms at the solid/liquid interface, where the flow is best computed using matched asymptotics. This leads to an explicit expression for a slip velocity, which includes the effect of the oscillations of the boundary. The Lagrangian mean flow is naturally separated into three contributions: an interior-driven Eulerian mean flow, a boundary-driven Eulerian mean flow and the Stokes drift. A scale analysis indicates that the latter two contributions can be neglected in devices much larger than the acoustic wavelength but need to be taken into account in smaller devices. A simple two-dimensional model of mean flow generation by surface acoustic waves is discussed as an illustration.

## 1. Introduction

The numerous applications of microfluidic technology, in biology and chemistry in particular, have stimulated a great deal of research about the physics of fluids at small scales (e.g. Nguyen & Wereley 2002; Squires & Quake 2005 or Tabeling 2006 for an introduction). Many microfluidic devices require to mix reacting solutions efficiently. As is well-recognized, this poses a challenge because the small values of the diffusivity of most chemicals make molecular diffusion impractically slow, while the low Reynolds numbers in typical devices preclude turbulent mixing. Much effort has therefore been devoted to the design of efficient micromixers, which typically rely on the formation of small-scale structures in the concentration fields to enhance the effect of molecular diffusion.

Although small-scale structures can be generated in steady flows, they develop much more rapidly in time-dependent flows through the process of chaotic advection (e.g. Ottino & Wiggins 2004; and other articles in the same journal issue). It is therefore highly desirable to develop methods for the stirring of microfluids that are flexible enough to generate time-dependent flows while interfering as little as possible with whatever chemical or biological processes might take place. One particularly attractive type of such methods is based on acoustic waves: (ultra)sound waves propagating in a fluid induce a time-averaged flow (or mean flow) through the nonlinear process known as acoustic streaming (e.g. Lighthill 1978*a*,*b*). Time-dependent mean flows can readily be obtained by varying the frequency of the waves and/or the location of the wave sources.

The potential of these methods has been demonstrated in a number of experiments (Moroney *et al.* 1991; Suri *et al.* 2002; Wixforth 2003; Guttenberg *et al.* 2004; Sritharan *et al.* 2006; Frommelt *et al.* 2008*a*; Du *et al.* 2009; Tan *et al.* 2009; Yeo & Friend 2009). Most of these share the same method of excitation of acoustic waves in the fluid, based on leaky surface acoustic waves (LSAWs). Surface acoustic waves (Rayleigh waves) are generated in a solid substrate by piezoelectric excitation. These waves propagate along the surface of the solid and, when this is in contact with a fluid above, radiate (leak) energy into the fluid in the form of acoustic waves. In turn, these waves generate a mean flow by acoustic streaming.

In this paper, we consider this physical mechanism and discuss its mathematical modelling. Several recent papers complement the experimental work with numerical simulations of the LSAWs and of the mean flow they generate (Köster 2007; Frommelt *et al.* 2008*b*; Tan *et al.* 2009; Antil *et al.* 2010). To compute the mean flow, these papers implement the formulation of Nyborg (1953, 1965) based on the time-averaged momentum equation. This formulation has the disadvantage of ignoring the implications of vorticity conservation: as other mechanisms of wave–mean flow interaction (e.g. Bühler 2009), acoustic streaming in a homogeneous fluid is strongly constrained by the fact that vorticity cannot be generated by purely inviscid processes. As a result, the mean flow forcing by acoustic waves depends entirely on dissipative processes. Because of the high frequency of the waves involved (typically in the range 10–200 MHz), these processes are very weak: indeed, the averaged momentum equation is completely dominated by the inviscid stress and pressure terms, which, however, balance and hence do not contribute to mean flow generation. To make this explicit, we use here the averaged vorticity equation in place of the momentum equation, following the original formulation of Eckart (1948) and Westervelt (1953). This formulation isolates the balance between the dissipation of the waves and of the mean flow which controls the streaming, and leads to the well-known observation that the streaming flow is independent of the value of the viscosity (or, more precisely, depends only on the ratio of the shear viscosity to the bulk viscosity). The vorticity formulation is particularly advantageous for the problem at hand because the mean flow is non-divergent up to negligible terms.

There are two forms of acoustic streaming whose relative contribution to mean flow generation depends on the size of the devices. *Interior streaming* is associated with the vorticity that appears when the irrotational acoustic waves dissipate in the fluid’s interior. *Boundary streaming*, on the other hand, is associated with a vortical part of the acoustic waves which is confined to thin boundary layers. To date, the boundary contribution to streaming in LSAW devices has received little attention (see Tan *et al.* 2009; however). In this paper, we consider both forms of streaming. We provide explicit expressions for the wave and mean flow structure inside the boundary layers and hence obtain the slip velocity that can be used as a boundary condition for the averaged vorticity equation. Because the boundaries of LSAW devices are oscillating, this requires extending the standard analysis of boundary streaming (e.g. Lighthill 1978*a*) in a manner analogous to Longuet-Higgins’ (1953) treatment of incompressible flows.

The main point of interest in the study of LSAW micromixers is the mean motion of fluid particles. This is governed by the Lagrangian mean flow which sums the Eulerian mean flow discussed above with the Stokes drift (e.g. Bühler 2009). The Stokes drift, which results from the nonlinearity of the advection equation, contributes directly to the mean advection; it also impacts the Eulerian mean flow through boundary conditions, since it is the total Lagrangian mean flow which satisfies no-slip boundary conditions (e.g. Bradley 1996).

The paper is organized as follows. In §2, we present the equations used to model the fluid, introduce their small-amplitude expansion, and define the mean flows to be evaluated. Section 3 discusses a solution procedure yielding the form of the acoustic waves while taking advantage of the small viscosity/high frequency of LSAW devices. The Eulerian and Lagrangian mean flows are considered in §4. The averaged vorticity equation governing the interior mean flow is derived, and the slip velocity is computed from the boundary-layer solution. The modelling of the solid supporting the LSAWs is described in §5: for simplicity, we restrict ourselves to a simple two-dimensional configuration and treat the solid as a linear elastic isotropic material. Section 6 presents results for this configuration with a choice of physical parameters that correspond to the experiments mentioned above. The paper concludes with a discussion in §7. Two appendices give technical details.

## 2. Formulation

The fluid is modelled by the compressible Navier–Stokes equations
2.1
and
2.2
for the fluid’s velocity ** u**, density

*ρ*and pressure

*p*(e.g. Landau & Lifschitz 1987). Here,

*μ*is the shear viscosity and

*μ*′=

*μ*

^{b}+

*μ*/3, where

*μ*

^{b}is the bulk viscosity. The shear and bulk viscosities are assumed constant; both contribute substantially to the dissipation of acoustic waves (e.g. Lighthill 1978

*b*). We ignore the thermal dissipation of acoustic waves as is appropriate for liquids, but not for gases. Equations (2.1) and (2.2) are complemented by an equation of state 2.3

Since the phenomenon of interest is usually characterized by small wave amplitudes, we can simplify (2.1)–(2.3) using a perturbation expansion based on a parameter *α*≪1 estimating the ratio of the density fluctuations to the background density. Thus, we introduce expansions of the form
2.4
and
2.5
where *ρ*_{0} is a constant, into the equations.

We assume a time-harmonic forcing with (angular) frequency *ω*. The *O*(*α*) solution, which we term the wave solution, can then be written in the form
2.6
where c.c. denotes the complex conjugate of the previous term. Similar expressions hold for *ρ*_{1} and *p*_{1}. The hatted fields are time-independent, complex fields that are computed in §3 below. Once they are determined, the *O*(*α*^{2}) fields generated nonlinearly can be obtained. Our focus is on the time-averaged response, and specifically on the Eulerian mean velocity
2.7
where the overbar denotes time average, and on the Stokes drift
2.8
where *ξ*_{1} is the linearized particle displacement, which satisfies ∂_{t}*ξ*_{1}=*u*_{1}. Both the Eulerian mean flow and Stokes drift contribute at the same order to the Lagrangian mean velocity
2.9
which gives the average motion of particles in the fluids to *O*(*α*^{2}) (e.g. Bühler 2009). Note that
2.10
where the divergence term vanishes only in special circumstances. Thus, is not necessarily the *O*(*α*^{2}) approximation to the mean mass flux .

## 3. Acoustic waves

We introduce the expansions (2.4) and (2.5) into (2.1) and (2.2), and collect the *O*(*α*) terms. Using the Helmholtz decomposition into irrotational and divergence-free flow
3.1
we obtain, after some manipulations, the damped wave equation
3.2
for the scalar potential *ϕ*_{1}, and the equation
3.3
for the vector potential (e.g. Nyborg 1965). Here, we have introduced the kinematic viscosities *ν*=*μ*/*ρ*_{0} and *ν*′=*μ*′/*ρ*_{0}, and the sound speed *c*^{2}=d*p*/d*ρ*|_{0}.

In the absence of interior forcing, the vorticity **∇**^{2}*ψ*_{1} is exponentially confined near the boundaries of the fluid, so that the wave solution in the interior of the fluid is irrotational. If the viscosity is small, the computation of the wave solution can therefore be much simplified using matched asymptotics. The condition for this to apply is that the boundary-layer thickness, given by
3.4
be much smaller than the inverse wavenumber of the acoustic waves, *κ*^{−1} say. In typical applications of LSAW-induced flows in water, *ω* is in the range 0.1–1 GHz and we estimate that *δ* is of the order of 40 nm; correspondingly, is of the order of 10^{−2}, so that an asymptotic treatment of the wave problem is well justified. The smallness of *δκ* also implies that the waves are only weakly dissipated in the fluid interior: solving equation (3.2) for plane waves gives their spatial damping rate
3.5
(Nyborg 1965). Thus, *δκ*≪1 clearly implies that *γ*/*κ*≪1, corresponding to weak interior dissipation. The viscous dissipation of the waves is, however, significant for the wave solution in devices whose size is of the order *γ*^{−1} or larger; it is also key to the mean flow generation. Note that the boundary layer is much thicker than the typical fluid-particle displacements associated with the LSAWs, which are of the order of a nanometre or less. Therefore, the condition *δκ*≪*α*, required to expand the equations in powers of *α* before expanding in *δ*, is met.

Taking advantage of the smallness of *δκ* gives a straightforward procedure for the computation of the wave solution. First the time-harmonic version of equation (3.2),
3.6
is solved to obtain the potential part of the wave field. It is not appropriate here to neglect the viscous term since the solution needs to be valid over distances that may be large compared with *γ*^{−1}. The boundary conditions, however, are those appropriate for an inviscid fluid: no normal flow across the boundary and, in the case of a dynamic boundary, continuity of normal stress and zero tangential stress (see §5 below). The result, then, is an approximation to the exact potential , with *O*(*δκ*) errors stemming from the approximate boundary conditions.

Once this approximation to is obtained, the leading-order, *O*(*δκ*) vector potential is found by solving equation (3.3) using a boundary-layer approach of the type pioneered by Rayleigh (1896, art. 352). This part of the solution ensures that the tangential velocity in the fluid also matches that of the boundary. Consistent with the linearization leading to equations (3.2) and (3.3), the boundary conditions are imposed at the undisturbed location of the boundary in the case of a dynamic boundary.

Let us consider a solid boundary at *z*=0 that is possibly oscillating with velocity The fluid velocity is written as
3.7
with a contribution of the vector potential that decreases exponentially as
3.8
tends to . The scaling (3.7) ensures that the leading-order contribution of the divergence-free flow is the *O*(1) tangential velocity needed to enforce no slip.

To leading order, the three components of the velocity in the boundary layer thus have the form
3.9
where we have written *Ψ*_{1}=(*Ψ*^{(1)},*Ψ*^{(2)},*Ψ*^{(3)}) and used *O*(*δ*) as a shorthand for *O*(*δκ*). Introducing equation (3.7) into equation (3.3), we obtain the leading-order equations
The solutions which decay as and ensure the continuity of the tangential velocity are
3.10
where we have introduced
3.11

Note that by equation (3.4) the viscous stress associated with the divergence-free flow, dominated by terms of the form , is *O*(*δ*) and hence negligible to leading order. This confirms the validity of imposing zero tangential stress at the boundary in the case of a dynamic boundary.

## 4. Mean flow

We now turn to the determination of the Lagrangian mean flow induced by the acoustic waves. The Stokes drift contribution to can be computed straightforwardly from the wave fields found in §3. The Eulerian flow , on the other hand, is determined by solving the time-averaged Navier–Stokes equations (2.1) and (2.2) to order *O*(*α*^{2}). We emphasize that the Eulerian mean flow depends in an essential manner on the presence of a viscous dissipation of the acoustic waves: for an inviscid fluid, the wave Reynolds stresses are exactly balanced by a pressure gradient and there is no mean flow forcing in the interior of the fluid (e.g. Lighthill 1978*b*). When the viscous effects are small (in the sense that *δκ*≪1), this balance holds approximately and dominates the averaged momentum equation. This equation is therefore not well-suited for the computation of the mean flow forcing, which is several orders of magnitude smaller than the Reynolds stress and pressure gradient. Instead, we follow Eckart (1948) and Westervelt (1953) and use the averaged vorticity equation. Because the acoustic waves are essentially irrotational in the fluid interior, this equation takes a very simple form, one that isolates the mean flow forcing and makes explicit that the mean flow response depends on the ratio *μ*^{b}/*μ* and not on *μ* and *μ*^{b} separately.

Since the acoustic waves have vorticity in the boundary layers, the mean flow generation there is controlled by a different balance than in the interior. An expression for the mean flow in the boundary layer is derived below; its limit away from the boundary layer provides the boundary condition that is needed for the interior in the form of a slip velocity.

### (a) Interior

We start by considering the mean flow equations away from the boundaries. Averaging equation (2.2) and using the important property leads to . Taking the divergence of equation (2.10) and averaging then gives
4.1
as the last term can be shown to vanish using . On the other hand, it follows from the definition of the Stokes drift and the irrotational nature of *u*_{1} that
using Einstein’s notation, equation (3.2) and . Integrating by parts leads to
4.2
As we are interested only in the mean flow to leading-order in *δ*, we can now focus on the Eulerian mean vorticity : its knowledge in the fluid interior, together with equation (4.2) and boundary conditions determine entirely.

We derive an equation for the Eulerian mean vorticity by first taking the curl of equation (2.1) to obtain
4.3
where ** ζ**=

**∇**×

**is the vorticity. We then introduce the expansions (2.4) and (2.5) and average. Since the**

*u**O*(

*α*) flow is irrotational in the interior, i.e.

*ζ*_{1}=

**∇**×

*u*_{1}=

**∇**

^{2}

*ψ*_{1}=0 there, the result is This can be further simplified into 4.4 using again that

**∇**×

*u*_{1}=0. This is the form obtained by Eckart (1948) and Westervelt (1953). A convenient alternative to equation (4.4) is obtained using that

*u*_{1}satisfies a wave equation up to

*O*(

*δ*

^{2}) terms and reads 4.5 This form, in which the mean-flow forcing is directly related to the acoustic energy flux , makes clear that depends on the viscosities only through their ratio

*ν*′/

*ν*(or equivalently

*μ*

^{b}/

*μ*).

In what follows, we use equations (4.2) and (4.5) to compute the Eulerian mean flow generated by LSAWs. We expect this formulation to be better conditioned than the solution of the averaged momentum equation: in the latter, the averaged wave stresses are balanced at *O*(1) by pressure gradients while much smaller, *O*(*δ*^{2}) terms involving viscous effects determine the Eulerian mean flow. The boundary conditions for equation (4.5) are found next by considering the boundary layers.

### (b) Boundary layers

Near the boundaries, the full *O*(*α*) velocity field of the form (3.9) needs to be taken into account. In the case of a solid boundary at *z*=0, the averaged *x*- and *y*-momentum equations should be solved to obtain the mean tangential velocity in the boundary layer and, by taking its limit as the boundary-layer coordinate , the slip velocity that serves as boundary condition for the interior equation (4.5). Let us consider the *x*-momentum equation, which reads
when only the dominant viscous term is retained. Now, the contribution to the left-hand side of this equation that is associated with the potential part of *u*_{1} is balanced by the pressure gradient up to *O*(*δ*^{2}) (as it is in the interior). Subtracting this part leads to
4.6
The boundary condition for this equation is obtained by averaging the no-slip condition on the (possibly moving) boundary and is
4.7
It is, therefore, convenient to solve, rather than equation (4.6), the equivalent equation for the Lagrangian mean flow, namely
4.8

This computation, taking advantage of the boundary-layer scaling (3.7), is carried out in appendix A. The result is
4.9
where all the functions of *z* are evaluated at *z*=0. Note that is *O*(1) throughout the boundary layer, in contrast with and which have (cancelling) *O*(*δ*^{−1}) contributions. Taking the limit in equation (4.9) gives the slip velocity for the interior solution,
4.10
The slip velocity in the *y*-direction is given by an analogous expression with (*x*,*y*) and (*U*,*V*) interchanged. The normal component of the Lagrangian-mean velocity vanishes to leading order in the boundary layer. Thus, the boundary conditions for the interior flow on a boundary *z*=0 reads
4.11

### (c) Complete solution

Let us summarize how the Eulerian and Lagrangian mean flows can be computed. Once the wave potential has been obtained, in general by solving for the coupled linear motion of the fluid and underlying solid, the Stokes drift away from the boundary layers can be directly evaluated from its definition (2.8). The Eulerian mean flow is computed by solving equations (4.2) and (4.5), with the boundary condition implied by equation (4.11) and similar on other boundaries. The Lagrangian mean flow is then computed as the sum of the Eulerian mean flow and Stokes drift.

It is convenient to decompose the Eulerian mean flow according to 4.12 Here is the interior-driven part, which satisfies equation (4.2) and (4.5) and no-slip boundary conditions, and is the boundary-driven part, which satisfies equation (4.2) and the homogeneous version of equation (4.5), namely 4.13 with boundary conditions of the form 4.14 so that equation (4.11) is satisfied. Note that is driven by a combination of boundary-layer and Stokes-drift effects, and that it is not necessarily tangential to the boundary.

We can use the results of the previous sections to estimate the order of magnitude of the various contributions to the Lagrangian mean flow. Taking *α*∼*ρ*_{1}/*ρ*_{0}∼*u*_{1}/*c* gives the estimates
4.15
for the Stokes drift (2.8), Eulerian mean flow in equation (4.5), and boundary-driven flow in equations (4.13) and (4.14). Here ℓ is an ‘outer scale’, that is, the scale over which quadratic correlations such as vary; it is fixed by the size of the wave source (e.g. the decay scale of LSAWs) or the smallest dimension of the fluid domain, whichever is the smallest. If ℓ is substantially larger than the wavelength, the interior-driven Eulerian mean flow dominates the other contributions to the Lagrangian mean flow, and can be used as a proxy for the Lagrangian mean flow; if ℓ is just a few wavelengths, however, the interior-driven, boundary-driven and Stokes contributions should all be computed to estimate the Lagrangian mean flow accurately. Both situations are likely to arise in LSAW devices: in the experiments of Frommelt *et al.* (2008*a*), for instance, ℓ is much larger than the wavelengths, while the channel device of Tan *et al.* (2009) is only a few wavelengths wide.

## 5. Fluid–solid coupling

We now examine the specifics of the generation of mean flows by LSAWs by considering the full coupling between the fluid and the underlying solid in a particular configuration. This configuration is motivated by the experimental work but taken as two-dimensional for simplicity: an isotropic linear elastic solid occupies the lower half plane *y*<0, and a compressible fluid occupies the first quadrant *x*>0, *y*>0. Surface acoustic waves in the solid are generated at some and are scattered and refracted at the solid–fluid interface, leading to the propagation of sound waves in the fluid which, in turn, generate a mean flow. See figure 1 for a schematic.

### (a) Solid

The solid is modelled as a dissipationless, isotropic, linear elastic solid. We write its governing equations using the velocity field instead of the more usual displacement field so as to make a complete parallel with the treatment of the fluid. Decomposing the velocity field in the solid as
5.1
where **∇**^{⊥}=(−∂_{y},∂_{x}), the potential and streamfunction satisfy the wave equations
5.2
where *a* and *b* are the longitudinal and transverse wave speeds of the solid. These are given in terms of the Lamé parameters *λ*^{s} and *μ*^{s} and density *ρ*^{s} of the solid by and . Here and in what follows, variables with the superscript s characterize the solid, while those with no superscripts continue to characterize the fluid.

We emphasize the difference between the spatial coordinates used to describe the motion in the fluid and in the solid: equations (2.1) and (2.2) for the fluid are written in Eulerian representation, with ** x** representing fixed positions in space; in contrast, equations (5.2) are written in Lagrangian representation, with

**labelling particles by means of their position in the undeformed solid. The distinction is crucial for the boundary conditions matching the motion in the fluid to that in the solid.**

*x*### (b) Boundary conditions

The wall bounding the liquid to the left is assumed rigid and fixed. Thus, we impose the condition
5.3
on the fluid velocity. At the interface between the solid and air (treated as vacuum), the tangential and normal components of the stress tensor, which satisfy
5.4
vanish: *T*^{s}=*N*^{s}=0 for *x*>0, *y*=0.

To write down the boundary conditions at the fluid–solid interface, we use the fluid displacements ** ξ**(

**,**

*x**t*) whose exact definition, ∂

_{t}

**+**

*ξ***⋅**

*u***∇**

**=**

*ξ***, can be approximated by the linearization ∂**

*u*_{t}

*ξ*_{1}=

*u*_{1}already used in §3. The continuity of the velocity field between fluid and solid is then written as 5.5 where

**′ is related to**

*x***according to**

*x***′=**

*x***+**

*x***(**

*ξ***′,**

*x**t*) and

**(**

*ξ***′,**

*x**t*) can be approximated by

*ξ*_{1}(

**,**

*x**t*). Similarly, the continuity of the tangential and normal stresses read 5.6 with

*T*

^{s}and

*N*

^{s}given in equation (5.4), and 5.7 As anticipated, the boundary-layer solution shows that equation (5.7) can be approximated as

*T*=

*O*(

*δ*) and

*N*=−

*p*+

*O*(

*δ*) so that the fluid is treated as inviscid in its interaction with the solid. Note that expanding equation (5.5) in powers of

*α*and averaging gives at

*O*(

*α*

^{2}) the condition (4.7) of vanishing Lagrangian mean velocity at the boundary. There is obviously no mean velocity in the solid.

## 6. Results

We present results for the wave fields and mean flows obtained in the configuration sketched in figure 1. The parameters, chosen in rough agreement with those of the experiments and numerical simulations of Frommelt *et al.* (2008*b*) and Köster (2007), are listed in table 1. The liquid is water at room temperature. The value used for its bulk viscosity differs from the value used in the references just mentioned and has been taken from recent experimental work by Dukhin & Goetz (2009).

Since the solid used in experiments is not isotropic and characterized by more than two elastic moduli, we have chosen the value of the Lamé parameters to obtain a LSAW phase speed that is close to the observed speed. Specifically, we have taken *λ*^{s}=*μ*^{s}=100 N m^{−3} corresponding to the longitudinal and transverse wave speeds reported in table 1. The dispersion relation of LSAWs can be written as
6.1
where *q*=*ω*/*k* is the (complex) phase speed (e.g. Tew 1992; Craster 1996). For the parameters in table 1 it has the solution *q*=4.27−5.4⋅10^{−2} i, corresponding to a wavenumber
6.2
and hence to LSAW wavelength and decay scale of about 29 μm and 720 μm, respectively. The speed of the LSAW is, therefore, *c*_{LSAW}=*ω*/*k*=4.3⋅10^{3} m s^{−1}.

### (a) Wave fields

To obtain the wave fields that result from the scattering of incident LSAWs, we solve the relevant system of three coupled Helmholtz equations (two in the solid, one in fluid) numerically. The simple geometry makes it possible to relate the velocities and stresses on the interface *y*=0 to the potentials , and on this interface. These relationships involve pseudodifferential operators that are best expressed using Fourier transforms in the *x*-direction. Using these, the problem can be reduced to one-dimensional pseudodifferential equations for , and which we solve using a pseudospectral discretization (e.g. Trefethen 2000). Details about the numerical procedure are given in appendix B.

We show in figure 2 the wave field in both the solid and the fluid. The figure displays only a small portion of the computational domain: the full computational domain is 5 mm long in the *x*-direction, with the forcing located around *x*=−1.25 mm (since we solve a one-dimensional problem, a grid in the *y*-direction is only needed for visualization). The scattering appears relatively simple and dominated by the LSAWs, although other types of modes are no doubt excited (see Craster 1996 for an analogous problem). The exponential decay of the LSAWs as *x* increases from zero is clearly visible, but the viscous damping of the acoustic waves is not because the damping length is much larger than the domain plotted.

As discussed above, the solution satisfies the continuity of the normal velocity between the fluid and solid, but not the continuity of the tangential velocity. The thickness of the boundary layer that forms to ensure the continuity of the tangential velocity is computed from equation (3.4) and found to be *δ*=46 nm. With a wavenumber in the fluid given by *κ*=*ω*/*c*=6.3×10^{5} m^{−1}, the dimensionless parameter estimating both the validity of the boundary-layer approach and the importance of interior dissipation is *δκ*=3×10^{−2}. We emphasize that the use of the analytic form of the solution in the boundary layer avoids the need for the exceedingly high resolution that would be needed for a fully resolved numerical computation.

Figure 3 displays the amplitude of the acoustic wave field in the fluid and provides a more detailed picture of the wave beam that is generated by the LSAW. The beam emanates from the corner (*x*,*y*)=(0,0), at an angle *θ* from the horizontal that can be computed from Snell’s Law as *θ*≈70^{°}. The figure shows a rather complicated interference pattern in the beam which results from the reflection on the rigid wall at *x*=0. In spite of the scattering in the solid and of the reflection in the liquid, the form of the waves on the interface *y*=0 is simple and well-described by a pure LSAW. To illustrate this, we compare in figure 3 the normal velocity with that predicted for a LSAW with complex wavenumber (6.2). The amplitude and phase of the LSAW have been fitted to match the numerical result. The agreement is excellent. Since the problem is linear, the amplitude of the interface displacements, of the order of 0.1 nm, is directly proportional to the strength of the forcing which we have chosen somewhat arbitrarily. We comment below on the magnitude of these displacements in connection with the amplitude of the mean flow generation.

### (b) Mean flows

We now consider the mean flows forced by the acoustic waves in the fluid. Since they are approximately non-divergent, the Stokes drift and Eulerian mean flows can be written in terms of streamfunctions , and , with
where **∇**^{⊥}=(−∂_{y},∂_{x}). The equations satisfied by *ψ*^{s}, and are readily obtained from equations (2.8), (4.5), (4.13) and (4.14). We have solved these equations using a finite-difference discretization of the bi-Laplacian that needs to be inverted to find and . The results depend strongly on the extent of the fluid domain in the *y*-direction. With the realistic values of the shear and bulk viscosities that we employ, the amplitude of the wave beam decreases over distances that are large compared with the typical size of experimental devices (the decay scale is estimated from equation (3.5) as *γ*^{−1}≈2 mm). As a result, treating the fluid domain as infinite in the *y*-direction would lead to an unrealistically strong Eulerian mean flow. We have therefore chosen to consider a bounded domain of size 1 mm in the *y*-direction. As we use the wave field computed in a semi-infinite domain, we neglect the reflected beam that appears on the upper boundary and whose amplitude is about 1/2 that of the main wave beam. The structure of the mean flow and its magnitude are not expected to be modified in an essential way by the reflected beam. In the *x*-direction, we continue to consider the domain as semi-infinite: the numerical computation requires to take a finite size, here 2.5 mm, but this has only little impact on the results.

Figure 4*a* shows the streamfunction of the interior-driven Eulerian flow. The structure is similar to that observed in experiments and previous numerical models: a strong clockwise vortex is established to the right of the acoustic wave beam (figure 4*b*), with a much weaker counterclockwise circulation to the left (figure 4*a*). The qualitative features of this structure are relatively insensitive to the size of the domain, but the strength of the circulation is not. Here, the mean velocities obtained are of the order of 10 mm s^{−1}. These velocities are large when compared with those reported in Köster (2007), even though the vertical displacement of the boundary that we have assumed is a factor 10 smaller than the displacements assumed in this paper for a similar set-up. It is not clear what the reason for the differences may be, although it may be related to the lack of a treatment of the boundary layers in Köster (2007). We note that our numerical results are consistent with the scaling (4.15): using our earlier estimate *cα*≈2×10^{−3} mm s^{−1} and ℓ*κ*∼3×10^{2} for ℓ=1 mm gives mm s^{−1}. Direct comparison with the velocities reported in experiments is delicate because the vertical displacements of the boundary are difficult to measure. It should also be borne in mind that lateral confinement, ignored in our two-dimensional set-up, is important for the amplitude of the mean velocity field, since it reduces ℓ.

Since ℓ*κ* is large, the interior-driven mean flow dominates the Stokes drift and the boundary-driven flow. This is confirmed by figure 4*b* which shows the sum . The largest velocities are reached near the boundary *y*=0 and are of the order of 10 μm s^{−1}, much smaller than the interior-driven velocities, as expected. The complete Lagrangian circulation is therefore very well approximated by the Eulerian circulation in figure 4*a*.

We emphasize that this is not a general feature of acoustic streaming in LSAW devices but one that depends on the geometry of the device: here, both the domain and the region of wave forcing provided by the LSAW are large compared with the acoustic wavelength, so that ℓ*κ*≫1. If the domain is smaller, then the three contributions to Lagrangian mean flow can matter. We illustrate this by considering a domain of 0.1 mm in the *y*-direction instead of the 1 mm used so far. The size of the domain is now comparable to that used in the experiments of Tan *et al.* (2009). We continue to use the unbounded form of the acoustic waves in spite of the strong reflections; the results should therefore be treated as illustrative rather than as proper predictions of the mean flows in the narrow domain. Figure 5 shows the streamfunctions and . The interior-driven velocities reduced by a factor of the order of 10^{3} when compared with those in the larger domain, and although they still dominate the interior-driven contributions, it is only be a factor of 10 or so. As a result, the Lagrangian mean flow is not accurately approximated by the interior-driven Eulerian flow only, in particular near the boundaries. This effect will be reinforced in domains that are smaller still, or that are confined in two or three directions.

## 7. Discussion

This paper is motivated by a series of recent experiments which have demonstrated the potential of streaming by LSAWs to generate mixing flows in microfluidic devices. It discusses the computation of the Lagrangian mean flow, which controls the trajectories of fluid particles, for general wave-forcing protocols. The problem is formulated using the averaged vorticity equation, following Eckart (1948) and Westervelt (1953). This formulation extracts the balance between wave dissipation and mean flow dissipation that is at the core of the streaming process by eliminating large, cancelling terms in the averaged momentum equation. Because of this, it appears preferable for numerical computations to the direct use of the momentum equation.

The formulation takes advantage of the weakness of the wave dissipation for realistic parameter values and employs a boundary-layer approach. In the interior of the fluid, the waves are irrotational and can be computed numerically by solving for a scalar potential. The waves have a significant vortical part only in thin boundary layers. Analytical expressions for the solution there are readily obtained in terms of the scalar potential, making it unnecessary to resolve the small-scale motion near the boundaries explicitly. The vortical part of the waves in the boundary layer contributes to the mean flow generation in the form of an effective slip velocity that serves as boundary condition for the interior mean flow. We derive an expression for this slip velocity which takes into account the oscillatory motion of the boundary that forces the acoustic waves in the fluid.

The Lagrangian mean flow is the sum of the Stokes drift and of the Eulerian mean flow. The latter is naturally separated into two contributions: an interior-driven one, and a boundary-driven one which is associated with a combination of the slip velocity and Stokes drift at the boundary. We emphasize that the Stokes drift and boundary-driven Eulerian flows can be as large as the interior-driven Eulerian flow in LSAW devices if the size of the devices is comparable to the acoustic wavelength.

Even when much smaller than the interior-driven flow, the boundary-driven flow can be important for a micromixer. As discussed by several authors, the no-slip boundary condition of viscous flows makes such flows inefficient mixers compared with slip flows because the scalar is trapped in sizable boundary layers near the walls (Lebedev & Turitsyn 2004; Salman & Haynes 2007). However, flows driven by acoustic streaming can be expected to behave essentially as slip flows as far as mixing is concerned, because the boundary streaming results in a finite slip velocity immediately outside an exceedingly thin boundary layer. This is a potentially valuable property of LSAW mixers which deserves future investigation.

In this paper, we have followed most earlier work on streaming by LSAWs in treating the nonlinearity parameter *α* as small. However, in experiments such as those of Tan *et al.* (2009) and in our computations, mean velocities of the order of several millimetres per second are reached, corresponding to Reynolds numbers of order one. It would therefore be useful to assess the importance of the inertial terms on the mean flow. This is possible using an expansion based on a large Strouhal number instead of small *α*. The developments in this case can be expected to be very similar to those presented in this paper, with the important difference that the inertial terms are retained in the mean flow equations (e.g. Riley 2001). When inertial effects are important, the mean flows can become unstable and unsteady, leading to improved mixing performances.

Another possible extension of our results concerns streaming in very thin domains. This is motivated by experiments in which LSAWs generate mean flows in drops that are confined between two plates separated by a narrow gap (Guttenberg *et al.* 2004; Frommelt *et al.* 2008*a*). The standard Hele–Shaw approximation (e.g. Batchelor 1967) can then be used to compute the interior-driven flow.

Also of interest would be the development of a ray-tracing approximation to the solution similar to that of Frommelt *et al.* (2008*b*). The scale separation required for such an approach is valid in the bulk of the fluid, but special attention needs to be paid to possible diffraction effects as occur, in particular, at the corner (*x*,*y*)=(0,0) of our model. The ray-tracing approach should also take into account the constraint of irrotationality of the acoustic waves in the fluid interior.

Finally, we note that the extreme thinness of the boundary layer (a few tens of nanometres) suggests that some non-negligible slip can occur at the fluid–solid interface. It would be interesting to analyse the possible consequences of such a slip by replacing the no-slip boundary condition used in this paper by a more sophisticated model such as the Navier boundary condition (e.g. Tabeling 2006).

## Acknowledgements

J.V. acknowledges the support of a Leverhulme Research Fellowship and the hospitality of the Courant Institute where part of this research was carried out. He thanks J. Cosgrove for stimulating discussions. O.B. gratefully acknowledges financial support under grant DMS-0604519 of the National Science Foundation (USA).

## Appendix A. Mean flow in the boundary layers

In this appendix, we solve the averaged momentum equation (4.8) in the boundary layer near *z*=0. We start by writing equation (3.7) as *u*_{1}(*z*)=*u*^{ϕ}(*z*)+*u*^{ψ}(*Z*), where *u*^{ϕ}=(*u*^{ϕ},*v*^{ϕ},*w*^{ϕ})=**∇***ϕ*_{1} is irrotational and *u*^{ψ}=(*u*^{ψ},*v*^{ψ},*δw*^{ψ})=*δ***∇**×*Ψ*_{1} is divergence free, and where we have omitted the dependence on (*x*,*y*,*t*) for simplicity. The corresponding displacement fields are written as *ξ*^{ϕ}=(*ξ*^{ϕ},*η*^{ϕ},*ζ*^{ϕ}) and *ξ*^{ψ}=(*ξ*^{ψ},*η*^{ψ},*δζ*^{ψ}). At this point, we treat all these fields as known exactly: it is only in the later stages of the derivation that their approximations up to *O*(*δ*) errors will be used. We will, however, use from the outset the approximations
A 1
which follow from equations (3.3), (3.4) and (3.8).

We now obtain the derivative that appears on the right-hand side of equation (4.8). This contains three terms, computed as
A 2
A 3
and
A 4
The computation is straightforward, if tedious: it uses extensively integration by parts in time to eliminate the displacements in favour of the velocities, as well as equation (A 1) and its analogue for *ξ*^{ψ} and *η*^{ψ} to eliminate second derivatives in *Z*. Note that the *O*(*δ*^{−3}) term in equation (A 4) implies that the leading-order approximations to *w*^{ϕ} and *u*^{ϕ} are insufficient to estimate equation (A 4) with the same *O*(*δ*^{−2}) accuracy as equations (A 2) and (A 3). However, as shown below, this term is cancelled in the right-hand side of equation (4.8) by an equal and opposite contribution to .

The other terms on the right-hand side of equation (4.8) are next obtained in the form A 5

Rewriting equation (4.8) as
A 6
using equations (3.4) and (3.8), and introducing the results (A 2)–(A 5) now leads to an explicit equation for the Lagrangian mean velocity . Since this equation applies to the boundary layer, the *z*-dependent terms of equation (A2)–(A5) are evaluated at *z*=0. After a number of simplifications, some involving the condition **∇**⋅*u*^{ψ}=0, this equation can be written as
A 7
As anticipated, the largest terms in equations (A 4) and (A5) have cancelled out, and knowledge of the leading-order approximation to the wave fields is sufficient to obtain with negligible, *O*(*δ*) errors.

The right-hand side of equation (A7) is made explicit by introducing the harmonic form of the wave solution. Specifically, we introduce , and , with the right-hand side given in equation (3.10), into equation (A7). Only one term involves the vertical velocity *w*^{ψ} (or rather the vertical displacement *ζ*^{ψ}). This is expressed in terms of *u*^{ψ} and *v*^{ψ} using the incompressibility condition ∂_{Z}*w*^{ψ}=−(∂_{x}*u*^{ψ}+∂_{y}*v*^{ψ}). The result is
A 8

Integrating this equation twice with respect to *Z*, imposing at *Z*=0 and boundedness as finally leads to equation (4.9). We have verified equation (4.9) in the classical problems of a plane wave and a standing wave above a flat wall (e.g. Lighthill 1978*a*).

## Appendix B. Numerical method

To find the flow in the fluid interior, we need to solve equation (3.6) for the potential, together with the time-harmonic version of equation (5.2), namely B 1

The boundary conditions are obtained from those in §5*b*. In particular, the stress tensor is given by ,

In numerical computations, we generate incident waves by imposing a non-zero on the boundary of the solid in a small interval around some *x*≪−1. We also handle the normal component of the boundary condition (5.3) by solving for the fluid in the whole upper half plane *y*>0 and imposing the symmetry .

The boundary-value problem determining the wave fields can be reduced to a one-dimensional problem using Fourier transforms to solve the interior equations (3.6) and (B 1) explicitly. Thus, we write
B 2
B 3and
B 4
Here, , and , with a choice of branches which ensures that, for , decay or radiation boundary conditions are satisfied in the solid as and in the fluid as . Note that we have defined *m*_{a} and *m*_{b} differently from *m*: this is because our main interest is for LSAWs which are evanescent in the solid and propagating in the liquid and thus, with our definitions, are characterized by *m*_{a}, *m*_{b} and *m* with much larger real parts than imaginary parts for the scales *k* that are primarily excited.

With the definitions (B2) it is straightforward to write down all fields in terms of , and , and hence in terms of (pseudodifferential transforms of) , and . Expressing the boundary conditions in this manner reduces the boundary-value problem to a set of one-dimensional pseudodifferential equations for , and . In principle, these equations could be solved analytically, using the Wiener–Hopf technique (Craster 1996). Here we solve them numerically, using a straightforward pseudospectral collocation method (Trefethen 2000). The infinite domain in *x* is replaced by a large periodic domain, and the three unknown fields , and are discretized on a regular grid. The pseudodifferential operators are then approximated by matrices which are computed using fast Fourier transforms.

- Received August 31, 2010.
- Accepted November 8, 2010.

- This journal is © 2010 The Royal Society