2.1 Boussinesq equations

In this manuscript, we focus specifically on the Boussinesq equations as an approximation of mass and momentum conservation when density fluctuations are suitably small. As in the typical case, we employ the Boussinesq model to simplify the governing equations, with the understanding that they are sufficiently accurate to describe oceanographic and atmospheric flows. Nevertheless, it is obviously the case that many density-varying flows of practical importance cannot be accurately described with the Boussinesq equations. The base flow sensitivity in these situations is a topic reserved for future research.

In the following section, let us assume that variables and operators marked by $(\,\tilde{\cdot }\,)$ are dimensional, and all others are non-dimensional. Let $\tilde{\boldsymbol{x}}:=\left[\begin{array}{@{}ccc@{}}\tilde{x} & {\tilde{y}} & \tilde{z}\end{array}\right]^{\text{T}}$ and $\tilde{t}$ respectively denote displacement and time, and let $\tilde{\unicode[STIX]{x1D70C}}(\tilde{\boldsymbol{x}},\tilde{t})$ , $\tilde{\boldsymbol{u}}(\tilde{\boldsymbol{x}},\tilde{t}):=\left[\begin{array}{@{}ccc@{}}\tilde{u} (\tilde{\boldsymbol{x}},\tilde{t}) & \tilde{v}(\tilde{\boldsymbol{x}},\tilde{t}) & \tilde{w}(\tilde{\boldsymbol{x}},\tilde{t})\end{array}\right]^{\text{T}}$ and $\tilde{p}(\tilde{\boldsymbol{x}},\tilde{t})$ respectively denote the density, velocity and pressure fields. Also, let $\tilde{\unicode[STIX]{x1D707}}$ denote the viscosity, $\tilde{\unicode[STIX]{x1D6FC}}$ denote the diffusivity of density (which may be, for instance, the effect of thermal conductivity or salinity diffusion), $\tilde{\boldsymbol{g}}$ denote the gravitational acceleration vector and $\tilde{\unicode[STIX]{x1D70C}}_{c}$ denote a constant density that is characteristic of the entire density field $\tilde{\unicode[STIX]{x1D70C}}(\tilde{\boldsymbol{x}},\tilde{t})$ . The dimensional Boussinesq equations we employ are

subject to $\tilde{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}=0$ and a problem-dependent set of boundary conditions. This is the form used in most descriptions of stratified flows.

There exist at least two sensible ways to non-dimensionalise these equations. In either case, we begin by selecting some characteristic length scale $\tilde{x}_{c}$ , usually determined by the size of some solid body of interest, and some characteristic velocity scale $\tilde{u} _{c}$ , usually given by a free-stream flow speed. The non-dimensional parameters associated with these scales are

and likewise for $x$ , $y$ , $z$ , $u$ , $v$ and $w$ . The associated operators are

The remaining terms are the density field and gravitational acceleration. If the flow in question is known to be linearly density stratified, then it is typically preferable to non-dimensionalise the density by the background density change $\unicode[STIX]{x0394}\tilde{\unicode[STIX]{x1D70C}}\geqslant 0$ that is encountered over a displacement of $\tilde{x}_{c}$ , or some other vertical scale, in the direction of gravity (see, e.g. de Stadler, Sarkar & Brucker Reference de Stadler, Sarkar and Brucker2010 and Ortiz et al. Reference Ortiz, Donnadieu and Chomaz2015). In this framework, the dimensional background density field is

If $\tilde{g}:=\Vert \tilde{\boldsymbol{g}}\Vert$ and we define the Brunt–Väisälä frequency (i.e. the characteristic buoyancy frequency) by

then the Froude number

indicates the ratio the buoyancy time scale to the convective time scale in the flow. The limit $Fr\rightarrow \infty$ indicates that the flow is unstratified or that gravity is absent.

In this approach, we would define the non-dimensional density

The advantage of this non-dimensionalisation is that the Froude number directly appears in the momentum equation, which is written as

where $\boldsymbol{e}_{g}$ is the unit vector in the direction of gravity and the dimensional $\tilde{g}$ has been replaced by $Fr^{-2}$ (see (6) in de Stadler et al. Reference de Stadler, Sarkar and Brucker2010). Therefore, the strength of the stratification, which is embodied by $Fr$ , is immediately apparent in the resulting non-dimensional Boussinesq equations.

For the purposes of our theoretical development, however, we wish to employ the Boussinesq equations in their general form, without requiring that the background density be linearly stratified. Therefore, we will employ an alternate non-dimensionalisation of density and gravity. This approach affords us broader applicability, at the cost that the equations may be less intuitive in the special case of stratified flows. Let us define

The Boussinesq operator $\mathscr{N}$ and associated partial differential equation then become

subject to $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ and a problem-dependent set of boundary conditions.

If we do suppose that the flow in question is stably linearly stratified, then we define $g:=\Vert \boldsymbol{g}\Vert$ and set $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}:=\unicode[STIX]{x0394}\tilde{\unicode[STIX]{x1D70C}}/\tilde{\unicode[STIX]{x1D70C}}_{c}\geqslant 0$ as the background density change over a displacement of one non-dimensional unit in the direction of gravity. The Froude number (2.6) is then

In this formulation, the choice of $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ and $g$ for a given $Fr$ is admittedly somewhat ambiguous. An examination of (2.10) reveals that $\unicode[STIX]{x1D70C}$ and $g$ can be freely scaled inversely proportionally to each other, with no effect on $\boldsymbol{u}$ . How, then, should the Froude number be set? In the stratification-based non-dimensionalisation (2.7), (2.8), some authors fix the density scale while varying gravity (e.g. de Stadler et al. Reference de Stadler, Sarkar and Brucker2010). We choose to do the opposite, which is to fix $g$ but vary the $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ scale. This approach is more physically intuitive, since in reality, the $Fr=\infty$ case is more accurately an absence of stratification than an absence of gravity. The choice of density scaling will have an effect on the inner product over infinitesimal fields to be used later (2.14).

To enable the use of linear stability theory, we will assume that there exists a steady state $\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D70C}_{0}(\boldsymbol{x}) & \boldsymbol{u}_{0}^{\text{T}}(\boldsymbol{x})\end{array}\right]^{\text{T}}$ , with $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{0}=0$ , such that

Next, we suppose that $\unicode[STIX]{x1D70C}^{\prime }(\boldsymbol{x},t)$ , $\boldsymbol{u}^{\prime }(\boldsymbol{x},t):=\left[\begin{array}{@{}ccc@{}}u^{\prime }(\boldsymbol{x},t) & v^{\prime }(\boldsymbol{x},t) & w^{\prime }(\boldsymbol{x},t)\end{array}\right]^{\text{T}}$ , and $p^{\prime }(\boldsymbol{x},t)$ are infinitesimal complex-valued perturbations in density, velocity and pressure from the steady-state solution. The linearisation of (2.10) about $\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D70C}_{0}(\boldsymbol{x}) & \boldsymbol{u}_{0}^{\text{T}}(\boldsymbol{x})\end{array}\right]^{\text{T}}$ is then given by the linear operator $\mathscr{L}$ and the associated partial differential equation

subject to $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }=0$ . The boundary conditions for (2.13) are typically the homogeneous versions of the conditions on $\mathscr{N}$  (2.10), since $\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D70C}^{\prime } & (\boldsymbol{u}_{0}+\boldsymbol{u}^{\prime })^{\text{T}}\end{array}\right]^{\text{T}}$ must satisfy the same boundary conditions as $\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D70C}_{0} & \boldsymbol{u}_{0}^{\text{T}}\end{array}\right]^{\text{T}}$ .

It is typically convenient to define an adjoint linear operator with respect to some inner product, so as to enable the computation of various types of sensitivity. For complex-valued density perturbations $\unicode[STIX]{x1D70C}^{\prime }$ , $\hat{\unicode[STIX]{x1D70C}}^{\prime }$ and velocity perturbations $\boldsymbol{u}^{\prime }$ , $\hat{\boldsymbol{u}}^{\prime }:=\left[\begin{array}{@{}ccc@{}}\hat{u} ^{\prime } & \hat{v}^{\prime } & {\hat{w}}^{\prime }\end{array}\right]^{\text{T}}$ in some control volume $\unicode[STIX]{x1D6FA}$ , along with the complex conjugation operation $(\bar{\cdot })$ , we choose the inner product

We remark that the weighting of the density and velocity terms in the above integral is arbitrary. Nevertheless, we do not introduce any additional weighting, for convenience.

The adjoint linear operator $\mathscr{L}^{\ast }$ is then defined as the operator – with an appropriate set of boundary conditions – that satisfies

For $\hat{\unicode[STIX]{x1D70C}}^{\prime }$ , $\hat{\boldsymbol{u}}^{\prime }$ and the adjoint linearised pressure $\hat{p}^{\prime }$ , this adjoint operator is given by

subject to $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{\prime }=0$ (see (21) in Ortiz et al. Reference Ortiz, Donnadieu and Chomaz2015). To clarify the notation, for two vector fields $\boldsymbol{u}$ and $\boldsymbol{v}$ given respectively in indicial notation by $u_{i}$ and $v_{i}$ , $((\unicode[STIX]{x1D735}\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{v})_{i}=(\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i})v_{j}$ . Given

and the outward normal vector $\boldsymbol{n}$ on the boundary $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ of the control volume $\unicode[STIX]{x1D6FA}$ , the boundary condition on $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ associated with $\mathscr{L}^{\ast }$ is

The associated adjoint linearised Boussinesq equations are simply

The adjoint is commonly derived using a bilinear concomitant (see Schmid & Henningson Reference Schmid and Henningson2001; Giannetti & Luchini Reference Giannetti and Luchini2007), but it can be verified perhaps more intuitively from (2.15) itself. In summary, it can be shown that

Integrating (2.20) over $\unicode[STIX]{x1D6FA}$ and applying the divergence theorem, we find that

Hence, (2.16) and (2.18) are sufficient conditions for the adjoint relation (2.15) to be satisfied.

2.2 Sensitivity to base flow modifications

The crux of linear stability theory is that the eigendecomposition of the linearised dynamics (2.13) reveals the stability characteristics of the steady-state flow. The Hartman–Grobman theorem states that the nonlinear (2.10) and linearised (2.13) dynamics will have the same stability at the fixed point $\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D70C}_{0} & \boldsymbol{u}_{0}^{\text{T}}\end{array}\right]^{\text{T}}$ , so long as none of the eigenvalues have a zero real part. We remind that the theorem makes no implications about the size of the domain of attraction in the nonlinear dynamics, nor does the linear stability analysis necessarily reveal finite-time transient growths stemming from the non-normality of $\mathscr{L}$ . The topics of nonlinear and non-normal stability are discussed in works by Schmid & Henningson (Reference Schmid and Henningson2001), Chomaz (Reference Chomaz2005) and Schmid & Brandt (Reference Schmid and Brandt2014), among others.

The sensitivity of the stability to base flow modifications, as presented by Marquet et al. (Reference Marquet, Sipp and Jacquin2008), provides an elegant closed-form solution for the eigenvalue perturbation resulting from an infinitesimal base flow perturbation. The direct and adjoint eigenmodes can be combined to construct a sensitivity field, which yields the eigenvalue perturbation via an inner product with the base flow perturbation. Whereas Marquet et al. presented the theory for constant-density flows, we now expand this theory to the Boussinesq equations, for which we assume that the density field is not constant, but is still close to a characteristic value. We will see that certain terms from Marquet et al.’s theory will remain, but additional terms stemming from the density field will appear as well.

In our notation, we will express the direct and adjoint eigendecompositions respectively as

The direct and adjoint eigenmodes are respectively

with the $\unicode[STIX]{x1D70C}$ subscript indicating the density component and the $\boldsymbol{u}$ subscript indicating the velocity component. Note that based on the definitions of the linear operators, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D753}_{\boldsymbol{u}}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D74D}_{\boldsymbol{u}}=0$ . As usual, the eigenvalues $\unicode[STIX]{x1D706}$ are the complex scalars that satisfy (2.22); if $\unicode[STIX]{x1D753}$ is the initial condition to a linearised flow, then the solution of the flow is given by $\text{e}^{\unicode[STIX]{x1D706}t}\unicode[STIX]{x1D753}$ . Through our notation, we imply that the above equations hold for all eigenvalue–eigenmode pairs, and we leave out eigenvalue–eigenmode indices for brevity.

As in previous works (e.g. Giannetti & Luchini Reference Giannetti and Luchini2007), the sensitivity analysis essentially computes what eigenvalue perturbation $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706}$ would result from an infinitesimal perturbation in the linearised dynamics $\unicode[STIX]{x1D6FF}\mathscr{L}$ . Without making any assumptions on $\unicode[STIX]{x1D6FF}\mathscr{L}$ besides its infinitesimal scale, a perturbation in (2.22a ) yields

In this form, we have not set $\unicode[STIX]{x1D6FF}\mathscr{L}$ yet, so we cannot compute $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D753}$ at this point. We can compute $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706}$ without computing $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D753}$ , however, by taking the inner product of (2.24) with $\unicode[STIX]{x1D74D}$ . Noting that

(2.24) conveniently reduces to $\left<\unicode[STIX]{x1D6FF}\mathscr{L}\unicode[STIX]{x1D753},\unicode[STIX]{x1D74D}\right>=\left<\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706}\unicode[STIX]{x1D753},\unicode[STIX]{x1D74D}\right>$ , that is,

This equation forms the basis of stability sensitivity analyses; all that remains is to choose $\unicode[STIX]{x1D6FF}\mathscr{L}$ .

In the spirit of Marquet et al. (Reference Marquet, Sipp and Jacquin2008), we will define $\unicode[STIX]{x1D6FF}\mathscr{L}$ to be the change in the linearised dynamics resulting from an infinitesimal change in the base flow. We will primarily consider two cases: a change in the base density field $\unicode[STIX]{x1D70C}_{0}(\boldsymbol{x})$ with a fixed base velocity field $\boldsymbol{u}_{0}(\boldsymbol{x})$ , and a change in $\boldsymbol{u}_{0}(\boldsymbol{x})$ with a fixed $\unicode[STIX]{x1D70C}_{0}(\boldsymbol{x})$ . In § 2.3, we provide some remarks that tie these two cases together.

In the former case, it is apparent from the linearised dynamics (2.13) that if the base flow density changes infinitesimally from $\unicode[STIX]{x1D70C}_{0}(\boldsymbol{x})$ to $\unicode[STIX]{x1D70C}_{0}(\boldsymbol{x})+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}(\boldsymbol{x})$ , then the corresponding change in the linearised dynamics is

Defining the inner product

for scalar fields, we follow (2.26) and find that the corresponding eigenvalue shift is

It is more convenient and elegant to express the above as an inner product of some quantity with $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ , so as to clarify the relation between $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706}$ . Using integration by parts, we obtain

where we have taken advantage of the fact that $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D753}_{\boldsymbol{u}}=0$ . Let us consider the physically realistic case that $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ is real valued, and assume the boundary condition that

on $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ . Then, we can define the partial sensitivity to base flow density modifications

so that the eigenvalue shift is conveniently expressed as

We specifically denote (2.32) as a partial sensitivity, because we have assumed that the base flow velocity field $\boldsymbol{u}_{0}$ is unchanged.

Next, we can perform the same analysis for the case where $\boldsymbol{u}_{0}$ undergoes an infinitesimal shift $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ , but $\unicode[STIX]{x1D70C}_{0}$ remains unchanged. In this case, the corresponding shift in the linearised dynamics (2.13) is

Following (2.26), the resulting eigenvalue shift is

Integration by parts on the last summand of the integrand yields

where we have once again used the fact that $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D753}_{\boldsymbol{u}}=0$ . Combining (2.35) and (2.36),

Thus, if we assume that $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ is real valued and that the boundary condition

is satisfied on $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ , then we can define the partial sensitivity to base flow velocity modifications

With this expression, we can employ the velocity inner product

to express the eigenvalue shift conveniently as

We note that the latter two summands in the numerator of (2.39) also appear in the constant-density base flow sensitivity (Marquet et al. Reference Marquet, Sipp and Jacquin2008), but the first summand appears specifically because of the density-varying nature of the dynamics in question.

As in Marquet et al. (Reference Marquet, Sipp and Jacquin2008), we can follow through the derivations of the partial sensitivities and determine the origins of each of the additive terms. These origins reveal the physical mechanisms underlying the partial sensitivities. First, the partial sensitivity to base flow density modifications $\unicode[STIX]{x2202}\unicode[STIX]{x1D706}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}$  (2.32) is ultimately derived from the $-\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ term in the density component of the operator perturbation $\unicode[STIX]{x1D6FF}\mathscr{L}$  (2.27). Thus, it represents the production of density perturbations by $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ .

Next, we can see that in the partial sensitivity to base flow velocity modifications (2.39), the additive term $-\bar{\unicode[STIX]{x1D713}}_{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D70C}}/\left<\unicode[STIX]{x1D753},\unicode[STIX]{x1D74D}\right>$ is derived from $-\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{\prime }$ in the density component of $\unicode[STIX]{x1D6FF}\mathscr{L}$  (2.34). Thus, it represents the transport of density perturbations by $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ , which we denote by

By the same mechanism as Marquet et al. (Reference Marquet, Sipp and Jacquin2008), the terms $-(\unicode[STIX]{x1D735}\unicode[STIX]{x1D753}_{\boldsymbol{u}})\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D74D}}_{\boldsymbol{u}}/\left<\unicode[STIX]{x1D753},\unicode[STIX]{x1D74D}\right>$ and $\unicode[STIX]{x1D753}_{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\bar{\unicode[STIX]{x1D74D}}_{\boldsymbol{u}}/\left<\unicode[STIX]{x1D753},\unicode[STIX]{x1D74D}\right>$ in (2.39) are respectively derived from $-\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }$ and $-\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ in (2.34). Therefore, they respectively represent the transport and production of velocity perturbations by $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ , which we denote by

These transport and production terms can be interpreted respectively as global analogues of convective and absolute instabilities (Marquet et al. Reference Marquet, Sipp and Jacquin2008). All the terms arising from transport and production mechanisms are summarised in table 1.

Table 1.

The additive terms in the partial sensitivities.

2.3 Remarks on total sensitivity

The partial sensitivity analysis in the prior section is most useful for understanding the effects of base flow density and velocity modifications, individually, on the flow stability. Conceptually, it is reasonable to imagine that the density and velocity modifications should be linked in practice. For instance, we may wish to determine to the change in the stability of a stratified bluff-body wake as the Froude number changes from $Fr$ to some $Fr+\unicode[STIX]{x1D6FF}Fr$ . In this case, both $\unicode[STIX]{x1D70C}_{0}$ and $\boldsymbol{u}_{0}$ would be expected to change.

If $\unicode[STIX]{x1D70C}_{0}$ and $\boldsymbol{u}_{0}$ change together, the shift in the linearised dynamics would be given by the sum of the right-hand sides of (2.27) and (2.34), that is,

It then follows that the eigenvalue perturbation is

which is in line with our interpretation of the partial sensitivities as partial derivatives.

If we are further able to express the base flow velocity modification $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ in terms of the base flow density modification $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ via the Jacobian equation

then (2.45) becomes

Hence, if we define the total sensitivity to base flow density modifications

then the eigenvalue perturbation from both $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ , as well as $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ induced by $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ , is expressed conveniently as

The difficulty in utilising the total sensitivity (2.48) is that the Jacobian $\unicode[STIX]{x2202}\boldsymbol{u}_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}$ is numerically very large and difficult to compute in its explicit form. In discretised space, its size would typically be the number of spatial dimensions times the square of the number of grid points or cells. Therefore, (2.48) is unlikely to be a tractable computation for two- or three-dimensional flows. For this reason, we will not pursue the total sensitivity in this manuscript besides to state it here, and we will base our analyses on the partial sensitivities instead. Should there exist an efficient computation of the total sensitivity, however, it would provide a valuable extension to the results presented in this manuscript.

3.1 Base flow and global modes

For the remainder of the manuscript, we discuss the $Re=30$ flow around a thin plate at a $90^{\circ }$ angle of attack. The plate is one non-dimensional unit high and 0.05 units thick, and is centred in a $200\times 200$ finite volume domain. An exception is made for the unstratified flow, for which the plate is centred in a domain that is 200 units wide and 30 units high. The full details of the computational methods are given in appendix A. We choose $Re=30$ because it is close to the critical value of 35 at which a supercritical Hopf bifurcation occurs for $Fr=\infty$ , as will be evident later. A Prandtl number of $Pr=7.19$ is chosen to match values approximately expected in oceanographic flows. The exact value of $Pr$ is not expected to make a large difference in our results (de Stadler et al. Reference de Stadler, Sarkar and Brucker2010).

The base flows at $Fr=1$ , 8, and $\infty$ are shown respectively in figures 1–3. The expected features can be seen in these figures. Most notably, lee waves can be seen in the stratified flows. Given a lengthy set of assumptions – that the stratified flow is two-dimensional, incompressible, steady and inviscid; that the upstream density and velocity are functions of only $y$ ; that there is no flow reversal; and that the Froude number is constant in the flow – the displacement $\unicode[STIX]{x1D702}$ of streamlines in lee waves is given by

(Turner Reference Turner1973, chap. 3.1). This essentially establishes that the lee wavelength is $2\unicode[STIX]{x03C0}Fr$ .

Figure 2.

The stratified flow at $Fr=8$ ; see the caption of figure 1.

Figure 3.

The unstratified flow from left to right around a vertical thin plate (green), at $Re=30$ .

At $Fr=1$ , the buoyancy effects are very strong, and notable lee waves of wavelength $2\unicode[STIX]{x03C0}$ are seen to emanate from the bluff body. Although the inertial effects of the flow are too large to support the complete blocking of the upstream flow by the bluff body (Yih Reference Yih1969), a significant upstream wake is nonetheless apparent. Furthermore, the recirculation region is shorter than in the $Fr=\infty$ case, since the fluid that is vertically deflected around the bluff body experiences a large restoring force back to the centre line. In the $Fr=8$ case, the lee waves have a larger wavelength of approximately $16\unicode[STIX]{x03C0}$ , and these waves are weaker in magnitude. In addition, the recirculation is longer than in the $Fr=1$ case because of the weaker buoyancy force. Finally, no density variations or lee waves are present when $Fr=\infty$ , as expected. Some additional discussion on the lee wavelength is given in § A.4.

Next, we employ the Arnoldi iteration to compute the leading eigenvalues and eigenmodes; see the extended discussion in § A.1 for numerical details. It is very computationally expensive to resolve multiple eigenvalues and eigenmodes for these open flows, since spatial effects extend to great distances, and buoyancy waves pervade through the entire domain. Therefore, we only report on the most converged eigenvalue and eigenmode for each Froude number. Based on the distribution of Arnoldi eigenvalues at the conclusion of the iteration, it is likely that the modes we report are in fact the least stable ones. The mode at $Fr=\infty$ is in excellent agreement with the least stable mode for the circular cylinder (e.g. Marquet et al. Reference Marquet, Sipp and Jacquin2008), and the mode at $Fr=1$ is consistent with the long-time behaviour of direct numerical simulations, as § 4.2 will show. As we will describe in more detail later, the mode we analyse for each Froude number is not a continuation of other modes we report at other Froude numbers. Therefore, we do not imply that the modes represent the same set of physics.

At $Fr=\infty$ , the Arnoldi iteration gives good convergence toward the least stable eigenvalue $\unicode[STIX]{x1D706}=-0.012\pm 0.642\text{i}$ and its complex conjugate. The corresponding direct and adjoint eigenmodes are shown in figure 4. As we discuss in greater detail in § A.2.1, the density components of these modes is $\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D70C}}=0$ , so the modes are identical to the constant-density modes, even allowing for non-zero density perturbations $\unicode[STIX]{x1D70C}^{\prime }$ in the linearised equations. The eigenmodes of unstratified flows around bluff bodies are well understood at this point and have been discussed in many works (e.g. Giannetti & Luchini Reference Giannetti and Luchini2007). Figure 4 depicts the common features. The direct eigenmode is localised in a long downstream region, with spatially oscillating structures. On the other hand, the adjoint eigenmode is large immediately downstream of the bluff body, but the remaining wake is actually upstream because of the reversed base flow advection in (2.16).

Figure 4.

The (a) direct and (b) adjoint velocity eigenmodes of the $Fr=\infty$ flow, corresponding to the least stable eigenvalue $\unicode[STIX]{x1D706}=-0.012+0.642\text{i}$ ; shown as the pointwise magnitude and streamlines of (a) $\text{Re}(\unicode[STIX]{x1D753}_{\boldsymbol{u}})$ and (b) $\text{Re}(\unicode[STIX]{x1D74D}_{\boldsymbol{u}})$ .

In contrast, the eigenmodes contain remarkably different features at $Fr=1$ . For this computation, the Arnoldi iteration resolves the eigenvalue $\unicode[STIX]{x1D706}=-0.001\pm 0.902\text{i}$ . The direct eigenmode is shown in figure 5, and the corresponding adjoint eigenmode is shown in figure 6. It is immediately apparent that large-scale structures (i.e. on the scale of the entire numerical mesh) are dominant in figures 5(a) and 6(a). These structures primarily consist of large vertical bands of density and vertical velocity; the density and vertical velocity are horizontally out of phase by a quarter period. The effects of the lee waves (figure 1) can also be seen in the direct eigenmodes. These large-scale structures are further analysed from a parallel flow perspective in § 4.1.

Figure 5.

The direct eigenmode at $Fr=1$ corresponding to $\unicode[STIX]{x1D706}=-0.001-0.902\text{i}$ , shown as the real part of density (blue: $\text{Re}(\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D70C}})<0$ ; red: $\text{Re}(\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D70C}})>0$ ) and streamlines of $\text{Re}(\unicode[STIX]{x1D753}_{\boldsymbol{u}})$ . (a) The domain $[-60,60]\times [-60,60]$ ; the right region rotates clockwise. (b) The domain $[-2,14]\times [-8,8]$ .

Figure 6.

The adjoint eigenmode at $Fr=1$ corresponding to figure 5, shown as the real part of density (blue: $\text{Re}(\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D70C}})<0$ ; red: $\text{Re}(\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D70C}})>0$ ) and streamlines of $\text{Re}(\unicode[STIX]{x1D74D}_{\boldsymbol{u}})$ . (a) The domain $[-60,60]\times [-60,60]$ ; the right region rotates clockwise. (b) The domain $[-7,2]\times [-4.5,4.5]$ .

A closer look at figures 5(b) and 6(b) reveals density and velocity structures associated with the downstream wake of the bluff body. These structures, however, are nowhere near as dominant as the downstream velocity structures for $Fr=\infty$ . In fact, the velocity streamlines of the direct and adjoint eigenmodes remain nearly vertical in the upstream and downstream wakes of the bluff body. In contrast, the eigenmode streamlines at $Fr=\infty$ (figure 4) form closed loops on the length scale of the bluff body.

Another important feature of the eigenmodes is the appearance of angled waves emanating from the bluff body. In the direct eigenmode (figure 5 a), an envelope downstream of the bluff body initially opens at an angle of $\unicode[STIX]{x1D703}=13^{\circ }$ from the horizontal direction. The angle increases to about $24^{\circ }$ at a distance of approximately $x=20$ downstream of the body. Wave crests can be seen inside this envelope. The envelope itself roughly demarcates the end of the lee waves in the base flow (figure 1). The adjoint eigenmode (figure 6 a) similarly contains an envelope – which is now upstream because of the reversed base flow direction in (2.16) – but the wave crests reside in the exterior of the envelope. This envelope initially opens at approximately $\unicode[STIX]{x1D703}=19^{\circ }$ before expanding to approximately $28^{\circ }$ further upstream.

It is well known that in an inviscid Boussinesq flow with a mass source or force oscillating at dimensional angular frequency $\tilde{\unicode[STIX]{x1D714}}$ in an open domain, gravity waves are emitted along beams at an angle of $\unicode[STIX]{x1D703}=\cos ^{-1}(\tilde{\unicode[STIX]{x1D714}}/{\tilde{N}})$ from the vertical (Voisin Reference Voisin1991; see also § A.4). This particular eigenmode has an angular frequency of $|\text{Im}(\unicode[STIX]{x1D706})|=0.901$ for a Brunt–Väisälä frequency of $N=1$ , with both quantities non-dimensionalised by convective time. Thus, the corresponding angle of wave transmission is predicted under this theory to be $\cos ^{-1}0.901=26^{\circ }$ . This angle is consistent with those observed in the direct and adjoint eigenmodes beyond the immediate wakes of the bluff body.

We remark that the use of a finite domain with artificial boundary conditions presents challenges in the determination of large-scale modal structures. In particular, the finite domain necessarily implies that any continua of eigenvalues and eigenmodes would be represented discretely instead. This may imply, for instance, that the horizontal wavelength is determined in part by the mesh size and the boundary conditions at the inlet and outlet. Nevertheless, the eigenmodes depicted in figures 5 and 6 still contain important near-wake features associated with the bluff body. We will see next that it is these wake features, and not the large-scale structures, that end up being important in the sensitivity calculations.

3.2 Sensitivity

4.1 Predictions from parallel flow theory

In § 3.1, the direct and adjoint eigenmodes for the $Fr=1$ flow were depicted, and their features were described. It was then shown in § 3.2 that the eigenvalue sensitivity regions only occupied a small region local to the bluff body, as opposed to the large-scale structures found in the eigenmodes. In this section, we use inviscid and viscous parallel flow theory to explain some of these features. The main objectives in this section are to show that the global eigenmodes are consistent with predictions from canonical simplified representations, and to demonstrate in § 4.1.2 how the localised sensitivity structures arise from global eigenmode structures. Since we are now working exclusively with a stratified flow, we will use the alternate non-dimensionalisation of density (2.7), (2.8) that is more appropriate for stratified flows.

4.2 Linearised simulations with modified base flow density

With the partial sensitivity theory developed in § 2.2, it is then natural to inquire what effects certain base flow modifications $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ and $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ would have on the linearised dynamics (2.13). Namely, if a certain $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ or $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ is expected to destabilise the flow by giving $\text{Re}(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706})>0$ via (2.33) or (2.41), we would like to see what the resulting instability is. To that end, we perform time-resolved numerical experiments using the linearised Boussinesq equations (2.13) with $Fr=1$ . The simulations start with initial conditions $\unicode[STIX]{x1D70C}^{\prime }(\boldsymbol{x},0)$ and $\boldsymbol{u}^{\prime }(\boldsymbol{x},0)$ randomised from a uniform distribution at every grid point. Although such a velocity perturbation field will have non-zero divergence, the incompressible flow solver will ensure that $\boldsymbol{u}^{\prime }(\boldsymbol{x},t)$ at all subsequent time steps will be divergence free.

The density and velocity perturbation fields at $t=835.65$ is shown in figure 15. It is evident, for instance, that figure 15(a) contains the same downstream density wake structure as the eigenmode shown in figure 5(b). On the other hand, figure 15(c) shows that $\boldsymbol{u}^{\prime }$ consists primarily of vertical motions, as the eigenmode also exhibits in figure 5(b). Some velocity distortions are also apparent immediately downstream of the plate. Whereas figure 15(c) shows that the mode streamlines are primarily oriented upwards, figure 15(b) shows that downstream of the plate, the streamlines deflect to the right near the bottom, and to the left near the top. This, again, is in agreement with the global eigenmode in figure 5(b). Such a resemblance between the time-resolved flow and the eigenmode does not come as a surprise, since this eigenmode is one of the least stable ones present for this base flow.

Figure 15.

A snapshot of a time-resolved linearised Boussinesq simulation using a random initial condition. Shown as (a)  $\unicode[STIX]{x1D70C}^{\prime }$ , (b)  $u^{\prime }$ and (c)  $v^{\prime }$ . Red and blue respectively indicate positive and negative values.

To investigate the partial sensitivity to base flow density modifications, we next run a simulation with the base flow density perturbed to $\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ , where $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}=500\text{Re}(\text{d}\unicode[STIX]{x1D706}/\text{d}\unicode[STIX]{x1D70C}_{0})$ . Although the multiplicative coefficient may look large, we remind that $\text{Re}(\text{d}\unicode[STIX]{x1D706}/\text{d}\unicode[STIX]{x1D70C}_{0})$ is very small to begin with; see figure 9. The size of $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ may be prohibitively large to prevent an accurate analysis as an infinitesimal base flow modification, but we choose this size for the ease of visual comparison to the unmodified base flow. For this experiment, the base flow velocity $\boldsymbol{u}_{0}$ and the initial conditions $\unicode[STIX]{x1D70C}^{\prime }(\boldsymbol{x},0)$ and $\boldsymbol{u}^{\prime }(\boldsymbol{x},0)$ remain unchanged.

Following the inner product equation relating $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706}$  (2.33), we would expect this base flow density perturbation to destabilise the flow. A snapshot of the linearised flow at a large $t$ is shown in figure 16. To highlight the slight destabilisation provided by this base flow modification, we emphasise that $\unicode[STIX]{x1D70C}^{\prime }$ is approximately an order of magnitude larger in figure 16(a) than in figure 15(a). The primary effect of this $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ in $\unicode[STIX]{x1D70C}^{\prime }$ is to create a slightly unstable upstream oscillation along the centreline. This density wave prominently forms a few units upstream of the bluff body, propagates toward the body, and is horizontally compressed and vertically stretched as it meets the body. Part of the wave then climbs around the body and extends into the recirculation region.

Figure 16.

A snapshot of a time-resolved linearised Boussinesq simulation using a random initial condition and a modified base flow density $\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ , where $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}=500\text{Re}(\text{d}\unicode[STIX]{x1D706}/\text{d}\unicode[STIX]{x1D70C}_{0})$ . Shown as (a)  $\unicode[STIX]{x1D70C}^{\prime }$ (with colours slightly exaggerated for clarity), (b)  $u^{\prime }$ and (c)  $v^{\prime }$ . Red and blue respectively indicate positive and negative values.

An effect of $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ can also be seen in $\boldsymbol{u}^{\prime }$ (figure 16 b,c). In particular, figure 16(c) shows the presence of a horizontal wave in $v^{\prime }$ that is not present in the absence of $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ (figure 15 c). This velocity wave begins just upstream of the body and extends a few units downstream, and also propagates in the positive $x$ direction. Since the linearised velocity equation has not changed, the change in the behaviour of $v^{\prime }$ must be understood as an effect of $\unicode[STIX]{x1D70C}^{\prime }$ . Specifically, the linearised dynamics (2.13) reveals that if $\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{0}$ points in the positive $y$ direction – e.g. upstream of the plate as in figure 9(a) – then a region with $v^{\prime }<0$ would see an increase in $\text{d}\unicode[STIX]{x1D70C}^{\prime }/\text{d}t$ via the density equation. The density perturbation $\unicode[STIX]{x1D70C}^{\prime }$ is then coupled back to $v^{\prime }$ via the momentum equation, such that an increase in $\unicode[STIX]{x1D70C}^{\prime }$ would cause a decrease in $\text{d}v^{\prime }/\text{d}t$ though the gravity term.

Figure 17.

A snapshot of a time-resolved linearised Boussinesq simulation using a random initial condition and a modified base flow velocity $\boldsymbol{u}_{0}+\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ , where $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}=350\text{Re}(\text{d}\unicode[STIX]{x1D706}/\text{d}\boldsymbol{u}_{0})$ . Shown as (a)  $\unicode[STIX]{x1D70C}^{\prime }$ , (b)  $u^{\prime }$ , and (c)  $v^{\prime }$ . Red and blue respectively indicate positive and negative values.

Figure 17 is a snapshot of the linearised simulation with the unperturbed base flow density $\unicode[STIX]{x1D70C}_{0}$ , but with a modified base flow velocity $\boldsymbol{u}_{0}+\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ , where $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}=350\text{Re}(\text{d}\unicode[STIX]{x1D706}/\text{d}\boldsymbol{u}_{0})$ . Once again, the initial conditions $\unicode[STIX]{x1D70C}^{\prime }(\boldsymbol{x},0)$ and $\boldsymbol{u}^{\prime }(\boldsymbol{x},0)$ are unchanged. Since the growth rate sensitivity $\text{Re}(\text{d}\unicode[STIX]{x1D706}/\text{d}\boldsymbol{u}_{0})$ only begins just upstream of the bluff body and primarily extends downstream (figure 9 d), $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ has little effect on the linearised dynamics upstream of the body. The velocity perturbation, however, compresses the recirculation bubble in both the horizontal and the vertical directions. The primarily vertical orientation of $\boldsymbol{u}^{\prime }$ is in rough agreement with the eigenmode shown in figure 5. As in the simulation with the modified base flow density (figure 16), however, this base flow velocity modification once again introduces waves that are coupled primarily between $\unicode[STIX]{x1D70C}^{\prime }$ and $v^{\prime }$ . The precise mechanisms by which $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$ modifies the linearised dynamics are difficult to analyse because of the varying vector directions of $\text{Re}(\text{d}\unicode[STIX]{x1D706}/\text{d}\boldsymbol{u}_{0})$ , as well as the three different terms through which $\unicode[STIX]{x1D6FF}\mathscr{L}$ depends on $\unicode[STIX]{x1D6FF}\boldsymbol{u}_{0}$  (2.34).

It is important to note that in these numerical experiments, a base flow modification would not only modify the eigenmode $\unicode[STIX]{x1D753}$ in question to $\unicode[STIX]{x1D753}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D753}$ (see (2.24)), but it would also stabilise or destabilise other eigenmodes. In the direct numerical simulations as they are, it is not possible to determine how much each eigenvalue and eigenmode shift is responsible for the altered dynamics in figures 16 and 17, especially since the base flow modification had to be large for the effects to be clearly visible. The distinction between these effects can be pursued further in future research.