## Abstract

An expansion in terms of the ratio *λ* of the characteristic crossflow velocity *U*_{∞} to jet velocity *U*_{j}, where *λ*=*U*_{∞}/*U*_{j}≪1, is used to obtain a representation of the basic three-dimensional steady flow in the nearfield of a transverse jet at large Reynolds numbers and to study its dominant instability. The inviscid vortex sheet analysis of Coelho and Hunt is extended so as to include asymptotic analysis of the viscous shear layers forming along the boundaries of the jet. These not only allow for continuity of the velocity components but also create vorticity whose advection induces an *O*(*λ*) axial flow in the direction of the jet. A uniformly valid solution is then constructed for use in a stability analysis that concentrates on the effect of crossflow upon the dominant mode of the free jet. Both the characteristic frequency and growth rate of this mode are found to increase with *λ*, in qualitative agreement with recent experimental observations.

## 1. Introduction

The transverse jet or jet in crossflow has been extensively studied (Kamotani & Greber 1972; Fric & Roshko 1994; Kelso *et al*. 1996; Smith & Mungal 1998; Megerian *et al*. 2007) owing to its widespread applications, particularly in propulsion systems (Holdeman 1993), and its superior mixing characteristics when compared with the free jet issuing into quiescent surroundings (Margason 1993; Broadwell & Breidenthal 1984; Karagozian 1986).

Transverse jets involve the complex interaction of a round free jet with a perpendicular crossflow (figure 1), generating a variety of different vortical structures, including nearfield shear layer vortices. These nearfield vortices, which arise due to shear layer instabilities in the upstream portion of the flow, undergo deformations that are thought (Kelso *et al*. 1996; Cortelezzi & Karagozian 2001) to be associated with the formation of the well-known counter-rotating vortex pair or CVP, a flow structure associated with enhanced mixing by the transverse jet in comparison with that by the free jet.

In this paper, we investigate the instability that gives rise to the nearfield vortices after developing an asymptotic solution for the three-dimensional base flow, valid for large values of the Reynolds number and small *λ*. The investigation of the instability of three-dimensional shear flows dates back to the seminal paper of Gregory *et al*. (1955) for boundary-layer flows near a surface, but overall knowledge of instabilities in three-dimensional flows is limited in scope. Such flows occur often in technology.

Blossey & Schmid (2002) were the first to perform a stability analysis of a transverse jet. They report some results of a global stability analysis that uses a time average of a direct numerical simulation of an unsteady transverse jet to determine a base flow. For the single case explored, *λ*=1/6, the stability analysis shows that growth rates of the shear layer modes are greater than those for the case of a free jet (*λ*=0). By use of modelling and asymptotic methods, we obtain here similar results for a range of *λ* using a steady base flow.

The stability problem with the original three-dimensional inviscid vortex-sheet solution of Coelho & Hunt (1989) for the base flow of a transverse jet is investigated by Alves *et al*. (2007) by means of a Fourier expansion in the circumferential direction. Growth rates are found to increase as *λ* increases, consistent with the experimental observations of Megerian *et al*. (2007), but a maximum growth rate could not be determined due to the discontinuous nature of the base flow model used. Also, numerical results in this study are restricted to the range *λ*<0.1.

In a separate paper by Alves *et al*. (2008), a local linear stability analysis is carried out by expanding in powers of *λ*. Thus, as *λ*→0, the established stability results for a free jet (as described, for instance, in the survey paper by Michalke 1984) are retrieved and serve as a basis of reference for the transverse jet. Because no exact steady solution to the Navier–Stokes equations exists for this three-dimensional flowfield, an approximate base flow for the transverse jet based on a modified version of the inviscid solution developed by Coelho & Hunt (1989) is used in the stability analysis. Alves *et al*. (2008) remove the discontinuities present in the inviscid solution due to use of a vortex sheet model by employing *ad hoc* modelling similar to that used for the free jet (Michalke 1984). The behaviour near the vortex sheets is assumed to vary as hyperbolic tangent functions. Reasonable agreement with experimental results is obtained, showing an increase in both growth rate and respective frequency as the crossflow to jet velocity ratio *λ* increases.

This paper presents a derivation of an asymptotic solution of the Navier–Stokes equations, valid for large values of the Reynolds number and small *λ*, for the transverse jet nearfield so that a more accurate description of the basic flow is available. This is done using a boundary-layer approximation for the jet shear layer and coupling the result with the inviscid solution inside and outside the jet. This procedure is applied at each order in *λ* within a perturbation expansion methodology. Furthermore, disturbance growth rates and related frequencies are computed through a local linear stability analysis that is also performed at each order in *λ*. This perturbation expansion approach allows us to study this complex flowfield in a reasonably simple manner. However, the analysis is limited to understanding the effect of the crossflow on flow structures present in the free jet.

## 2. Inviscid base flow model

The inviscid transverse jet base flow originally derived by Coelho (1988) and corrected by Alves *et al*. (2007) can be simplified and written as(2.1)where (*U*, *V*, *W*) denote velocity components in the (*x*, *r*, *θ*) directions, respectively, and at each order in *λ* are given by(2.2)

A vortex sheet exists and a discontinuity is present in the above solution at *r*=1. The flow for *r*>1 is the same as that past a circular cylinder up to *O*(*λ*), but the flow is fully three-dimensional everywhere at *O*(*λ*^{2}). The order *λ*^{2} distortion of the vortex sheet from a circular cylinder will be ignored because our analysis of the basic flow will be carried out to *O*(*λ*) only.

## 3. Viscous base flow model

We will concentrate on how the crossflow interacts with the jet downstream of the exit of the jet. Thus, we will ignore interaction between the crossflow and the plane surface bounding it by assuming the surface to be stress free, thereby eliminating the horseshoe vortices. Experimental results by Kelso & Smits (1995) indicate that the crossflow has a minor effect on the almost uniform core velocity profile of the exiting jet for values of *λ* considered here, say *λ*<0.33. It is also assumed that the crossflow does not exhibit large scale separation as would be the case if the jet were a solid cylinder instead of a fluid one. This assumption is consistent with the experimental observations of Fric & Roshko (1994).

### (a) Boundary-layer approximation

We consider the flow of an incompressible fluid whose conservation equations in cylindrical coordinates (*x*, *r* and *θ*) are given by the continuity equation and the steady Navier–Stokes equations. All spatial variables are made non-dimensional with the jet exit radius *R*_{0}, all velocity components are made non-dimensional with the jet exit velocity *U*_{j} and the relevant jet Reynolds number is given by *Re*≡*U*_{j}*R*_{0}/*ν*.

We seek a solution to the mean flow in the limit where *λ*≡*U*_{∞}/*U*_{j}≪1 by employing the expansion(3.1)

In Coelho and Hunt's analysis, *V*_{0}=*U*_{1}=*P*_{1}=0 and *P*_{0}=const. In the current analysis, *V*_{0} is the nonparallel component in the viscous shear layer at lowest order, *U*_{1} is found to be necessary to achieve a solution at *O*(*λ*) and *P*_{1} is retained momentarily but will be taken later to be zero. The governing equations and various matching conditions near *r*=1 for each velocity component are now discussed at each order in the analysis.

#### (i) Zeroth-order problem

At the lowest order, we have the well-known mixing-layer problem of classical boundary-layer theory (Schlichting 1986), assuming that the Reynolds number is large so that the boundary-layer thickness is much less than *R*_{0}. The governing equations can be obtained by substituting expansion (3.1) into the original equations and retaining terms of order *λ*^{0}. We consider axisymmetric flow at the lowest order.

We now assume a thin viscous shear layer of non-dimensional thickness *δ*≪1 and define the radial coordinate representative of the shear layer region as(3.2)

In order to be consistent with the continuity equation at *O*(*λ*^{0}), we expand as(3.3)and, to establish the boundary-layer balance in the Navier–Stokes equations at *O*(*λ*^{0}), we define(3.4)

The assumptions above lead to and . Hence, the governing equations at *O*(*λ*^{0}) become(3.5)

The boundary-layer problem at the next order in *δ*, which would yield governing equations for and , is purposefully ignored here since it turns out that these two velocity components are identically zero for the case of a planar mixing layer (Ting 1959; Klemp & Acrivos 1972).

We also note that *λ* and *δ* can be varied independently. Taking *λ*→0 for any arbitrary but small value of *δ* yields the well-known high Reynolds number free jet solution. On the other hand, taking *δ*→0 for any arbitrary but small value of *λ* yields our modified version of the inviscid solution from Coelho & Hunt (1989). Our focus is on cases where *λ*>*δ*.

#### (ii) First-order problem

Now we proceed to the derivation of the viscous shear layer equations at order *λ* by first analysing the matching conditions at *r*=1 for the velocity components based on the inviscid solutions (2.2). These conditions imply that there is no velocity boundary layer for *V*_{1}, although there is a velocity boundary layer for *W*_{1}. For this reason, we need to write the matching conditions in the shear layer coordinate and so we scale the radial and azimuthal velocities as(3.6)to obtain the matching conditions for based on the *O*(*λ*) version of the inviscid solutions (2.2)(3.7)

By substituting expansions (3.1), as well as relation (3.2), into the original governing equations and collecting the terms of order *λ*, we obtain the first-order governing equations. From the radial momentum equation we find that . If we take our viscous solution to be *U*_{1}=0 on the basis that there is no inviscid solution to match at the edges of the shear layer, we find from the streamwise momentum equation that . Because *P*_{1}=0 on an inviscid basis, we have a contradiction. We therefore take *U*_{1}≠0 to obtain, respectively, the continuity and streamwise and azimuthal momentum equations(3.8)(3.9)(3.10)where we have set *P*_{1}=0 and ignored terms of order *δ* and higher.

At *O*(*λδ*), we need to consider only the azimuthal momentum equation(3.11)in order to achieve a continuous first derivative near *r*=1.

The stability analysis described later does not require information regarding the mean flow at *O*(*λ*^{2}), and so we limit the expansion to *O*(*λ*).

### (b) Similarity solution

The lowest order problem defined by equations (3.5) has a well-known similarity solution (Schlichting 1986) involving the similarity variable *η*, where(3.12)

We will show how this similarity solution can be extended so as to describe the flow in the shear layers at *O*(*λ*). However, its use is likely to be unsuccessful at higher orders of *λ*.

#### (i) Zeroth-order problem

At this order, the similarity solution yields the following forms for *U*_{0} and :(3.13)where a prime denotes a derivative with respect to *η*.

The solution for *f*_{0} is well established (Schlichting 1986). Our solution is in excellent agreement with that of Lock (1951).

#### (ii) First-order problem

Equation (3.10) is decoupled from the other equations arising at this order and so can be solved separately. Equations (3.8) and (3.9) on the other hand are coupled and have to be solved together. The order *λ* velocity components are defined in terms of similarity functions as(3.14)which leads to rewriting equations (3.8)–(3.10) as(3.15)(3.16)whose boundary conditions, obtained from the matching conditions (3.7), are(3.17)

The procedure used to obtain the zeroth-order similarity solution, which involves obtaining asymptotic solutions for large |*η*|, marching inwards numerically and matching the solutions at *η*=0, is employed here again. We first look for asymptotic solutions for equations (3.15) and (3.16). In the limit *η*→−∞, where *f*_{0}≃*η*+*a*_{1}, we have the following asymptotes:(3.18)

In the limit *η*→+∞, where *f*_{0}≃*a*_{3}, we have the following asymptotes:(3.19)

Equation (3.16, left) is a first-order ordinary differential equation, but we established two boundary conditions for *h*_{1} in (3.17). Hence, the system seems overdetermined. The numerical procedure described and employed previously was used to solve not only for *f*_{1} but for *h*_{1} as well. Hence, both the boundary conditions were employed. However, we can only enforce continuity of *f*_{1}, and *h*_{1} at *η*=0. This means that only three out of the four constants above could be determined in the searching procedure. Because equation (3.15) is a second-order differential equation, constants *b*_{1} and *b*_{3} were maintained in our calculations. Furthermore, a few numerical experiments showed that the value obtained for *b*_{4} was independent of the value chosen for *b*_{2}. On the other hand, no solution for these constants could be obtained if we specified *b*_{4}. In order to satisfy the inviscid boundary condition inside and outside the jet, we set *b*_{2}=0 and obtained the values for the other constants from the searching procedure. The values for the unknown constants were found to be *b*_{1}=151.74, *b*_{3}=312.18 and *b*_{4}=4.5029. It seems to us that *b*_{4} is associated with the displacement of the zero streamline at order *λ*.

Finally, we consider the *O*(*δ*) azimuthal velocity component given in equation (3.11) and define(3.20)which leads to rewriting equation (3.11) as(3.21)whose boundary conditions, obtained from the matching conditions (3.7), are(3.22)

Once again we look for asymptotic solutions to be used in our numerical procedure. In the limits *η*→∓∞, the asymptotic solutions are given by(3.23)

The unknown constants in this analysis, obtained by enforcing continuity of *g*_{1,1} and at *η*=0, were found to be *b*_{5}=0.17100 and *b*_{6}=0.66945.

Figure 2 shows the first-order similarity solutions for *f*_{1}, *g*_{1,0}, *h*_{1} and *g*_{1,1} as functions of *η*. The result for *f*_{1} is of special interest because it was found to be zero in the inviscid analysis of Coelho & Hunt (1989). It arises due to radial advection of the vorticity associated with *U*_{0} and, as a result, its amplitude is independent of viscosity. As indicated in figure 2, *f*_{1} has the nature of a jet-like motion in the axial direction along the shear layer, but *U*_{1} can have either sign depending on *θ* due to the cos[*θ*] term in equation (3.14).

We note in passing that the fully numerical solution obtained by Alves (2006) contains a 2*π*-periodic component for *U*(*r*, *θ*), which corroborates the present asymptotic solution.

### (c) Uniformly valid asymptotic solution

We here use the method of matched asymptotic expansions to combine the inviscid and boundary-layer solutions presented previously into one uniformly valid viscous solution by adding both solutions and subtracting the common terms. The reader is referred to Van Dyke (1964) and Cole (1968) for more details.

The streamwise boundary-layer velocity component *U*_{0} already coincides with the inviscid solution outside of the shear layer, whereas *U*_{1} is zero there. However, that is not the case for the other velocity components for which the inviscid solution is more complicated and a correction term must be added (Alves 2006). Special care is needed when composing the uniformly valid viscous solution for *V*_{1} because we need to take into account the effect of *b*_{4}. If the procedure used for *W*_{1}(*r*) is repeated for *V*_{1}(*r*), the entire solution will be displaced outwards. This is the reason why we associated *b*_{4} with a boundary-layer displacement effect. Hence, we replace *η* by *η*−*b*_{4} in solution (3.6) in order to eliminate this problem. It is important to note that the same needs to be done with *U*_{1} since the governing equations for both of these components are coupled.

Figure 3 shows the uniformly valid solution for *W*_{1}(*r*) and *V*_{1}(*r*), together with the inviscid and boundary-layer solutions, as a function of the radial coordinate *r* for *Re*=1000 and at *x*=1. Although not very clear on the scales shown in this figure, the second-order derivative of both functions *V*_{1} and *W*_{1} is discontinuous because terms of *O*(*δ*^{2}) have been ignored in their respective expansions.

## 4. Local linear stability analysis

We consider the flow of an incompressible fluid whose conservation equations in cylindrical coordinates (*x*, *r* and *θ*) are given by the continuity equation and the unsteady Navier–Stokes equations. The next step in the analysis is to split all dependent variables into mean plus disturbance values, where the mean velocity is assumed to satisfy the governing equations in an asymptotic sense. By subtracting the mean governing equations from the original set and neglecting the nonlinear terms, we obtain the disturbance equations(4.1)(4.2)

Because the instability mechanism of interest is essentially inviscid, we consider in this work inviscid disturbances only, i.e. (*Re*→∞) in equation (4.2), although viscous effects may be of importance for low enough Reynolds numbers as shown by Morris (1976) for the free jet problem. At the large Reynolds numbers envisioned here, the viscous correction to the growth rate is very small.

As discussed in §1, the free jet is sensitive to local instabilities that are convectively unstable. Hence, we expect the same to be true for the transverse jet, at least for small values of *λ*, and introduce spatially evolving disturbances that can be separated into the Fourier modes as follows:(4.3)where c.c. denotes the complex conjugate and . Since we assumed that disturbances grow downstream, the wavenumber *α* is taken to be a complex number and the frequency *ω* a real number. Owing to the normal mode assumption, the streamwise development of our base flow is neglected and *x* becomes a parameter in our analysis, as commonly done in locally parallel linear stability analysis (Michalke 1971, 1984; Michalke & Hermann 1982).

Applying relations (4.3) to the disturbance continuity equation (4.1) yields(4.4)whereas the linearized disturbance Euler equations, with streamwise derivatives of the basic flow neglected, yield(4.5)(4.6)and(4.7)

### (a) Perturbation expansion approach

We now take advantage of the fact that our base flow was developed by assuming the crossflow to jet velocity ratio *λ* to be a small parameter and expand our dependent variables in terms of *λ* in such a way as to obtain the free jet disturbance equations in the limit *λ*→0. The disturbance quantities are then expanded as(4.8)where the azimuthal mode number at the lowest order has been set to zero (*m*=0) since previous studies have shown that the axisymmetric mode is the most unstable one in the near field of transverse jets (Megerian *et al*. 2007; Alves *et al*. 2008). Also, the axisymmetric mode is the most unstable mode for the free jet upstream of the end of the potential core.

We also expand the wavenumber *α* in the same way,(4.9)and set *ω*=ω_{0}, which is the frequency corresponding to *α*_{0}.

For our *O*(*λ*^{2}) expansion to be accurate, we need to have *λ*^{2}<0.1 or *λ*<0.33. The solution procedure, which is demonstrated in this section, was applied using the symbolic computation capabilities of the software Mathematica (Wolfram 1999).

#### (i) Zeroth-order problem

Substituting expansions (3.1), (4.8) and (4.9) into the disturbance equations (4.5)–(4.7), collecting the terms of order *λ*^{0} and neglecting *V*_{0} in analogy with the ‘quasi-parallel’ analysis of two-dimensional flows, a single equation for the zeroth-order pressure disturbance *P*_{0} can be written in the form(4.10)which is subject to the following boundary conditions:(4.11)

The solution procedure at this order is standard and consists of obtaining asymptotic solutions as *r*→0 and *r*→∞, using them to march towards an intermediate location and then matching the solutions.

#### (ii) First-order problem

Substituting expansions (3.1), (4.8) and (4.9) into the disturbance equations (4.5)–(4.7) and collecting the terms of order *λ*^{1}, we are able to obtain the set of governing equations at this order. The *θ* dependence of *U*_{1}, *V*_{1} and *W*_{1} is known from the base flow relations (3.14), allowing us to write(4.12)and so we can examine the *θ* dependence of the forcing terms on the r.h.s. of the *O*(*λ*) equations. Note that the zeroth-order disturbances do not depend on *θ*. As a result, the complex wavenumber correction *α* does not interact with the crossflow terms of our base flow or, in other words, the forcing terms involving the basic flow are non-resonant. This fact leads us to take *α*_{1}=0.

The *θ* dependence of the first-order disturbances can now be taken as(4.13)allowing the *O*(*λ*) equations to be simplified further, where forcing terms contain only resonant terms. These equations can be combined into one single equation for the first-order pressure disturbance . However, in doing so, 's second derivative in *r* would appear in one of the nonhomogeneous terms in the resulting equation. Since this derivative is discontinuous, as mentioned earlier, this procedure is avoided in order to minimize numerical error in the solution. Nevertheless, and can still be eliminated from the system of equations without having to take any second-order *r* derivatives of the first-order base flow components. In doing so, we obtain(4.14)and(4.15)where the nonhomogeneous forcing terms are given by(4.16)and(4.17)

The first-order pressure disturbance boundary conditions are(4.18)and the first-order radial velocity disturbance boundary conditions are obtained by evaluating equation (4.15) at *r*→0 and *r*→∞. Since the nonhomogeneous term *Λ*^{2} vanishes in both limits, we have in these limits(4.19)

One can easily see that the problem given by equations (4.14)–(4.19) has a non-trivial solution, which will be carried on to the second-order problem. Hence, we need to solve for the first-order disturbances. The solution procedure for these equations is very similar to the one used for the zeroth-order problem and further details can be found elsewhere (Alves 2006).

#### (iii) Second-order problem

Substituting expansions (3.1), (4.8) and (4.9) into the disturbance equations (4.5)–(4.7) and collecting the terms of order *λ*^{2}, we are able to obtain the set of governing equations at this order. Since the *θ* dependence of *U*_{2}, *V*_{2} and *W*_{2} can be inferred from the inviscid base flow relations (2.1), and the *θ* dependence of the zeroth- and first-order disturbances and respective base flows were given in previous subsections, we can examine the *θ* dependence of the forcing terms on the r.h.s. of the *O*(*λ*^{2}) equations. Similar to the first-order problem, the wavenumber correction is multiplied by terms that depend on *r* only. Hence, resonance will only exist if there are crossflow-related forcing terms in the governing equations that depend only on *r* as well. This allows us to define(4.20)in order to focus on the resonating part of the solution. Similar to the first-order problem, it is possible to combine the resulting equations into one single equation for the second-order pressure disturbance(4.21)where the nonhomogeneous forcing terms are given by(4.22)(4.23)and(4.24)

Equations (4.21) and (4.24) can be used to solve for the second-order correction for the complex wavenumber *α*_{2} as a function of the known lowest order complex wavenumber *α*_{0} and real frequency *ω*_{0}. In order to do so, we require the homogeneous solution of equation (4.21) to be orthogonal to the forcing term of the self-adjoint form of the same equation. This leads to the solvability condition(4.25)

## 5. Results and discussion

Figure 4*a*,*b*, respectively, shows the wavenumber and growth rate corrections for the free jet axisymmetric mode as a function of the frequency *ω*_{0} for different values of *x*. As happened with the linear stability analysis of the inviscid base flow (Alves *et al*. 2007) and of the continuous (tanh) base flow model (Alves *et al*. 2008), we observe divergent behaviour when analysing the results obtained in the limit of low frequency. The reader is referred to the aforementioned papers for a more detailed discussion. These corrections are different from the ones obtained with the hyperbolic tangent model (Alves *et al*. 2008) in several ways: (i) the low-frequency region (*ω*_{0}/*π*<0.4) has a more pronounced variation, (ii) the corrections decrease significantly in the high-frequency end, and (iii) these corrections are less sensitive to variations in *x*, as shown in results not displayed here.

Figure 5*a* shows the impact of the behaviour just discussed on the transverse jet growth rate as a function of the frequency *ω*_{0} for different values of *λ*. For the highest value of *λ*, figure 5*a* indicates an increase in the value of the frequency associated with the most unstable disturbance, which is consistent with the result obtained with the tanh model and experimental results. This figure shows one major difference between the two base flow models: the range of unstable frequencies for the uniformly valid asymptotic solution (UVAS) decreases as the crossflow becomes stronger.

In order to demonstrate the impact of the base flow component *U*_{1} on the instability, we show in figure 5*b* the spatial growth rate of the axisymmetric mode as a function of frequency for different velocity ratios in the absence of *U*_{1}. A significant qualitative change can be seen in this figure as compared with figure 5*a*. In figure 5*b*, the maximum growth rate decreases as does the frequency associated with it when *λ* increases. We therefore conclude that the axial flow *U*_{1}(*r*, *θ*) has a destabilizing effect, although other factors also undoubtedly contribute to the increased growth rate.

## 6. Conclusions

We have demonstrated that the vortex sheet solution of Coelho & Hunt (1989) for a steady, inviscid transverse jet can be improved upon by means of asymptotic analysis at large *Re* so as to allow for continuous variation of the velocity in shear regions forming along the boundary of the jet. It was found that the vorticity created by viscosity at *O*(*λ*) induces an axial flow at this order. Similarity solutions for flow in the shear layer were obtained to *O*(*λ*), and a UVAS was obtained by combining these solutions with the inviscid solution of Coelho & Hunt (1989). This three-dimensional base flow was then used in a linear stability analysis for a disturbance corresponding to the most unstable mode existing for a free jet at large Reynolds numbers. As *λ* increases from zero, we found that both the growth rate and characteristic frequency increase, in qualitative agreement with the experimental results of Megerian *et al*. (2007). Moreover, the range of frequency for which instability occurs is reduced so that the spectrum is more peaked about the frequency corresponding to the most unstable disturbance.

## Acknowledgments

The senior author (R.E.K.) was studying at Imperial College in 1960 on a Fulbright Fellowship when he became acquainted with Trevor Stuart who was then at the Aerodynamics Division of the NPL. He transferred to the NPL for the second year of the fellowship and returned as a post-doctoral researcher following further doctoral study at MIT. He is deeply indebted to Trevor Stuart for this stimulating introduction to the field of flow instabilities and for his support and influence since then. Both authors would like to thank Dr Wu for the time put into improving the overall presentation of the paper and his constructive criticism of our work. This research was supported by the National Science Foundation under grants CTS-0200999 and CTS-0457413. L.S. de B.A. would like to acknowledge the financial support from CAPES/Brazil and UCLA Graduate Division.

## Footnotes

One contribution of 10 to a Theme Issue ‘Theoretical fluid dynamics in the twenty-first century’.

- © 2008 The Royal Society