On the streaming in a microfluidic Kundt's tube

Abstract We derive an analytical solution for the acoustic streaming inside a rigid tube resulting from a pseudo-standing wave field, generated by two counterpropagating travelling waves. We solve the second-order axisymmetric problem that follows from the perturbation expansion of the governing equations. In the process, we impose no restriction on the diameter of the tube with respect to the thickness of the viscous boundary layer and acoustic wavelength. The derived solution is then used to study the evolution of streaming patterns inside the tube when geometrical and material parameters are varied. We show how the Schlichting streaming torus at the wall bounds the Rayleigh streaming near the axis of the tube. Decreasing the ratio ($\varXi$) of the tube radius to the viscous boundary layer thickness gradually expands the Schlichting streaming, suppressing the Rayleigh streaming. Considering the average mass transport velocity, the Rayleigh streaming vanishes at the critical ratio $\varXi _{{S}}^{{M}}=6.2$. The critical ratio is independent of fluid properties in the limit of large acoustic wavelength relative to the radius of the tube ($\varLambda$). When $\varLambda$ decreases towards unity, large-scale Eckart-like streaming develops near the axis, superseding the Rayleigh streaming, while the Schlichting streaming remains at the wall. In addition, we demonstrate the relevance of the compressibility of the streaming flow and of the full inclusion of the spatial variation of the Reynolds stress that acts as the streaming source. The study is especially relevant for microfluidic systems, wherein the viscous boundary layer can reach significant thicknesses.

assumed a two-dimensional geometry extending infinitely in the out-of-plane direction. More than half a century later, Westervelt (1953) and Nyborg (1953) extended Rayleigh's approach to compressible fluid at the first order of perturbation expansion. A different case was considered by Schlichting (1932), who solved the incompressible streaming problem near the wall. The so-called Schlichting streaming is known to appear between the wall and the Rayleigh streaming in the bulk of a fluid. Its direction of rotation is opposite that of the Rayleigh streaming.
The two-dimensional theory was extended by Hamilton, Ilinskii & Zabolotskaya (2003a) to a channel of an arbitrary width relative to the viscous boundary layer thickness, extending infinitely in the out-of-plane direction. They showed that the oppositely directed Schlichting streaming vortices, between the Rayleigh vortices and the wall, could become dominant as the thickness of the viscous boundary layer was increased. More recently, Doinikov, Thibault & Marmottant (2017) considered a fluid bounded between an elastic solid wall and a reflector, posing no restrictions on the width of the channel. The standing wave in their case was generated via two counterpropagating leaky waves, originating from the solid wall.
The first to extend Rayleigh's approach to cylindrical geometry were Schuster & Matz (1940). They derived a very concise expression for the streaming velocity and analysed the magnitude of the streaming. (The expression for the streaming velocity by Schuster & Matz (1940) is later used for comparison, and is given explicitly in appendix G.) The latter was found to scale with the square of the pressure amplitude inside the tube, which was confirmed experimentally. The solution of Schuster & Matz (1940), although very convenient, relies on many assumptions: the viscous boundary layer is assumed small compared to the tube radius; the acoustic wavelength is assumed large with respect to the tube radius; the wave is assumed undamped along the tube axis. All of this considerably restricts the applicability of their solution. In more recent years, others have also analysed various aspects of streaming in a tube (Qi, Johnson & Harris 1995;Menguy & Gilbert 2000;Bailliet et al. 2001;Hamilton, Ilinskii & Zabolotskaya 2003b).
Recently, Baltean-Carlès et al. (2019) analysed the contributions of individual parts of the streaming source terms, and how they are affected by the finite length of the tube relative to the diameter. They discovered that the standard exclusion of terms is invalid when the length of the tube becomes comparable to its diameter.
Here, we derive a general analytical solution for the acoustic streaming that results from a pseudo-standing wave field in a tube of an arbitrary diameter. The acoustic field is generated by two counterpropagating decaying travelling waves, which is a common approach for generating acoustic fields in standing surface acoustic wave (SSAW) systems (Devendran et al. 2016). We use an approach similar to that of Doinikov et al. (2017), but consider a cylindrical geometry and a rigid wall. We also extend the formulation by considering the complete spatial variation of the Reynolds stress as a source term for the streaming, and by considering also the irrotational component of the streaming velocity (i.e. assuming that the streaming flow is compressible). First, we derive the dispersion equation to determine the wavenumber in the fluid. Then we solve the first-and the second-order problems that follow from the perturbation expansion. The solution is valid inside as well as outside the viscous boundary layer, which is of unrestricted thickness. It therefore covers also very thin capillaries, where the viscous boundary layer is comparable to the radius. The solution is then used to analyse the evolution of streaming vortices with respect to changes in the fluid viscosity, tube diameter and excitation frequency. In addition, we analyse the importance of the irrotational part of the streaming velocity and of individual parts of the spatial variation of the Reynolds stress (Lighthill 1978) that acts as a source term in the streaming equations. Some of those parts are often assumed to be negligible (e.g. Rayleigh 1884;Schlichting 1932;Lighthill 1978;Doinikov et al. 2017). The impact of different assumptions on the streaming with respect to the increasing viscous boundary layer thickness relative to the tube diameter is investigated.

Problem statement and assumptions
We assume that the fluid, initially at rest, is situated in an infinite rigid tube with an inner diameter of 2a. The geometry of our problem is parametrized in cylindrical coordinates, as depicted in figure 1, and is symmetric with respect to the z-axis oriented along the centre of the tube. In the fluid, there are two counterpropagating spatially-decaying harmonic travelling waves along the axis of the tube, which form a pseudo-standing wave when superimposed. The acoustic field resembles a standing wave near the plane of symmetry (z = 0), approaching a travelling-wave-like behaviour far away from the plane of symmetry. We neglect thermal effects, and the motion of a barotropic compressible viscous fluid is therefore governed by the Navier-Stokes equations To constrain the problem, the no-slip boundary condition is imposed at the wall. By applying a regular perturbation technique (Hamilton & Blackstock 1998, p. 281), the problem can be solved in successive steps of increasing order in terms of the small Mach (2.4) with the amplitude of the fluid velocity denoted by v a , and the speed of sound in the fluid by c f . The perturbed variables can then be written in the form of a series, namely = 0 + ε 1 + ε 2 2 + · · · , where the subscript denotes the order in the perturbation expansion. In our formulation, we will use 0 = 0 , 1 = ε 1 , 2 = ε 2 2 , and solve the associated first-and second-order problems. Following our assumptions, v 0 = 0.

Fluid motion at the first order
The perturbation expansion of the governing equations to the first order results in the momentum equations and the continuity equation with equilibrium density ρ 0 . The first-order pressure and density are related through the linear equation of state Using the Helmholtz decomposition v 1 = ∇ϕ 1 + ∇ × ψ 1 (3.4) (e.g. Blackstock 2001, p. 76), the first-order velocity field can be separated into a sum of the gradient of a scalar velocity potential ϕ 1 and the curl of a divergence-free vector velocity potential ψ 1 . First-order fields are assumed to have a harmonic time-dependence, i.e. 1 (x, t) = Re[˜ 1 (x)e −iωt ] with angular frequency ω, where Re[ ] denotes the real part of , and˜ 1 (x) is the complex amplitude of the corresponding first-order field 1 (x, t). Using this assumption and inserting (3.4) in (3.1), (3.2), (3.3) leads to the set of first-order potential equations with the wavenumber in an unbounded viscous fluid (3.7) (Blackstock 2001, p. 305), and the viscous wavenumber where δ represents the thickness of the viscous boundary layer and is computed as We assume that there is a forward and a backward travelling wave inside the fluid, which form a pseudo-standing wave along the z-axis when superimposed. The pseudo-standing wave potentials therefore follow as a sum of potentials of forward and backward travelling waves, namely (3.11) with superscripts + and − identifying the forward and the backward travelling waves, respectively. Based on the axial symmetry of the problem, we assume the following form of potentials: with the basis vector e θ of unit length, unknown functions F(r) and G(r), and an unknown complex wavenumber k.
To find the functions F(r) and G(r), we insert the potentials of the forward travelling wave, (3.12) and (3.14), into (3.5) and (3.6), which leads to The unknown functions F(r) and G(r) are the solutions to the Sturm-Liouville equations that follow from the r-dependent part of (3.16) and (3.17), and can be expressed as with constants A 1 , A 2 , B 1 , B 2 , and nth-order Bessel functions J n and Y n of the first and of the second kind, respectively. Because the velocity has to be finite at r = 0, A 2 and B 2 have to be zero, which yields (3.23) 3.1. Dispersion relation and the unknown constants at the first order The first-order velocity field of a forward travelling wave can be written in a potential form as We substitute (3.12) and (3.14) into (3.24), which leads to the velocity field of a forward travelling wave, with basis vectors e r and e z of the cylindrical coordinate system. In order to find the dispersion relation and the unknown constants A 1 and B 1 , we apply the no-slip boundary condition to the fluid at the wall of the tube, namely where | r=a denotes the evaluation of at r = a. Substituting the solutions (3.22) and (3.23) into (3.25) and using the condition (3.26) leads to a system of two equations, (3.28) The system of (3.27) and (3.28) has a non-trivial solution only if the determinant of the system equals zero, i.e.
This leads to the dispersion relation which is used to numerically find the real and the imaginary parts of the wavenumber k, defined as k R = Re[k] and k I = Im[k], respectively. The equations leading to both parts of the wavenumber are The wavenumber is later calculated for water, oil, glycerol and air (see § 5 and appendix A).
We can now use one of the equations, e.g. (3.27), and express one of the unknown constants in terms of the second, which is a fitting parameter determined by the pressure profile of the wave. In our case, we choose A 1 as a fitting parameter, and B 1 then follows as (3.33)

First-order velocity field
The velocity field of the backward travelling wave is obtained by substituting (3.13) and (3.15) into (3.4): (3.34) The total velocity field is then obtained by summing (3.25) and (3.34): Substituting the solutions (3.22) and (3.23) for F(r) and G(r), respectively, into (3.35) yields the total first-order velocity: 3.3. Constants in terms of the pressure amplitude For easier physical interpretation of the constants A 1 and B 1 , we express them here in terms of the pressure amplitude p a at z = 0. Using the equation of state (3.3) in combination with the continuity equation (3.2) leads to the following relation: Applying the solution for the first-order velocity (3.36) to the pressure-velocity relation (3.37) yields the pressure field where the fraction inside the brackets is the amplitude of the pressure field p 1 ; that is, Now, the constant A 1 can be expressed directly from (3.39) as The remaining constant follows from combining (3.40) and (3.33): (3.41) Even though we consider only the lowest mode of propagation, the wave is not plane, but, as evident in (3.38), contains r-dependency due to the wave attenuation in the viscous boundary layer. A similar result was also obtained by Qi et al. (1995), but their study is restricted to wide channels (a δ).

Steady fluid motion at the second order
The steady second-order equations are obtained by time-averaging the equations that follow from the second-order perturbation expansion. The set of the so-called streaming equations comprises the continuity equation and the momentum equations where the body force on the right-hand side is the spatial variation of the Reynolds stress (Lighthill 1978). The time-average is denoted by with the first-order oscillation period T. The streaming velocity can be decomposed as where ϕ 2 is the scalar velocity potential and Ψ 2 is the divergence-free vector potential. Substituting (4.4) in (4.2) and applying the curl to the resulting equation leads to In the process of simplifying the terms on the right-hand side of (4.5), we will use the first-order velocity field (3.35) and the following identities: where k R and k I are respectively the real and imaginary part of the wavenumber, and * denotes the complex conjugate of . Equation (4.5) can now be reformulated into (4.8) with the functions R R (r), R I (r) resulting from the first term in the angle brackets on the right-hand side of (4.5), and the functions E R (r), E I (r) originating from the second term. We use separate functions for easier analysis of the individual contributions at a later stage, because the dropping of the second term in the angle brackets on the right-hand side of (4.5) is a very common assumption (e.g. Rayleigh 1884;Schlichting 1932;Lighthill 1978;Doinikov et al. 2017). The four functions are given explicitly in appendix C.
The geometry of the problem suggests (4.9) Substituting (4.9) into (4.8) and eliminating the basis vector yields the following expression: where we applied the operator Δ θ , defined as Based on the right-hand side of (4.10), we assume with the unknown functions S R (r) and S I (r). Substituting (4.12) into (4.10) gives two ordinary fourth-order differential equations, which we must solve in order to find the stream function Ψ 2 (r, z). In the process, it will be convenient to define the linear differential operator and use it to transform (4.13) and (4.14) into and respectively.
To solve (4.16) and (4.17), we will make use of the linear finite Hankel transform of order one (introduced by Sneddon 1946Sneddon , 1972, which is defined as for a function f (r) integrable on the interval 0 r b. The corresponding inverse transform follows as with¯ denoting the transformed expression, and provided that ξ i with i = 1, 2, 3, . . . are the positive roots of the transcendental equation The transformation of Δ 1 that is applied to a function f (r) gives (Sneddon 1946(Sneddon , 1972) Applying the finite Hankel transform defined in (4.18) to the differential (4.16), limiting the domain of transformation with the radius of the tube a, and exploiting the property specified in (4.21) yields The transformed unknown functionS R (ξ i ) can be expressed from (4.22) as (4.23) with the constants C 1 and C 2 defined as The function S R (r) can be obtained by applying the inverse transform (4.19)-(4.23), i.e. (4.26) The second unknown function, S I (r), can be obtained by applying an analogous procedure to the differential equation (4.17). The resulting expression can be written as where the constants are defined as Returning to the decomposition of the streaming velocity (4.4) and substituting it into (4.1) leads to Using the first-order velocity field (3.35), the first-order pressure field (3.38), the equation of state (3.3) and the identities with the source terms M R (r) and M I (r) given in appendix D. We assume, based on (4.33), that the solution has the following form: The unknown functions U R (r) and U I (r) can be obtained by solving two differential equations that follow from substituting (4.34) into (4.33), namely and The general solutions to the homogeneous equations related to (4.35) and (4.36), denoted by the superscript H , follow as and respectively, with the unknown constants C 5 , C 6 , C 7 , C 8 . The particular solutions to (4.35) and (4.36), denoted by the superscript P , are obtained by the method of variation of parameters (e.g. Boyce, DiPrima & Meade 2017), and follow as and respectively. The general solutions to (4.35) and (4.36) are sums of solutions to the related homogeneous equations, U H R (r) and U H I (r), and of particular solutions, U P R (r) and U P I (r). Therefore, the general solutions are given as and where we applied C 6 = 0 and C 8 = 0, since the solution has to remain finite at r = 0.

Streaming velocity
The streaming velocity v 2 from (4.4) can be expressed in terms of the functions S R (r), where the derivatives dS R (r)/dr, dS I (r)/dr, dU R (r)/dr, dU I (r)/dr follow from (4.26), (4.27), (4.41), (4.42), respectively, and are given in appendix E.
The expression (4.43) is a time-averaged Eulerian streaming velocity, which has often been the final result of classical studies (e.g. Rayleigh 1884;Schlichting 1932;Schuster & Matz 1940). However, the observable behaviour of the fluid undergoing streaming, e.g. by particle image velocimetry (Wiklund, Green & Ohlin 2012), is better represented by the time-averaged Lagrangian velocity (Westervelt 1953) or by the average mass transport velocity (Nyborg 1965). The time-averaged Lagrangian velocity is defined as given in appendix F. The average mass transport velocity follows as and is divergence-free, which follows from (4.1). Using (3.3), (3.38) and (3.36), together with (4.46), leads to 4.2. Unknown constants at the second order To find the constants C 1 , C 2 , C 3 , C 4 , C 5 and C 7 appearing in (4.26), (4.27), (4.41) and (4.42), we apply the no-slip boundary condition at the wall. The condition is imposed by fixing the time-averaged Lagrangian streaming velocity (4.44) to zero at r = a. However, the assumption of a rigid wall leads to zero first-order velocity at r = a, and consequently to v SD = 0 at r = a. Therefore, the no-slip boundary condition can be imposed by setting the Eulerian streaming velocity (4.43) to zero at r = a.
This leads to four equations, namely from the r-component of the streaming velocity, and ρ 0 η from the z-component.
Since (4.26) and (4.27) contain a factor of J 1 (ξ i r) that evaluates to zero for every ξ i at r = a, S R (r) and S I (r) are also zero at r = a by definition. The constants C 2 and C 4 then follow trivially from their definitions (4.25) and (4.29), respectively: From (4.49) and (4.51), we obtain Exploiting, again, that S R (r) and S I (r) are zero at r = a, together with boundary conditions (4.50) and (4.51), yields the last two constants in the following form:

Streaming analysis
Here, we analyse the acoustic streaming for multiple combinations of material, geometrical and wave properties in a tube excited with two opposing decaying travelling waves that form a pseudo-standing wave. The analysis is performed in , at which the streaming velocity field is sufficiently converged. The limit is determined through convergence studies, information on which is given in appendix B.
In our analysis, we use water, oil, glycerol and air. The material properties of all four fluids are given in table 1.

Definitions
To improve understanding of the fluid behaviour at the first order, we show, in figure 2, the real and the imaginary part of the wavenumber. In addition, the first-order pressure profile is given in figure 3, which indicates the symmetry of the pseudo-standing pressure wave with respect to the z = 0 plane. The amplitude decay relates to the imaginary part (k I ) of the wavenumber, while the change of the wavelength with tube radius relates to the real part (k R ) of the wavenumber. For the liquids, we will use a pressure amplitude of 100 kPa, which corresponds to the usual pressure amplitudes in acoustofluidic systems (Barnkob et al. 2010). In the case of air, the usual pressure amplitude is much lower, up to 1 kPa (Schuster & Matz 1940;Imani & Robert 2015), which is the value that will be used for the results presented in appendix A.
The first velocity node is positioned at z = 0, which is also a plane of symmetry of our pseudo-standing pressure wave (as indicated in figure 3). Since the fluid is viscous and bounded, the actual wavelength is defined as λ R = 2π/k R . The velocity node therefore repeats with a period of λ R /2. The velocity antinode is in this case positioned at z = λ R /4, and repeats with the same period as the node. Even though the fluids in our analysis are viscous (i.e. the imaginary part of the wavenumber k is nonzero), we will use the acoustic wavelength of an unbounded ideal fluid λ, defined as λ = c f /f , for a more generalized analysis. Re Re The second-order solution is also symmetric about the z = 0 plane, which follows from the analysis of the z-dependency of the streaming velocity fields (4.43) and (4.47). The streaming patterns will therefore be analysed on the interval 0 ≤ z ≤ λ. When a large-scale behaviour is analysed, the interval of analysis will be increased to 0 ≤ z ≤ 4λ.
In continuation, we will distinguish between two types of streaming vortices, namely Rayleigh streaming and Schlichting streaming. The Rayleigh streaming vortex denotes the flow from the velocity node towards the antinode near the axis of the tube, and the oppositely directed flow near the wall. The Schlichting streaming vortex represents the flow in the direction opposite to that of the Rayleigh streaming, and appears between the Rayleigh streaming and the wall. The two types of streaming are also shown in figure 4. The streaming vortices are analysed on a single rz-plane, as the problem is axisymmetric. However, when three-dimensional geometry is considered these vortices extend through the whole range of θ, namely 0 θ < 2π, and are of toroidal shape. The two types of First-order pressure along z-axis at r = a/2 Oil (Ξ = 22.5) Oil (Ξ = 10) Oil (Ξ = 7.5) Figure 3. The first-order pressure in oil along the z-axis, at r = a/2, for f = 100 kHz and p a = 100 kPa. The dimensionless ratio Ξ = a/δ represents the varying tube radius relative to the viscous boundary layer thickness of δ = 12 μm. For each radius of the tube, we plot the pressure in the range of −5λ R < z < 5λ R , with λ R = 2π/k R .  streaming vortices are boundary-driven, and are normally confined between a velocity nodal plane and antinodal plane of an acoustic field. However, we will also encounter streaming driven by the decaying travelling wave grazing the wall (Nyborg 1965, p. 274). This type of large-scale (with respect to λ) streaming will be called Eckart-like streaming, after Eckart (1948), who theoretically investigated the streaming driven by the attenuation of a travelling wave in the bulk of a fluid, extending over several acoustic wavelengths in the direction of the wave propagation. This happens, for example, when the waveguide of an acoustic beam is wider than the acoustic wavelength (Boluriaan & Morris 2003). We will often refer to the ratio of a tube radius to the thickness of the viscous boundary layer, i.e. the quantity Ξ = a δ . (5.1) The critical ratio Ξ S , as indicated in figure 5, denotes the transition between the Rayleigh-plus-Schlichting streaming regime to the Schlichting-only streaming regime. It is defined as the value of Ξ at which the axial component of the streaming velocity changes direction, at z = λ R /8. We will also use another dimensionless parameter, the ratio of the The critical ratio Λ E signifies the transition between the Rayleigh-plus-Schlichting streaming regime to the Eckart-plus-Schlichting regime, and is defined as the value of Λ at which the axial component of the streaming velocity changes direction, at z = 3λ R /8, while Ξ is sufficiently above Ξ S so that the transition associated with Ξ S does not interfere. The specified critical values of Ξ and Λ, depicted in figure 5, are problem-specific parameters, and their dependency on the material properties and on the excitation frequency will be investigated in the following sections. The superscript M indicates that the ratio is determined via the average mass transport velocity (4.47), while critical ratios without the superscript are determined via the Eulerian streaming velocity (4.43). Figure 6 shows the evolution of the average mass transport velocity patterns inside an oil-filled tube with respect to the decreasing tube radius at a constant frequency of f = 100 kHz, which corresponds to a viscous boundary layer thickness of δ = 12 μm and an acoustic wavelength of λ = 14.45 mm. We observe a qualitative change in the behaviour.  Figure 7 gives the Eulerian streaming velocity patterns corresponding to the same cases as presented in figure 6. We observe that the periodicity of the streaming pattern with respect to the z-axis, and within the first few vortices, is very distinct when only a single type of streaming is present (e.g. in figure 6a,e, f ), compared to the situations where Rayleigh and Schlichting streaming vortices coexist (e.g. in figure 6b-d). The first-order acoustic fields are periodic in z, but their amplitudes are exponential functions of z. Then, at the second order, the time-averaging of the products of first-order fields leads to the decoupling of the sinusoidal ), for f = 100 kHz and p a = 100 kPa. Also shown for comparison are the z-components of the Eulerian streaming velocity, v z 2 , the Lagrangian streaming velocity, v Lz 2 , the average mass transport velocity, v Mz 2 , and a reference velocity (evaluated at z = λ/8) from Schuster & Matz (1940), v Sz 2 , given as (G1) in appendix G. and exponential z-dependencies from products into sums. This results in the observed non-periodic behaviour along the z-axis at the second order.

The evolution of streaming patterns
The streaming patterns in figure 6 are plotted for one ideal wavelength (λ) along the z-axis. However, we can observe that the streaming patterns and the actual wavelength (λ R ) shrink along the z-axis with the decrease of Ξ . This follows from the increase of k R with the decrease of Ξ , which is shown in figure 2.
In figures 8 and 9, we look at the axial and radial velocity profiles, respectively, corresponding to the cases analysed in figures 6 and 7. The superscripts z and r denote ), for f = 100 kHz and p a = 100 kPa. Also shown for comparison are the r-components of the Eulerian streaming velocity, v r 2 , the Lagrangian streaming velocity, v Lr 2 , the average mass transport velocity, v Mr 2 , and a reference velocity from Schuster & Matz (1940), v Sr 2 , given as (G1) in appendix G. the z-and r-components of the velocity, respectively (e.g. v z 2 = v 2 · e z ). We analyse profiles of four different velocities: v z 2 and v r 2 relate to the Eulerian streaming velocity v 2 ; v Lz 2 and v Lr 2 relate to the Lagrangian streaming velocity v 2 L ; v Mz 2 and v Mr 2 relate to the average mass transport velocity v 2 M ; and we also plot v Sz 2 and v Sr 2 , which correspond to the model by Schuster & Matz (1940), given as (G1) in appendix G. The axial velocity profiles, in figure 8, are plotted through the radius of the tube, at z = λ R /8, i.e. at the middle of the first streaming vortex along the z-axis, starting from z = 0. We observe that the magnitudes of v z 2 , v Lz 2 and v Mz 2 are largest when one type of streaming vortex is 911 A28-23 dominant, i.e. at Ξ = 22.5 (a = 270 μm) and at Ξ = 3.33 (a = 40 μm). The profiles are, outside the viscous boundary layer, well aligned with v Sz 2 at Ξ = 22.5, but the difference becomes more and more significant as the ratio Ξ decreases, which is reasonable since Schuster & Matz (1940) assumed Ξ 1. After reaching a certain critical ratio Ξ (referred to later as Ξ S ), the change of the direction of the z-component of the streaming velocity at r = 0 indicates the qualitative change in the behaviour which is not captured by the model of Schuster & Matz (1940). Namely, the Schlichting streaming vortex spreads throughout the whole radius of the tube, as seen in figure 8(e, f ).
The analysis of the critical ratios Ξ M S and Ξ S at f = 100 kHz, in figures 10 and 11, shows that the ratios are invariant with respect to the fluid properties, as they are the same for water, oil and glycerol. Figures 10 and 11 show the evolution of the z-components of the Eulerian streaming velocity and average mass transport velocity, respectively, at r = 0. Both evaluations are done at z = λ R /8, (i.e. the middle of the first streaming vortex along the z-axis). The ratios Ξ M S and Ξ S are defined as the values of Ξ at which the Eulerian 911 A28-24 streaming velocity and the average mass transport velocity, respectively, in the centre of the tube change sign, as this indicates the disappearance of the Rayleigh streaming. The critical ratios resulting from the discussed analysis are Ξ M S = 6.2 and Ξ S = 5.6. In addition, figure 11 shows the values for the streaming velocity magnitude using the model of Schuster & Matz (1940). The equation for the velocity magnitude is given as (G2) in appendix G. We observe that the simplified model (Schuster & Matz 1940) matches well with our solution when the radius of the tube is large relative to the thickness of the viscous boundary layer, i.e. for Ξ 1. However, for Ξ 20 the simplified model (Schuster & Matz 1940) is not valid anymore.

Effect of the acoustic wavelength on the streaming
The analysis of the streaming patterns at higher frequencies and for the same range of Ξ reveals no significant differences in the streaming patterns, assuming Λ = λ/a 1 and consequently λ/δ 1. However, for glycerol, the tube radius corresponding to Ξ = 22.5 (the largest value of Ξ analysed in figure 6) can approach the wavelength at frequencies that are attainable in microfluidic devices. Keeping Ξ constant at 22.5 and reducing Λ through manipulation of the frequency and radius leads to the change of the z-component of the streaming velocity at z = 3λ R /8, as depicted in figure 12. This eventually leads to the change of the direction of the analysed velocity, which signals the appearance of large-scale Eckart-like streaming, at Λ = Λ E . The critical ratios based on the average mass transport velocity and the Eulerian streaming velocity are Λ M E = 3.1 and Λ E = 3.5, respectively. In addition, we observe that the streaming velocity is constant with respect to Λ, for Λ 20. Additional studies reveal that Λ M E and Λ E are not invariant with respect to Ξ , and increase as Ξ decreases.
The transition from the Rayleigh-plus-Schlichting streaming regime at Λ > Λ M E to the Eckart-plus-Schlichting regime at Λ < Λ M E is also studied through the evolution of the average mass transport velocity patterns in figure 13. The ratio Ξ is kept constant at 22.5, while the ratio Λ is varied through adjustment of the frequency and the radius of the tube. The resulting patterns are qualitatively different from the patterns in oil shown of a travelling wave which grazes the wall (Nyborg 1965, p. 274). Since the acoustic wave propagates through the whole cross-section of the tube, the described behaviour differs from the regular Eckart streaming (Eckart 1948), where the acoustic source is smaller than the tube cross-section.

Contribution of the compressibility of the streaming flow
Streaming flow is sometimes assumed incompressible in similar situations (e.g. Schuster & Matz 1940;Doinikov et al. 2017), which means that the second-order continuity equation (4.1), namely ∇ · v 2 = −1/ρ 0 ∇ · ρ 1 v 1 , is simplified to ∇ · v 2 = 0. Nevertheless, the validity of this assumption is not clear, and the effect of compressibility at the second order, as remarked by Menguy & Gilbert (2000), still poses one of the main unresolved issues in the field of acoustic streaming. Here, we use the solution for v 2 that results if the aforementioned continuity equation (4.1) is replaced by ∇ · v 2 = 0 (this renders the irrotational part of the streaming velocity, namely ∇ϕ 2 , zero, and changes the remaining constants C 1 and C 3 ).
Since the concept of the average mass transport velocity does not arise naturally in this case, i.e. it is not divergence-free by definition, we plot only the Eulerian streaming velocity patterns, in figure 15. Comparison with figures 6 and 7 shows that the incompressible Eulerian streaming velocity patterns in figure 15 better resemble Critical ratio Ξ S assuming ∇ · v 2 = 0 Water Oil Glycerol Figure 16. The z-component of the Eulerian streaming velocity, v z 2 , using the divergence-free formulation, i.e. ∇ · v 2 = 0, at r = 0 and z = λ R /8, depending on the normalized tube radius Ξ = a/δ, for f = 100 kHz and p a = 100 kPa. The critical ratio Ξ S is defined as the value of Ξ at which v z 2 changes direction. The velocity magnitudes (S. & M.) calculated with the model of Schuster & Matz (1940) are given for reference (see (G2) in appendix G). the average mass transport velocity patterns from figure 6 than the corresponding Eulerian streaming velocity patterns from figure 7. Additionally, in the transition phase, when Rayleigh streaming is gradually disappearing (figure 15b-d), the incompressible patterns preserve the periodicity with respect to the z-axis better than the patterns of the compressible solution.
With the modified solution for the incompressible streaming flow, we analyse the critical ratio Ξ S . The results, in figure 16, indicate that the compressibility slightly decreases Ξ S , as the incompressible solution yields Ξ S = 5.8, compared to the value Ξ S = 5.6 of the compressible solution (figure 11).

Contribution of individual streaming source terms
We repeat the average mass transport velocity pattern analysis from figure 6, but this time we set E R (r) = E I (r) = 0, i.e. we omit the usually neglected source term in the streaming equations (4.2) (see e.g. Rayleigh 1884;Schlichting 1932;Lighthill 1978;Doinikov et al. 2017). The simplified solution gives the results shown in figure 17. We observe that the average mass transport velocity patterns differ from those given in figure 6 for the (compressible) average mass transport velocity patterns. Furthermore, the patterns closely resemble the incompressible Eulerian streaming velocity patterns from figure 15.
In figure 18, we analogously repeat the analysis of Ξ M S from figure 11. We observe that the transition from the Rayleigh-plus-Schlichting streaming to the Schlichting streaming occurs at lower value of Ξ when the source terms E R (r) and E I (r) are neglected. Specifically, the critical normalized radius Ξ M S is reduced from 6.2 to 5.8. Interestingly, this value corresponds to the Eulerian-streaming-velocity-based Ξ S of the incompressible streaming flow solution (figure 15).

Conclusions
We derived an analytical solution for the motion of a viscous fluid at the first and at the second order in a tube of infinite length containing a pseudo-standing wave composed of two counterpropagating decaying travelling waves along the axis of the tube. The solution 911 A28-29 for evaluation of arbitrary streaming fields is available in the supplementary material. In our analysis, we observed two types of streaming vortices, namely Rayleigh streaming (fluid flow from the velocity node to the antinode near the axis of the tube) and the oppositely directed Schlichting streaming. The Rayleigh streaming is suppressed by the Schlichting streaming when the ratio of the tube radius to the thickness of the viscous boundary layer (i.e. Ξ = a/δ) is below the critical value of Ξ M S = 6.2, considering the average mass transport velocity, for Λ = λ/a 1. If the Eulerian streaming velocity is considered, the transition occurs at Ξ S = 5.6. The critical radii Ξ M S and Ξ S at which the transitions between different streaming regimes occur are invariant with respect to the material parameters, as long as we are in the limit of Λ 1. Similar behaviour was also predicted for a two-dimensional channel by Hamilton et al. (2003a), who evaluated Ξ M S at 5.7, but with the radius replaced by the half-width of their channel.
The vortices appear to be very distinct when a single vortex type, either Rayleigh or Schlichting, is dominant, e.g. figure 6(a, f ). In the intermediate region, where Ξ is decreasing towards Ξ S , we observe non-periodic behaviour along the z-axis (e.g. figure 6c,d), which can be attributed to the decaying nature of the pseudo-standing wave. The amplitude of the acoustic streaming was found to be dependent on the normalized radius Ξ . The amplitude is lowest when Ξ corresponds to the regime where neither Rayleigh nor Schlichting streaming dominates (i.e. 20 Ξ > Ξ S ). However, for Ξ 20, the streaming velocity magnitude at the centre of the tube matches relatively well to the solution of Schuster & Matz (1940).
When Ξ > Ξ S , decreasing Λ can result in a transition from Rayleigh streaming vortices to a large-scale Eckart-like streaming (figure 13). The critical ratio Λ at which this transition occurs was evaluated for the case of glycerol at a constant Ξ of 22.5, and stands at Λ M E = 3.5 for the average mass transport velocity, and at Ξ E = 3.1 for the Eulerian streaming velocity (figure 12).
The investigation of the contribution of the compressibility of the streaming flow indicated that the resulting Eulerian streaming velocity patterns (figure 15) better resemble the average mass transport velocity patterns of the compressible solution (figure 6) than the patterns of the compressible Eulerian streaming velocity ( figure 7). Compressibility appears to decrease the critical ratio Ξ S , since the incompressible solution yields Ξ S = 5.8, whereas the compressible solution yields Ξ S = 5.6. The analysis of the effects of the individual streaming source terms revealed that all terms originating from the spatial variation of the Reynolds stress have to be considered to correctly predict the streaming patterns. Neglecting the commonly-neglected the source term, namely v 1 ∇ · v 1 = 0 (for example, Rayleigh (1884), Schlichting (1932), Lighthill (1978) and Doinikov et al. (2017) neglect this term), leads to the decrease of the critical ratio Ξ M S (from 6.2 to 5.8), but still offers a good approximation. Interestingly, the average-mass-transport-velocity-based Ξ M S = 5.8 of the solution with partially neglected source terms matches the Eulerian-streaming-velocity-based Ξ S of the incompressible solution, and the corresponding streaming patterns match as well. First-order pressure along z-axis at r = a/2 Air (Ξ = 22.5) Air (Ξ = 10) Air (Ξ = 7.5) Figure 20. The first-order pressure in air along the z-axis, at r = a/2, for f = 100 kHz and p a = 1 kPa. The dimensionless ratio Ξ = a/δ represents the varying tube radius relative to the viscous boundary layer thickness of δ = 7.13 μm. For each radius of the tube, we plot the pressure in the range of −5λ R < z < 5λ R , with λ R = 2π/k R .  Figure 21. The z-components of the Eulerian streaming velocity, v z 2 , and of the average mass transport velocity, v Mz 2 , in air, at r = 0 and z = λ R /8, depending on the normalized tube radius Ξ = a/δ, for f = 100 kHz and p a = 1 kPa. The critical ratios Ξ S and Ξ M S are defined as the values of Ξ at which v z 2 and v Mz 2 , respectively, change direction. The velocity magnitude (S. & M.) calculated with the model of Schuster & Matz (1940) is given for reference (see (G2) in appendix G). and similarly for (4.27), The differentiation of (4.41) with respect to r leads to and similarly for (4.42),