Introduction to hydrodynamic stability

This notebook has been written in Mathematica by 

Mark J. McCready
Professor and Chair of Chemical Engineering
University of Notre Dame
Notre Dame IN 46556
USA

Mark.J.McCready.1@nd.edu
http://www.nd.edu/~mjm/


It is copyrighted to the extent allowed by whatever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy that this notice remain visible to other users. 
There is no charge for copying and dissemination.  

Version:  5/14/99
A more recent version of this notebook may be available at
http :  // www . nd . edu/~mjm/linstab . nb

Summary

This notebook is intended to give a first introduction to hydrodynamic instability.  This includes both the physical concepts and several useful mathematical manipulations.

There are three parts.  The first discusses the concept of linear instability theory and uses a simple wave equation to demonstrate the linearization and calculation of temporal and spatial growth.  

The second part derives the stability relation for a two-layer inviscid flow, the Kelvin-Helmoltz instability.

The third part shows how to derive the basic equation of hydrodynamic stability for Newtonian fluids, the Orr-Sommerfeld equation.

Mathematica aside

In Mathematica, it is convenient to give all expressions a "name".  I try to pick ones that are consistent with what is being done (but sometimes "temp" is used).  This assignment is done with an "=" sign.  To make an equation, a "==" is used.  This distinction is very useful in computer algebra and is employed in all of the packages with which I am familiar .

Input notation

I would enter the dynamic boundary condition from a key pad as

\(dc = 
      \(-\ γ\)\ D[η[x, t], {x, 2}]\  - 
		
        \(ρ\_1\) 
          \((\(-D[ϕ\_1[x, y, t], t]\) - U1\ D[ϕ\_1[x, y, t], x] - 
              g\ η[x, t])\) + 
			
        \(ρ\_2\) 
          \((\(-D[ϕ\_2[x, y, t], t]\) - 
					
              U2\ D[ϕ\_2[x, y, t], x] - g\ η[x, t])\)\)

In Standard form, which is a little shorter, you would need the typeset window to make this practical.

\(dc = 
      \(-\((γ\ ∂\_{x, 2}η[x, t])\)\) - 
        ρ\_1\ 
          \((\(-∂\_t ϕ\_1[x, y, t]\) - 
              U1\ ∂\_x ϕ\_1[x, y, t] - g\ η[x, t])\) + 
        ρ\_2\ 
          \((\(-∂\_t ϕ\_2[x, y, t]\) - 
              U2\ ∂\_x ϕ\_2[x, y, t] - g\ η[x, t])\)\)

Back to start of calculations.

Finally it may be more easily read in Traditional form as shown here.  Note that this particular cell is inactive and will not be evaluated.  You can change this, if you really want to, by selecting the cell, going into the "preferences" under the Edit menu andunde Cell Options, Evaluation Options, set Evaluatable to True.

dc=-(γ (((∂^2) η(x, t))/(∂(x^2))))-(Subscript[ρ, 1]) (-((∂(Subscript[ϕ, 1])[x, y, t])/(∂t))-U1 ((∂(Subscript[ϕ, 1])[x, y, t])/(∂x))-g η(x, t))+(Subscript[ρ, 2]) (-((∂(Subscript[ϕ, 2])[x, y, t])/(∂t))-U2 ((∂(Subscript[ϕ, 2])[x, y, t])/(∂x))-g η(x, t))

Back to start of calculations.

Note that you can convert any "cryptic" input expression by selecting the cell, going up to the cell menu and selecting "convert to traditional form".  This is "shift curly t" on the keypad.

How to use this notebook

The best way to use this notebook is to open it in Mathematica and work through the examples and make changes in parameters or procedural steps to explore these problems.  On line help is available in Mathematica 3 so that a definition and in most cases examples for any unknown command can be obtained.  If you do not have a license for Mathematica, you can download MathReader (http://www.wolfram.com/mathreader/) free of charge.  It does not let you change anything or run the calculations, but it does allow full access to the notebook.  

Other notebooks that cover a range of fluid dynamics problems are available at:
http://www.nd.edu/~mjm/

These can also be used to explore multifluid flows and learn more about using Mathematica.

Motivation: Importance of waves

Waves form as the result of wind on natural bodies of water and contribute to increased gas transfer and droplet production.  Despite the high Reynolds numbers involved, these originate as a hydrodynamic instability that requires a fully-viscous solution of the linearized Navier-Stokes equations to characterize.  You may also observe waves on the windshield of   your car during a rain storm or in the car wash or on the wing of an airplane -- perhaps as the deicing fluid is sheared of.  These examples involve much lower liquid Reynolds numbers and  (with some additional simplifications) can be analyzed with wave equations similar to the example used in the first part below.  

A big area of interest in interfacial waves are those that occur in gas-liquid (i.e., two-phase flow) pipelines and process equipment.  In these situations they greatly increase the pressure drop and interface transport rates, and can lead to flow regime transitions.  Some examples of these kinds of waves are available as still pictures and movies at http://www.nd.edu/~mjm/waves.descrip.html.

The general idea of flow instability

It has been observed in nature, that the steady state solutions for different systems can become unstable to infinitesimal disturbances which should be expected to always be present, (the ground is always vibrating, buildings breath and bend, etc. ...) and possibly because of molecular motions.  A common example is the formation of waves on bodies of water owing to the action of wind.  The "Taylor- Couette Flow" instability is a popular laboratory instabilility that arises due to centrifugal force, and Rayleigh-Benard convection, which arises because of density differences is important both in nature and in laboratories.

Each of these instabilities has a precise, although not necessarily well understood, physical mechanism.  The common feature of an instability is that infinitesimal velocity or density perturbances are amplified (by the base flow or global forces) and thus grow to finite size.  Growth of distrubances could be algebraic or exponential.  Typical analysis (such as those shown below) assume an exponential growth because it is expected that this would overwhelm any algebraic growth.  However, algebraic analyses have been used in some situations where exponential models did not match data.  It is not clear that these have matched any better, but this discussion is beyond the point of this introductory module.

Infinitesimal perturbations are expected to be in the form of noise.  The question is, how to represent this when we want to model instability.  Fortunately, the noise is infinitesimal, which means its amplitude is small compared to any length scale such as its wavelength.  This allows the nonlinear governing equations to be well-approximated by linearized versions.  The linear equations are amenable to a Fourier mode  analysis that can be used to represent any noice signal as a linear combination of independent modes.   
 
If we assume exponential growth (because is the strongest possible and what is observed in nature), and if the growth is in time (which could be how it occurs)  an equation for the amplitude, a, of some disturbance, a is  a = a0 Exp[ ω t], where a0 is the initial amplitude of the disturbance and ω is the temporal growth rate.  Even is a0 is of molecular dimensions (i.e. 10^-10 cm) , and the growth rate is a very reasonable (for common systems),  ω = 1 /s.  We find that it would take only 23 seconds for the disturbance to reach an amplitude of 1 cm!  We thus expect a linearly unstable flow to show some evidence of growing disturbances, unless the residence time is short.

Back to conclusions

Here is the time calculation

N[Log[(1/((1/(10^10))))]]

Note that if the amplitude is growing exponentially, at some point nonlinear processes will become important.  Nonlinear analysis is beyond the scope of this module.

Analysis of instability

To do mathematical analysis of an instability, we need to chose a "basestate" that is the base flow in the absence of an instability.  This could be 0 velocity or it could be a falling film with no waves or a stratified flow with no waves, etc..  

Since we expect that linear equations should govern the initial growth of instability, we will linearize the complete governing equations around the base flow.  This is done by taking the baseflow, say u0 and allowing a small perturbation, say ε.  Thus the complete velocity field would be u = u0 + ε u1.  Note this this is for a one dimensional problem.  For higher dimensions, you would have order ε components in the other directions even if the base flow in that direction was 0.

The magnitude of ε is small (<<1) because it is, for example, the amplitude to wavelength of the noise signal.  

We proceed by substituting u = u0 + ε u1into the governing equations and the boundary conditions and collecting powers of ε.  Because we desire the analysis to be valid for any arbitrary ε, we can separate the system into powers of ε.  
The (ε^0) equations will be the equations for the base state and should be identically 0.  
The (ε^1) equations, should give the behavior of the very small amplitude disturbance and will contain u0, which we know and u1, the distrubance that we wish to study.  These equations are necessarily linear (we just linearized them with this procedure).  Any higher powers of ε will be ignored and saved for when we want to do nonlinear analysis.

To determine the response of the equation to an arbitrary noise signal, we choose a mode that represents the kind of disturbance that we expect to see.  If the domain is fixed and there is no flow through, (e.g., a solid beam), we might expect to use fixed spatially periodic modes that grow in time.  For waves on water, we would use (traveling) spatially and temporally periodic disturbances that could grow in space and/or time.     

Since the system is linear, we can examine the response of any separate mode without worrying about the effect of other modes. This linearity allows us to decompose an arbitrary disturbance into a sum of Fourier modes; each mode will ultimately satisfy the equations and boundary conditions.  By scanning the entire frequency or wavenumber range, we can be sure that we understand the effect of any initial (infintesimal) disturbance.

Back to conclusions

Example: Waves on a falling liquid film

Problem set up

For a system of traveling waves on a liquid film, the disturbance might be conveniently represented by:

 u1= a Exp[i ( k x - ω t) ] ,
 
which is periodic in time with a circular frequency, ω [rad/s] and travels in the positive x direction with a speed, Re[ω/k].  Since we are considering that this disturbance might grow (the whole point of this exercise) we will allow for the possibility that k or ω is complex.  Then this disturbance could growth spatially, Im [k] <0 or temporally, Im[ω]>0 depending on the circumstance of the flow. 

For a simple, but nontrivial, example, we choose a (generalized) form of the Kuramoto-Sivshinsky equation which is derived for a falling film of water by (e.g., H. -C. Chang, Chemical Eng. Sci. 1986; Alekseenko et al., AIChE J, 1985) as:

∂(h/∂)t  + c0 ∂(h/∂)x  + α h ∂(h/∂)x + β  (((d^2) h)/∂)(x^2) + γ (((d^3) h)/∂)(x^3) + σ (((d^4) h)/∂)(x^4) == 0

where h is the deviation of the surface height from some steady value, t is time and x is the flow direction.
This is a "wave" equation in which all of the complexities of the flow have been reduced to a  single equation for one dependent variable, the liquid height, h.

In this equation, all of the variables are real.

Clearly a constant value of h is a solution to this equation.  We will choose
h=h0 and begin by linearizing the equation around h =h0.

Back to conclusions

The linearized equation

We choose, h = h0 + ε h1and substitute into the wave equation.  

This first cell is inactive and is formatted in "traditional form" to show the substitution clearly.  Traditional form is not a good input form for Mathematica

temp1=((∂(h0 + ε h1(x, t)))/(∂t))+c0 ((∂(h0 + ε h1(x, t)))/(∂x))+β (((∂^2)(h0 + ε h1(x, t)))/(∂(x^2)))+γ (((∂^3)(h0 + ε h1(x, t)))/(∂(x^3)))+σ (((∂^4)(h0 + ε h1(x, t)))/(∂(x^4)))+α ((∂(h0 + ε h1(x, t)))/(∂x)) (h0 + ε h1(x, t))

Here is the start of the calculation where we use the conventional "standard form"

\(temp1 = 
      ∂\_t\((h0 + ε\ h1[x, t])\) + 
        c0\ ∂\_x\((h0 + ε\ h1[x, t])\) + 
        β\ ∂\_{x, 2}\((h0 + ε\ h1[x, t])\) + 
        γ\ ∂\_{x, 3}\((h0 + ε\ h1[x, t])\) + 
        σ\ ∂\_{x, 4}\((h0 + ε\ h1[x, t])\) + 
        α\ ∂\_x\((h0 + ε\ h1[x, t])\)\ 
          \((h0 + ε\ h1[x, t])\)\)

Get the separate powers of ε

temp2 = Collect[Expand[temp1], ε]

The powers of ε give:

Coefficient[temp2, ε, 0]

Which is the base state solution; any constant value of h0 is a solution.

stabeq = Coefficient[temp2, (ε^1)]

Which is the linearized equation.  Note that for a base state of h0, the nonlinear term contributes h0 α ∂(h1/∂)x.

Back to conclusions

This is the nonlinear term that is neglected.

Coefficient[temp2, (ε^2)]

Traveling wave mode expansion

Now we substitute a traveling wave disturbance into the linearized equation.  Note that the substitution is done in a way that works for any arbitrary derivative.

temp3=stabeq/.{Derivative[a1_, a2_][h1][x, t]:>(∂_({x, a1}, {t, a2}))(
                a Exp[I (k x - ω t)])}

We now divide out the original disturbance, which is possible since we used an exponential form.

temp4=Cancel[Apart[Expand[(temp3/(a Exp[I (k x-ω t)]))]]]

Back to conclusions

The result is the complete dispersion relation for waves in this system.  The relation between speed and frequency and wavenumber is given in this relation as is the linear growth.

Temporal growth

Consider first, spatially periodic disturbances that grow in time. This is the easy case (only 1 time derivative and it is a first derivative) and the one that most people naturally do.  Of course waves on a falling film seem to grow with distance (!) 
 For this case, k is real and ω is complex.

We solve the above equation for ω and then get a relation for ω.  (It takes this format to do it)

omega = 
      ω /. 
        Solve[temp4 == 0, ω][[1
          ]]

Expand[omega]

At this point we can see that certain terms contribute to the imaginary part of ω and others only to the real part.  Let's look at the real part which is the frequency and relates to the wave speed.

The wave speed is ω/k so

speedx = Cancel[Apart[(omega/k)]]

We might as well use the canned program for Real and Imaginary variables.

Needs["Algebra`ReIm`"]

We need some definitions of this form to make Mathematica know which variables are real.

x /:  Im[x]  =  0;
t /:  Im[t]  =  0;
h /:  Im[h]  =  0;
h0 /:  Im[h0]  =  0;
α /:  Im[α]  =  0;
β /:  Im[β]  =  0;
c0 /:  Im[c0]  =  0;
k/: Im[k] = 0;
γ /:  Im[γ]  =  0;
σ /:  Im[σ]  =  0;

We can get the real part by:

realspeed  =  Re[speedx]

Plot[
      realspeed /. {h0 -> 1.5, α -> 1, c0 -> 0, 
          β -> 3, γ -> -3, σ -> 4}, {
        k, 0, 1},  AxesLabel -> {"k", " growthrate"}];

We see that the wave speed is given by a base parameter, c0, and some kind of speed change due to the film thickness, h0 α.  The wave speed is dependent upon the wavenumber and this is given by the term, -γ (k^2).

The imaginary part of omega gives the growth rate

Im[omega]

For β and σ positive, which they should be (σ is a surface tension parameter), we have the first term trying to cause the instability, (fluid inertia) and the second one countering it.

It is instructive to plot the growth rate.  We see that at low wavenumber (long waves) the growth rate is increasing.  However at higher wavenumber the stabilization becomes stronger and eventually the waves become stable .

Plot[
      Im[omega /. {h0 -> 1.5, α -> 1, c0 -> 0, 
            β -> 3, γ -> -3, 
            σ -> 4}], {k, 0, 1},  
      AxesLabel -> {"k", " growthrate"}];

So that there is a region for k < ~ .85 where the waves are unstable.  The fastest growing wave is at k =~ 0.6.  

Now for a linear system, we can never produce a mode that is not initially excited.  Thus the only waves that grow are those that were initially present.  If the initial disturbance is noise, then all modes are present and we would expect to see the fastest growing mode dominate.  However, this will not usually be a pure mode but a band of waves in the vicinity of the peak. 

If the initial disturbance does not include modes near the peak, then others may dominate.

Spatial growth

We could have also considered spatially growing waves, which are most likely to be the physically relevant case for a falling film, by defining ω as real and then allowing for complex k.  Even for this simple equation, the analysis is more complicated because the spatial derivatives are of higher order than the time derivatives.

temp5 = Expand[temp4]

kays=Solve[temp5 == 0,k];

Mathematica knows how to solve quartic equations analytically, so it can get the 4 roots.

Here is the first one:

Table[
      N[kays[[1]] /. {h0 -> 1.5, 
            α -> 1, c0 -> 0, β -> 3, 
            γ -> -3, σ -> 4}], {ω, 0, 
        5,  .5}]

This root, [[1]], corresponds to a wave that is going upstream.  This would not happen.

Here is the second one:

Table[
      N[kays[[2]] /. {h0 -> 1.5, 
            α -> 1, c0 -> 0, β -> 3, 
            γ -> -3, σ -> 4}], {ω, 0, 
        5,  .5}]

This roo,t [[2]], also corresponds to a wave that is going upstream.  This also would not happen.

Here is the third one:

Table[
      N[kays[[3]] /. {h0 -> 1.5, 
            α -> 1, c0 -> 0, β -> 3, 
            γ -> -3, σ -> 4}], {ω, 0, 20, 
        1}]

This root, [[3]], seems to be more and more unstable as the frequency increases, this is not physical for this system.  However, it travels down stream and so we cannot ignore it.  When an "unphysical mode" appears, this is cause for concern that the derivation has neglected something important or that our solution procedure has flaws in it.  In this case, I think that the equation, which has been "generalized" to make it more interesting to study, is being used with parameters that are outside where it behaves well.

A more interesting consequence of this analysis is that you might find it very difficult to integrate this equation numerically!  

To proceed, consider both, [[3]], and [[4]].

Let's plot the two parts to see what happens.
Here we extract the 3rd root.

k3 = k /. kays[[3]];

Short[k3]

We can plot the real wave speed as  ω/k.

\(\(speed3 = 
      Plot[ω\ /
          Re[N[k3 /. {h0 -> 1.5, α -> 1, c0 -> 0, 
                  β -> 2, γ -> \(-3\), 
                  σ -> 8}]], 
	{ω,  .001, 6}, 
        PlotStyle -> \ Dashing[{ .04,  .02}]]; \)\)

The spatial growth rate is -Im[k].

\(\(growth3 = 
      Plot[\(-\ Im[
            N[k3 /. {h0 -> 1.5, α -> 1, c0 -> 0, 
                  β -> 2, γ -> \(-3\), 
                  σ -> 8}]]\), 
	{ω,  .001, 6}, 
        PlotStyle -> \ Dashing[{ .04,  .02}]]; \)\)

Do this again for the 4th root.

k4 = k /. kays[[4]];

First the speed

\(\(speed4 = 
      Plot[ω\ /
          Re[N[k4 /. {h0 -> 1.5, α -> 1, c0 -> 0, 
                  β -> 2, γ -> \(-3\), 
                  σ -> 8}]], 
	{ω,  .001, 6}]; \)\)

Then the growth

\(\(growth4 = 
      Plot[\(-\ Im[
            N[k4 /. {h0 -> 1.5, α -> 1, c0 -> 0, 
                  β -> 2, γ -> \(-3\), 
                  σ -> 8}]]\), 
	{ω,  .001, 6}]; \)\)

Here we plot the two wave speeds:

Show[speed3, speed4, 
      AxesLabel -> {"ω", "wave speed"}];

Here we plot the two growth curves.

twogrowth = 
      Show[growth3, growth4, AxesLabel -> {"ω", "growth"}];

The lower mode, which starts out dashing and then becomes solid, is behaving more like a real wave than the top one.  We do not see how a wave can have unbounded growth for very short waves when surface tension is present.  The speed is better also as it more closely matches the qualitative behavior of the temporal analysis.  

We will leave the issue of the "unphysical" mode because we can do nothing else with it right now (Although you could vary some of the parameters to see how it behaves.) but in a real problem, you would have to find out if arises because of the limitations of the derviation or the solution -- or if the instability is really there!

Relation between spatial and temporal growth

The spatial and temporal growth rates can be related by a Gaster Transformation (M. Gaster, J. Fluid Mech, 14, p222, 1962).  This is done by using the group velocity, 
∂ω/∂k, to convert temporal to spatial growth, 

spatial growth  =  temporal growth/ (∂ω/∂k)

Here we go with this calculation.

omega

The group velocity is defined to be real.

groupvel = Re[(∂_k) omega]

Here is a plot that shows the spatial growth obtained from the temporal growth.

\(\(pseudospatial\  = 
      ParametricPlot[{
          Re[omega /. {h0 -> 1.5, α -> 1, c0 -> 0, 
                β -> 2, γ -> \(-3\), 
                σ -> 8}], 
          \((Im[omega]/groupvel)\) /. {h0 -> 1.5, α -> 1, 
              c0 -> 0, β -> 2, γ -> \(-3\), 
              σ -> 8}}, 
	{k, 0, 1}, 
        AxesLabel -> {"frequency", "converted spatial growth"}]; 
    \)\)

Now we show the real spatial and the converted spatial growths. They agree well for positive growth but not well as the amplitude of the negative growth gets larger.

Generally, it is possible to make this conversion for regions close to neutral stability.

\(\(Show[twogrowth, pseudospatial, 
      PlotRange -> 
		{{0, 4}, {\(- .25\), \  .1}}]; \)\)

Back to conclusions

Why can we use a complex representation for a real disturbance ?

Back to conclusions

An important question could be that the disturbance that we used was not a real number.   To make the disturbance a real number, we can use a0 Exp[I(k x - ω t)] + c.c. or
a0 Exp[I(kx - ωt)] + a0cc Exp[-I(kx - ω t)].  If this is substituted into the equation we will get twice as many terms.  One set will have Exp[I(kx - ωt)] as a factor and the other set will have Exp[-I(kx - ωt)] as a factor.  For this expression to always be true (i.e. for any x and t ), these sets must each be equal to 0 separately.  The analysis would then be identical to the one done above.  Note that the appearance of an i in an equation denotes a π/2 phase shift (i.e., Sine versus Cosine).  

Note that if we are doing a nonlinear problem, we would have to keep track of the complex conjugate part as it may contain additional information.

Kelvin-Helmholtz Instability for finite depth

Problem set up

A good reference for this section is R. L. Panton, Incompressible flow, Wiley, 1984

The generation of water waves by wind has puzzled and fascinated scientists for centuries.  It is a classic stability problem that, as it turns out, was not really completely understood until about the 1970's.  

This problem was attacked very early in the study of hydrodyanmic stability for the idealized case of inviscid flows.  This solutions is known as the Kelvin-Helmholtz instability and it does not predict the observed behavior of water wave generation!
It does not work for any two layer system where the viscosities of the two fluids are very different, no matter how high the Reynolds numbers are.  Still it is a useful instability for study as it provides a nice physical mechanism to an easily understood situation.

Back to conclusions

Consider a stratified horizontal flow of two inviscid fluids.  The top fluid is "2" and has a velocity of U2.   The bottom is "1" and has a velocity of U1.  The interface is at y = 0 and the top wall is at y = h2, the bottom wall is at y = -h1.    
 
 Since each phase is inviscid, Laplace's equation for the velocity potential can be written,
 
 ∇2Φ1 = 0  	and	 	∇2Φ2 = 0.  
 
 where
 
 U1 = ∇Φ1 	 and 		U2 = ∇Φ2.

Laplace equation solution

First we should find solutions for Φ in each phase that fits the necessary boundary conditions.  This would be that there is no flow through the top and bottom walls.  
The other boundary conditions will be introduced along the way.
 
The waves will be traveling in the x direction, grow and oscillate in time and must decay away from the interface in y.    

For the lower phase, why don't we guess a form as

Φ1[x,y,t] = U1 x + ϕ1[y] Exp[I k(x - c t)]

ϕ1[y] = B1  Cosh[k (y+h1)]

this gives the base flow plus the disturbance.   

For the upper phase, how about

Φ2[x,y,t] = U2 x + ϕ2[y] Exp[I k(x - c t)]

ϕ2[y] = A2  Cosh[k (y-h2)]

In these expressions, k is wavenumber, c is wave speed, x is the direction of travel and t is time.  

Let's see how well we have guessed!

For the upper phase

(Subscript[Φ, 1]) = 
      U1 x + B1 Cosh[k (h1 + y)] Exp[I k (x - c t)]

Check in Laplace's equation and it works

(first see it in traditional form)

(((∂^2) (Subscript[Φ, 1]))/(∂(x^2)))+(((∂^2) (Subscript[Φ, 1]))/(∂(y^2)))

\(∂\_{x, 2}Φ\_1 + 
      ∂\_{y, 2}Φ\_1\)

Back to conclusions

Now check the lower boundary condition, it works also

(∂_y) (Subscript[Φ, 1]) /. y -> -h1

now try the upper phase.

(Subscript[Φ, 2]) = 
      U2 x + A2 Cosh[k (y - h2)] Exp[I k (x - c t)]

test in LaPlaces' equation,

\(∂\_{x, 2}Φ\_2 + 
      ∂\_{y, 2}Φ\_2\)

Test the boundary condition,

(∂_y) (Subscript[Φ, 2]) /. y -> h2

OK so far, so good.

Kinematic boundary condition

Whenever there is an interface between fluids that has waves on it, we need a relation between the velocity field and the motion of the interface.  This can formulated by using the substantial derivative of the surface position function, η(x,t) and equating this to the normal velocity of the fluid at the interface.

((Dη (x, t))/Dt)  = normal velocity at the interface = ((∂ ϕ (x, y, t))/∂)y  (note that this potential is just for the perturbation velocity because there is no average normal velocity)

this becomes

∂η/∂t + U ∂η/∂x + W ∂η/∂z = ∂ ϕ/∂y  @ y= η. 

In this example we will not consider the transverse direction,  so W = 0 and ∂ ()/∂z = 0.  In its present form, the second term on the left is nonlinear.  We just need a linearized version for this linear stability problem.  This is

 ∂η/∂t + U1 ∂η/∂x  = ∂ ϕ/∂y     @ y= η,  

because the perturbation term,  ∂ϕ/∂x ∂η/∂x is order a2, which can be ignored.

There is still one problem remaining.  It would be exceedingly inconvenient to evaluate this condition at y = η, since this changes with x and t.  It turns out that this problem can be solved by domain perturbation.  We would like to evaluate these equations at y = 0 and we need the real values of the variables at this point.  We can get them by expanding the velocities in a Taylor series.  The surface variable is not a function of y and thus does not require any expanding. 

U(y=h) = U(y=0) + a ∂U/∂y(y=0) + H.O.T.

v(y=h)  =  v(y=0) + a ∂v/∂y (y=0) + H.O.T 

The correction terms are of order a and would not enter in this linear problem.  Further, for this inviscid flow, there is no velocity gradient in the y direction so there is no correction for the flow direction velocity.  

As a result of this discussion, we have this "kinematic" condition for both phases related to the interfacial shape.     

 ∂η/∂t + U1 ∂η/∂x  = ∂ ϕ1/∂y     @ y= 0.   
 
  ∂η/∂t + U2 ∂η/∂x  = ∂ ϕ2/∂y     @ y= 0.

Back to conclusions

Here is the kinematic condition, first in traditional form,

kc1=-((∂η(x, t))/(∂t))-U1 ((∂η(x, t))/(∂x))+((∂(Subscript[ϕ, 1])(x, y, t))/(∂y))

kc1 = 
      -(∂_t) η[x, t] - U1 (∂_x) η[x, t] + 
        (∂_y) (Subscript[ϕ, 1])[x, y, t]

Here we impose the kinematic condition for phase 1 by  substituting for ϕ1 and η.  Note that \(\(η\&^\)\) is just the wave amplitude.

temp1=kc1/.{\(η[x, t] -> η\&^\ Exp[I\ k\ \((x - c\ t)\)]\),(Subscript[ϕ, 1])[x, y, t] -> 
              B1 Cosh[k (h1 + y)] Exp[I k (x - c t)],(Subscript[ϕ, 2])[x, y, t] -> 
              A2 Cosh[k (y - h2)] Exp[I k (x - c t)],Derivative[a1_, a2_][η][x, t]:>\(∂\_\({x, a1}, {t, a2}\)\((
                  η\&^\ Exp[I\ k\ \((x - c\ t)\)])\)\),(Subscript[ϕ, 1]^((a1_, a2_, a3_)))[x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  B1 Cosh[k (h1 + y)] Exp[I k (x - c t)])
                  ,(Subscript[ϕ, 2]^((a1_, a2_, a3_)))[x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  A2 Cosh[k (y - h2)] Exp[I k (x - c t)])
                  }

temp2 = 
      Cancel[Expand[(temp1/Exp)[I k (x - c t)]]] /. y -> 0

BEE1 = 
      B1 /. Solve[temp2 == 0, B1][[1
          ]]

So now we know B1 in terms of \(\(η\&^\)\). 
Now do the second kinematic condition

Here it is is in traditional form

kc2=-((∂η(x, t))/(∂t))-U2 ((∂η(x, t))/(∂x))+((∂(Subscript[ϕ, 2])(x, y, t))/(∂y))

kc2 = 
      -(∂_t) η[x, t] - U2 (∂_x) η[x, t] + 
        (∂_y) (Subscript[ϕ, 2])[x, y, t]

temp3=kc2/.{\(η[x, t] -> η\&^\ Exp[I\ k\ \((x - c\ t)\)]\),(Subscript[ϕ, 1])[x, y, t] -> 
              B1 Cosh[k (h1 + y)] Exp[I k (x - c t)],(Subscript[ϕ, 2])[x, y, t] -> 
              A2 Cosh[k (y - h2)] Exp[I k (x - c t)],Derivative[a1_, a2_][η][x, t]:>\(∂\_\({x, a1}, {t, a2}\)\((
                  η\&^\ Exp[I\ k\ \((x - c\ t)\)])\)\),(Subscript[ϕ, 1]^((a1_, a2_, a3_)))[x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  B1 Cosh[k (h1 + y)] Exp[I k (x - c t)])
                  ,(Subscript[ϕ, 2]^((a1_, a2_, a3_)))[x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  A2 Cosh[k (y - h2)] Exp[I k (x - c t)])
                  }

temp4 = 
      Cancel[Expand[(temp3/Exp)[I k (x - c t)]]] /. y -> 0

Now we can get A2

AY2 = 
      A2 /. Solve[% == 0, A2][[1
          ]]

Dynamic Boundary condition

The pressure will not be constant along the interface, if waves are present.  We can find out the variation by applying the Bernoulli equation for each phase at the interface.  The pressure must be the same on each side of the interface, except for the effect of durface tension.  This is a "dynamic condition" (pressure changes with velocity) and is referred to in this way.  

We can write the Bernoulli equation for each phase as:

∂ϕ1/∂t + 1/2 (∇ϕ)2 + P1/ρ1 + g η(x,t)  = 0  

∂ϕ2/∂t + 1/2 (∇ϕ)2 + P2/ρ2 + g η(x,t)  = 0      

We will linearize the (∇ϕ)2 term.  The difference in the pressure can be written as
(P2-P1) = - γ (((∂^2) η (x, t))/(∂x2))

where γ is the surface tension coefficient.

Here is the dynamic boundary condition.

First see it in traditional form

dc=-(γ (((∂^2) η(x, t))/(∂(x^2))))-(Subscript[ρ, 1]) (-((∂(Subscript[ϕ, 1])[x, y, t])/(∂t))-U1 ((∂(Subscript[ϕ, 1])[x, y, t])/(∂x))-g η(x, t))+(Subscript[ρ, 2]) (-((∂(Subscript[ϕ, 2])[x, y, t])/(∂t))-U2 ((∂(Subscript[ϕ, 2])[x, y, t])/(∂x))-g η(x, t))

\(dc = 
      \(-\((γ\ ∂\_{x, 2}η[x, t])\)\) - 
        ρ\_1\ 
          \((\(-∂\_t ϕ\_1[x, y, t]\) - 
              U1\ ∂\_x ϕ\_1[x, y, t] - g\ η[x, t])\) + 
        ρ\_2\ 
          \((\(-∂\_t ϕ\_2[x, y, t]\) - 
              U2\ ∂\_x ϕ\_2[x, y, t] - g\ η[x, t])\)\)

Substitute all that we know

temp6=dc/.{\(η[x, t] -> η\&^\ Exp[I\ k\ \((x - c\ t)\)]\),(Subscript[ϕ, 1])[x, y, t] -> 
              B1 Cosh[k (h1 + y)] Exp[I k (x - c t)],(Subscript[ϕ, 2])[x, y, t] -> 
              A2 Cosh[k (y - h2)] Exp[I k (x - c t)],Derivative[a1_, a2_][η][x, t]:>\(∂\_\({x, a1}, {t, a2}\)\((
                  η\&^\ Exp[I\ k\ \((x - c\ t)\)])\)\),(Subscript[ϕ, 1]^((a1_, a2_, a3_)))[x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  B1 Cosh[k (h1 + y)] Exp[I k (x - c t)])
                  ,(Subscript[ϕ, 2]^((a1_, a2_, a3_)))[x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  A2 Cosh[k (y - h2)] Exp[I k (x - c t)])
                  }

temp7 = 
      Cancel[Expand[(temp6/Exp)[I k (x - c t)]]] /. y -> 0

Substitute for A2 and B1

This command substitutes for A2 from above

temp8 = temp7 /. A2 -> AY2

This command substitutes for B1 from above

temp9 = temp8 /. B1 -> BEE1

temp10 = Expand[temp9]

Dispersion relation and growth expression

By substituting, we have obtained a solution for c in terms of known variables.  Note that \(\(η\&^\)\), which is arbitrary, is multiplying every term linearly and will cancel.

temp11 = Collect[temp10, c]

This can be written more concisely as

temp12 =  FullSimplify[%]

This is the dispersion relation for the wave speed.  
The solutions are:

cees = Solve[temp12 == 0, c];

Short[cees]

Which should be a set of upstream and downstream traveling waves.  
c1 will be the upstream wave and will not appear in a flow system
c2 will be the downstream wave that we would like to study the behavior of

c1 = c /. cees[[1]]

c2 = c /. cees[[2]]

Example calculations

We are using cgs units for the air-water system.  The wind speed is 300 cm/s and the layers are deep.   The system is stable.

c2numb = 
      c2 /. {(Subscript[ρ, 2]) ->  .00125, (Subscript[ρ, 1]) -> 1.020, 
          g -> 980, γ -> 74, U1 -> 1, U2 -> 300, 
          h1 -> 10, h2 -> 10}

kh1 = 
      Plot[Re[N[c2numb]], {k,  .01, 30}, 
        AxesLabel -> {"wavenumber", "wavespeed"}];

This is the typical shape of a water wave dispersion curve.  Gravity waves (long waves) get faster as they get longer.  In contrast, capillary waves get faster as they get shorter.  (These waves are linearly stable.)

Here is the growth curve which is always 0.  This means that the theory predicts that no wave modes will grow.  Note that for real fluids, stable waves have negative growth rates.

kh2 = 
      Plot[Im[N[c2numb]], {k,  .01, 30}, 
        AxesLabel -> {"wavenumber", "wave growth"}];

Now suppose the depth is not as large

c2numb2 = 
      c2 /. {(Subscript[ρ, 2]) ->  .00125, (Subscript[ρ, 1]) -> 1.020, 
          g -> 980, γ -> 74, U1 -> 1, U2 -> 300, 
          h1 ->  .2, h2 ->  .2}

kh3 = 
      Plot[Re[N[c2numb2]], {k,  .01, 30}, 
        AxesLabel -> {"wavenumber", "wavespeed"}];

We might as well compare the wave speeds for two different depths.  If the depth is too small, the long waves do not speed up.

Show[kh1, kh3];

Now the growth curve which is again always 0.

kh4 = 
      Plot[Im[N[c2numb2]], {k,  .01, 30}, 
        AxesLabel -> {"wavenumber", "wave growth"}];

Now, let's turn up the gas speed

c2numb3 = 
      c2 /. {(Subscript[ρ, 2]) ->  .00125, (Subscript[ρ, 1]) -> 1.020, 
          g -> 980, γ -> 74, U1 -> 1, U2 -> 600, 
          h1 ->  .2, h2 ->  .2}

\(\(kh5 = 
      Plot[Re[N[c2numb3]], {k,  .01, 30}, 
        AxesLabel -> {"wavenumber", "wave speed"}, 
		
        PlotRange -> All]; \)\)

We see that the behavior of the speed in the low wavenumber region is different.  This is what the model predicts when the waves become unstable.  This does not seem correct and is not predicted by more accurate models.

Now the growth curve is not always 0.  Wavenumbers less than about 5/cm are predicted to be unstable.

kh6 = 
      Plot[Im[N[c2numb3]], {k,  .01, 30}, 
        AxesLabel -> {"wavenumber", "wave growth"}];

You would find that the velocity difference will need to be about 660 cm/s for large depths but less if the depths are lowered and the speed is changed such as the example above.

Now make both layers deep.  We need to make the velocity higher.

c2numb4 = 
      c2 /. {(Subscript[ρ, 2]) ->  .00125, (Subscript[ρ, 1]) -> 1.020, 
          g -> 980, γ -> 74, U1 -> 1, U2 -> 670, 
          h1 -> 10, h2 -> 10}

kh7 = 
      Plot[Re[N[c2numb4]], {k,  .01, 30}, 
        AxesLabel -> {"wavenumber", "wave speed"}];

We again see that the behavior of the speed in the wavenumber region where waves are unstable is discontinuous and then constant.  This is not correct as discussed below.

Now the growth curve is bounded away from 0 wavenumber.

kh8 = 
      Plot[Im[N[c2numb4]], {k,  .01, 30}, 
        AxesLabel -> {"wavenumber", "wave growth"}];

Which shows the temporal grow rate for an interfacial disturbance.  
This calculation is easy and neat, however there is one big problem, it does not match the linear behavior for the air water system.  The predicted point of instability is too high; the real value would be more like 3 m/s.  Furthermore, note the speed of unstable waves.  They are predicted to travel at the speed ((Subscript[ρ, 1]) u1+(Subscript[ρ, 2]) u2)/((Subscript[ρ, 1])+(Subscript[ρ, 2])), which is not correct.

Now check the upstream wave:

c1numb1 = 
      c1 /. {(Subscript[ρ, 2]) ->  .00125, (Subscript[ρ, 1]) -> 1.020, 
          g -> 980, γ -> 74, U1 -> 1, U2 -> 670, 
          h1 -> 10, h2 -> 10}

kh9 = 
      Plot[Re[c1numb1], {k,  .01, 20}, 
        AxesLabel -> {"wavenumber", "wave speed"}];

The wave speed is now usually negative

kh10 = 
      Plot[Im[c1numb1], {k,  .01, 20}, 
        AxesLabel -> {"wavenumber", "wave growth"}];

The upstream wave does have a negative growth region that corresponds to the wavenumbers where the downstream wave is predicted to be unstable.

Examination of the physics

No surface tension

It is worthwhile to get some additional physical feel for this instability, even if it is not applicable to a real air-water system.  First, what if  there is no surface tension.

c2numb5 = 
      c2 /. {(Subscript[ρ, 2]) ->  .00125, (Subscript[ρ, 1]) -> 1.020, 
          g -> 980, γ -> 0, U1 -> 1, U2 -> 670, 
          h1 -> 10, h2 -> 10}

kh11 = 
      Plot[Re[c2numb5], {k,  .01, 20}, 
        AxesLabel -> {"wavenumber", "wave speed"}, 
        PlotRange -> All];

The wave speed looks peculiar.  Without surface tension to stabilize the short waves, these are unstable.  If they are unstable, the wave speed is just the density weighted average of the upper and lower depths.

kh12 = 
      Plot[Im[c2numb5], {k,  .01, 20}, 
        AxesLabel -> {"wavenumber", "wave growth"}];

Without surface tension, all waves shorter than some cutoff are unstable.

So surface tension stabilizes the short waves, which is reasonable.

What is the point of neutral stability and what is the mechanism?

Let's consider the simplest case, both layers deep, air-water ( (Subscript[ρ, 2]) <<< (Subscript[ρ, 1])) and that the velocity, U1 is much less than U2.

nuet1 = c2 /. {h2 -> h1, U1 -> 0}

We know (that is I know from looking at the expression that this can be done), by definition, for h1->Infinity ...

nuet2 = nuet1 /. {Coth[h1 k] -> 1}

For instability, we need the expression inside the root to be negative.  Let's look at it.

nuet3 = Apart[nuet2]

Grab the part inside the radical,

nuet4 =  Numerator[nuet3[[2]]][[1]]

Do a little more manipulation on it,

nuet5  = 
      Collect[Collect[ExpandAll[nuet4/(Subscript[ρ, 1])/4], (Subscript[ρ, 2])], k]

Now we see that for each power of k, how important the (Subscript[ρ, 2])terms are.  We choose the important ones by making the others 0.  This is inelegant, but the easiest way to do it.

neutral1  = 
      nuet5 /. { γ (Subscript[ρ, 2]) -> 0,   ((Subscript[ρ, 2])^2) -> 0}

If this quantity is greater than 0, the flow is unstable.  The first term is the restoring force of surface tension, the third term is the restoring force of gravity.  The second term is the action of the gas flow which is causing deformation because of the low pressure region at the wave crest and high pressure regions in the wave trough.

Thus the mechanism for instability is the pressure caused by a Bernoulli effect becoming stronger than the restoring forces of gravity and / or surface tension. 

For flows involving real fluids this effect is definitely present.  However, this mechanism is not as efficient as another mechanism (which involves a shift in the pressure minimum, see Jurman, Deutsch and McCready, JFM, 238, 178-219, 1992) and hence is not responsible for the initial formation of waves in air-water flows and other sytems where there is a significant viscosity difference between the phases.  However this general mechanism may be important for growth of large amplitude waves.

Orr-Sommerfeld Equation

A good reference for this section is R. L. Panton, Incompressible flow, Wiley, 1984

Here we derive the Orr-Sommerfeld equation which is a 4th order ODE that describes the growth on infinitesimal periodic distrubances that are governed by the Navier-Stokes equations.

The linearized x-momentum equations for a nearly-parallel flow are
x direction (first in traditional form)

xmom=((∂p(x, y, t))/(∂x))-(((((∂^2) u(x, y, t))/(∂(x^2)))+(((∂^2) u(x, y, t))/(∂(y^2))))/Re)+((∂u(x, y, t))/(∂t))+((∂u(x, y, t))/(∂x)) u0(y)+((∂u0(y))/(∂y)) v(x, y, t)

\(xmom = 
      ∂\_x p[x, y, t] - 
        \(∂\_{x, 2}u[x, y, t] + 
            ∂\_{y, 2}u[x, y, t]\)\/Re + ∂\_t u[x, y, t] + 
        ∂\_x u[x, y, t]\ u0[y] + 
        ∂\_y u0[y]\ v[x, y, t]\)

Back to conclusions

y direction (first in traditional form)

ymom=((∂p(x, y, t))/(∂y))-(((((∂^2) v(x, y, t))/(∂(x^2)))+(((∂^2) v(x, y, t))/(∂(y^2))))/Re)+((∂v(x, y, t))/(∂t))+((∂v(x, y, t))/(∂x)) u0(y)

\(ymom = 
      ∂\_y p[x, y, t] - 
        \(∂\_{x, 2}v[x, y, t] + 
            ∂\_{y, 2}v[x, y, t]\)\/Re + ∂\_t v[x, y, t] + 
        ∂\_x v[x, y, t]\ u0[y]\)

The continuity equation is

cont=((∂u(x, y, t))/(∂x))+((∂v(x, y, t))/(∂y))

cont = (∂_x) u[x, y, t] + (∂_y) v[x, y, t]

We need to expand these in terms of normal modes where the variable will be assumed to have the form ξ[x,y,t]=\(ξ\&^\)[y] Exp[I(α x -α c t)] based on the idea that we need functions that can describe spatially and temporally varying disturbances

Back to conclusions

os1=xmom/.{u[x, y, t] -> 
              uh[y] Exp[I (α x - α c t)],v[x, y, t] -> 
              vh[y] Exp[I (α x - α c t)],p[x, y, t] -> 
              ph[y] Exp[I (α x - α c t)],Derivative[a1_, a2_, a3_][u][x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  uh[y] Exp[I (α x - α c t)]),Derivative[a1_, a2_, a3_][v][x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  vh[y] Exp[I (α x - α c t)]),Derivative[a1_, a2_, a3_][p][x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  ph[y] Exp[I (α x - α c t)])}

os2 = 
      Cancel[Expand[(os1/Exp)[I (α x - α c t)]]]

Which is the result for the x equation

os3=ymom/.{u[x, y, t] -> 
              uh[y] Exp[I (α x - α c t)],v[x, y, t] -> 
              vh[y] Exp[I (α x - α c t)],p[x, y, t] -> 
              ph[y] Exp[I (α x - α c t)],Derivative[a1_, a2_, a3_][u][x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  uh[y] Exp[I (α x - α c t)]),Derivative[a1_, a2_, a3_][v][x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  vh[y] Exp[I (α x - α c t)]),Derivative[a1_, a2_, a3_][p][x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  ph[y] Exp[I (α x - α c t)])}

os4 = 
      Cancel[Expand[(os3/Exp)[I (α x - α c t)]]]

which is the result for the y equation

os5=cont/.{u[x, y, t] -> 
              uh[y] Exp[I (α x - α c t)],v[x, y, t] -> 
              vh[y] Exp[I (α x - α c t)],p[x, y, t] -> 
              ph[y] Exp[I (α x - α c t)],Derivative[a1_, a2_, a3_][u][x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  uh[y] Exp[I (α x - α c t)]),Derivative[a1_, a2_, a3_][v][x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  vh[y] Exp[I (α x - α c t)]),Derivative[a1_, a2_, a3_][p][x, y, t]:>(∂_({x, a1}, {y, a2}, {t, a3}))(
                  ph[y] Exp[I (α x - α c t)])}

os6 = 
      Cancel[Expand[(os5/Exp)[I (α x - α c t)]]]

which is the continuity equation

Now we use the continuity equation and the disturbance stream function definition
which is uh[y]=∂ψ[y]/∂y and vh[y]=-I αψ[y].

Back to conclusions

os7=os2/.{uh[y] -> (∂_y) ψ[y],vh[y] -> -I α ψ[y],Derivative[a1_][uh][y]:>\(∂\_{y, a1}\((∂\_y ψ[y])\)\),Derivative[a1_][vh][y]:>\(∂\_{y, a1}\((\(-I\)\ α\ ψ[y])\)\)}

os8=os4/.{uh[y] -> (∂_y) ψ[y],vh[y] -> -I α ψ[y],Derivative[a1_][uh][y]:>\(∂\_{y, a1}\((∂\_y ψ[y])\)\),Derivative[a1_][vh][y]:>\(∂\_{y, a1}\((\(-I\)\ α\ ψ[y])\)\)}

Now combine the x and y equations to eliminate pressure,

os9 = (∂_y) os7 - I α os8

Simplify[Cancel[Expand[(os9/(I α))]]]

This, with little effort, is the Orr-Sommerfeld equation.

Conclusions

1.  We first believe or observe (or learn from a combination of mathematical analysis and observation) that fluid dynamical instabilities arise when the flow is such that very small amplitude disturbances are amplified by the main flow.   (click    here to read again)

2.  Because of this infinitesimal origin of instabilities, the problems are amenable to linear analysis.  That is, the nonlinear governing equations can be linearized before they are solved.  This makes the problems much more tractable.   (click   here to read again) 

3.  Analysis of instability requires us to choose the mathematical form of the disturbance.  For waves on a liquid interface this would be a spatially and temporally periodic form.  We recognize that waves could grow in space and or time and they can travel in space.   It is important to use a form that  can be used as the basis of a Fourier Series so that any initial shape of "noise" could be represented.  Recall that because the equations are linear, the modes are independent.  

4.  A complex exponential is a good form for a periodic disturbance as it can be divided back out of any derivation (for linear equations) even after arbitrary derivatives have been taken.  We see that this introduction of an imaginary part of an expression does not lead to anything unphysical.  The presence of an  i just denotes a phase shift.  

5.  We find that waves can grow in either space or time and there is a relation between these,  as recognized by  Gaster. 

6.  It is important to understand (and not later forget!) that the linearization  that was done means that the solutions are not valid once the wave amplitude becomes large enough that the quadratic terms in amplitude can no longer be neglected.   

7.  The Kelvin-Helmholtz example shows the manipulations necessary to deal with separate governing equations for each fluid (because the fluid properties are different), and boundary conditions.  

8.  The Orr-Sommerfeld equation is the general starting point for stability analysis of Newtonian fluids with constant properties.  The derivation from the linearized Navier-Stokes equations using the definition of a disturbance stream function (obtained from the continuity equation) and a general  traveling, periodic disturbance  is shown.