The topic of hydrodynamic instability has already been introduced in Demonstration 9. The model considered there was taken to represent large-scale, latitudinally sheared flow, with the space coordinate representing latitude, but the model could be applied just as well to a smaller scale (or laboratory scale) flow in which there was vertical shear and the space coordinate was taken to be height. (It would then, of course, be appropriate to set $\beta$ equal to zero.) We now consider a modification to the latter model, in which the velocity is a function of height $z$, the new ingredient being that we allow the fluid to be stably stratified in density, so that buoyancy restoring forces act. One of the subtleties of this problem is that the effect of the stratification is in some cases destabilising as it allows new modes of instability over those in the unstratified problem. (It has already been remarked how the restoring effect of the Rossby-wave propagation mechanism acts to destabilize an unstratified shear flow.)

Again the stability of the flow is examined by solving an eigenvalue problem, from which the growth rate of disturbances to the flow may be predicted.

1: The model

We consider an inviscid shear flow in the $x$-direction, with the velocity varying across the flow and having components $( u , v , w ) = ( U ( z ) , 0 , 0 )$. In addition it is assumed that there is a background density distribution $\rho = \rho_{0} (z)$.

We consider the evolution of disturbances with vertical velocity

$$w ( x , z, t ) = {\rm Re}[ \hat w ( z ) \exp ( i k ( x - c t ) ) ] \eqno (11.1a)$$

and horizontal velocity

$$u ( x , z ,t ) = {\rm Re}[ \hat u ( z ) \exp ( i k ( x - c t ) ) \eqno (11.1b)$$

where ${\rm Re}(.)$ denotes taking the real part, $k$ is the wavenumber of the disturbance in the $x$-directions, and $c$ represents a complex phase speed in the $x$-direction ( $k c$ being the corresponding complex frequency). The functions $\hat w$ and $\hat u$ are related by $\hat u = (i / k ) d \hat w / d z$. When substituted into the linearised equations of motion, this gives the equation

$$( U - c ) \bigl \{ { d^{2} \hat w \over d z^{2} } - k^{2} \h... ...2} } \hat w + {J N^{2} \hat w \over U - c} = 0 . \eqno (11.2)$$

For simplicity, the disturbance has been assumed to be independent of the $y$-coordinate, though relaxing this assumption simply results in some minor changes to the above equation. We have neglected the effects of density variation in the inertial terms (the Boussinesq approximation) and defined the buoyancy frequency $N$ by

$$N^{2} (z) = -{ g \over \rho_0 } { d \rho_{0} \over d z}$$

The factor $J$ has been introduced into (11.2) as a convenient way to vary the strength of the stratification. Varying $J$ is equivalent to varying $g$ or $d \rho_0 / d z$.

For the purposes of this demonstration we shall assume that the flow is confined between rigid walls at $z = \pm L$. The appropriate boundary conditions on $\hat w$, corresponding to no flow normal to the walls, are then $\hat w = 0$ at $z = \pm L$. With this boundary condition (11.2) defines an eigenvalue problem. It is convenient to consider the wavenumber $k$ as given (and real and positive) and then seek the complex phase speed $c = c_{r} + i c_{i}$, with $c_{r}$ and $c_{i}$ real, as an eigenvalue.

Recall from Demonstration 10 that the growth rate of a disturbance of the form (11.1) is equal to $k c_{i}$, because $\exp \{ i k ( x - c t ) \} = \exp \{ i k ( x - c_{r} t ) \} \exp \{ k c_{i} t ) \}$, whose absolute magnitude is $\exp \{ k c_{i} t ) \}$. The problem of determining the stability therefore amounts to determining whether there is an eigenvalue $c$ with a positive imaginary part.

The eigenvalue problem is solved by writing (11.2) in finite-difference form and then finding the eigenvalues of the resulting matrix.

2: Running the demonstration

On entering this demonstration you must:

(i) Specify a basic flow $U ( z )$ by selecting Get U. You may draw the graph of the basic flow using the mouse, or you may give the functional form of $U ( z )$. On a small computer it is impractical to take the number $M$ of interior grid points more than about 100, and so the results will be meaningful only for fairly simple profiles representing shear layers or jets. For instance, the shear layer profile $U = \tanh ( z )$ works well with the default value $M = 19$.

(ii) The corresponding vorticity gradient $d^{2} U / d z^{2}$ is now calculated. Since this involves repeated numerical differentiation it tends to give a vorticity gradient that is very noisy when the input came from the mouse - unless you have extraordinarily steady hands! You may now apply a smoothing to $U ( z )$ and hence to the vorticity gradient, by selecting Smooth (repeatedly if necessary).

(iii) The buoyancy frequency $N$ (or equivalently the profile of density gradient) may now be specified, again using the mouse or via a functional form.

(iii) Now select Calculate. By default the programme takes what it thinks to be a reasonable value of the wavenumber $k$, and eventually a plot in the complex plane of the eigenvalues $c$ will be shown. It is those eigenvalues that appear in the upper half-plane ( $c_{i} = {\rm Re}( c ) > 0$ ) that are of physical interest, demonstrating the existence of disturbances that can grow exponentially from arbitrarily small beginnings. (Their complex-conjugate reflections represent the same disturbances with time running backwards.)

3: Notes

The method used to find the eigenvalues is fairly robust, but as in Demonstration 10 seems to fail to converge for some choices of basic flow, particularly when the number of points in the finite-difference representation is large, or there are several points over which the velocity is exactly constant. Again as in Demonstration 10, there are a large set of eigenvalues of (11.2) that have zero imaginary part. In this case, it is found that some of them are shifted off the real axis in the finite-differenced problem and therefore appear as spurious growing or decaying modes. The only reliable way to discover whether a particular calculated eigenvalue corresponds to a true solution of the continuous problem is to repeat the calculation with a larger number $M$ of grid points and see whether it persists as a calculated solution. You might like to note that a single calculation with $M = 19$ takes about $10$ seconds, with $M=39$ about $70$ seconds and with $M = 79$ about $8$ minutes.

4: Menu Options

Calculate: starts the eigenvalue calculation.

Get U(z): defines the velocity profile $U ( y )$ selecting either the Mouse or Equation options.

Get N2(z): defines the density profile $N^{2} ( z )$ selecting either the Mouse or Equation options.

Submenu: leads to:

Misc: allowing the number $M$ of points in the grid, the number of Cases to be considered and the width of the domain $L$ to be altered.

J values: specifies the value of the constant $J$ for each case.

k values: specifies the value of the wavenumber $k$ for each case.

Smoothing: changes the parameters $\mu$ and $n$ involved in the definition of the smoothing. $\mu$ controls the strength of a binomial filter and $n$ controls the number of times that it is applied each time that Smooth is selected.

Polish: this option is under development and not yet available.

Print: dumps the screen to the printer.

Smooth: applies smoothing to the flow profile $U ( z )$ and correspondingly to the vorticity gradient $d^{2} U / d z^{2}$.

Flip: swaps between a screen showing the complex phase speeds $c$ and one showing the complex frequencies $k c$. The latter screen also shows the pattern of disturbance velocities in the $(x,z)$ plane as a contour plot of $\psi$. The $z$-axis is taken to be horizontal and the $x$-axis is vertical. A plot of the flow speed $U$, the vorticity gradient $d^{2} U / d z^{2}$ and the density profile $N^2$ is shown in the same configuration for easy comparison.

E-func.: calculates the eigenfunction corresponding to a specified eigenvalue. When this option is selected the screen will first flip to that showing the phase speeds rather than the frequencies. An eigenvalue is selected by using the arrow keys to move amongst the various complex phase speeds depicted on the plot. The < $\gets$ >and < $\to$ >keys may be used to move within a given Case. The < $\uparrow$ >and < $\downarrow$ >keys may be used to move from Case to Case , if Cases is greater than 1. The eigenvalue for which the eigenfunction is to be calculated is specified by pressing the <ENTER>key. The eigenfunction is displayed by graphs of the real and imaginary parts of the complex Fourier coefficients $\hat u$ and $\hat w$, plotted against $z$. Also shown is the Reynolds stress or momentum flux $\overline { u w } = {\rm Re}( \hat u^{*} \hat w )$.

Exit: leaves this demonstration and returns to the main menu.

5: Suggested experiments

(i) Investigate the effect of simple stable stratification (e.g. constant positive $N^{2}$) on the stability of a shear layer (e.g. $U(z) = \tanh (z)$). Look at a sequence of cases where the stratification increases from zero (e.g. change the value of $J$). You should find that in most cases the stratification inhibits the growth rate of the instability and that for large enough stratification there is no instability. There is a famous theorem due to Miles and Howard that instability is possible only if the Richardson number $Ri = J N^{2} / ( d U / d z )^{2}$ is less than ${1 \over 4}$ somewhere in the flow. Check your results against this theorem (but recall the warnings in §3).

(ii) Now consider a profile which is stable with no stratification (e.g. $U (z) = z$ or $U ( z ) = \sinh ( z )$. Then add stratification. You should be able to find unstable modes in cases where the stratification has a minimum in the centre of the flow domain ( $J N^{2} = 4 \tanh ^{2} ( z )$ works well). This mode of instability may be understood as involving two gravity waves centred on either side of the minimum in $J N^{2}$ and with phase propagation relative to the flow such that the effect of the shear is to bring them to rest relative to each other. This interpretation is particularly clear in a model where there are discrete jumps in the density field, as studied by Taylor, but it seems to be useful in more general cases.

(iii) Returning the shear layer, consider the case where the stratification is strong in a narrow region in the centre of the layer (e.g. $U(z) = \tanh (z)$ and $N^{2} ( z )$ proportional to $\exp ( - (z/a)^{2})$, with $a = 10$ as a good starting choice). You will find that for small values of $J$ the modes with the fastest growth rate have zero phase speed (i.e. they propagate at the same speed as the flow at the centre of the layer). This is generally referred to as a Kelvin-Helmholz mode and is essentially the same as that in the unstratified case, with some small modification. For larger $J$ there are instead two modes with about the same growth rate and non-zero phase speed (one with positive, the other with negative, relative to the centre of the layer). Investigation of the eigenfunctions will show that each of the modes is confined to one side of the layer. These are referred to as Holmboe modes and may be understood as interacting gravity and Rossby waves - the gravity waves centred in the region of strong stratification and the Rossby waves centred on the stronger vorticity gradients at the edge of the shear layer. Enthusiasts might also investigate the changes in unstable modes as the thickness $a$ of the stratified layer changes.

(iv) The plot of complex phase speeds again shows a circle passing through the points on the real axis corresponding to the maximum and minimum values of $U ( z )$, and symmetric about the real axis. For this problem there is a mathematical theorem which states that any growing modes must have phase speeds lying in the upper semicircle. You should be able to construct simple examples of non-growing modes that have phase speeds which do not lie in the circle.