Linear stability of long waves in two-layer channel flow

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, if you use it, this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version:  5/14/99
More recent versions of this notebook should be available at the web site:
http://www.nd.edu/~mjm/longwave.nb

Summary

This notebook solves the longwave stability problem for two-layer, pressure driven channel flow.  Through this example, the method for using a regular perturbation technique to solve a differential eigenvalue problem, with boundary conditions, is demonstrated.  

References: 
C. -S Yih (1967) "Instability due to viscosity stratification", J. Fluid Mech., 27 pp337-352.
S. G. Yiantsios and B. G. Higgins (1988) "Linear stability of plane Poiseullie flow of two-superposed fluids", Phys. Fluids, 31 pp3225-3238.

The 0 order and first order terms for the wave velocity for a pressure driven channel flow are obtained from a long wave expansion.  All of the necessary manipulations are shown for this direct perturbation solution method.

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;
Importance of long wave stability on flow regime transitions.

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 importance of flow regime is demonstrated in Demonstration of the effect of flow regime on pressure drop in multifluid flows.  For regime transitions where formation of waves are important, stratified to slug or stratified to annular, it is the growth of long waves that have been implicated in the transition.   

This notebook also shows how to use regular perturbation theory to solve an eigenvalue problem that involves boundary conditions.

[BACK] to "How to"

When a gas and a liquid are forced to flow together inside a pipe, there are at least 7 different geometrical configurations, or flow regimes, that are observed to occur. The regime depends on the fluid properties, the size of the conduit and the flow rates of each of the phases. The flow regime can also depend on the configuration of the inlet; the flow regime may take some distance to develop and it can change with distance as (perhaps) the pressure, which affects the gas density, changes. For fixed fluid properties and conduit, the flow rates are the independent variables that when adjusted will often lead to changes in the flow regime.

Air-water flow in a 1.27 cm diameter pipe oriented horizontally. 

These videos are from "Two Phase Flow Regimes in Reduced Gravity", NASA Lewis Research Center Motion Picture Directory 1704.

 By J. B. McQuillen, R. Vernon and A. E. Dukler. 
 
The video was taken at 400 fps and the projection is at 29.97 fps

Earth gravity

Bubbly flow(.mov movie)

Superficial gas velocity = .16 m/s,  Superficial liquid velocity = .90 m/s.

In this example of bubbly flow, the liquid flow rate is high enough to break up the gas into bubbles, but it is not high enough to cause the bubbles to become mixed well within the liquid phase. The Froude number, 

F = ((g d)/((Subscript[U, L])^2)) = .15

where g is gravitation acceleration, d is the pipe diameter and (Subscript[U, L]) is liquid superficial velocity, shows that while inertia is slightly larger than gravity forces, the experiments show that there is not enough mixing to scatter the bubbles.  

If the pipe were oriented vertically, the phase orientation would be symmetric, but there would likely be "slip" between the phases and the gas would not move at the same speed as the liquid.  

Annular Flow (.mov movie)

Superficial gas velocity = 7.4 m/s, Superficial liquid velocity = .08 m/s.

In annular flow, the liquid coats the walls.  However, because of gravity, the liquid distribution is not symmetric.  There is much more liquid on the bottom of the pipe than the top.  The velocity of the gas is large enough to cause waves to form in the liquid and also to atomize some liquid.  The maximum possible wave amplitudes scale, for liquid layers that are not too thin, as roughly the liquid thickness.  

Slug flow (.mov movie)

Superficial gas velocity = .17 m/s, Superficial liquid velocity = .08 m/s.

The slug regime, is characterized by the presence of liquid rich slugs that span the entire channel or pipe diameter.  These travel at a speed that is a substantial fraction of the gas velocity and they occur intermittently.  Slugs cause large pressure and liquid flow rate fluctuations.  The movie shows the approach of first a large wave and later a long slug.  Other movies of slugs would show much more gas entrainment and a flow that looks much more violent.  The length to diameter ratio of slugs varies greatly with flow rates, pipe diameter and fluid properties.  If the diameter is very large, F can always be large and slug flow, where the entire diameter is bridged, will not form.  Instead roll waves waves, which are breaking traveling waves, will be seen.  Liquid may or may not coat the entire pipe because there will be substantial atomization.  


Here are the same flow rates if gravity is reduced to an insignificant level. 

Bubbly flow

Superficial gas velocity = .13 m/s, Superficial liquid velocity = .89 m/s.

In this example of bubbly flow, there is no gravity so that there is no buoyancy force on the gas bubbles.  Thus they mix freely within the liquid.  While there is definitely a continuous phase (liquid) and a dispersed phase (gas), this is close to the idealized homogenous or dispersed flow that will be used for comparison in the examples below.  

Annular flow

Superficial gas velocity = 8.1 m/s, Superficial liquid velocity = .08 m/s.

Now that gravity is removed, the liquid distribution is uniform around the pipe.  Large disturbance waves still occur but they are seen as "ring - like" disturbances.  The absence of gravity also increases the amplitude to film depth ratio of traveling waves because there is no liquid drainage from the wave caused by gravity. 

Slug Flow

Superficial gas velocity = .16 m/s, Superficial liquid velocity = .08 m/s.

In the absence of gravity, the liquid distribution is uniform and the slugs are now liquid "trapped" between traveling "Taylor Bubbles".  This flow will not experience large pressure fluctuations and the flow rate fluctuations occur only on the size of the bubbles.  This is close to the idealized slug regime that is considered in the calculations below. 

The calculations below show directly that for the easily calculable case of laminar flow, the flow regime greatly influences the pressure - drop flow rate relation.  Thus, if design or operation of a device requires accurate knowledge of flow rate and pressure drop, there is a need to know the flow regime.  The rates of heat and mass transfer are also often important in process equipment and the movies suggest that these will also depend significantly on the flow regime!

The base state flow

This notebook addresses the case of Poiseulle flow so that the walls are fixed.
The upper fluid is 1 and the lower is 2, the velocity at the interface, y=0, is used to normalize the velocity so that we have

u1  = 1 + a1 y + b1 (y^2);

u2  = 1 + a2 y + b2 (y^2);

The constants are given here but we will substitute these as 
      needed later

a1 = ((m - (n^2))/((n^2) + n));

b1 = -((m + n)/((n^2) + n));

a2 = (a1/m);

b2 = (b1/m);

The m's and n's are the viscosity ratio, m= ((Subscript[μ, 2])/(Subscript[μ, 1])) and the thickness ratio, n=((Subscript[d, 2])/(Subscript[d, 1])).

\(\(profile1 = 
      Plot[u1 //. {a1 -> \((m - n^2)\)/\((n^2 + n)\), 
				a2 -> a1/m, \ 
            b1 -> \(-\((m + n)\)\)/\((n^2 + n)\), 
				b2 -> b1/m, 
            m -> 20, n -> 1\/2}, {y, 0, 1}, PlotRange -> All, 
        
		DisplayFunction -> Identity]; \)\)

\(\(profile2 = 
      Plot[u2 //. {a1 -> \((m - n^2)\)/\((n^2 + n)\), 
				a2 -> a1/m, \ 
            b1 -> \(-\((m + n)\)\)/\((n^2 + n)\), 
				b2 -> b1/m, 
            m -> 20, n -> 1\/2}, {y, \(-\(1\/2\)\), 0}, 
        PlotRange -> All]; \)\)

\(\(vprofile = 
      Show[profile1, profile2, 
			DisplayFunction -> $DisplayFunction]; 
    \)\)

This looks OK for the case of the lower fluid (phase 2) being more viscous.

The stability equations for the disturbance flow

The governing equations will be the Orr-Sommerfeld equation for each phase.

The equations for the disturbance stream functions in both phases are the Orr-Sommerfeld equation.

Upper phase

os1=Derivative[4][ϕ][y]-2 (α^2) ϕ''[y]+(α^4) ϕ[y]-I α R (((Subscript[U, 1])[y] - c) (ϕ''[y]-(α^2) ϕ[y])-((Subscript[U, 1]))''[y] ϕ[y])

lower phase

os2=Derivative[4][ψ][y]-2 (α^2) ψ''[y]+(α^4) ψ[y]-(1/m)(I α R r (((Subscript[U, 2])[y] - c) (ψ''[y]-(α^2) ψ[y])-((Subscript[U, 2]))''[y] ψ[y]))

Some narration about boundary conditions

Here is a brief discussion of the boundary conditions.

First the velocity is 0 on the walls
***   	ϕ[0] = ϕ'[0] = 0,
***   	ψ[-n] = ψ'[-n] = 0

Next the tangential velocity at the interface
  
 Note that this says that u1(η) = u2(η).  Thus we must expand the velocity to include the effect of the average velocity
  
 ***	 ϕ'[0] + η \(U1\&_\)'[y]  = ψ'[0] + η \(U2\&_\)'[y]
  
 this is
 
η is given from the Kinematic condition D[η,t] = (Subscript[v, y])[0]

substituting gives:

***   	ϕ'[0]-ψ'[0] = ϕ[0](a2-a1)/(c-\(U\&_\)[0])

Next the normal velocity at the interface 

d[ψ,x] = d[ϕ,x]  which gives simply

***  	ψ[0] = ϕ[0]

Next, the tangential stress is continuous across the interface

(Subscript[μ, 1]) (∂u1/∂y + ∂v1/∂x)  =  (Subscript[μ, 2]) (∂u2/∂y + ∂v2/∂x)
			
which gives
			
*** 	ϕ''[0] + (α^2) ϕ[0] = m (ψ''[0] + (α^2) ψ[0])
			
Finally, we have the continuity of Normal stress.  This is comprised of, in dimensional terms,

(p - ρ g a - μ dv/dy ) = σ ∂2 η/∂x2 
		
which says that the pressure change differs by the surface tension jump.  
		
We substitute for the pressure from the governing equation and make some other simplifications and get
		
		-i α R   ((c-\(U\&_\)[0]) ϕ'[0] + a1 ϕ[0]) -      ϕ'''[0] + 3 (α^2) ϕ'[0] +
		 i α R r ((c-\(U\&_\)[0]) ψ'[0] + a2 ψ[0]) + m( ψ'''[0] - 2 (α^2) ψ'[0]) = 
		 		        I α R (F + (α^2) S) ϕ[0]/(c-\(U\&_\)[0])

Expansions for the stream functions and wave speed

The problem will now be solved using the long wave expansion, i.e. α-->0 so that that we can write the stream functions and the wave speed as perturbation expansions in α.

\(ϕexp = 
      ϕ\_0[y] + α\ ϕ\_1[y] + α\^2\ ϕ\_2[y]; 

    ψexp = 
      ψ\_0[y] + α\ ψ\_1[y] + α\^2\ ψ\_2[y]; 

    cexp = c0 + α\ c1 + α\^2\ c2; \)

The procedure is to expand the governing equations and the boundary conditions in the expansion variables.  These can then be collected in powers of α.  Each separate power of α must be independent which allows a sequential solution to the problem.  Higher orders depend on known solutions to lower orders.

Expand the upper phase os-equ

os11=os1/.{Derivative[a1_][ϕ][y]:>\(∂\_{y, a1}ϕexp\),c -> cexp,ϕ[y] -> ϕexp,((Subscript[U, 1]))''[y]->D[u1, {y, 2}],(Subscript[U, 1])[y] -> u1}

os12 = Expand[os11];

os13 = Collect[os12, α]

os1O0 = Coefficient[os13, α, 0]

os1O1 = Simplify[Coefficient[os13, α, 1]]

os1O2 = Simplify[Coefficient[os13, α, 2]]

Expand the lower phase os-equ

os21=os2/.{Derivative[a1_][ψ][y]:>\(∂\_{y, a1}ψexp\),c -> cexp,ψ[y] -> ψexp,((Subscript[U, 2]))''[y]->D[u2, {y, 2}],(Subscript[U, 2])[y] -> u2}

os22 = Expand[os21];

os23 = Collect[os22, α];

os2O0 = Coefficient[os23, α, 0]

os2O1 = Simplify[Coefficient[os23, α, 1]]

os2O2 = Coefficient[os23, α, 2]

expand the boundary conditions

Each of these boundary conditions are written so that the expression is equal to 0.

Here is, no tangential velocity at top wall, bc1

bc1=ϕ'[y]

bc11=bc1/.{Derivative[a1_][ϕ][y]:>\(∂\_{y, a1}ϕexp\),Derivative[a1_][ψ][y]:>\(∂\_{y, a1}ψexp\),c -> cexp,ϕ[y] -> ϕexp,((Subscript[U, 1]))''[y]->2 b1,(Subscript[U, 0]) -> 1}

bc12 = bc11 /. y -> 1

bc13 = Expand[bc12]

bc1O0 = Coefficient[bc13, α, 0]

bc1O1 = Coefficient[bc13, α, 1]

bc1O2 = Coefficient[bc13, α, 2]

Here is, no normal velocity at top wall, bc2

bc2 = ϕ[y]

bc21=bc2/.{Derivative[a1_][ϕ][y]:>\(∂\_{y, a1}ϕexp\),Derivative[a1_][ψ][y]:>\(∂\_{y, a1}ψexp\),c -> cexp,ϕ[y] -> ϕexp,((Subscript[U, 1]))''[y]->2 b1,(Subscript[U, 0]) -> 1}

bc22 = bc21 /. y -> 1

bc23 = Expand[bc22]

bc2O0 = Coefficient[bc23, α, 0]

bc2O1 = Coefficient[bc23, α, 1]

bc2O2 = Coefficient[bc23, α, 2]

Here is, the tangential velocity match, bc3

bc3=ϕ'[y]-ψ'[y]-((ϕ[y] (a2 - a1))/(c - (Subscript[U, 0]) ))

bc31=bc3/.{Derivative[a1_][ϕ][y]:>\(∂\_{y, a1}ϕexp\),Derivative[a1_][ψ][y]:>\(∂\_{y, a1}ψexp\),c -> cexp,ϕ[y] -> ϕexp,((Subscript[U, 1]))''[y]->2 b1,(Subscript[U, 0]) -> 1}

bc32 = bc31 /. y -> 0

bc33 = Expand[bc32] ;

bc34 = Simplify[Series[bc33, {α, 0, 3}]]

bc3O0 = Coefficient[bc34, α, 0]

bc3O1 = Coefficient[bc34, α, 1]

bc3O2 = Coefficient[bc34, α, 2]

Here is the normal velocity match,bc4

bc4 = ψ[y] - ϕ[y]

bc41=bc4/.{Derivative[a1_][ϕ][y]:>\(∂\_{y, a1}ϕexp\),Derivative[a1_][ψ][y]:>\(∂\_{y, a1}ψexp\),c -> cexp,ϕ[y] -> ϕexp,ψ[y] -> ψexp,((Subscript[U, 1]))''[y]->2 b1,(Subscript[U, 0]) -> 1}

bc42 = bc41 /. y -> 0

bc43 = Expand[bc42];

bc4O0 = Coefficient[bc43, α, 0]

bc4O1 = Coefficient[bc43, α, 1]

bc4O2 = Coefficient[bc43, α, 2]

Here is the interfacial shear stress match,bc5

bc5=ϕ''[y]+(α^2) ϕ[y]-m (ψ''[y]+(α^2) ψ[y])

bc51=bc5/.{Derivative[a1_][ϕ][y]:>\(∂\_{y, a1}ϕexp\),Derivative[a1_][ψ][y]:>\(∂\_{y, a1}ψexp\),c -> cexp,ϕ[y] -> ϕexp,ψ[y] -> ψexp,((Subscript[U, 1]))''[y]->2 b1,(Subscript[U, 0]) -> 1}

bc52 = bc51 /. y -> 0

bc53 = Expand[bc52];

bc54 = Collect[bc53, α];

bc5O0 = Coefficient[bc54, α, 0]

bc5O1 = Coefficient[bc54, α, 1]

bc5O2 = Coefficient[bc54, α, 2]

Here is the pressure boundary condition,bc6

bc6=-I α R ((c - (Subscript[U, 0])) ϕ'[y]+a1 ϕ[y])  -     Derivative[3][ϕ][y]+3 (α^2) ϕ'[y]+I α R r ((c - (Subscript[U, 0])) ψ'[y]+a2 ψ[y])+m (Derivative[3][ψ][y]-3 (α^2) ψ'[y])-((I α R (F + (α^2) S) ϕ[y])/(c - 
              (Subscript[U, 0])))

bc61=bc6/.{Derivative[a1_][ϕ][y]:>\(∂\_{y, a1}ϕexp\),Derivative[a1_][ψ][y]:>\(∂\_{y, a1}ψexp\),c -> cexp,ϕ[y] -> ϕexp,ψ[y] -> ψexp,((Subscript[U, 1]))''[y]->2 b1,(Subscript[U, 0]) -> 1}

bc62 = bc61 /. y -> 0

bc63 = Expand[bc62];

bc64 = Series[bc63, {α, 0, 3}]

bc6O0 = Coefficient[bc64, α, 0]

bc6O1 = Coefficient[bc64, α, 1]

bc6O2 = Coefficient[bc64, α, 2]

Here is no tangential velocity at bottom wall, bc7

bc7=ψ'[y]

bc71=bc7/.{Derivative[a1_][ϕ][y]:>\(∂\_{y, a1}ϕexp\),Derivative[a1_][ψ][y]:>\(∂\_{y, a1}ψexp\),c -> cexp,ϕ[y] -> ϕexp,((Subscript[U, 1]))''[y]->2 b1,(Subscript[U, 0]) -> 1}

bc72 = bc71 /. y -> -n

bc7O0 = Coefficient[bc72, α, 0]

bc7O1 = Coefficient[bc72, α, 1]

bc7O2 = Coefficient[bc72, α, 2]

Here is, no normal velocity at bottom wall, bc8

bc8 = ψ[y]

bc81=bc8/.{Derivative[a1_][ϕ][y]:>\(∂\_{y, a1}ϕexp\),Derivative[a1_][ψ][y]:>\(∂\_{y, a1}ψexp\),c -> cexp,ϕ[y] -> ϕexp,ψ[y] -> ψexp,((Subscript[U, 1]))''[y]->2 b1,(Subscript[U, 0]) -> 1}

bc82 = bc81 /. y -> -n

bc8O0 = Coefficient[bc82, α, 0]

bc8O1 = Coefficient[bc82, α, 1]

bc8O2 = Coefficient[bc82, α, 2]

Here is the solution to the 0 order problem

Here are the equations for the 0 order problem

os1O0

os2O0

bc1O0

bc2O0

bc3O0

bc4O0

bc5O0

bc6O0

bc7O0

bc8O0

Here we solve for the eigen functions

We solve as much as we can with one command.  We leave out the boundary condition with the eigen value in it and treat it later.

\(zerofuncs = 
      DSolve[{os1O0 == 0, os2O0 == 0, 
		bc1O0 == 0, bc2O0 == 0, 
          bc4O0 == 0, 
		bc5O0 == 0, bc6O0 == 0, bc7O0 == 0, bc8O0 == 0}, 
        
	{ϕ\_0[y], ψ\_0[y]\ }, y]\)

Let's choose (Subscript[c, 8])so as to make this simple as it is the arbitrary constant for this problem.  This makes the two stream functions unity at the interface.

zerofuncs1 = 
      Simplify[zerofuncs /. 
          C[8] -> -6 (m - (n^2))/2/m/(n^2)/(n + 1)]

Here is a check with eq 8 of Y&H

phi0 = (Subscript[ϕ, 0])[y] /. zerofuncs1[[1]]

Simplify[Coefficient[Expand[phi0], y, 3]]

Simplify[Coefficient[Expand[phi0], y, 2]]

Simplify[Coefficient[Expand[phi0], y, 1]]

Simplify[Coefficient[Expand[phi0], y, 0]]

psi0 = (Subscript[ψ, 0])[y] /. zerofuncs1[[1]]

Simplify[Coefficient[Expand[psi0], y, 3]]

Simplify[Coefficient[Expand[psi0], y, 2]]

Simplify[Coefficient[Expand[psi0], y, 1]]

Simplify[Coefficient[Expand[psi0], y, 0]]

Here we get the 0 order eigenvalue

bc3O0

\(eig01 = 
      bc3O0 //. {ϕ\_0[0] -> \((phi0 /. y -> 0)\)\ , \ 
          \(ϕ\_0'\)[0] -> \((D[phi0, y] /. y -> 0)\), 

          		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          ψ\_0[0] -> \((psi0 /. y -> 0)\)\ , \ 
          \(ψ\_0'\)[0] -> \((D[psi0, y] /. y -> 0)\), 
		
          a1 -> \((m - n^2)\)/\((n^2 + n)\), 
				a2 -> a1/m, \ 
          b1 -> \(-\((m + n)\)\)/\((n^2 + n)\)}\)

eig02 = eig01 /. c0 -> cx + 1

eig03 = Solve[eig02 == 0, cx]

c0eig = 1 + FullSimplify[cx /. eig03[[1]]]

Here we can check with Y&H eq. (9)

\(c00\  = 
      1 + \ \((2\ \((m - n^2)\)\ \((m - 1)\)\ \((n^3 + n^2)\))\)/
	
          \((\((n^2 + n)\) 
              \((n^4 + 4\ n^3\ m\  + \ 6\ n^2\ m + \ 4\ n\ m\  + \ m^2)\))
            \)\)

It works!!

check = Simplify[c00 - c0eig]

Here is the solution to the first order problem

Here are the equations for the first order problem

os1O1

os2O1

bc1O1

bc2O1

bc3O1

bc4O1

bc5O1

bc6O1

bc7O1

bc8O1

The solution to the first order problem

We are fortunate again that c1 appears only in bc3.  This condition will be saved for last.  
We will solve the two OS equations and boundary conditions in three sets because they are so big.

Start with the upper phase OS eq, + top wall BC's

\(eqs12 = 
      Simplify[{os1O1 == 0, \ bc1O1 == 0, bc2O1 == 0} /. {
            ϕ\_0[0] -> \((phi0 /. y -> 0)\)\ , \ \ \ \ 
            \(ϕ\_0'\)[0] -> \((D[phi0, y] /. y -> 0)\), 
			
            \(\(ϕ\_0'\)'\)[0] -> \((D[phi0, {y, 2}] /. y -> 0)\), 

            			\(\(\(ϕ\_0'\)'\)'\)[0] -> 
              \((D[phi0, {y, 3}] /. y -> 0)\), 
			ϕ\_0[y] -> phi0, 
            D[\ ϕ\_0[y], {y, a1_}] :> D[phi0, {y, a1}], 
		\ \ \ 
            ψ\_0[0] -> \((psi0 /. y -> 0)\)\ , \ 
            \(ψ\_0'\)[0] -> \((D[psi0, y] /. y -> 0)\), 
			
            \(\(ψ\_0'\)'\)[0] -> \((D[psi0, {y, 2}] /. y -> 0)\), 

            			\(\(\(ψ\_0'\)'\)'\)[0] -> 
              \((D[psi0, {y, 3}] /. y -> 0)\), 
			\ 
            ψ\_0[y] -> psi0\ , 
            D[\ ψ\_0[y], {y, a2_}] :> D[psi0, {y, a2}], 
			
            c0 -> c0eig}]\)

ans12phi = DSolve[eqs12, (Subscript[ϕ, 1])[y], y];

phi1temp1 = ans12phi /. {C[3] -> c3f, C[4] -> c4f};

phi1temp2 = (Subscript[ϕ, 1])[y] /. phi1temp1[[1]]

Here is the lower phase OS equation + bottom wall BC's

\(eqs12a = 
      Simplify[{os2O1 == 0, \ bc7O1 == 0, bc8O1 == 0} /. {
            ϕ\_0[0] -> \((phi0 /. y -> 0)\)\ , \ \ \ \ 
            \(ϕ\_0'\)[0] -> \((D[phi0, y] /. y -> 0)\), 
			
            \(\(ϕ\_0'\)'\)[0] -> \((D[phi0, {y, 2}] /. y -> 0)\), 

            			\(\(\(ϕ\_0'\)'\)'\)[0] -> 
              \((D[phi0, {y, 3}] /. y -> 0)\), 
			ϕ\_0[y] -> phi0, 
            D[\ ϕ\_0[y], {y, a1_}] :> D[phi0, {y, a1}], 
		\ \ \ 
            ψ\_0[0] -> \((psi0 /. y -> 0)\)\ , \ 
            \(ψ\_0'\)[0] -> \((D[psi0, y] /. y -> 0)\), 
			
            \(\(ψ\_0'\)'\)[0] -> \((D[psi0, {y, 2}] /. y -> 0)\), 

            			\(\(\(ψ\_0'\)'\)'\)[0] -> 
              \((D[psi0, {y, 3}] /. y -> 0)\), 
			\ 
            ψ\_0[y] -> psi0\ , 
            D[\ ψ\_0[y], {y, a2_}] :> D[psi0, {y, a2}], 
			
            c0 -> c0eig}]\)

ans12psi = DSolve[eqs12a, (Subscript[ψ, 1])[y], y];

psi1temp2 = (Subscript[ψ, 1])[y] /. ans12psi[[1]];

Short[psi1temp2]

Here are the remaining boundary conditions, except for the one that has c1 in it

\(\(eqsO1a = 
      \(({bc4O1 == 0, bc5O1 == 0, bc6O1 == 0})\) /. {
          ϕ\_0[0] -> \((phi0 /. y -> 0)\)\ , \ \ \ \ 
          \(ϕ\_0'\)[0] -> \((D[phi0, y] /. y -> 0)\), 
			
          \(\(ϕ\_0'\)'\)[0] -> \((D[phi0, {y, 2}] /. y -> 0)\), 
			
          \(\(\(ϕ\_0'\)'\)'\)[0] -> \((D[phi0, {y, 3}] /. y -> 0)\), 

          			ϕ\_0[y] -> phi0, 
          D[\ ϕ\_0[y], {y, a1_}] :> D[phi0, {y, a1}], 
		\ \ \ 
          ψ\_0[0] -> \((psi0 /. y -> 0)\)\ , \ 
          \(ψ\_0'\)[0] -> \((D[psi0, y] /. y -> 0)\), 
			
          \(\(ψ\_0'\)'\)[0] -> \((D[psi0, {y, 2}] /. y -> 0)\), 
			
          \(\(\(ψ\_0'\)'\)'\)[0] -> \((D[psi0, {y, 3}] /. y -> 0)\), 

          			\ ψ\_0[y] -> psi0\ , 
          D[\ ψ\_0[y], {y, a2_}] :> D[psi0, {y, a2}], 
          ϕ\_1[0] -> \((phi1temp2 /. y -> 0)\)\ , \ \ \ \ 
          \(ϕ\_1'\)[0] -> \((D[phi1temp2, y] /. y -> 0)\), 
			
          \(\(ϕ\_1'\)'\)[0] -> \((D[phi1temp2, {y, 2}] /. y -> 0)\), 

          			\(\(\(ϕ\_1'\)'\)'\)[0] -> 
            \((D[phi1temp2, {y, 3}] /. y -> 0)\), 
			
          ϕ\_1[y] -> phi1temp2, 
          D[ϕ\_1[y], {y, a1_}] :> D[phi1temp2, {y, a1}], 
			
          ψ\_1[0] -> \((psi1temp2 /. y -> 0)\)\ , \ \ \ \ 
          \(ψ\_1'\)[0] -> \((D[psi1temp2, y] /. y -> 0)\), 
			
          \(\(ψ\_1'\)'\)[0] -> \((D[psi1temp2, {y, 2}] /. y -> 0)\), 

          			\(\(\(ψ\_1'\)'\)'\)[0] -> 
            \((D[psi1temp2, {y, 3}] /. y -> 0)\), 
			
          ψ\_1[y] -> psi1temp2, 
          D[ψ\_1[y], {y, a1_}] :> D[psi1temp2, {y, a1}], 
			
          c0 -> c0eig}; \)\)

Short[eqsO1a]

O1bcstemp = Solve[eqsO1a, {C[4], c3f, c4f}];

Short[O1bcstemp]

Finally we work on bc3 which contains c1

\(\(O1bc3 = 
      bc3O1 /. {ϕ\_0[0] -> \((phi0 /. y -> 0)\)\ , \ \ \ \ 
          \(ϕ\_0'\)[0] -> \((D[phi0, y] /. y -> 0)\), 
			
          \(\(ϕ\_0'\)'\)[0] -> \((D[phi0, {y, 2}] /. y -> 0)\), 
			
          \(\(\(ϕ\_0'\)'\)'\)[0] -> \((D[phi0, {y, 3}] /. y -> 0)\), 

          			ϕ\_0[y] -> phi0, 
          D[\ ϕ\_0[y], {y, a1_}] :> D[phi0, {y, a1}], 
		\ \ \ 
          ψ\_0[0] -> \((psi0 /. y -> 0)\)\ , \ 
          \(ψ\_0'\)[0] -> \((D[psi0, y] /. y -> 0)\), 
			
          \(\(ψ\_0'\)'\)[0] -> \((D[psi0, {y, 2}] /. y -> 0)\), 
			
          \(\(\(ψ\_0'\)'\)'\)[0] -> \((D[psi0, {y, 3}] /. y -> 0)\), 

          			\ ψ\_0[y] -> psi0\ , 
          D[\ ψ\_0[y], {y, a2_}] :> D[psi0, {y, a2}], 
          ϕ\_1[0] -> \((phi1temp2 /. y -> 0)\)\ , \ \ \ \ 
          \(ϕ\_1'\)[0] -> \((D[phi1temp2, y] /. y -> 0)\), 
			
          \(\(ϕ\_1'\)'\)[0] -> \((D[phi1temp2, {y, 2}] /. y -> 0)\), 

          			\(\(\(ϕ\_1'\)'\)'\)[0] -> 
            \((D[phi1temp2, {y, 3}] /. y -> 0)\), 
			
          ϕ\_1[y] -> phi1temp2, 
          D[ϕ\_1[y], {y, a1_}] :> D[phi1temp2, {y, a1}], 
			
          ψ\_1[0] -> \((psi1temp2 /. y -> 0)\)\ , \ \ \ \ 
          \(ψ\_1'\)[0] -> \((D[psi1temp2, y] /. y -> 0)\), 
			
          \(\(ψ\_1'\)'\)[0] -> \((D[psi1temp2, {y, 2}] /. y -> 0)\), 

          			\(\(\(ψ\_1'\)'\)'\)[0] -> 
            \((D[psi1temp2, {y, 3}] /. y -> 0)\), 
			
          ψ\_1[y] -> psi1temp2, 
          D[ψ\_1[y], {y, a1_}] :> D[psi1temp2, {y, a1}], 
			
          c0 -> c0eig}; \)\)

Short[O1bc3]

bc3O1t2 = Solve[O1bc3 == 0, c1];

Short[bc3O1t2]

bc3O1t3 = bc3O1t2 /. O1bcstemp[[1]];

Here is the answer for c1

\(\(c1eig1 = 
      bc3O1t3 //. {a1 -> \((m - n^2)\)/\((n^2 + n)\), 
				a2 -> a1/m, \ 
          b1 -> \(-\((m + n)\)\)/\((n^2 + n)\), 
				b2 -> b1/m}; \)\)

Short[c1eig1]

Finally if we get the expression by itself and simplify it, the remaining constant, (Subscript[c, 3]) drops out.

see1 = Simplify[c1 /. c1eig1[[1]]]

Here is a check against the answers of Yiantsios and Higgins

We need can see that the real part of c1 is 0.

Simplify[N[see1 /. {m -> 2, n ->  .5, r -> 1, R -> 1, F -> 0}]]

Simplify[N[see1 /. {m ->  .4, n -> 4, r -> 2, R -> 10, F -> 5}]]

Let's try to get figure 2.  To produce the ordinate we choose R=1 and take Im[ ] to get rid of the i.

\(\(plot1 = 
      Plot[Im[N[see1 /. {n ->  .5, r -> 1, 
						R -> 1, F -> 0}]], {
          m, 0, 10}, PlotStyle -> {Dashing[{ .03,  .02}]}]; \)\)

plot2 = 
      Plot[Im[N[see1 /. {n -> 1, r -> 1, R -> 1, F -> 0}]], {m, 0, 10}];

plot3 = 
      Plot[Im[N[see1 /. {n -> 2, r -> 1, R -> 1, F -> 0}]], {m, 0, 10}, 
        PlotStyle -> {Dashing[{ .04,  .02,  .02,  .02}]}];

We seem to have perfect agreement. with fig 2a.

Show[{plot1, plot2, plot3}, PlotRange -> {- .01,  .02}];

If we try to do figure 2b there is a problem.  Y&H state in the text that the flow is neutrally stable if m=√n.  Fortunately they tell us why.  Looking at equation 2 we see that if m = (n^2) then (Subscript[a, 1]) = (Subscript[a, 2]) = 0!.  If this is the case then by equation 5c, we need (Subscript[c, 1]) = 0 so that \(\(ϕ\^'\)\) - \(\(ψ\^'\)\) can be finite at y=0.

Thus there is an error in figure 2b.  The curve should be m = (n^2), not what is plotted!!  You can check this further if you need to.