c Cambridge University Press 2010

499

Eﬀect of compressibility on the global stability of axisymmetric wake ﬂows P. M E L I G A1 †, D. S I P P1 1

AND

J.- M. C H O M A Z2

ONERA/DAFE, 8 rue des Vertugadins, 92190 Meudon, France

2

LadHyX, CNRS-Ecole Polytechnique, 91128 Palaiseau, France

(Received 10 August 2009; revised 11 May 2010; accepted 12 May 2010; ﬁrst published online 19 August 2010)

We study the linear dynamics of global eigenmodes in compressible axisymmetric wake ﬂows, up to the high subsonic regime. We consider both an afterbody ﬂow at zero angle of attack and a sphere, and ﬁnd that the sequence of bifurcations destabilizing the axisymmetric steady ﬂow is independent of the Mach number and reminiscent of that documented in the incompressible wake past a sphere and a disk (Natarajan & Acrivos, J. Fluid Mech., vol. 254, 1993, p. 323), hence suggesting that the onset of unsteadiness in this class of ﬂows results from a global instability. We determine the boundary separating the stable and unstable domains in the (M, Re) plane, and show that an increase in the Mach number yields a stabilization of the afterbody ﬂow, but a destabilization of the sphere ﬂow. These compressible eﬀects are further investigated by means of adjoint-based sensitivity analyses relying on the computation of gradients or sensitivity functions. Using this theoretical formalism, we show that they do not act through speciﬁc compressibility eﬀects at the disturbance level but mainly through implicit base ﬂow modiﬁcations, an eﬀect that had not been taken into consideration by previous studies based on prescribed parallel base ﬂow proﬁles. We propose a physical interpretation for the observed compressible eﬀects, based on the competition between advection and production of disturbances, and provide evidence linking the stabilizing/destabilizing eﬀect observed when varying the Mach number to a strengthening/weakening of the disturbance advection mechanism. We show, in particular, that the destabilizing eﬀect of compressibility observed in the case of the sphere results from a signiﬁcant increase of the backﬂow velocity in the whole recirculating bubble, which opposes the downstream advection of disturbances. Key words: compressible ﬂows, instability, wakes/jets

1. Introduction Axisymmetric wakes have been studied experimentally and numerically for diﬀerent geometries of revolution (Mair 1965; Achenbach 1974; Fuchs, Mercker & Michel 1979). It has generally been acknowledged that this class of ﬂows is dominated by an instability of the helical mode, resulting in the low-frequency shedding of large-scale coherent structures. In the context of afterbody ﬂows, this vortex shedding is detrimental to the engineering application, as it may increase drag and cause ﬂow-induced vibrations, resulting in ﬂuctuating dynamic loads whose magnitude can be critical during the transonic phase of ﬂight. We consider here a compressible † Email address for correspondence: [email protected]ﬂ.ch

500

P. Meliga, D. Sipp and J.-M. Chomaz

afterbody ﬂow at moderate Reynolds number and at a Mach number M = 0.5, a parameter setting which may be of practical interest for the low-density ﬂows encountered in the stratosphere by high-altitude rockets and re-entry vehicles. In this range of Reynolds numbers, the vortex-shedding activity has been linked to an instability of helical modes of azimuthal wavenumbers m = ±1 (Siegel & Fasel 2001; Seidel et al. 2008). For the higher Reynolds numbers found at lower altitudes, vortex shedding persists as a coherent large-scale phenomenon superimposed on the turbulent ﬂow ﬁeld (Achenbach 1972; Taneda 1978; Depr´es, Reijasse & Dussauge 2004), which suggests that the present results rigorously derived at moderate Reynolds numbers can carry over as a ﬁrst step towards the turbulent case. We consider two model geometries of blunt and bluﬀ bodies, namely an axisymmetric blunt-based afterbody modelling an ideal rocket shape and a sphere. An eigenvalue solution of any compressible stability problem depends on the base ﬂow proﬁles and on the Mach number, the base ﬂow itself being an implicit function of the Mach number. When the latter is varied, compressible eﬀects are thus simultaneously at work at the perturbation and the base ﬂow levels, since the Mach number is modiﬁed both in the disturbance and the base ﬂow equations. The eﬀect of compressibility on the stability of shear ﬂows has been so far addressed in the framework of the local theory of parallel ﬂows (Michalke 1971; Pavithran & Redekopp 1989, amongst others). The most widely used approximation is to prescribe analytical base ﬂow proﬁles satisfying the inviscid Navier–Stokes equations. In this case, the base ﬂow is independent of the Mach number, meaning that this approach is relevant only to assessment of the eﬀect of compressibility at the perturbation level. Another approximation is to add the Crocco–Busemann relation derived from the steady boundary-layer equations (Schlichting 1978) to include a Mach squared compressible correction in the density and temperature proﬁles (Jackson & Grosch 1990; Jendoubi & Strykowski 1994). This allows recovery of the eﬀect of compressibility on the base ﬂow thermodynamic variables. Still, it fails to consider its eﬀect on the base ﬂow velocity proﬁles, which turns out to be of crucial importance for the non-parallel wake ﬂows considered here, as will be shown in the following. The present study uses the framework of the global stability of non-parallel ﬂows, i.e. ﬂows that are inhomogeneous in both the cross-stream and the streamwise directions (Jackson 1987), for which the base ﬂow underlying the stability analysis is ﬁrst computed as a solution of the compressible Navier–Stokes equations, without any other approximation. In the incompressible regime, Natarajan & Acrivos (1993) have shown that the axisymmetric solution prevailing at low Reynolds numbers is generically unstable to a steady and a time-periodic global mode. The onset of unsteadiness in the real ﬂow ultimately results from the leading-order nonlinear interaction of these concomitantly unstable modes (Fabre, Auguste & Magnaudet 2008; Meliga, Chomaz & Sipp 2009a): the sequence of bifurcation undergone by the real ﬂow then involves the destabilization of the steady three-dimensional branch by an antisymmetric or a symmetric perturbation made of the superposition of two counter-rotating oscillating modes with diﬀerent phases. This competition between modes is not within the scope of the present study, which deals instead with the eﬀect of compressibility on the oscillating mode, believed to trigger the onset of the periodic regime at higher Reynolds numbers (Ormi`eres & Provansal 1999; Tomboulides & Orszag 2000). Increased computer capacities now make it possible to apply this approach to the compressible regime, where a number of computational issues arise (see Crouch, Garbaruk & Magidov 2007 or Crouch et al. 2009 on the shock-induced buﬀet over a two-dimensional airfoil, Br`es & Colonius 2008 on the ﬂow over an open

Eﬀect of compressibility on the global stability of axisymmetric wakes

501

l D Flow

D

Figure 1. Schematic of the conﬁgurations under study. (a) Slender body of revolution of diameter D and total length l = 9.8D. (b) Sphere of diameter D.

cavity, Robinet 2007 on shock wave/boundary-layer interactions or Mack, Schmid & Sesterhenn 2008 on the ﬂow around a swept parabolic body). We develop a gradient-based sensitivity formalism for the study of compressible non-parallel ﬂows. Compressible eﬀects are viewed as an input/output sensitivity problem relying on the evaluation of the gradient of an eigenvalue (output) with respect to small modiﬁcations of the Mach number (input). Emphasis is put on the role of the base ﬂow in the perturbation dynamics, as we generalize to compressible ﬂows the sensitivity formalism originally formulated for parallel ﬂows by Bottaro, Corbett & Luchini (2003) and Hwang & Choi (2006), and recently extended to spatially developing ﬂows by Marquet, Sipp & Jacquin (2008). By investigating how the growth rate of an unstable mode is aﬀected by changes in the shape of the base ﬂow proﬁles, this analysis is appropriate to theoretical investigation of the mechanisms leading the instability. We use here adjoint methods to compute the gradients of interests by solving only once the state and adjoint problems. As will be shown, such an approach requires a relatively ‘low’ computational cost and allows us to provide physical interpretations for the observed eﬀects by splitting a given gradient into the sum of production, streamwise advection and cross-stream advection terms. The paper is organized as follows. The ﬂow conﬁguration and numerical method are presented in § 2. The base ﬂow and disturbance equations are solved in § 3, where we investigate the impact of compressibility on the bifurcating modes. Section 4 presents a brief summary of the adjoint-based sensitivity formalism encompassing small modiﬁcations of the Mach number. The observed compressibility eﬀects are ultimately discussed in § 4.3, where we show that the interpretations admitted up to now are not relevant to the case of wake ﬂows. A physical interpretation is proposed in Section 5, in terms of a competition between the advection and production of disturbances. 2. Flow conﬁguration and theoretical formulation We investigate the stability of the axisymmetric ﬂow developing in the compressible regime past two model geometries of revolution, namely the afterbody and the sphere shown in ﬁgure 1. The afterbody models a rocket shape, with a blunt trailing edge of diameter D placed into a uniform ﬂow at zero angle of attack (Mair 1965; Weickgenannt & Monkewitz 2000), and is identical to that experimentally investigated ` (2004), with a total length l = 9.8D and an ellipsoidal by Sevilla & Mart´ınez-Bazan nose of aspect ratio 3 : 1. The problem is formulated using a standard cylindrical coordinate system (r, θ, z) of axis Γa , whose origin is taken at the centre of the base for the afterbody, and at the centre of the body for the sphere.

502

P. Meliga, D. Sipp and J.-M. Chomaz

2.1. Governing equations The ﬂuid is a non-homogeneous compressible perfect gas with constant speciﬁc heat cp , thermal conductivity κ and dynamic viscosity µ, related by a unit Prandtl number. The ﬂuid motion is described by the state vector q = (, u, Θ, p)T , where the superscript T designates the transpose, is the density, Θ is the temperature, p is the pressure and u = (u, v, w)T is the three-dimensional velocity ﬁeld with u, v and w its radial, azimuthal and streamwise components, respectively. The state vector q obeys the unsteady compressible Navier–Stokes equations, thus leading to a set of six nonlinear equations (continuity, momentum, internal energy and perfect gas) formulated in non-conservative variables as ∂t + ∇ · u + u · ∇ = 0, 1 1 ∇p − ∇ · τ (u) = 0, ∂t u + ∇u · u + 2 γM Re M2 γ τ (u) : d(u) − ∇2 Θ = 0, ∂t Θ + u · ∇Θ + (γ − 1)p∇ · u − γ (γ − 1) Re P rRe p − Θ = 0,

(2.1a) (2.1b) (2.1c) (2.1d)

with d(u) and τ (u) the strain and stress tensors deﬁned as 2 1 (∇u + ∇u T ), τ (u) = − (∇ · u )I + 2d(u). (2.2) 2 3 Note that a diﬀerent set of equations can be used, if the internal energy equation is replaced by its total energy or entropy counterpart. Equations (2.1) have been made non-dimensional using the body diameter D and the upstream quantities W∞ , ∞ , Θ∞ and p∞ as respective velocity, density, temperature and pressure scales, and the Reynolds and Mach numbers are deﬁned as d(u) =

Re =

W∞ ∞ DW∞ , , M= µ γ R g Θ∞

(2.3)

with Rg the ideal gas constant. For the afterbody ﬂow, the Reynolds number for transition to turbulence in the developing boundary layer is Re 12 000, as estimated from Weickgenannt & Monkewitz (2000). The Reynolds numbers prevailing here being such that Re < 1500, we assume that the boundary layer thus remains laminar up to the trailing edge. This assumption also holds for the sphere, for which the transition occurs at Reynolds numbers Re & 800 (Tomboulides & Orszag 2000). 2.2. Numerical method From now on, all governing equations are given as formal relations between diﬀerential operators. Equations (2.1) are thus conveniently written as B(q)∂t q + M(q, G) = 0,

(2.4)

where B and M are diﬀerential operators and G is a set of relevant control parameters made here of the Reynolds and Mach numbers only (in particular, the angle of attack remains zero throughout the study). In the following, one must distinguish between the complete form of these operators, deﬁned for the state vector q = (, u, Θ, p)T , and their reduced forms deﬁned for the state vector q = (, u, Θ)T , which can be straightforwardly deduced by replacing the pressure terms by their expressions arising from the perfect gas state equation. The complete form is more suitable to the presentation of the theoretical framework, whereas the reduced form is used in the

Eﬀect of compressibility on the global stability of axisymmetric wakes Γext

503

r∞ + ls

r∞ Γin

Γout

4. 5 2 Γw

z –∞ – ls

z–∞

–12.3

Γa 0

5. 25

z∞

z∞ + ls

Figure 2. Schematic of the computational domain Ω for the afterbody ﬂow: z−∞ , z∞ and r∞ denote the location of the physical inlet, outlet and lateral boundaries. This physical domain is padded with a sponge zone of width l s , shown as the light grey shaded area. The inner solid lines delimit regions characterized by diﬀerent vertex densities. The dark grey shaded area corresponds to the near-wake domain Ωin used to normalize the eigenmodes.

numerics as it requires smaller computational resources. To ease the reading, we omit voluntarily the diﬀerence between both forms, the choice of the relevant one being clear from the context. The complete form of all operators is detailed in Appendix D. The choice of the boundary conditions is crucial in compressible ﬂows. In order to apply appropriate far-ﬁeld conditions, the body is enclosed into two concentric cylinders deﬁned as r 6 r∞ + ls

r 6 r∞ and z−∞ 6 z 6 z∞ and z−∞ − ls 6 z 6 z∞ + ls

(inner cylinder), (outer cylinder).

(2.5a) (2.5b)

The inner enclosing cylinder corresponds to the footprint of the computational domain that would have been used for an incompressible ﬂow, whereas the outer cylinder deﬁnes the location of the inlet, outlet and external boundaries (denoted Γin , Γout and Γext , respectively) in the numerics. In the domain enclosed between the cylinders, shown as the light grey shaded area in ﬁgure 2, all ﬂuctuations are progressively damped to negligible levels through artiﬁcial dissipation, as the Reynolds number is smoothly decreased from its value deﬁned in (2.3) to the small value Res = 0.1 at the boundary of the computational domain. The purpose of such sponge regions is to minimize numerical-box-size eﬀects by gradually attenuating all vortical and acoustic ﬂuctuations before they reach the boundary of the domain (Colonius 2004). The Reynolds number in all equations should therefore be replaced by a ‘computational’ deﬁned by Re(r, z) = Re in the inner cylinder, and Reynolds number Re z) = Re + (Res − Re)ζ (z, z∞ ) Re(r, z) = Re + (Res − Re)ζ (z, z−∞ ) Re(r,

∞ , z))ζ (r, r∞ ) = Re(r ∞ , z) + (Res − Re(r Re

if

r 6 r∞

and z > z∞ ,

(2.6a)

if

r 6 r∞

and z < z−∞ ,

(2.6b)

if

r > r∞ ,

(2.6c)

where ζ is the function deﬁned by ζ (a, b) =

π |a − b| 1 1 , + tanh α tan − + π 2 2 2 ls

(2.7)

504

P. Meliga, D. Sipp and J.-M. Chomaz

along with α = 4. In addition to this artiﬁcial damping, numerical dissipation in the sponge zones is increased by progressive grid stretching. The governing equations are then solved using a uniform free-stream ﬂow condition u = (0, 0, 1)T , , Θ = 1

on Γin ∪ Γext ∪ Γout .

(2.8)

We enforce additional no-slip adiabatic wall conditions at the body wall u = 0, ∂n Θ = 0

on Γb ,

(2.9)

where ∂/∂n is the derivative normal to the surface. We use the FreeFem++ software (http://www.freefem.org) to generate a twodimensional triangulation of the azimuthal plane θ = 0 with the Delaunay–Voronoi algorithm. The mesh reﬁnement is controlled by the vertex densities imposed on both external and internal boundaries. Regions where the mesh density varies are depicted by the solid lines in ﬁgure 2. All equations are numerically solved by a ﬁnite-element method using the mesh M1 shown in Appendix A to provide the most accurate results, built with z−∞ = −100, z∞ = 300, r∞ = 25 and ls = 200 and made of 662 816 triangles. A set of equations is ﬁrst multiplied by r to avoid the singularity on the r = 0 axis. The associated variational formulation is then derived and spatially discretized onto a basis of Arnold–Brezzi–Fortin MINI elements (Matsumoto & Kawahara 2000), with four-node P1b elements for the velocity components and three-node P1 elements for the density and temperature. The sparse matrices resulting from the projection of these variational formulations onto the basis of ﬁnite elements are built with the FreeFem++ software. 3. Global stability analysis The stability analysis relies on the existence of a steady solution about which perturbations are superimposed. The total ﬂow ﬁeld is split into a steady axisymmetric base ﬂow Q = (ρ, U, 0, W, T , P )T and a three-dimensional perturbation q = (ρ , u , v , w , T , p )T of small amplitude . Unless speciﬁed otherwise, we present here results pertaining only to the afterbody conﬁguration, the results obtained for the sphere being very similar. 3.1. Base ﬂow calculations The base ﬂow Q is solution of the steady axisymmetric form of the nonlinear system (2.4), written formally as (3.1) M0 ( Q, G) = 0, where M0 is the axisymmetric form of operator M. The base ﬂow satisﬁes boundary conditions (2.8) and (2.9) along with the additional condition U = 0, ∂r W = ∂r ρ = ∂r T = 0 on Γa obtained for axisymmetric solutions from mass, momentum and internal energy conservation as r → 0. The base ﬂow is obtained using an iterative Newton method (Barkley, Gomes & Henderson 2002): starting ¯ this method involves iterations of the guess value Q¯ + δ Q, ¯ from a guess value Q, ¯ where δ Q is the solution of a simple linear problem. At each step, a matrix inversion is performed using the UMFPACK library, which consists in a sparse direct LU solver (Davis & Duﬀ 1997; Davis 2004). The process is carried out until the L2 -norm of the residual of the governing equations for Q¯ becomes smaller than 10−12 . In the low-Mach-numbers limit, the ﬂow quantities are expanded as power series in γ M 2 (Nichols, Schmid & Riley 2007) and the initial guess is obtained by continuation from the incompressible solution computed using the solver presented in Meliga

Eﬀect of compressibility on the global stability of axisymmetric wakes –0.3

505

1.1

2

r 1 0 –10

–8

–6

–4

0

–2

2

4

6

8

z Figure 3. Afterbody ﬂow: spatial distribution of the base ﬂow streamwise velocity, Re = 998.5; M = 0.5.

et al. (2009a). For Mach numbers M > 0.3, the initial guess is simply chosen as a solution of the compressible equations computed for a lower value of the Mach number. Since we do not use the governing equations under their conservative form, the numerical method cannot easily account for the presence of shock waves in the computational domain, as this would require the use of mesh reﬁnement techniques to fully resolve the viscous structure of the shock. Consequently, the local Mach √ number Ml = MU/ T must remain smaller than one everywhere in the ﬂow, and the free-stream Mach number can therefore be increased up to M ∼ 0.7 for the present computations. The dynamics is from now on exempliﬁed by setting M = 0.5. Figure 3 shows iso-contours of the base ﬂow streamwise velocity W computed in the high subsonic regime (Re = 998.5, M = 0.5). The solid line is the streamline linking the separation point to the stagnation point on the r = 0 axis, and deﬁnes the separatrix delimiting the recirculation bubble developing in the lee of the afterbody, whose length is approximately 2.5 diameters. The negative values of the streamwise velocity close to the axis reach 30 % of the free-stream velocity. 3.2. Eigenvalue calculations All perturbations are sought under the form of normal modes ˆ z)e(σ +iω)t+imθ + c.c., q = q(r,

(3.2)

where qˆ is the eigenmode, the so-called global mode, for which both the cross-stream and streamwise directions are eigendirections. The azimuthal wavenumber of the global mode is m, its growth rate and pulsation are σ and ω, respectively. Substituting (3.2) into (2.4) and retaining only terms of order yields a system of equations governing the normal mode under the form of a generalized eigenvalue problem for ˆ λ = σ + iω and q: (3.3) λB( Q)qˆ + Am ( Q, G)qˆ = 0, with Am the complex linearized evolution operator obtained from A = ∂M/∂q by replacing the θ derivatives by im, whose expression is detailed in Appendix D. The perturbation satisﬁes homogeneous boundary conditions linearized from the Navier–Stokes conditions, and the additional conditions at the axis Γa depend on the azimuthal wavenumber m. For the m = ±1 modes discussed throughout this ˆ = ρˆ = Tˆ = 0, ∂r uˆ = ∂r vˆ = 0. This eigenproblem study, we use the speciﬁc condition w is solved using the ‘implicitly restarted Arnoldi method’ of the ARPACK library (http://www.caam.rice.edu/software/ARPACK) based upon a shift and invert strategy (Ehrenstein & Gallaire 2005). To normalize the m = ±1 global modes, we impose ﬁrst the phase of the radial ˆ 1) is real positive. To normalize the velocity to be zero at r = 0 and z = 1, i.e. u(0,

506

P. Meliga, D. Sipp and J.-M. Chomaz (a)

–1.6

0.03

2

r 0 –20

–16

–12

–8

–4

0

4

8

(b)

12

16

–3.6

20 3.6

2

r 0 –20

–16

–12

–8

–4

0

4

8

12

16

20

z Figure 4. Afterbody ﬂow: spatial distribution of streamwise velocity for the leading global modes, M = 0.5. The black hue corresponds to vanishing magnitudes of perturbation. (a) Stationary mode 1 at threshold of the ﬁrst instability (Re1 = 483.5). (b) Oscillating mode 2 at the threshold of the second instability (Re2 = 998.5). Only the real part is shown.

mode amplitude, we introduce the near-wake domain Ωin deﬁned as z ∈ [−12.3, 5.25 ] z ∈ [−2.5, 5.25 ]

and r < 2 and r < 2

(afterbody), (sphere),

(3.4a) (3.4b)

shown as the dark grey shaded area in ﬁgure 2. We also use the inner product ˆ rdΩ, where aˆ and bˆ belong to Cn , dΩ is the surface element on the ˆ a · b Ω computational domain Ω, and · refers to the canonical Hermitian scalar product in Cn . The eigenmode is then normalized so that qˆ · Bqˆ rdΩ = 1. (3.5) Ωin

This normalization choice has no physical eﬀect but eases the comparison between results obtained on diﬀerent meshes when convergence tests are carried out. For incompressible ﬂows, this choice has a simple physical interpretation, as the integrand reads simply ˆ 2, qˆ · Bqˆ = u (3.6) so that condition (3.5) imposes the kinematic energy of the perturbation to be unity in Ωin . In the present compressible case, this inner product is convenient for the numerics but is not physically motivated, as the integrand does not represent any meaningful physical quantity, either the total energy or the total enthalpy of the perturbation. Still, we insist that the compressible eﬀects discussed in the following are intrinsic and do not depend on the choice of the inner product. For all values of the Mach number in the range M < 0.7 prevailing here, the sequence of bifurcations undergone by the axisymmetric solution is identical to that previously documented in the incompressible regime for spheres and disks (Natarajan & Acrivos 1993). When the Reynolds number is increased, the axisymmetric base ﬂow is ﬁrst destabilized at Re1 (M) by a stationary mode 1 (ω = 0) whose eigenvector chosen as ˆ 1 , Tˆ1 )T is real using the present normalization. For the afterbody qˆ 1 = (ρˆ 1 , uˆ 1 , iˆv1 , w ﬂow at M = 0.5, we ﬁnd the critical Reynolds number Re1 = 483.5 (Re1 = 212.5 for

Eﬀect of compressibility on the global stability of axisymmetric wakes ( a)

507

( b) 285

1100

U

U 280

Re 1000 275

S 900

S 270

0

0.2

0.4

M

0.6

0

0.2

0.4

0.6

M

Figure 5. Oscillating mode 2: boundary separating the unstable domain (shaded area labelled U) from the stable domain (area labelled S) in the (M, Re) plane. (a) Afterbody. (b) Sphere.

ˆ 1 at the sphere at M = 0.5). Figure 4(a) shows the streamwise velocity disturbances w threshold for the afterbody, which extend far downstream of the body. The velocity perturbation is negative in ﬁgure 4(a), meaning that the total ﬂow slows down in this azimuthal plane. The azimuthal wavenumber of this mode being m = 1, the streamwise velocity perturbation is opposite on the other side of the revolution axis, where the total ﬂow speeds up, hence inducing an oﬀ-axis displacement of the wake, as in the case of a sphere at zero Mach number (Johnson & Patel 1999). Owing to compressibility, similar eﬀects exist for the temperature and density perturbations, although at a lower level of magnitude. Namely we ﬁnd that the low-velocity wake region is hotter and lighter than the base ﬂow, whereas the high-velocity wake region is cooler and heavier (not shown here for conciseness). When the Reynolds number is increased further above Re1 , the axisymmetric solution is destabilized at the second threshold value Re2 (M) by a pair of oscillating m = 1 modes of frequency ω = ±ω2 , whose eigenvectors are complex conjugates. In the following, only the mode of frequency ω2 , named mode 2, will be discussed, ˆ 2 , Tˆ2 )T . For the afterbody ﬂow at its eigenvector being denoted qˆ 2 = (ρˆ 2 , uˆ 2 , iˆv2 , w M = 0.5, we ﬁnd a critical Reynolds number Re2 = 998.5 and a frequency ω2 = 0.40 corresponding to a Strouhal number St = ω2 D/(2πU∞ ) = 0.063 (resp. Re2 = 275.2, ω2 = 0.66 and St = 0.11 for the sphere at M = 0.5). Figure 4(b) shows the real parts ˆ 2r at threshold for the afterbody: it exhibits of the streamwise velocity disturbances w positive and negative velocity perturbations alternating downstream of the body, in a regular, periodic way that deﬁnes a spatial period of approximately 12 diameters. The imaginary part (not shown here) displays a similar structure, but are approximately in spatial quadrature with extrema located close to the nodes of the real parts. This mode therefore corresponds to the development of a spiral in the wake of the body, which rotates in time at the frequency ω2 . 3.3. Impact of compressibility We investigate from now on the eﬀect of the free-stream Mach number on the ﬂow stability by focusing on the oscillating mode 2, which is expected to dominate the ﬂow dynamics at suﬃciently large Reynolds numbers. The subscript 2 is therefore systematically omitted to ease the notation. We show in ﬁgure 5(a) the boundary of the stability domain for mode 2 in the (M, Re) plane: the ﬂow is unstable (resp. stable) for combinations of parameters

508

P. Meliga, D. Sipp and J.-M. Chomaz

located in the shaded region labelled U (resp. in the white region labelled S). The values for M = 0 arise from the resolution of the incompressible stability problem on a mesh made only of the inner cylinder deﬁned by (2.5a), the numerical method being derived from Meliga et al. (2009a) and Meliga, Chomaz & Sipp (2009b). Owing to the choice of the reference scales, the low-Mach-number limit agrees with the incompressible results without supplemental rescaling. Compressibility has a stabilizing eﬀect as increasing the Mach number is seen to yield a moderate increase in the critical Reynolds number, by approximately 17 % (from Re2 = 909.1 at M = 0 to 1061.1 at M = 0.7). Still, it must not be inferred from the present study that compressibility systematically acts as a stabilizing mechanism for compressible wake ﬂows. Indeed, we present in ﬁgure 5(b) the boundary of the stability domain for the sphere: interestingly, increasing the Mach number from small values now yields a destabilizing eﬀect, as the critical Reynolds number decreases by approximately 3 % (from Re2 = 280.8 at M = 0 to 273.8 at M 0.63). In contrast, the stabilizing eﬀect already documented for the afterbody ﬂow is ultimately retrieved for high subsonic Mach numbers. Although such variations may seem small, it will be shown in the following section that they are the consequence of competitive signiﬁcant eﬀects that are simultaneously at work. 3.4. Eﬀect of the baroclinic torque It has been suggested by Soteriou & Ghoniem (1995) that the diﬀerence in the stability properties of homogeneous and non-homogeneous shear layers results from the existence of a baroclinic torque. The main idea is that a baroclinic torque arising from misaligned gradients of base ﬂow density and pressure perturbations, reading 1 1 ∇pˆ × ∇ , (3.7) Γ = 2 γM ρ can act as a source term for the vorticity perturbations. For the parameter setting prevailing here, the existence of this torque has been used in a previous paper by the authors to interpret the stabilizing eﬀect of compressibility observed for model parallel wakes (Meliga, Sipp & Chomaz 2008). The main idea, originally formulated by Nichols et al. (2007), is that the baroclinic torque, which arises from the shear layer undulation, induces a further deformation which is out of phase with the total shear layer displacement, thus decreasing the temporal growth of the absolutely unstable mode. To assess the eﬀect of the baroclinic torque on the global stability of the fully non-parallel wake ﬂows considered here, we generalize the idea initially introduced by Lesshaﬀt & Huerre (2007) in the context of parallel hot jets. These authors proposed solving a modiﬁed dispersion relation in which the linearized momentum equations are artiﬁcially forced so as to cancel the baroclinic torque. In this case, the torque has only one non-trivial component owing to the parallel assumption and to the axisymmetry of the relevant disturbances. In contrast, the torque is fully three-dimensional here, owing to the non-parallelism of the ﬂow and to the non-axisymmetry of the relevant disturbances. In practice, we solve the forced compressible stability problem ˆ (λ + δΓ λ)B( Q)qˆ + Am ( Q, G)qˆ = ρ Sq,

(3.8)

where δΓ λ is the eigenvalue variation resulting from the addition of the forcing term S, appropriately deﬁned by 1 1 ˆ p ∇ . (3.9) Sqˆ = − γ M2 ρ

Eﬀect of compressibility on the global stability of axisymmetric wakes (a)

509

( b)

0.006

0.0006

0.03

0.015

0.004

0.0004

0.02

0.01

δΓ σ

||Γ ||max

0.002

0.0002

0.005

0.01

0 0

0.2

0.4

M

0.6

0 0

0.2

0.4

0.6

M

Figure 6. Eﬀect of the baroclinic torque on the oscillating mode 2 at the threshold of instability, i.e. for parameter settings along the neutral curve shown in ﬁgure 5. (a) Afterbody. (b) Sphere. The solid curve stands for the growth rate variation δΓ σ induced by an appropriate forcing deﬁned by (3.9), whose aim is to cancel the eﬀect of the baroclinic torque. The dashed curve represents the maximum magnitude of the baroclinic torque.

ˆ = − Γ , meaning that this forcing exactly It can be checked that ∇ × (Sq) counterbalances the baroclinic term when one recasts the momentum equations into their vorticity counterpart. When considering the growth rate variation δΓ σ , a positive value means that cancelling the baroclinic torque has a destabilizing eﬀect, i.e. the baroclinic eﬀects are stabilizing. On the contrary, a negative value means that cancelling the baroclinic torque has a stabilizing eﬀect, i.e. that baroclinic eﬀects are destabilizing. The solid lines in ﬁgures 6(a) and 6(b) represent the variation δΓ λ computed along the neutral curves shown in ﬁgure 5. For both conﬁgurations, we ﬁnd positive variations, meaning that the torque has a stabilizing eﬀect, a result consistent with the results arising from the local analysis. Moreover, this eﬀect increases with the Mach number, as both curves are monotonically increasing. This is due to the increasing magnitude of the torque, as depicted by the dashed lines in ﬁgure 6, which present the maximum amplitude of the baroclinic torque Γ max as a function of the Mach number. Still, such results question the usual interpretation of compressibility in terms of baroclinic eﬀects, as the eﬀect of the torque is identical for both the afterbody and the sphere, the diﬀerence being only in order of magnitude. In particular, this approach does not explain the destabilizing eﬀect observed in the case of the sphere. 4. Sensitivity analysis to a modiﬁcation of the Mach number We investigate how the stability of the global mode, taken at threshold of instability, is aﬀected by a small but ﬁnite modiﬁcation of the Mach number of magnitude δM. A given eigenvalue λ is explicitly a function of the base ﬂow Q and of the Mach number. The base ﬂow itself being an implicit function of the Mach number, the eigenvalue can be written formally as λ = λ( Q(M), M). To understand the complex compressible eﬀects discussed in § 3.3, one must therefore keep in mind that when the Mach number varies from M to M + δM, compressible eﬀects are simultaneously at work at two diﬀerent levels: (i) at the perturbation level, as a result of the change in the Mach number in the disturbance equations, (ii) at the base ﬂow level, as a result of the change in the Mach number in the base ﬂow equations. From this point of view,

510

P. Meliga, D. Sipp and J.-M. Chomaz

the baroclinic torque deﬁned in (3.7) can be written formally as Γ = Γ ( Q(M), M), meaning that the introduction of the forcing term (3.9) acts both at the disturbance and the base ﬂow levels. The eﬀect at the disturbance level corresponds to a ﬁctitious ﬂow for which the base ﬂow would be artiﬁcially frozen. Actually, it corresponds to the eﬀect investigated up to now in the framework of the local stability of parallel ﬂows. Indeed, in this approach, the base ﬂow is not a solution of the governing equations. It is prescribed under the form of analytical proﬁles independent of the Mach number and satisfying the inviscid equations: see, for instance, the model velocity proﬁles introduced by Monkewitz & Sohn (1988) to study the stability of hot jets, and recently used to assess the eﬀect of compressibility on the stability of jets and wakes (Lesshaﬀt & Huerre 2007; Meliga et al. 2008). In contrast, the eﬀect at the base ﬂow level corresponds to a ﬁctitious ﬂow for which the Mach number would remain constant in the disturbance equations whereas the base ﬂow would adapt the variations of compressibility, an eﬀect which lies outside the scope of most local analyses, for which a compressible Mach squared correction can be included only in the base density and temperature proﬁles by adding in the Crocco–Busemann relation to the governing equations (Pavithran & Redekopp 1989; Jackson & Grosch 1990). Although such an approach has proved fruitful in providing quantitative results for weakly non-parallel shear ﬂows such as mixing layers and jets, it will be shown here that the change in the base momentum proﬁles is exactly the key to interpretation of the stability of compressible, non-parallel wake ﬂows. Provided the Mach number modiﬁcation is small enough for the linear assumption to hold, it is possible to analyse separately both eﬀects, which are in the end likely to aﬀect the disturbance dynamics. The eigenvalue variation induced by a change in the Mach number δM can therefore be written as ∂λ ∂ Q ∂λ (4.1) δM + δM = δM λ + δ Q λ, δλ = ∂M ∂ Q ∂M with δM λ the variation induced at the perturbation level, and δ Q λ the variation induced by the implicit modiﬁcation of the base ﬂow. The expression of each speciﬁc variation can be derived by carrying out direct stability calculations, as developed in Appendix B. However, such an approach is extremely computationally intensive and time-consuming owing to the tremendous number of degrees of freedom involved. We rather use here a more systematic technique relying on sensitivity analyses, whose aim is to compute the gradients of the eigenvalue with respect to the Mach number and to the base ﬂow variables. This analysis relies on the computation of an adjoint global mode qˆ † = (ρˆ † , uˆ † , Tˆ † , pˆ † )T , i.e. a Lagrange multiplier for the global mode, obtained as the solution of an eigenvalue problem. Such an approach is classically used in ﬂow control and optimization problems (Luchini & Bottaro 1998; Gunzburger 1999; Corbett & Bottaro 2000, 2001). The detailed calculations are postponed to Appendix C, and we only mention here that the adjoint global mode is normalized with respect to the global mode so that qˆ † · Bqˆ r dΩ = 1. (4.2) Ω

The speciﬁc localization of the adjoint global modes resulting from the convective nonnormality of the linearized evolution operator (Chomaz 2005; Marquet et al. 2009) is illustrated for the afterbody ﬂow in ﬁgure 7, where one observes high amplitudes within the recirculating bubble and steadily decreasing amplitudes upstream of the

Eﬀect of compressibility on the global stability of axisymmetric wakes 8

–12

(a)

511

2

r 0 –20

–16

–12

–8

–4

0

4

8

12

16

–23

(b)

20 23

2

r 0 –20

–16

–12

–8

–4

0

4

8

12

16

20

z Figure 7. Same as ﬁgure 4 for the oscillating adjoint global modes at threshold of instability.

body, as classically reported for incompressible wakes (Meliga et al. 2009b). The correctness and accuracy of the adjoint method is assessed in Appendix B, where we present a comparison between the values of the variations δσ/δM arising from the adjoint formalism and that arising from direct stability calculations. 4.1. Compressible eﬀects at the disturbance level A change of the Mach number in the disturbance equations induces a change in the pressure gradient exerted on the perturbation and in the amount of energy it exchanges with the base ﬂow. In practice, when the Mach number is increased from M to M + δM, the ﬂow stability is altered as a forcing term δ fˆ proportional to the perturbation quantities occurs in the right-hand side of the disturbance equations: T M 2 ˆ ˆ 2γ (γ − 1) ˆ + τ (u) ˆ : d(u)) ˆ ∇p, δM. (4.3) (τ (U) : d(u) δ f = 0, γ M3 Re The resulting eigenvalue variation δM λ is computed here assuming that the base ﬂow is artiﬁcially kept constant. If δM is small, the forcing (4.3) acts as a weak perturbation of the linearized evolution operator under the form of momentum and internal energy sources in the perturbation equations (Giannetti & Luchini 2003, 2007). Following these authors, the variation is simply given by projection of the forcing onto the adjoint global mode: qˆ † · δ fˆ r dΩ, (4.4) δM λ = Ω

the corresponding growth rate variation being obtained by retaining only the real part of (4.4). 4.2. Compressible eﬀects at the base ﬂow level A change to the Mach number in the base ﬂow equations induces a change in the pressure gradient exerted on the base ﬂow and in the viscous dissipation of its energy. In practice, the ﬂow stability is altered, as the base ﬂow on which disturbances develop is modiﬁed. Provided δM is small, the base ﬂow solution of (3.1) for the new Mach number M + δM can indeed be approximated by Q + δ Q, where δ Q is the base ﬂow

512

P. Meliga, D. Sipp and J.-M. Chomaz

modiﬁcation computed as the solution of the linear problem T 2 M A0 ( Q, M)δ Q = 0, ∇P , 2γ (γ − 1) τ (U) : d(U) δM. γ M3 Re

(4.5)

The resulting variation δ Q λ is computed here assuming that the Mach number remains constant in the disturbance equations. To this end, we generalize now to compressible ﬂows the sensitivity analysis to base ﬂow modiﬁcations, as originally formulated for parallel ﬂows by Bottaro et al. (2003) and Hwang & Choi (2006) and recently generalized to spatially developing ﬂows by Marquet et al. (2008), whose aim is to compute the gradients of the eigenvalue with respect to the base ﬂow variables. All results are from now on discussed in terms of the physically relevant conservative variables, as the modiﬁcation of base ﬂow density also integrates eﬀects of the modiﬁcation of momentum and internal energy when using non-conservative variables. To this end, the base ﬂow is recast into conservative variables, namely Q c = (ρ, ρU, ρT , P )T . The corresponding base ﬂow modiﬁcation can be simply expressed in terms of δ Q as (4.6) δ Q c = Hδ Q, with H the matrix mapping from non-conservative onto conservative perturbation quantities. Introducing the complex ﬁelds ∇ρ λ , ∇ρU λ , ∇ρT λ and ∇P λ deﬁning the sensitivity of the eigenvalue to a modiﬁcation of the base ﬂow density, momentum, internal energy and pressure, the eigenvalue variation can be written as (∇ρ λ · δρ + ∇ρU λ · δ(ρU) + ∇ρT λ · δ(ρT ) + ∇P λ · δP ) r dΩ, (4.7) δ Qλ = Ω

the corresponding sensitivities for the growth rate being obtained by retaining only the real parts of these complex ﬁelds. The sensitivity of the eigenvalue to base ﬂow modiﬁcations is deﬁned by the ﬁeld ∇ Q c λ = (∇ρ λ, ∇ρU λ , ∇ρT λ , ∇P λ)T , so that δ Qλ = ∇ Q c λ · δ Q c r dΩ. (4.8) Ω

Using the normalization condition (4.2), we obtain † ∂ † qˆ † , λBqˆ + Am qˆ ∇ Q c λ = −H ∂Q

(

)

(4.9)

where the † superscript denotes the adjoint of the preceding operator (see Appendix C). 4.3. Application to the observed compressible eﬀects For both conﬁgurations, we have computed the variations δ Q σ/δM and δM σ/δM along the neutral curve of the oscillating modes 2 presented in ﬁgure 5. Results are shown as the dashed and dash-dotted lines in ﬁgure 8, where the overall variation δσ/δM computed from (4.1) is also reported as the solid line. Positive values (resp. negative values) mean that an increase in the Mach number induces an increase (resp. decrease) in the growth rate, and therefore correspond to a destabilizing (resp. stabilizing) eﬀect. The variation δM σ/δM is negative for both conﬁgurations, meaning that the compressible eﬀects at the disturbance level are always stabilizing. This result is consistent with the idea arising from local analyses that an increase in the Mach number prevents the upstream propagation of disturbance waves and thus yields a stabilization of shear ﬂows by promoting a convective instability

Eﬀect of compressibility on the global stability of axisymmetric wakes (a)

(b)

0 δM σ/δM

–0.02

0.08

δQ σ/δM δσ/δM

0.04

δσ/δM

513

–0.04 δQ σ/δM

–0.06

δσ/δM

0 δM σ/δM

–0.04

–0.08 –0.08 –0.10 0

0.2

0.4

M

0.6

0

0.2

0.4

0.6

M

Figure 8. Variation δσ/δM corresponding to a small modiﬁcation of the Mach number, computed along the neutral curves shown in ﬁgure 5 for the oscillating mode 2 at the threshold of instability. (a) Afterbody. (b) Sphere. The solid curve stands for the overall variation δσ/δM. The dashed and dash-dotted curves represent the variations δ Q σ/δM and δM σ/δM corresponding respectively to the base ﬂow variation ∂ Q/∂M, and to the modiﬁcation of the Mach number in the disturbance equations.

(Pavithran & Redekopp 1989). In contrast, the compressible eﬀects at the base ﬂow level are stabilizing for the afterbody (δ Q σ/δM 6 0) but destabilizing for the sphere (δ Q σ/δM > 0). The key idea here is thus that the compressible eﬀects discussed in § 3.3 are triggered by the modiﬁcation of the base ﬂow proﬁles on which disturbances develop. For the afterbody, the overall stabilizing eﬀect of compressibility is indeed due to the domination of the variation δ Q σ/δM, which represents 90 % of the total variation in the range of Mach numbers considered. Both contributions are signiﬁcant for the sphere, hence explaining the small variations of the Reynolds numbers mentioned previously on this conﬁguration. Still, the overall destabilizing eﬀect again results from the variation δ Q σ/δM being dominant for M 6 0.63 (this limit value being in excellent agreement with the one above which the eﬀect of compressibility is reversed and starts being stabilizing, as seen from ﬁgure 5b). 5. Physical interpretation 5.1. Eﬀect of the advection and production mechanisms For a base ﬂow modiﬁcation δ Q, it is possible to interpret the eigenvalue variation in terms of a competition between an advection mechanism and a production mechanism. In the local theory, this distinction has been formalized via the concepts of convective and absolute instability: the ﬂow is said to be locally convectively unstable if the advection of disturbances by the base ﬂow dominates over their production, and locally absolutely unstable when production is strong enough to sustain the ﬂush of the base ﬂow. For incompressible ﬂows, Marquet et al. (2008) have shown that it is straightforward to split the sensitivity function and to identify contributions accounting for the advection and production of disturbances. The case of compressible ﬂows is more involved, as the perturbation may exchange energy with the base ﬂow in diﬀerent ways. To identify advection and production terms, we linearize the governing equations, ﬁrst expressed using an integral formulation and conservative

514

P. Meliga, D. Sipp and J.-M. Chomaz

variables. The physical origin of all terms in (3.3) then naturally arises when turning back into non-conservative variables. For instance, the nonlinear term u · ∇u in the momentum equation (2.1b) corresponds to the advection of momentum by the ﬂow. Its linearization gives rise to two classes of terms in the linearized momentum ˆ · ∇U + ρU · ∇uˆ is an advection term corresponding to the advection equation: (a) ρU ˆ + ρ uˆ by the base ﬂow, (b) ρ uˆ · ∇U is a production of the momentum disturbance ρU term corresponding to the reciprocal advection of the base ﬂow momentum ρU by the perturbation. It is thus possible to gather all advection terms into the single advection operator Cm accounting for the advection of the perturbation (see Appendix D for a detailed expression). All other terms are production terms accounting for the production of disturbances through the advection of the base ﬂow quantities by the perturbation and through the sink/source terms of the governing equations. It is now possible to split the eigenvalue variation into δ Q λ = δ Q,A λ + δ Q,P λ, where δ Q,A λ is the variation arising from the change in the advection terms and δ Q,P λ is the variation arising from the change in the production terms. Physically, a positive value of δ Q,A λ indicates a destabilization of the eigenmode owing to a weakening of the disturbances advection. Similarly, a positive value of δ Q,P λ indicates a destabilization owing to an increase in the production of disturbances. These terms are computed respectively as δ Q,A λ = Ω

∇ Q c ,A λ · δ Q c r dΩ,

δ Q,P λ = Ω

∇ Q c ,P λ · δ Q c r dΩ,

(5.1)

where ∇ Q c ,A λ and ∇ Q c ,P λ are the advection and production sensitivity functions, computed by isolating the contribution of the advection and production terms in the sensitivity functions (4.9). We obtain simply ∇ Q c ,A λ = −H†

∂ λBqˆ + Cm qˆ ∂Q

(

)

†

qˆ † ,

∇ Q c ,P λ = ∇ Q c λ −∇ Q c ,A λ.

(5.2)

Figure 9(a) presents the values of δ Q,A σ/δM (solid line) and δ Q,P σ/δM (dash-dotted line) computed from (5.1) as functions of the Mach number, at the critical Reynolds number. The total variation δ Q σ/δM shown in ﬁgure 8 is also reported as the dashed line. The contribution of the production terms is negative, meaning that increasing the Mach number is stabilizing by weakening the production of disturbances for both the afterbody and the sphere. This contrasts with the contribution of the advection terms, which is found to be strikingly diﬀerent: it is negative for the afterbody, for which the downstream advection of disturbances is strengthened and thus yields a simultaneously stabilizing eﬀect when the Mach number increases; however, it is positive for the sphere, for which the downstream advection strongly weakens and yields a competitive destabilizing eﬀect. For both the afterbody and the sphere, the eﬀect of advection dominates over that of production. When the Mach number is varied, the change in the advection of disturbances resulting from the modiﬁcation of the base ﬂow proﬁles is therefore the leading mechanism, hence explaining why the compressible eﬀects discussed in § 3.3 diﬀer for both conﬁgurations. As an attempt to further analyse this mechanism, we have computed separately the contribution of the four components of the state vector to the variation δ Q,A σ ,

Eﬀect of compressibility on the global stability of axisymmetric wakes (a) 0

(b) δQ,P σ/δM

0.4

δQ,A σ/δM

–0.02

δQ σ/δM

515

δQ,A σ/δM

–0.04

δQ σ/δM

0

δQ σ/δM

–0.06

0.2

δQ,P σ/δM

–0.2

–0.08 –0.10

0

0.2

0.4

0.6

–0.4

0

M

0.2

0.4

0.6

M

Figure 9. Growth rate variation corresponding to a small modiﬁcation of the Mach number, computed for the oscillating mode 2 at the threshold of instability. (a) Afterbody. (b) Sphere. The dashed line represents the variation δ Q σ/δM. The solid line (resp. dash-dotted line) represents the contribution of the variation δ Q,A σ/δM owing to the modiﬁcation of the advection mechanism (resp. the variation δ Q,P σ/δM owing to the modiﬁcation of the production mechanism), so that δ Q σ/δM is the sum of these two contributions.

deﬁned as

δρ,A λ =

Ω

∇ρ,A λ · δρ r dΩ, ∇ρT ,A λ · δ(ρT ) r dΩ,

δρT ,A λ = Ω

⎫ ⎪ ∇ρU,A λ · δ(ρU) r dΩ, ⎪ ⎬

δρU,A λ = Ω δP ,A λ = ∇P ,A λ · δP r dΩ,

⎪ ⎪ ⎭

(5.3)

Ω

so that δ Q,A λ is the sum of these four contributions. Physically, δρU,A σ corresponds to the growth rate variation induced by the modiﬁcation of the base ﬂow momentum in the terms already identiﬁed as advection terms, i.e. the variation that would be computed in a ﬁctitious ﬂow for which only the momentum components would be allowed to vary, all other components being kept artiﬁcially ﬁxed. Of course, for real developing ﬂows such as those considered here, the Mach number acts by modifying all components of the base ﬂow, meaning that the modiﬁcations of density, momentum, internal energy and pressure cannot be prescribed individually but are connected to one another through relation (4.5). Such a decomposition is therefore qualitative and is used only as a means to gain insight at the mechanisms at work by estimating the importance of each individual base ﬂow component in the complex eﬀect observed. Results of decomposition (5.3) are given in table 1, where the dominant contribution is systematically displayed in a grey shaded cell, so as to ease the reading. We ﬁnd that when the Mach number is varied, the induced growth rate variation is entirely dominated by the contribution of momentum, as the density, energy and pressure modiﬁcations contribute very little to the overall variations. Such results have important physical interpretations in terms of base ﬂow calculations, as they suggest that small errors in the computation of ρ and T have a limited impact on the stability of the oscillating mode 2. They also indicate that the compressible correction of the base density and temperature proﬁles that can be included in the local approach by use of the Crocco–Busemann relation is not relevant to the case

516

P. Meliga, D. Sipp and J.-M. Chomaz

Afterbody Sphere

Total

δρ

δ(ρU)

δ(ρT )

δP

−4.4 × 10−2 3.6 × 10−1

2.4 × 10−4 −2.3 × 10−3

−4.5 × 10−2 3.6 × 10−1

−6.4 × 10−7 −7.3 × 10−4

0 0

Table 1. Oscillating mode 2 at the threshold of instability, M = 0.5: growth rate variation δ Q,A σ/δM induced by the modiﬁcation of the advection mechanism. The variations δρ,A σ/δM, δρU,A σ/δM, δρT ,A σ/δM and δP ,A σ/δM are obtained by evaluating the individual variations arising from the modiﬁcation of density, momentum, internal energy and pressure, so that the overall variation δ Q,A σ/δM is the sum of these four contributions. δ(ρU)⊥

Total Afterbody Sphere

−2

−4.5 × 10 3.6 × 10−1

δ(ρU) −2

−8.1 × 10 −1.0 × 10−1

3.7 × 10−2 4.6 × 10−1

Table 2. Oscillating mode 2 at the threshold of instability, M = 0.5: growth rate variation δσ/δM induced by the modiﬁcation of the cross-stream and streamwise momentum components.

of wake ﬂows, for which the correct compressible eﬀects are entirely triggered by the momentum components. The dominant variation δρU,A σ/δM can itself be split into its contributions related to cross-stream and streamwise momentum, as reported in table 2 using the subscripts ⊥ and , respectively. For both the afterbody and the sphere, the modiﬁcation of the cross-stream momentum component is stabilizing, whereas that of the streamwise component is destabilizing. Still, if cross-stream momentum induces variations of same order of magnitude in both cases, the streamwise component induces a variation larger by one order of magnitude for the sphere. This competition between cross-stream and streamwise advection precisely explains the opposite compressible eﬀects reported for the two geometries. The afterbody ﬂow is indeed stabilized because the modiﬁcation of the cross-stream advection of disturbances is dominant, whereas the sphere ﬂow is destabilized since the modiﬁcation of the streamwise advection dominates. 5.2. Discussion Regions in space responsible for the destabilizing weakening of the advection mechanism triggered by the streamwise momentum component may be identiﬁed by plotting the spatial distribution of the integrand ∇ρU,A σ · ∂(ρW )/∂M (ﬁgures 10a and 10b), whose integration over space yields the variation δρU,A σ/δM. The colour lookup table has been set up so that the dark grey hue indicates vanishing contributions, hence showing that the magnitude is almost zero everywhere in the ﬂow, except in the core of the recirculating bubble. The latter region exhibits a complex alternation of regions contributing either to a stabilization (negative values) or to a destabilization (positive values) of the ﬂow. The diﬀerence between the afterbody and the sphere is limited to the magnitude of the integrand, which is about 2–3 times larger for the sphere, and to the existence of a strongly stabilizing region located very close to the afterbody base. This diﬀerence does not arise from a diﬀerent level of sensitivity, as the spatial distributions of the gradient ∇ρU,A σ presented in ﬁgures 10(c) and 10(d) are remarkably similar, but from the base ﬂow modiﬁcation itself. We present in ﬁgures 10(e) and 10(f ) the spatial distribution of the modiﬁcation of streamwise

517

Eﬀect of compressibility on the global stability of axisymmetric wakes (a)

–2

1

(b)

2

–5

1

5

r 0

0 0

1

2

(c)

–50

1

3

0

1

2

50 (d) 1

–30

3 30

r 0

0 0

1

2

(e)

–0.2

1

3

0

1

2

0.2 ( f ) 1

–0.3

3 0.3

r 0

0 0

1

2

z

3

0

1

2

3

z

Figure 10. Oscillating mode 2 at the threshold of instability, M = 0.5. (a) Spatial distribution of the momentum integrand ∇ρU,A σ · ∂(ρW )/∂M for the afterbody. The integration over space of this ﬂow ﬁeld yields the variation δρU,A σ/δM involved in the stabilizing/destabilizing eﬀect of compressibility. The grey hue corresponds to vanishing magnitudes of the integrand. (b) Same as (a) for the sphere. (c) Spatial distribution of the gradient ∇ρU,A σ for the afterbody. (d) Same as (c) for the sphere. (e) Spatial distribution of the streamwise momentum variation ∂(ρW )/∂M for the afterbody. (f ) Same as (e) for the sphere.

momentum ∂(ρW )/∂M for both conﬁgurations. We ﬁnd negative values in the recirculating bubble of the afterbody ﬂow, in a region limited to the internal periphery of the separation line. An increase in the Mach number therefore induces a moderate increase in the backﬂow velocity, a result consistent with the fact that the modiﬁcation of the streamwise momentum component weakens the advection of the perturbations. Results are somewhat similar for the sphere, but the negative values found in the recirculating bubble are seen to be larger by one order of magnitude, so that the backﬂow velocity now strongly increases with the Mach number. Such a diﬀerence may be explained by the blockage eﬀect induced by both geometries, which is somewhat limited for the high-aspect ratio afterbody but large for the sphere. In return, the corresponding destabilizing eﬀect is expected to be larger for the sphere, consistent with the results discussed from table 2. A similar analysis can be carried out for the cross-stream momentum component, which shows that the stabilizing eﬀect is due to an increase in the cross-stream velocity in the recirculating bubble (not shown here for conciseness), the latter being responsible for the strengthening of the cross-stream advection mechanism. In closing this section, it should be noted that the momentum variations found along the separation lines, namely ∂(ρW )/∂M 6 0 and ∂(ρU )/∂M > 0, also mean that the recirculation bubble of the base ﬂow extends as the Mach number increases, as illustrated in ﬁgure 11 where we present the evolution of the recirculation length. This may be understood by recalling that the low-pressure levels prevailing in the recirculation bubble limit its spatial extension. When increasing the Mach number, this

518

P. Meliga, D. Sipp and J.-M. Chomaz (a)

(b)

2.7

2.9

2.7 2.5 2.5

Lr

2.3

2.3

2.1 0

0.2

0.4

M

0.6

0

0.2

0.4

0.6

M

Figure 11. Recirculation length Lr of the base ﬂow computed as a function of the Mach number, the Reynolds number being kept constant and equal to its critical value at M = 0.5: (a) afterbody, Re2 = 998.5; (b) sphere, Re2 = 275.2.

eﬀect is relaxed as 1/γ M 2 decreases, hence explaining the increase in the recirculation length. This mechanism persists at high Reynolds numbers, as Merz, Page & Przirembel (1978) earlier reported a similar increase in the recirculation length of the turbulent mean ﬂow developing past an axisymmetric afterbody. The results presented here are therefore believed to be valid also for turbulent ﬂows, at least qualitatively. The importance of the recirculation length on the stability of compressible ﬂows has been brieﬂy discussed by Bouhadji & Braza (2003), who carried out direct numerical simulations of the wake developing past a two-dimensional NACA 0012 wing at zero angle of attack and at a Reynolds number of 10 000. These authors report that the ﬂow is steady at M = 0 but that vortex shedding can be triggered by simply increasing the Mach number, an eﬀect that they have attributed to the increase of the recirculation length. Our results question this interpretation since we show that the ﬂow stability properties are diﬀerent for the afterbody and the sphere, even though the recirculation length increases in both cases. As a result, the variation of the recirculation length must not be seen as the leading mechanism involved in these complex compressible eﬀects, but rather as a systematical consequence of the base ﬂow modiﬁcation itself. In return, the ﬂow can be either stabilized or destabilized, depending on the competition between the strenghs of the stabilizing cross-stream and destabilizing streamwise advection mechanisms, as sketched in ﬁgure 12. 6. Conclusion A theoretical framework for the study of global modes in compressible ﬂows has been developed and applied to strongly non-parallel axisymmetric wake ﬂows in the high subsonic regime. The base ﬂow, stability and adjoint stability equations have been derived and numerically solved for an axisymmetric blunt-based afterbody modelling an ideal rocket shape and a sphere. A consistent sequence has been found for the destabilization of the steady, axisymmetric solution, that does not depend on the value of the Mach number. A ﬁrst instability occurs for a stationary global mode (named mode 1) of azimuthal wavenumber m = 1, and a second instability occurs at a larger Reynolds number, for an oscillating global mode (named mode 2), of the same azimuthal wavenumber m = 1. The similarity between this sequence of destabilization

Eﬀect of compressibility on the global stability of axisymmetric wakes (a)

519

(b)

Figure 12. Eﬀect of the Mach number on the advection of disturbances for the (a) afterbody and (b) sphere. The solid and dashed lines delimit the recirculation bubbles for Mach numbers M and M + δM, respectively. The arrows indicate the direction in which advection of disturbances is strengthened, the dominant contribution being shown as the black arrow.

and the one known from the incompressible ﬂow past a disk and a sphere gives credence to the interpretation of the large-scale oscillation observed in this class of ﬂows at large Reynolds numbers in terms of a global instability triggered by the destabilization of the oscillating mode 2. The boundaries separating the stable and unstable domains in the (M, Re) plane have been determined, and it has been shown that increasing the Mach number has a stabilizing eﬀect on mode 2 in the case of the blunt-based afterbody. As an attempt to generalize such results to other conﬁgurations of axisymmetric wakes, we have computed the same stability boundary for the sphere conﬁguration, and have shown that a similar increase in the Mach number now yields a destabilization of mode 2. This complex stabilizing/destabilizing eﬀect has been further investigated using adjoint-based sensitivity analyses, aimed at predicting the variations of the eigenvalue with the Mach number. The key point here is that compressible eﬀects are simultaneously at work at the base ﬂow and at the perturbation level. The eigenvalue variation arising from both contributions has been derived analytically, and we have shown that the observed compressible eﬀects are triggered by a modiﬁcation of the base ﬂow proﬁles. Using this adjoint-based theoretical formalism, a physical interpretation has been proposed, based on a competition between the terms of production, cross-stream advection and streamwise advection. The proposed decomposition provides evidence linking the dominant compressible eﬀect to the modiﬁcation of the streamwise momentum proﬁles. If blockage eﬀects are large, an increase in the Mach number strongly strengthens the backﬂow velocity in the recirculating bubble and makes it able to oppose the downstream advection of disturbances, hence explaining the overall destabilizing eﬀect found for the sphere. The authors acknowledge ﬁnancial support of CNES (the French Space Agency) within the framework of the research and technology programme Aerodynamics of Nozzles and Afterbodies. Appendix A. Sensitivity results to mesh spacing Seven diﬀerent meshes, denoted M1 to M7 have been used to assess convergence in the numerical results. These meshes, detailed in table 3, exhibit various spatial extents and vertex densities, as well as various sizes for the sponge zones. Results presented in this appendix correspond to the ﬁnest mesh M1 . A comparison of the results obtained with the meshes M1 to M7 is provided in table 4 and shows that a good convergence is already achieved for the coarser mesh M6 , as results are identical down to the third digit.

520

P. Meliga, D. Sipp and J.-M. Chomaz

M1 M2 M3 M4 M5 M6 M7

z−∞

z∞

r∞

ls

−100 −100 −70 −100 −100 −100 −100

200 150 200 200 200 200 200

25 25 25 20 25 25 25

100 100 100 100 70 100 100

nt

DoF0

DoFm

662 816 631 559 653 569 642 370 644 067 441 227 662 816

2 664 620 2 539 022 2 627 322 2 582 908 2 589 454 1 774 726 2 664 620

3 662 183 3 489 557 3 610 937 3 549 820 3 558 851 2 439 021 3 662 183

Table 3. Properties of the meshes as a function of the parameters z−∞ , z∞ , r∞ and l s , corresponding to the location of the physical inlet, outlet and lateral boundaries, and to the size of the sponge zone: nt is the number of triangles, DoF0 is the number of degrees of freedom for axisymmetric state vectors used in the base ﬂow calculations, and DoFm is the number of degrees of freedom for three-dimensional state vectors used in the global stability calculations. Meshes M1 and M2 have the same vertex densities but diﬀer in the location of the outlet boundary. In the same way, M1 and M3 diﬀer in the location of the inlet boundary, while M1 and M4 diﬀer in the location of the lateral boundary and M1 and M5 diﬀer in the size of the sponge zone. M1 and M6 have the same spatial extent but M6 is built with lower vertex densities. M1 and M7 are identical but we use a diﬀerent damping function in the sponge zone, deﬁned by (2.7) along with α = 3. σ1 M1 M2 M3 M4 M5 M6 M7

7.0 × 10−5 −7.7 × 10−5 −1.1 × 10−4 −3.9 × 10−5 −6.2 × 10−5 4.6 × 10−5 7.0 × 10−5

σ2

ω2

1.0 × 10−5 −1.4 × 10−4 −1.9 × 10−4 −1.4 × 10−4 −1.0 × 10−4 2.8 × 10−5 1.0 × 10−5

0.3936 0.3931 0.3930 0.3931 0.3931 0.3935 0.3936

Table 4. Dependence of the eigenvalues on the diﬀerent meshes characterized in table A. The eigenvalue σ1 corresponding to the stationary mode 1 is computed at the ﬁrst instability threshold (ReA = 212.6 − M = 0.1), and the eigenvalue σ2 + iω2 corresponding to the oscillating mode 2 is computed at the second instability threshold (Re2 = 280.6 − M = 0.1).

Appendix B. Validation of the adjoint-based gradients This appendix aims to assess the accuracy of the adjoint method presented in this study. We recall that the variation of a given eigenvalue δλ is expressed as δλ = δ Q λ + δM λ where δ Q λ and δM λ are respectively the variation arising from the modiﬁcation of the base ﬂow (the Mach number in the perturbation equations being kept constant) and that arising from the modiﬁcation of the Mach number in the perturbation equations (the base ﬂow being kept constant). Only the afterbody conﬁguration is considered here. Since the sensitivity analysis is linear in essence, the variation of the eigenvalue computed using the adjointbased gradients should agree with direct eigenvalue calculation carried out in the limit δM → 0. Considering the oscillating global mode 2 at the threshold of the second instability, i.e. Re2 = 998.5 and M = 0.5, we compute the linear estimation of the growth rate for diﬀerent values of δM, thanks to expressions (4.4)–(4.8). These variations are then computed exactly by carrying out the following eigenvalue calculations: assume Q 1 is the base ﬂow solution at the Mach number M. For each value of δM, we compute ﬁrst the base ﬂow Q 2 , solution of the nonlinear equations

Eﬀect of compressibility on the global stability of axisymmetric wakes

521

0.001 δM σ

0

δσ –0.001

δQ σ

δσ

–0.002

–0.003 10–4

10–3

δM

10–2

10–1

Figure 13. Oscillating mode 2 at the threshold of instability, Re2 = 998.5 − M = 0.5: variation of the growth rate as a function of the modiﬁcation of the Mach number δM. The solid, dashed and dash-dotted lines stand for the variations δσ , δ Q σ and δM σ obtained from the sensitivity analysis. The dark grey, light grey and white circle symbols stand for the nonlinear results obtained from direct stability calculations.

(3.1) for the Mach number M + δM. We solve then numerically the stability problems δλ −→ (λ + δλ)B( Q 2 )qˆ + Am ( Q 2 , M + δM)qˆ = 0, δ Q λ −→ (λ + δ Q λ)B( Q 2 )qˆ + Am ( Q 2 , M)qˆ = 0, δM λ −→ (λ + δM λ)B( Q 1 )qˆ + Am ( Q 1 , M + δM)qˆ = 0,

(B 1a) (B 1b) (B 1c)

the associated growth rate variations being obtained simply by retaining the real parts. Figure 13(a) depicts the growth rate variations computed as functions of the amplitude δM. The dark grey symbols (resp. light grey and white circle symbols) stand for the exact nonlinear variation δσ (resp. δQ σ and δM σ ) obtained by direct stability calculations. The corresponding linear estimations arising from the sensitivity analysis are presented as the solid, dashed and dash-dotted curves, respectively. For small amplitudes δM < 10−3 , the relative diﬀerence is not measurable and results are superposed, indicating that the base ﬂow modiﬁcation owing to the increase in M is linear in this range. These results validate the sensitivity analysis and in particular the accuracy of the sensitivity functions computed in this study. For larger amplitudes, we observe small discrepancies, as the decrease in the growth rate is slightly larger if computed by stability calculations. This means that the true nonlinear stabilizing eﬀect of the Mach number is slightly larger than that estimated by the sensitivity analysis. Still, the variations obtained up to δM = 0.05 are very well approximated by the linear estimation, as the maximum relative diﬀerence is about 6 %. Appendix C. Derivation of the sensitivity functions to base ﬂow modiﬁcations The eigenvalue variation δ Q λ, simply noted δλ here to ease the notation, is investigated with respect to the base ﬂow modiﬁcation δ Q, the Mach number being kept constant. The variations are such that δλ = δσ + iδω = ∇ Q λ , δ Q, ˆ where we use from now on aˆ , b = aˆ · bˆ r dΩ for compact notation. Ω

(C 1)

522

P. Meliga, D. Sipp and J.-M. Chomaz

ˆ λ} In the present formalism, the base ﬂow Q is the control variable, the eigenpair {q, is the state variable and eigenproblem (3.3) is the state equation, i.e. the constraint to be satisﬁed. We introduce a Lagrange multiplier qˆ † (also known as ‘adjoint’ or ‘co-state’ variable) for the state variable, now referred to as the adjoint perturbation, and deﬁne the functional ˆ λ) = λ − qˆ † , λB( Q)qˆ + Am ( Q)q. ˆ L( Q, qˆ † , q,

(C 2)

The gradient, with respect to any variable s, is deﬁned as L(s + δs) − L(s) ∂L δs = lim . (C 3) →0 ∂s We assume that the state equation is satisﬁed for any arbitrary base ﬂow modiﬁcation, so that the gradient of the functional with respect to the adjoint variable is zero. It can be checked that the gradient with respect to the state variable is zero, provided we deﬁne qˆ † as the solution of the adjoint eigenvalue problem λ∗ B† ( Q) qˆ † + A†m ( Q) qˆ † = 0,

(C 4) †

A†m

along with the normalization condition (4.2). In (C 4), B and are the adjoint of operators B and Am , obtained by integrating by parts the disturbance equations (Schmid & Henningson 2001). We obtain simply B† = B since B is a diagonal operator, whereas the complete expression of operator A†m can be found in Appendix D. The boundary conditions to be fulﬁlled by the adjoint perturbations are such that all boundary terms arising during the integration by part vanish, which imposes conditions identical to that of the global modes. Eigenproblem (C 4) is solved via the Arnoldi method, and adjoint global modes are then normalized with respect to the direct global modes, according to (4.2). Since this adjoint problem is formulated for continuous operators, the spatial discretization of operators Am and A†m leads to discrete operators that are not Hermitian conjugates. Consequently, we check a posteriori that the adjoint eigenvalues are complex conjugate with the direct eigenvalues and that a bi-orthogonality relation is satisﬁed for the 10 leading global modes (i.e. that the scalar product of one of the 10 leading adjoint modes with any of the 10 leading direct global modes is less than 10−6 , except when the direct and adjoint modes correspond to complex conjugate eigenvalues), and conclude that our numerical procedure accurately estimates the direct and adjoint global modes of the compressible problem. The eigenvalue variation now reads δλ =

∂L δ Q. ∂Q

(C 5)

The gradient of the functional with respect to the base ﬂow can be expressed as

∂L ∂ ˆ Q δ Q = −qˆ † , (λB( Q)qˆ + Am ( Q)q)δ ∂Q ∂Q † ∂ † ˆ (λB( Q)qˆ + Am ( Q)q) qˆ , δ Q , (C 6) =− ∂Q so that the sensitivity function ∇ Q λ is given by † ∂ ˆ qˆ † . ∇ Qλ = − (λB( Q)qˆ + Am ( Q)q) ∂Q

(C 7)

Eﬀect of compressibility on the global stability of axisymmetric wakes

523

Because we use non-conservative variables in the numerics, expression (C 7) corresponds to the sensitivity function ∇ Q λ = (∇ρ λ, ∇U λ , ∇T λ , ∇P λ)T , where ∇ρ λ, ∇U λ, ∇T λ and ∇P λ deﬁne the sensitivity of the eigenvalue to a small modiﬁcation of the base ﬂow density, velocity, temperature and pressure, such that δλ = (∇ρ λ · δρ + ∇U λ · δU) + ∇T λ · δT + ∇P λ · δP ) r dΩ. (C 8) Ω

To derive the sensitivity functions in term of the conservative variables, as deﬁned by (4.9), we simply substitute δ Q by its conservative counterpart Hδ Q into (C 6), since both relations (4.7) and (C 8) are to be simultaneously satisﬁed. In closing this section, it should be noted that such an approach is very similar to that used in optimization problems, where one enforces the stationarity of a Lagrangian as a means to minimize a given functional under speciﬁc constraint. We would like to insist that no such stationarity is enforced here, and that the functional is only used as a means to compute the diﬀerent gradients of interest.

Appendix D. Detailed expression of the diﬀerential operators All operators given here pertain to the complete state vector q = (ρ, u, T , p)T . The reduced form of these operators, i.e. that pertaining to the state vector q = (ρ, u, T )T used in the numerics, can be straightforwardly obtained by replacing the pressure terms by their expressions arising from the perfect gas state equation. The solid symbols • are used to clarify the action of a selected number of diﬀerential operators. I being the identity operator, the non-zero terms of operators B, Am and Cm describing the evolution of the global modes are B11 = 1, B22 = ρI, B33 = ρ, Am11 = U · ∇ + ∇ · U, Am12 = ∇ρ· + ρ∇, Am21 = ∇U · U, Am22 = ρ∇[ • ] · U + ρ∇U · [ • ] −

1 ∇ · τ [ • ], Re

1 ∇, γ M2 = U · ∇T ,

Am24 = Am31

Am32 = ρ∇T · + (γ − 1)P ∇ · − γ (γ − 1) γ ∇2 , P rRe = (γ − 1)∇ · U, = −T , = −ρ, = 1,

Am33 = ρU · ∇ − Am34 Am41 Am43 Am44

M2 (τ (U) : d [ • ] + τ [ • ]: d(U)) , Re

524

P. Meliga, D. Sipp and J.-M. Chomaz Cm11 Cm21 Cm22 Cm31 Cm33

= U · ∇ + ∇ · U, = ∇U · U, = ρ∇[ • ] · U, = U · ∇T , = ρU · ∇.

Similarly, the non-zero terms of the adjoint operators A†m and C†m are A†m11 = −U · ∇, A†m12 = (∇U · U)·, A†m13 = U · ∇T , A†m14 = −T , A†m21 = −ρ∇, 1 ∇ · τ [ • ], Re M2 ∇ · ([ • ]τ (U)) , = ρ∇T − (γ − 1)∇(P [ • ]) + 2γ (γ − 1) Re γ ∇2 , = −ρU · ∇ − P rRe = −ρ, 1 =− ∇·, γ M2 = (γ − 1)∇ · U,

A†m22 = −ρ∇[ • ] · U + ρ∇U T · − A†m23 A†m33 A†m34 A†m42 A†m43

A†m44 = 1 , C†m11 = −U · ∇, C†m12 = (∇U · U)·, C†m13 = U · ∇T , C†m22 = −ρ∇[ • ] · U, C†m33 = −ρU · ∇. REFERENCES Achenbach, E. 1972 Experiments on the ﬂow past spheres at very high Reynolds numbers. J. Fluid Mech. 54, 565–575. Achenbach, E. 1974 Vortex shedding from spheres. J. Fluid Mech. 62, 209–221. Barkley, D., Gomes, M. G. M. & Henderson, R. D. 2002 Three-dimensional instability in ﬂow over a backward-facing step. J. Fluid Mech. 473, 167–190. Bottaro, A., Corbett, P. & Luchini, P. 2003 The eﬀect of base ﬂow variation on ﬂow stability. J. Fluid Mech. 476, 293–302. Bouhadji, A. & Braza, M. 2003 Physical analysis by numerical simulation of organised modes and shock-vortex interaction in transonic ﬂows around an aerofoil. Part 1. Mach number eﬀect. J. Comput. Fluids 32, 1233–1260. Bre` s, G. A. & Colonius, T. 2008 Three-dimensional instabilities in compressible ﬂow over open cavities. J. Fluid Mech. 599, 309–339. Chomaz, J.-M. 2005 Global instabilities in spatially developing ﬂows: non-normality and nonlinearity. Annu. Rev. Fluid. Mech. 37, 357–392.

Eﬀect of compressibility on the global stability of axisymmetric wakes

525

Colonius, T. 2004 Modeling artiﬁcial boundary conditions for compressible ﬂow. Annu. Rev. Fluid Mech. 36, 315–345. Corbett, P. & Bottaro, A. 2000 Optimal perturbations for boundary layers subject to streamwise pressure gradient. Phys. Fluids 12, 120–130. Corbett, P. & Bottaro, A. 2001 Optimal control of nonmodal disturbances in boundary layers. Theor. Comput. Fluid Dyn. 15, 65–81. Crouch, J. D., Garbaruk, A. & Magidov, D. 2007 Predicting the onset of ﬂow unsteadiness based on global instability. J. Comput. Phys. 224, 924–940. Crouch, J. D., Garbaruk, A., Magidov, D. & Travin, A. 2009 Origin of transonic buﬀet on aerofoils. J. Fluid Mech. 628, 357–369. Davis, T. A. 2004 A column pre-ordering strategy for the unsymmetric-pattern multifrontal method. ACM Trans. Math. Softw. 30 (2), 165–195. Davis, T. A. & Duff, I. S. 1997 An unsymmetric-pattern multifrontal method for sparse LU factorization. SIAM J. Matrix Anal. Appl. 18 (1), 140–158. Depre´ s, D., Reijasse, P. & Dussauge, J.-P. 2004 Analysis of unsteadiness in afterbody transonic ﬂows. AIAA J. 42 (12), 2541–2550. Ehrenstein, U. & Gallaire, F. 2005 On two-dimensional temporal modes in spatially evolving open ﬂows: the ﬂat-plate boundary layer. J. Fluid Mech. 536, 209–218. Fabre, D., Auguste, F. & Magnaudet, J. 2008 Bifurcations and symmetry breaking in the wake of axisymmetric bodies. Phys. Fluids 20 (5), 051702. Fuchs, H. V., Mercker, E. & Michel, U. 1979 Large-scale coherent structures in the wake of axisymmetric bodies. J. Fluid Mech. 93, 185–207. Giannetti, F. & Luchini, P. 2003 Receptivity of the circular cylinders ﬁrst instability. In Proceedings of 5th European Fluid Mechanics Conference, Toulouse, France. Giannetti, F. & Luchini, P. 2007 Structural sensitivity of the ﬁrst instability of the cylinder wake. J. Fluid Mech. 581, 167–197. Gunzburger, M. D. 1999 Sensitivities, adjoints and ﬂow optimization. Intl J. Numer. Math. Fluids 31 (1), 53–78. Hwang, Y. & Choi, H. 2006 Control of absolute instability by basic-ﬂow modiﬁcation in parallel wake at low Reynolds number. J. Fluid Mech. 560, 465–475. Jackson, C. P. 1987 A ﬁnite-element study of the onset of vortex shedding in ﬂow past variously shaped bodies. J. Fluid Mech. 182, 23–45. Jackson, T. L. & Grosch, C. E. 1990 Absolute/convective instabilities and the convective Mach number in a compressible mixing layer. Phys. Fluids A 2, 949–954. Jendoubi, S. & Strykowski, P. J. 1994 Absolute and convective instability of axisymmetric jets with external ﬂow. Phys. Fluids 6 (9), 3000–3009. Johnson, T. A. & Patel, V. C. 1999 Flow past a sphere up to a Reynolds number of 300. J. Fluid Mech. 378, 19–70. Lesshafft, L. & Huerre, P. 2007 Linear impulse response in hot round jets. Phys. Fluids 19 (2), 024102. Luchini, P. & Bottaro, A. 1998 G¨ ortler vortices: a backward-in-time approach to the receptivity problem. J. Fluid Mech. 363, 1–23. Mack, C. J., Schmid, P. J. & Sesterhenn, J. L. 2008 Global stability of swept ﬂow around a parabolic body: connecting attachment-line and crossﬂow modes. J. Fluid Mech. 611, 205–214. Mair, W. A. 1965 The eﬀect of a rear-mounted disc on the drag of a blunt-based body of revolution. Aeronaut. Q. 16, 350–360. Marquet, O., Lombardi, M., Chomaz, J.-M., Sipp, D. & Jacquin, L. 2009 Direct and adjoint global modes of a recirculation bubble: lift-up and convective non-normalities. J. Fluid Mech. 622, 1–21. Marquet, O., Sipp, D. & Jacquin, L. 2008 Sensitivity analysis and passive control of the cylinder ﬂow. J. Fluid Mech. 615, 221–252. Matsumoto, J. & Kawahara, M. 2000 Stable shape identiﬁcation for ﬂuid–structure interaction problem using MINI element. J. Appl. Mech. 3, 263–274. Meliga, P., Chomaz, J.-M. & Sipp, D. 2009a Global mode interaction and pattern selection in the wake of a disk: a weakly nonlinear expansion. J. Fluid. Mech. 633, 159–189.

526

P. Meliga, D. Sipp and J.-M. Chomaz

Meliga, P., Chomaz, J.-M. & Sipp, D. 2009b Unsteadiness in the wake of disks and spheres: instability, receptivity and control using direct and adjoint global stability analyses. J. Fluid. Struc. 25, 601–616. Meliga, P., Sipp, D. & Chomaz, J.-M. 2008 Absolute instability in axisymmetric wakes: compressible and density variation eﬀects. J. Fluid Mech. 600, 373–401. Merz, R. A., Page, R. H. & Przirembel, C. E. G. 1978 Subsonic axisymmetric near-wake studies. AIAA J. 16, 656–662. Michalke, A. 1971 Instabilit¨ at eines kompressiblen runden Freistrahls unter Ber¨ ucksichtigung des Einﬂusses der Strahlgrenzschichtdicke. Z. Plugwiss. (19), 319–328 (English translation in NASA TM 75190, 1977). Monkewitz, P. A. & Sohn, K. D. 1988 Absolute instability in hot jets. AIAA J. 26 (8), 911–916. Natarajan, R. & Acrivos, A. 1993 The instability of the steady ﬂow past spheres and disks. J. Fluid Mech. 254, 323–344. Nichols, J. W., Schmid, P. J. & Riley, J. J. 2007 Self-sustained oscillations in variable-density round jets. J. Fluid Mech. 582, 341–376. Ormie` res, D. & Provansal, M. 1999 Transition to turbulence in the wake of a sphere. Phys. Rev. Lett. 83, 80–83. Pavithran, S. & Redekopp, L. G. 1989 The absolute-convective transition in subsonic mixing layers. Phys. Fluids A 1 (10), 1736–1739. Robinet, J.-C. 2007 Bifurcations in shock-wave/laminar-boundary layer interaction: global instability approach. J. Fluid Mech. 579, 85–112. Schlichting, H. 1978 Boundary Layer Theory, 7th edn. McGraw-Hill. Schmid, P. J. & Henningson, D. S 2001 Stability and Transition in Shear Flows. Springer. Seidel, J., Siegel, S., Jeans, T., Aradag, S., Cohen, K. & McLaughlin, T. 2008 Analysis of an axisymmetric bluﬀ body wake using Fourier transform and POD. AIAA Paper 2008-552. ` C. 2004 Vortex shedding in high Reynolds number axisymmetric Sevilla, A. & Mart´ınez-Bazan, bluﬀ-body wakes: local linear instability and global bleed control. Phys. Fluids 16 (9), 3460– 3469. Siegel, S. G. & Fasel, H. F. 2001 Eﬀect of forcing on the wake drag of an axisymmetric bluﬀ body. AIAA Paper 2001-0736. Siegel, S., Seidel, J., Cohen, K., Aradag, S. & McLaughlin, T. 2008 Open loop transient forcing of an axisymmetric bluﬀ body wake. AIAA Paper 2008-595. Soteriou, M. C. & Ghoniem, A. F. 1995 Eﬀects of the free-stream density ratio on free and forced spatially developing shear layers. Phys. Fluids 7 (8), 2036–2051. Taneda, S. 1978 Visual observations of the ﬂow past a sphere at Reynolds numbers between 104 and 106 . J. Fluid Mech. 85, 187–192. Tomboulides, A. G. & Orszag, S. A. 2000 Numerical investigation of transitional and weak turbulent ﬂow past a sphere. J. Fluid Mech. 416, 45–73. Weickgenannt, A. & Monkewitz, P. A. 2000 Control of vortex shedding in an axisymmetric bluﬀ body wake. Eur. J. Mech. B Fluids 19, 789–812.