An introduction to the Kelvin-Helmholtz Instability

We will now explore a mathematical demonstration of the Kelvin-Helmholtz instability in its most basic form. Let’s begin by considering the simple flow schematic shown below, two horizontal, parallel, infinite streams of fluid, each uniformly moving in opposite directions; fluid one is blue, fluid two is green. In a fuller description of the KH instability, we could consider the fluids as being of different densities, but we are going to keep things extremely simple here, and imagine both streams are of the same density, and the only difference is in the velocity. The line dividing these two fluids can be considered to be a “vortex sheet“, as it is a discontinuity in the fluid velocity.

Our task is now to consider the stability of this very simple system. For our purposes here, we will consider the most basic possible analysis, which is to assume that the fluid is both inviscid and incompressible. The analysis we are going to conduct, is to consider the equations of motion that describe this system, we will then linearize the system, and consider the response of this linear system to infinitely small perturbations; if the system grows in response to an infinitely small perturbation, it is unstable. To determine if this system is stable, we will undertake the following steps:

  1. Establish the equations of motion.
  2. Linearize these equations of motion.
  3. Propose an ansatz (a model solution) and substitute it into our linearized equations
  4. Obtain an ordinary differential equation from re-arranging our linearized equations
  5. Impose appropriate boundary conditions to reframe the problem as an eigenvalue problem
  6. Solve the eigenvalue problem to determine the stability of the system

 

To signpost what we’ll be doing, here it is represented diagramatically.

 

Establish the equations of motion.

The governing equations of motion for fluids are the conservation of mass, momentum, and energy. At this early stage we can omit the conservation of energy statement, as our assumption of incompressibility means we don’t need it to form a solution. We will also remove a number of terms by assuming that the flow is inviscid and two-dimensional. For compactness, we will write our initial conservation statements using vector notation, we could make this more compact still using tensor notation, but as that is unfamiliar to many people when starting out, we’ll stick with vector form. If you are unfamiliar with the vector form, click the little “expand” option to see them in standard differential form.

We are going to rework these equations into a number of different forms, so to help keep track of them, I will colour code them. When we arrive at key results with these equations in different form, I will use the same colours.

Conservation of Mass in Vector Form (2D, Incompressible)

    \[ \centering \nabla \cdot\mathbf{u} = 0 \]

Expand to see in non-vector form

    \[ \centering \frac {\partial u} {\partial x} + \frac{\partial v}{\partial y} = 0 \]

Conservation of X-Momentum in Vector Form (2D, Incompressible, Inviscid)

    \[ \centering \rho ( \frac{\partial u}{\partial t}+\nabla u \cdot \mathbf{u})=-\frac{dp}{dx} \]

Expand to see in non-vector form

    \[ \centering \rho ( \frac{\partial u}{\partial t}+ u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y}) =-\frac{dp}{dx} \]

Conservation of Y-Momentum in Vector Form (2D, Incompressible, Inviscid)

    \[ \centering \rho ( \frac{\partial v}{\partial t}+\nabla v \cdot \mathbf{u} )=-\frac{dp}{dy} \]

Expand to see in non-vector form

    \[ \centering \rho ( \frac{\partial v}{\partial t}+ u\frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}) =-\frac{dp}{dy} \]

Here u and v are our velocity components in the x and y direction, and \mathbf{u} is a vector of these velocities \mathbf{u}=[u,v].

And that’s it, these are our equations of motion!

 

Linearize the equations of motion.

 

As a first step towards linearization, we now apply the Reynolds Decomposition to these equations, meaning we decompose all variables (q=[u,v,p]) into a temporal mean \bar{q} and fluctuating q' component.

    \begin{equation*} q=\bar{q}+q' \end{equation*}

First, we apply this to the continuity equation:

Conservation of Mass – Reynolds Decomposition

    \begin{equation*} \nabla \cdot \mathbf{u} = 0 \nonumber \end{equation*}

See details of substitution

When we are not used to manipulating these equations, it is easiest to expand this into its full form:

    \begin{equation*} \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} = 0 \nonumber \end{equation*}

We then substitute in our Reynolds-decomposed variables:

    \begin{equation*} \frac{\partial \bar{u}}{\partial x}+\frac{\partial u'}{\partial x}+\frac{\partial \bar{v}}{\partial y}+\frac{\partial v'}{\partial y} = 0 \\ \nonumber \end{equation*}

Now, to simplify this further, we consider our initial problem description, which is that we have two infinite parallel streams, each uniform within its region..

  • If the flows are parallel, they can have no vertical velocity component, so \frac{\partial \bar{v}}{\partial y}=0.
  • Since each flow is uniform within its region, they can have no velocity gradient in x, so \frac{\partial \bar{u}}{\partial x}=0

 

 

Thus our equation reduces to:

    \[ \centering \boxed{\frac {\partial u'} {\partial x} + \frac{\partial v'}{\partial y} = 0} \]

 

Conservation of x-momentum – Reynolds Decomposition

We now apply the same decomposition to the x-momentum equation. The decomposition for the time derivative and pressure gradient are straightforward, so these can be done first:

    \begin{equation*} \rho ( \frac{\partial u}{\partial t}+\nabla u \cdot \mathbf{u})=-\frac{dp}{dx} \\ \nonumber \end{equation*}

Becomes:

    \begin{equation*} \frac{\partial \bar{u}}{\partial t}+\frac{\partial u'}{\partial t}+\nabla u \cdot \mathbf{u}=\frac{1}{\rho}(-\frac{d\bar{p}}{dx}-\frac{dp'}{dx}) \end{equation*}

 

Now, we expand the term \nabla u \cdot \mathbf{u} and perform our Reynolds decomposition and substitution:

See details of substitution

    \begin{equation*} \nabla u \cdot \mathbf{u}=u \frac{\partial{u}}{\partial x} + v \frac{\partial{u}}{\partial y} \end{equation*}

And now substitute, using the product rule:

    \begin{eqnarray*} u \frac{\partial{u}}{\partial x}=\bar{u}\frac{\partial{\bar{u}}}{\partial x}+\bar{u}\frac{\partial{u'}}{\partial x}+u'\frac{\partial{\bar{u}}}{\partial x}+u'\frac{\partial{u'}}{\partial x} \\ \nonumber v \frac{\partial{u}}{\partial y}=\bar{v}\frac{\partial{\bar{u}}}{\partial y}+\bar{v}\frac{\partial{u'}}{\partial y}+v'\frac{\partial{\bar{u}}}{\partial y}+v'\frac{\partial{u'}}{\partial y} \end{eqnarray*}

 

Now, if we were to substitute these two into our x-momentum equation, we would end up with an additional 8 terms, which is going to end up with a very small fontsize if we want to fit it on one page. Luckily, there are two steps we can take to simplify here. The first is a repeat of the simplifications made when considering the continuity equation: our flows are parallel and infinite, and uniform within their region. Thus once again, \bar{v}= \frac{\partial \bar{u}}{\partial x} = \frac{\partial \bar{u}}{\partial y} = 0. Once these terms are set to zero, our eight terms reduce to three:

    \begin{eqnarray*} u \frac{\partial{u}}{\partial x}=\bar{u}\frac{\partial{u'}}{\partial x}+u'\frac{\partial{u'}}{\partial x} \\ \nonumber v \frac{\partial{u}}{\partial y}=v'\frac{\partial{u'}}{\partial y} \end{eqnarray*}

 

At this stage of a section called “linearization”, we have not yet actually linearized anything! We now introduce our first linearizing assumption: our fluctuations are small relative to our mean, sufficiently small that the product of fluctuations can be neglected. In the above equation, two of our three terms appear to involve products of fluctuations (u'\frac{\partial{u'}}{\partial x}, v'\frac{\partial{u'}}{\partial y}). However, strictly speaking, these are actually the products of fluctuations and the derivative of fluctuations, which is not quite the same thing; even if the fluctuations are small, their derivatives may not be. However, we can perform a little re-arrangement to get the result we need:

Expand to see details of re-arranging

 

Let’s start with something that we know is equal to zero, i.e. the product of two fluctuations \mathbf{u'}u'. If we take the divergence of this product:

    \begin{equation*} \nabla \cdot (\mathbf{u'}u')= u'\nabla \cdot\mathbf{u'}+\nabla u' \cdot \mathbf{u'} \end{equation*}

We have already shown from continuity that \nabla \cdot \mathbf{u'}=0, and the term (\mathbf{u'}u') on the LHS is the product of two fluctuations, so is also assumed to be zero under our linear approximation, therefore \nabla u' \cdot \mathbf{u'}=0, demonstrating in this case the products of the fluctuating components with their derivatives are also zero.

 

We set u'\frac{\partial{u'}}{\partial x}, v'\frac{\partial{u'}}{\partial y} to zero, and are left only with a single term:

    \begin{equation*} \nabla u \cdot \mathbf{u}=\bar{u}\frac{\partial{u'}}{\partial x} \end{equation*}

 

Thus our original equation becomes:

 

    \begin{equation*} \frac{\partial \bar{u}}{\partial t}+\frac{\partial u'}{\partial t}+\bar{u}\frac{\partial{u'}}{\partial x}=\frac{1}{\rho}(-\frac{d\bar{p}}{dx}-\frac{dp'}{dx}) \end{equation*}

 

 

 

We remove a few remaining terms by considering that the temporal mean of a time derivative must be zero by definition, so \frac{\partial \bar{u}}{\partial t}=0, and under the assumptions of inviscid, infinite, uniform flow, there can be no mean pressure gradient: \frac{d\bar{p}}{dx}=0.

 

It is worth noting at this point that we are removing these terms based on arguments that are intended to be more intuitive than rigorous. Later, when we revisit stability analysis for flows with less simplifying assumptions, we will be more rigorous in our justification for how these equations are simplified.

 

Thus concludes our substitution of the Reynolds decomposition, and our linearization of the x-momentum equation:

 

    \[ \centering \boxed{\frac{\partial u'}{\partial t}+\bar{u}\frac{\partial{u'}}{\partial x}=\frac{1}{\rho}(-\frac{dp'}{dx})} \]

 

 

 

Conservation of y-momentum – Reynolds Decomposition

The procedure for the y-momentum linearization is identical to the x-momentum equation, so you can click here to see a few intermediate steps, or just skip straight to the result:

Expand to see a few extra steps for y-momentum substitution

    \begin{equation*} \rho ( \frac{\partial v}{\partial t}+\nabla\cdot(v \mathbf{u} ))=-\frac{dp}{dy} \end{equation*}

 

    \begin{equation*} \frac{\partial \bar{v}}{\partial t}+\frac{\partial v'}{\partial t}+v\nabla \cdot\mathbf{u}+\nabla v \cdot \mathbf{u}=\frac{1}{\rho}(-\frac{dp'}{dy}) \end{equation*}

    \[ \centering \boxed{\frac{\partial v'}{\partial t}+\bar{u}\frac{\partial{v'}}{\partial x}=\frac{1}{\rho}(-\frac{dp'}{dy})} \]

 

 

Linearized Incompressible Inviscid 2D Equations of Motion

Try saying that three times fast. Here they are:

    \[ \frac {\partial u'} {\partial x} + \frac{\partial v'}{\partial y} = 0 \]

    \[ \frac{\partial u'}{\partial t}+\bar{u}\frac{\partial{u'}}{\partial x}=\frac{1}{\rho}(-\frac{dp'}{dx}) \]

    \[ \frac{\partial v'}{\partial t}+\bar{u}\frac{\partial{v'}}{\partial x}=\frac{1}{\rho}(-\frac{dp'}{dy}) \]

 

Propose an ansatz

We now have a system of linear equations, but just because we have a set of equations, does not mean that we know how to solve them. To work towards a solution, we will start by proposing an ansatz, i.e. a guess of what form the solution might take. There are a few different ways we could justify our approach here, but for now, let us do a post hoc justification based on empiricism; we have already seen from the experiments in the videos above that our instabilities tend to result in wavelike shapes, so our proposed solution should be in the form of a wave. When we are considering solutions in the form of waves, we can describe these solutions in terms of normal modes. In a system that undergoes some kind of oscillatory behaviour:

A normal mode is the motion in which all parts of the system move sinusoidally with the same frequency and with a fixed phase relation.

So, we have two key assumptions thus far:

  • We have a linear system (i.e. a set of linear equations)
  • We are assuming a normal-mode ansatz

Together, these two assumptions allow us to use the method of normal modes, which has two useful outcomes for us:

  • Each normal mode satisfies the linear system, and can thus be treated independently – the solutions are independent.
  • The summation of all the individual normal modes represents the complete motion of the system.

These statements are really true for any linear system, but if we are trying to find the motion of a system that has oscillatory behaviour, wavelike motion is the correct way to describe it, and thus normal modes are the appropriate way to break it down.

Essentially, we are going to consider the system that we began with above, experiencing perturbations that look like this:

 

Mathematically, we do this by proposing the following ansatz:

    \begin{equation*} q'=\hat{q}(y)e^{i(kx-\omega t)} \end{equation*}

Here q'=[u',v',p'], i.e. we are saying that we will use this ansatz for all of the variables in our linearized equations. \hat{q}(y) is a function only of y; the ansatz we have proposed suggests that the waves we are interested in are travelling in the x-direction, not the y-direction. If the perturbations are all in the form of waves, then it’s appropriate the consider the x-direction in terms of periodic functions, but we have made no such assumption in the y-direction. So we have an unknown function \hat{q}(y) to describe the behaviour of our perturbations as we change our position in y. Working out what this function is will come later.

Now we come to the exponential term e^{i(kx-\omega t)}. This is the actual normal mode part, the part that dictates all parts of the system move sinusoidally with the same frequency and with a fixed phase relation. If necessary, revise Euler’s formula as a means of representing periodic functions.

Now that we have our linearized equations and our ansatz, the next step is to substitute the ansatz into these linearized equations

Ansatz Substitution

So, just so we have them clearly written out, let’s explicitly state all three of our primary variables in terms of our ansatz:

    \begin{eqnarray*} u'=\hat{u}(y)e^{i(kx-\omega t)} \\ \nonumber v'=\hat{v}(y)e^{i(kx-\omega t)} \\ \nonumber p'=\hat{p}(y)e^{i(kx-\omega t)} \\ \nonumber \end{eqnarray*}

If we look back at our linearized equations, what we actually need are the derivatives of these variables with respect to space (x,y) and time (t), in terms of our ansatz. When evaluating these derivatives, consider that the exponential term is a function of (x,t), while the shape function term is a function only of y. Given this separation of variables, the derivatives can be determined in a single operation without the need for any application of the product rule:

Expand to see full list of derivatives in terms of normal-mode ansatz

    \begin{eqnarray*} \frac{\partial u'}{\partial t}=-i \omega \hat{u}(y)e^{i(kx-\omega t)} \\ \nonumber \frac{\partial u'}{\partial x}=i k \hat{u}(y)e^{i(kx-\omega t)} \\ \nonumber \frac{\partial v'}{\partial t}=-i \omega \hat{v}(y)e^{i(kx-\omega t)} \\ \nonumber \frac{\partial v'}{\partial x}=i k \hat{v}(y)e^{i(kx-\omega t)} \\ \nonumber \frac{\partial v'}{\partial y}= \frac{\partial \hat{v}}{\partial y} e^{i(kx-\omega t)} \\ \nonumber \frac{\partial p'}{\partial x}=i k \hat{p}(y)e^{i(kx-\omega t)} \\ \nonumber \frac{\partial p'}{\partial y}= \frac{\partial \hat{p}}{\partial y} e^{i(kx-\omega t)} \\ \nonumber \end{eqnarray*}

 

Now it is simply a matter of substituting these terms into our linearized equations of motion.

 

Conservation of Mass – Ansatz Substitution

    \[ \frac {\partial u'} {\partial x} + \frac{\partial v'}{\partial y} = 0 \]

Two straightforward substitutions produces:

    \begin{equation*} ik\hat{u} + \frac{\partial \hat{v}}{\partial y} = 0 \end{equation*}

For compactness we will now use the notation \hat{v}_y=\frac{\partial \hat{v}}{\partial y}

    \[ ik\hat{u} + \hat{v}_y = 0 \]

 

Note that here and for the substitutions that follow, the exponential component e^{i(kx-\omega t)} is removed from both sides of the equation by division.

Conservation of x-Momentum- Ansatz Substitution

    \[ \frac{\partial u'}{\partial t}+\bar{u}\frac{\partial{u'}}{\partial x}=\frac{1}{\rho}(-\frac{dp'}{dx}) \]

Again, we need only make three straightforward substitutions here:

    \[ -i \omega \hat{u} +\bar{u} i k \hat{u} =-\frac{1}{\rho}ik\hat{p} \]

Conservation of y-Momentum- Ansatz Substitution

    \[ \frac{\partial v'}{\partial t}+\bar{u}\frac{\partial{v'}}{\partial x}=\frac{1}{\rho}(-\frac{dp'}{dy}) \]

And once again, straightforward substitution produces:

 

    \[ -i \omega \hat{v} +\bar{u} ik\hat{v}=-\frac{1}{\rho}\hat{p}_y \]

 

Re-arrange ansatz form of linearized equations into a single ordinary differential equation

So, we have our linearized conservation of mass, x-momentum and y-momentum, expressed in terms of our normal-mode ansatz:

    \[ ik\hat{u} + \hat{v}_y = 0 \]

    \[ -i \omega \hat{u} +\bar{u} i k \hat{u} =-\frac{1}{\rho}ik\hat{p} \]

    \[ -i \omega \hat{v} +\bar{u} ik\hat{v} =-\frac{1}{\rho}\hat{p}_y \]

Remembering that our goal here is to find an equation that we can actually solve, our next step is to try to reduce this set of equations into a single equation we can solve. At this point, we need to decide which variable we wish to look for in our solution. Though there is no “correct” answer, the most common approach is to seek a solution in terms of the pressure perturbation. So our goal here is an equation in terms of only pressure (and its derivatives). We begin by re-arranging our conservation of mass expression, to produce an expression for \hat{u}:

    \[ ik\hat{u} + \hat{v}_y = 0 \rightarrow \hat{u} = - \frac{\hat{v}_y}{ik} \]

We now substitute this into our x-momentum equation:

    \begin{equation*} \frac{\omega}{k}\hat{v}_y-\bar{u}\hat{v}_y=-\frac{1}{\rho}ik\hat{p} \end{equation*}

We now need a way to substitute this into our y-momentum equation, which may not at first be obvious; the y-momentum equation is in terms \hat{v} and \hat{p}_y, while the equation we’ve just obtained here is in terms of \hat{v}_y and \hat{p}. One equation concerns pressure and the derivative of velocity, the other concerns velocity and the derivative of pressure. Remember here that our goal is to eliminate all non-pressure terms from the equation, and the best way to do this is to take the y-derivative of the y-momentum equation:

    \begin{eqnarray*} LHS: \frac{\partial}{\partial y}(-i \omega \hat{v} +\bar{u} ik\hat{v})= -i \omega \hat{v}_y +\bar{u} ik\hat{v}_y \\ \nonumber RHS: \frac{\partial}{\partial y}(-\frac{1}{\rho}\hat{p})_y= -\frac{1}{\rho}\hat{p}_{yy} \\ \nonumber \end{eqnarray*}

Yielding:

    \begin{equation*} -i \omega \hat{v}_y +\bar{u} ik\hat{v}_y = -\frac{1}{\rho}\hat{p}_{yy} \end{equation*}

This can be re-arranged as:

    \begin{equation*} \hat{v}_y= -\frac{1}{\rho i (\bar{u}k-\omega)}\hat{p}_{yy} \end{equation*}

Which is now ready for substitution into the equation we obtained from continuity and x-momentum:

Expand to see details of substitution

    \begin{equation*} \frac{\omega}{k}\hat{v}_y-\bar{u}\hat{v}_y=-\frac{1}{\rho}ik\hat{p} \end{equation*}

    \begin{equation*} \frac{\omega}{k} (-\frac{1}{\rho i (\bar{u}k-\omega})\hat{p}_{yy} +\frac{\bar{u}}{\rho i (\bar{u}k-\omega}\hat{p}_{yy}=-\frac{1}{\rho}ik\hat{p} \end{equation*}

    \begin{equation*} \hat{p}_{yy}( \frac{\bar{u}}{i (\bar{u}k-\omega}-\frac{\omega}{k i (\bar{u}k-\omega}) +ik\hat{p}=0 \end{equation*}

 

Multiply by ik

    \begin{equation*} \hat{p}_{yy}( \frac{\bar{u}k-\omega}{\bar{u}k-\omega}) - k^2\hat{p}=0 \end{equation*}

Which simplifies to:

 

    \begin{equation*} \hat{p}_{yy} - k^2\hat{p}=0 \end{equation*}

And because math is beautiful, we end up with

 

    \[ \boxed{\hat{p}_{yy} - k^2\hat{p}=0} \]

 

This single expression represents the conservation of mass, momentum and energy, for our linearized system, under the assumption of our normal mode ansatz. And this is the kind of equation that can we solve.

 

Impose appropriate boundary conditions to reframe the problem as an eigenvalue problem

 

What we have now is a homogeneous second-order differential equation for the pressure perturbation pressure, in terms only of the wavenumber k. The solution to this simple ODE is well known:

\hat{p}=K_1 e^{-ky}+K_2 e^{ky}

Expand to see explanation of ODE solution if required

From the original ODE we need to have a function such that \frac{\partial f}{\partial x}=k f. The simplest function we can propose that satisfies this is an exponential function, in this case \hat{p}=e^{ky}. Except of course, \hat{p}=e^{-ky} works just as well as a solution. And we could have any multiplicative constant out the front of our exponential, and it would still be a solution.

Thus \hat{p}=K_1 e^{-ky}+K_2 e^{ky}, where K_1 and K_2 are arbitrary constants.

Our job is now to determine the values of these multiplicative constants K_1 and K_2. We do this by considering the boundary conditions of the problem.

Boundary Conditions – Farfield

The first, and most intuitively straightforward boundary condition, is that if we are dealing with a small perturbation at the interface between the two fluids, this perturbation needs to decay away to zero as we move infinitely far from the interface. We express this mathematically as:

    \begin{equation*} \hat{p}\rightarrow 0 \quad as \quad y\rightarrow \pm \infty \end{equation*}

Thus far we have largely been able to ignore that our fluid is divided into two discrete regions, separated by an interface between them, but we now need to start considering these two regions individually. Considering the upper fluid (fluid one), as we move to y \rightarrow \infty:

    \begin{equation*} \hat{p}_1 \rightarrow 0 \quad as \quad y\rightarrow \infty \end{equation*}

    \begin{equation*} \hat{p}_1=0=K_1 e^{-k \infty}+K_2 e^{k \infty} \end{equation*}

Thus for our expression for \hat{p}_1, K_2=0 and:

    \begin{equation*} \hat{p}_1=K_1 e^{-ky} \end{equation*}

 

Likewise:

    \begin{equation*} \hat{p}_2 \rightarrow 0 \quad as \quad y\rightarrow -\infty \end{equation*}

    \begin{equation*} \hat{p}_2=0=K_1 e^{k \infty}+K_2 e^{-k \infty} \end{equation*}

Thus for our expression for \hat{p}_2, K_1=0 and

    \begin{equation*} \hat{p}_2=K_2 e^{ky} \end{equation*}

Boundary Conditions – Pressure Matching

For our next boundary condition, we consider the interface between the two fluid streams, at y=0. We have defined the interface between our two fluids as an infinitely thin vortex sheet (a plane over which velocity is discontinuous). Pressure, however, cannot be discontinuous (a shock wave is close, but even then in a real fluid is not a discontinuity), as this would imply infinite acceleration of the fluid across the sheet. Thus our next boundary condition is our pressure matching boundary condition, which is to say:

    \begin{equation*} \hat{p}_1=\hat{p}_2 \quad at \quad y=0 \end{equation*}

    \begin{equation*} K_1 e^{-k 0}=K_2 e^{k 0} \end{equation*}

This has the seemingly trivial result of

    \begin{equation*} K_1=K_2 \end{equation*}

It’s a necessary step, but we aren’t finished yet, we have another boundary condition to go.

 

Boundary Conditions – Kinematic Condition

We still need one more boundary condition, and the one we will impose is called the kinematic condition, so called because it considers the displacement of the interface between the two fluids, rather than the forces acting upon the interface. We introduce a new variable \zeta, which is the displacement of the vortex sheet from its mean position. As this displacement is due to the perturbations we have been considering, the same normal-mode ansatz applies here:

    \begin{equation*} \zeta =\hat{\zeta}e^{i(kx-\omega t)} \end{equation*}

The kinematic condition we will impose is that the displacement defined on either side of the vortex sheet must be identical, i.e. \zeta_1=\zeta_2. We then just need a link between the displacement of the sheet and our other equations of motion.

We find this link by noting that the rate of change of the displacement of the vortex sheet must equal the vertical velocity of the disturbance at either side of the sheet, which can be written as:

    \begin{equation*} \frac{D \zeta}{Dt} =v_{y} \end{equation*}

D here refers to the material derivative, i.e. \frac{D}{Dt}=\frac{\partial}{\partial t}+\mathbf{u} \cdot \nabla

If we recall our linearized y-momentum equation:

    \[ \frac{\partial v'}{\partial t}+\bar{u}\frac{\partial{v'}}{\partial x}=\frac{1}{\rho}(-\frac{dp'}{dy}) \]

 

We can rewrite the term on the LHS in terms of the material derivative:

 

    \[ \frac{D v_{y}}{Dt}=\frac{1}{\rho}(-\frac{dp'}{dy}) \]

So then, we can write our linearized y-momentum equation in terms of the displacement of the vortex sheet:

    \[ \frac{D^2 \zeta}{Dt^2}=\frac{1}{\rho}(-\frac{dp'}{dy}) \]

Now we must work out how to express the \frac{D^2 \zeta}{Dt^2} term:

Click to see expansion of material derivative term

    \begin{equation*} \zeta =\hat{\zeta}e^{i(kx-\omega t)} \end{equation*}

 

    \begin{equation*} \frac{D \zeta}{Dt} =(-i \omega \hat{\zeta}+\bar{u}ik\hat{\zeta})e^{i(kx-\omega t)} \end{equation*}

 

    \begin{equation*} \frac{D^2 \zeta}{Dt^2} =(i^2 \omega^2 \hat{\zeta} -\bar{u}i^2\omega k \hat{\zeta}+\bar{u}^2i^2k^2\hat{\zeta}-\bar{u}i^2\omega k \hat{\zeta})e^{i(kx-\omega t)} \end{equation*}

 

    \begin{equation*} \frac{D^2 \zeta}{Dt^2} =(- \omega^2 +\bar{u} \omega k -\bar{u}^2k^2\+\bar{u}\omega k )\hat{\zeta}e^{i(kx-\omega t)} \end{equation*}

 

    \begin{equation*} \frac{D^2 \zeta}{Dt^2} =(- \omega^2 +2\bar{u} \omega k -\bar{u}^2k^2 )\hat{\zeta}e^{i(kx-\omega t)} \end{equation*}

 

    \begin{equation*} \frac{D^2 \zeta}{Dt^2} =-(\bar{u}k-\omega)^2 \hat{\zeta}e^{i(kx-\omega t)} \end{equation*}

 

We can then substitute this back into our linearized y-momentum equation, dropping the exponential term from both sides.

    \begin{equation*} (\bar{u}k-\omega)^2 \hat{\zeta}=\frac{1}{\rho}\hat{p}_{y} \end{equation*}

We re-arrange this to isolate our displacement term on the LHS:

    \begin{equation*} \hat{\zeta}=\frac{1}{\rho((\bar{u}k-\omega)^2 )}\hat{p}_{y} \end{equation*}

And we can now enforce our condition that displacements on either side of the vortex sheet must match (\hat{\zeta}_1 = \hat{\zeta}_2) in terms of the pressures \hat{p}_1 and \hat{p}_2, which we already have expressions for.

Recalling that:

    \begin{equation*} \hat{p}_1=K_1 e^{-ky} \quad and \quad \hat{p}_2=K_2 e^{ky} \end{equation*}

We can take the derivative of both of these with respect to the y-direction:

    \begin{equation*} \hat{p_{y}}_1=-kK_1 e^{-ky} \quad and \quad \hat{p_{y}}_2=kK_2 e^{ky} \end{equation*}

This allows us now to specify K_1 and K_2, by considering that the vortex sheet is located at y=0:

    \begin{eqnarray*} K_1= -(\bar{u_1}k-\omega)^2 \frac{\hat{\zeta}}{k} \\ \nonumber K_2= (\bar{u_2}k-\omega)^2 \frac{\hat{\zeta}}{k} \\ \nonumber \end{eqnarray*}

We then recall what we determined from matching our pressures at the interface, that K_1=K_2, which means we can state:

    \begin{equation*} -(\bar{u_1}k-\omega)^2 = (\bar{u_2}k-\omega)^2 \end{equation*}

This equation is now in the form of a quadratic eigenvalue problem; our goal here is to find the roots of this quadratic for \omega. 

Solve the eigenvalue problem to determine the stability of the system

The solution of our system is to find the frequency \omega for a given wavenumber k and velocities \bar{u_1} and \bar{u_2}.

Expand to see the details of finding the roots of the quadratic

    \begin{equation*} -(\bar{u_1}^2k^2-2\bar{u_1} k \omega+\omega^2) = (\bar{u_2}^2k^2-2\bar{u_2} k \omega+\omega^2) \end{equation*}

    \begin{equation*} 2 \omega^2 -2\omega k (\bar{u_1}+\bar{u_2}) +k^2(\bar{u_1}^2+\bar{u_2}^2) \end{equation*}

Finding the roots using the quadratic formula:

 

And here, at last, is our result! It’s a relationship between frequency and wavenumber, so we term this a dispersion relation.

    \[ \boxed{\omega = k (\bar{u_1}+\bar{u_2}) \pm \sqrt{k^2( \bar{u_1}^2-\bar{u_2}^2)}} \]

The consequences of this expression are not immediately obvious or intuitive if we’re not used to looking at stability theory, so let’s drill down a bit. Looking at our solution, assuming that we provide a real-valued wavenumber k, we can see that one part of the equation will always produce a real-valued result (k (\bar{u_1}+\bar{u_2})), the other will always produce an imaginary valued result (\sqrt{k^2( \bar{u_1}^2-\bar{u_2}^2)}), so our frequency is going to be a complex number. How can you have a complex frequency? We need to go back to our original discussions of our normal mode ansatz:

    \begin{equation*} q'=\hat{q}(y)e^{i(kx-\omega t)} \end{equation*}

We can separate our time and space components and write the equation as follows:

    \begin{equation*} q'=\hat{q}(y) e^{ikx} e^{-i\omega t} \end{equation*}

Let’s just ignore the e^{ikx} component for now, and focus on the e^{-i\omega t}, as it is \omega we have just determined will be complex. We can again use the properties of exponents to separate, but this time between the real and imaginary components:

    \begin{equation*} e^{-i (\omega_r+\omega_i) t} = e^{-i\omega_r t} e^{-i\omega_i t} \end{equation*}

Given that the ansatz we used had an exponent of i \omega t, the real-valued frequency ends up imaginary, and the imaginary valued frequency ends up real once substituted into our normal mode ansatz! Again, if you are comfortable with the use of complex exponentials to represent growing/decaying periodic functions, skip ahead, otherwise expand this section for some further discussion:

Expand to see a quick exploration of complex frequency

So, our overall behaviour with respect to frequency is the product of the real-valued frequency and imaginary valued frequency. What do they look like in isolation?

Let’s start with a real-valued frequency \omega = 1. Keep in mind once we put this in with our ansatz, what we end up with is y=e^{-it}. That produces a pure oscillation like so:

What about an imaginary valued frequency \omega = -0.1i? We put that into our ansatz and we get y=e^{-0.1t}, a purely decaying wave:

Likewise then, for \omega = +0.1i our ansatz produces y=e^{0.1t}, and we get an exponentially growing wave:

If we then consider complex valued frequencies, such as \omega = 1-0.1i or \omega = 1+0.1i, we get a superposition of this oscillation with decay or growth respectively:

 

 

 

 

 

 

 

So, what does our dispersion relation tell us? For any flow with shear, i.e. any flow where U_1 \neq U_2, for any real wavenumber k (so any physical wave), there will be a pair of complex-valued solutions to our dispersion relation. The real component will give us the frequency of the perturbation, but the critically important part lies in the imaginary valued component: one root is associated with an exponential decay, but the other is associated with exponential growth. So we arrive finally at our critical result:

For an inviscid, incompressible planar mixing layer, in the presence of any shear at all, no matter how small, an infinitely small perturbation will grow infinitely large in time. In other words, all flows of this type are unstable to perturbations at all wavenumbers. This is the result we have been seeking, a mathematical demonstration of why the structures we observe in mixing layers arises.

 

It is worth keeping in mind of course that we have considered a massively simplified problem. Later in the course we will drop some of these simplifications, and consider something that better approximates a real flow. Here our result says that no matter what the wavenumber of the disturbance, it will grow exponentially in time. This turns out to be a consequence of our vortex-sheet approximation; we have assumed that there is a discontinuity in velocity between the two fluids, which does not occur in real flows. Once we use a more realistic velocity profile, we will find that the flow is still very unstable, but very short waves may be damped, rather than growing exponentially. One perhaps surprising result, is that this instability does not depend on viscosity! Even people with some experience in fluid mechanics tend to assume it must be the action of viscous forces that causes this instability to form, but we see here that there is no requirement of viscosity for the instability to arise:

 

To finish with, here is an animation of our result, setting U_1=1, U_2=-1,  and k=2, we can watch our perturbations grow exponentially in time:

 

 

 

 

Now, would we expect a real flow to behave in this way? No, at this point we need to remember once again what we are considering here, which is a linearized problem. The linearization step was that we assumed our fluctuations were small relative to our mean quantities, and clearly that assumption is rapidly violated as the perturbation grows exponentially. We might expect our linear theory to do well in the initial stages, then do increasingly worse as the perturbation grows.

So, this brings our introduction to the KH instability to an end. Let’s recap what we have learned.

  • By linearizing our governing equations and proposing an appropriate ansatz, we can determine whether perturbations in fluid systems will grow or decay in time.
  • The Kelvin-Helmholtz instability renders all mixing-layers unstable to any perturbation for any amount of shear, if the velocity profile is discontinuous.
  • The Kelvin-Helmholtz instability mechanism is an inviscid mechanism.
  • We must remember the limitations of our linear model – it is only intended to predict the initial response of the system to small perturbations.

 

 

If any step of this was unclear, please let me know. The purpose of this section is to be a basic introduction to the material; no question is too basic, no question is “stupid”.