Laminar channel flow over a (long)wavy surface

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/pressure.shear.fluctuations.nb

Summary

This notebook attempts to provide insight into how different models for linear stability will predict the formation of waves and by implication, the transition from stratified flow to another regime.  It is expected for two fluids, where one is much more viscous and dense, and the other has a Reynolds number large compared to unity, that the pressure shear variations over a linear wave strongly control the wave growth.  See Hanratty (1983) in Waves on fluid interfaces, p221-259. Academic Press, Kuru, Sangalli, Uphold and McCready (1995) IJMF, 21, p733-753.

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 have used a mix of input notation to show how it can be done.

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])\)\)

Finally it may be more easily read in Traditional form as

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))

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.

Motivation for this notebook

A number of studies have suggested that regime transitions from stratified to slug and/or annular flow occur when then stratified flow becomes unstable to long waves.  The plot above shows a number of different models, most based on linear stability theory, to predict these flow regime transitions.  They do not agree.  Further, the integral equation models, Lin & Hanratty, Barnea, Crowley et al., use an integrated form of the momentum equations -- as opposed to the differential equations, and these results do not agree (compare orange versus green lines).  This lack of agreement is troubling and there are two major issues.  

First we need to determine the inherent error in using differential as opposed to an integral form of the momentum equations.  We cannot expect these two to be the same.

The second issue concerns the effect of turbulence.  While for a laminar flow the differential equations can be solved for all conditions with the appropriate numerical or possibly analytical approach, a turbulent flow requires a model for the turbulence which is necessarily of somewhat limited accuracy.  (Note that it is effectively of unknown accuracy because turbulence models have been "tuned" for various single phase flows but not for the region near the interface of two-phase flows).  Given this question of accuracy, it might be reasonable to consider using an integral model for the turbulent case.  However, it is still important to know the inherent accuracy of integral models.

In this notebook, we attempt to get some understanding of both of these issues by solving the much simpler problem of flow over a wavy surface and looking at the fluctuations in pressure and shear stress that result.  We consider only the long wave region, H/λ <<1 (amplitude/H and amplitude/λ are also small).  Hanratty (1983) (in Waves on fluid interfaces, p221-259. Academic Press) emphasizes how important shear and pressure variations are in the wave formation problem for gas-liquid systems.  The pressure and shear stress components are be compared for different classes of problems to see the quantitative and qualitiative differences.  The integral momentum model of Lin and Hanratty (1986) IJMF is solved along with the differential laminar equations and a differential model that uses a polynomial model.

The wave is chosen to be long because it is the formation and growth of the longest waves that could cause regime transitions.  Further, the integral equations make sense only when the waves are long.

We hope that a comparison for laminar flow between the integral and differential models gives an estimate of the inherent accuracy of the integral model at predicing linear stability.  Further, the comparison of the polynomial differential and laminar differential models should give some idea how much difference turbulence might be in the stability problem.

Description and physical situation

This notebook solves the problem of laminar flow over a wavy solid surface.  The linearized equations are used.  A long wave expansion perturbation technique is employed.  The major point of this notebook is to calculate the pressure and shear stress components at the wave surface.  These provide insight into the wave growth process.  A comparison of the longwave perturbation and boundary-layer approximations are given.

The total average channel height is H and we nondimensionalize with it so that the top wall is at y=1.  Also, we use the average velocity, \(TraditionalForm\`U\&_\), to nondimensionalize the velocities.

How to navigate this notebook

It should be possible to run the "Derivations" section which will fill in answers for all of the stress and pressure variations.  The polynomial profile section will take several minutes to run on a Mac G3.  You could also make changes or correct my mistakes by altering these sections.  Alternatively, you can just use my answers by clicking on the "Results" section.

Preview of conclusions

A comparison of the laminar differential and integral models suggests that the first terms match exactly which for the stability problem control the wave speed.  However there is about a 20% difference between the pressure at second order and the integral model does not produce a second term in the shear stress.  The second term controls the wave growth and while the shear stress seems to be less important than pressure for wave growth, we expect at least this 20% difference to show up.  Thus we must expect a significant difference between differential and integral stability analyses.

The comparison of the laminar and polynomial differential models shows consistent significant differences.  Thus we can expect a difference in the stability boundaries.  This is, of course, expected.  The problem is that we can't tell exactly how to normalize the effects.

Derivations of pressure and shear variations using all of the models

Pressure and shear stress variations from the linearized differential  equation for a laminar flow.

Base state channel flow

We want to solve the problem of a laminar flow in a channel in which one wall is wavy.  We will locate the origin at the average location of the wavy wall and the channel top at y=h.  The average velocity profile will be parabolic, and we will be doing a perturbation solution, so first find the parabolic profile solution.  The governing equation is 0 = -(1/ρ) ∂(p/∂)x +  
        ν (∂^2) (u/(((∂y))^2))  with no slip boundary conditions at y = 0 and h.  

If we nondimensionalize with a length, h and the average velocity, U, we could solve this directly and get the parabolic profile as u (y)  =  6  y  +  6 ((y )^2) for u(y) nondimensionalized with the average velocity.  

However, lets solve this using the stream function, (Subscript[ψ, 0]) (y) . 

For (Subscript[ψ, 0]) (y) defined using:

	u(y) =  ((∂(Subscript[ψ, 0]) (y))/∂)y and v(y) = - ((∂(Subscript[ψ, 0]) (y))/(∂x ))
	
where u (y) is the flow direction velocity, we would substitute for u(y) with ((∂(Subscript[ψ, 0]) (y))/∂)y and then take a derivative with y to get rid of the pressure.  The result is: (((∂^4) (Subscript[ψ, 0]) (y))/∂)(y^4) = 0.

The form of (Subscript[ψ, 0])(y) is clearly a cubic polynomial.

ψ0  = a3 (y^3) + a2 (y^2) + a1 y + a0

The value of a0 is arbitrary because the governing equations are for the velocity, not the stream function, so they always involve derivatives of ∂(Subscript[ψ, 0]) (y).  

Let's fit the boundary conditions.

At the bottom (wavy) wall, u(0)=0

bc1 = ((∂_y) ψ0 /. y -> 0)  == 0

At the top wall, u(1)=0

bc2 = ((∂_y) ψ0 /. y -> 1)  == 0

Lets solve these,

ans1 = Solve[{bc1, bc2}, {a1, a2}]

We can substitute:

psi0temp1 = ψ0 /. ans1[[1]]

We might as well choose (Subscript[ψ, 0]) (0) = 0

psi0temp2 = psi0temp1 /. a0 -> 0

To get the last constant, we need to make the average flow rate match the real flow rate which is unity for the nondimensional problem.

temp1 = Integrate[D[psi0temp2, y], {y, 0, 1}]

We then can get a3 by

temp2 = Solve[temp1 == 1, a3]

Thus the stream function for the average problem is

psi0ans = psi0temp2 /. temp2[[1]]

Check to make sure that the velocity field is correct

uave = D[psi0ans, y]

Plot[uave, {y, 0, 1}];

One more check for the average velocity being correct,

(∫_0^1) uave \[DifferentialD]y

The wall shear stress is

lamavestress = 1/R D[uave, y] /. y -> 0

The pressure drop is

1/R  D[uave, {y, 2}]

This concludes the problem for the average flow.

Expansion of the equations for the perturbation problem

Momentum equation to get flow field, from Orr-Sommerfeld

Here is the governing equation, which is the steady Orr-Sommerfeld equation with an eigenvalue of 0.  The Orr-Sommerfeld equation arises for a periodic flow that is governed by the linearzed Navier-Stokes equations.  If this equation is unfamiliar, see the NoteBook: Introduction to linear stability for a derivation of the Orr-Sommerfeld equation.

goveq=((I (α^3) ϕ[y])/R)-(α^2) ϕ[y] u[y]-((2 I α ϕ''[y])/R)+u[y] ϕ''[y]-ϕ[y] u''[y]+((I Derivative[4][ϕ][y])/(α R))

We can expandϕ(y) in powers of α as

ϕexpand = 
      (Subscript[ϕ, 0])[y] + α (Subscript[ϕ, 1])[y] + (α^2) (Subscript[ϕ, 2])[y] + 
        (α^3) (Subscript[ϕ, 3])[y] + (α^4) (Subscript[ϕ, 4])[y] + 
        (α^5) (Subscript[ϕ, 5])[y] + (α^6) (Subscript[ϕ, 6])[y] + 
        (α^7) (Subscript[ϕ, 7])[y]

Now substitute:

expand1=goveq/.{Derivative[4][ϕ][y]:>\(∂\_{y, 4}ϕexpand\),ϕ''[y]:>\(∂\_{y, 2}ϕexpand\),ϕ[y] ->  ϕexpand, u[y] -> uave, u''[y] :> D[uave, {y, 2}]	}

expand2 = ExpandAll[α expand1];

expand3 =  Collect[expand2, α];

Here are the expanded equations

Here is the 0 order equation

eq0 = Coefficient[expand3, α, 0]

Here is the 1st order equation

eq1 = Coefficient[expand3, α, 1]

Here is the 2nd order equation

eq2 = Coefficient[expand3, α, 2]

Here is the 3rd order equation

eq3 = Coefficient[expand3, α, 3]

Here is the 4th order equation

eq4 = Coefficient[expand3, α, 4]

Momentum equation to get pressure

The equation for pressure can be obtained from the x  momentum equation a couple of steps before it is combimed with the y-momentum equation to get the Orr-Sommerfeld equation.  Note that we will then have to set the wave speed , c, to 0.

presseq0=I α press[y]-I α c ϕ'[y]+(((α^2) ϕ'[y])/R)+I α u[y] ϕ'[y]-I α ϕ[y] u'[y]-((Derivative[3][ϕ][y])/R)

presseq1 = presseq0 /. c -> 0

presseq2 = -Expand[(presseq1/(I α))] + press[y]

Now expand this in terms of α as well.

pexpand1=presseq2/.{Derivative[4][ϕ][y]:>\(∂\_{y, 4}ϕexpand\),Derivative[3][ϕ][y]:>\(∂\_{y, 3}ϕexpand\),\(\(ϕ'\)[y] :> ∂\_{y, 1}ϕexpand\),ϕ''[y]:>\(∂\_{y, 2}ϕexpand\),ϕ[y] ->  ϕexpand, u[y] -> uave, u''[y] :> D[uave, {y, 2}]	}

pexpand2 = ExpandAll[ pexpand1];

pexpand3 =  Collect[pexpand2, α];

Here are the expanded pressure equations

Here is the 0 order equation

peq0 = Coefficient[pexpand3, α, -1]

Here is the 1st order equation

peq1 = Coefficient[pexpand3, α, 0]

Here is the 2nd order equation

peq2 = Coefficient[pexpand3, α, 1]

Here is the 3rd order equation

peq3 = Coefficient[pexpand3, α, 2]

Here is the 4th order equation

peq4 = Coefficient[pexpand3, α, 3]

We also need the interfacial shear

The shear stress at the wall is given by, 1/R (∂v/∂x+ ∂u/∂y).  We will need to expand this aroung y=0, (domain perturbation) because the wall is not at 0, it is a y=  η and we don't want to have to evaluate the boundary condition there all of the time.  

This should be

τ(0) = 1/R ((α^2) ϕ+∂ϕ/∂y + ∂\(TraditionalForm\`U\&_\)/∂y) @ y=0

Here we solve the 0 order problem

Here is the first perturbation term:

eq0

eq0t1 = DSolve[eq0 == 0, (Subscript[ϕ, 0])[y], y]

eq0t2 = 
      (Subscript[ϕ, 0])[y] /. eq0t1[[1]]

The first boundary condition is that the top wall normal velocity is 0

bc01 = eq0t2 /. {y -> 1}

The next boundary condition is that the top wall tangential velocity is 0

bc02 = (∂_y) eq0t2 /. y -> 1

At the wavy wall, the normal velocity is 0

bc03 = eq0t2 /. y -> 0

At the wavy wall, the tangential velocity is 0.  Note that the wall is not always at 0 so we need to use domain perturbation to fix this.  It involves expanding the tangential velocity around 0.

The result is that u(η)=0= u'(0) + ∂\(TraditionalForm\`\(\(U\&_\)(0)\)/∂y\)

\(bc04a = \((∂\_y eq0t2 + ∂\_{y, 1}uave)\)\)

bc04 = bc04a /. y -> 0

consts0 = 
      Solve[{bc01 == 0, bc02 == 0, bc03 == 0, bc04 == 0}, {C[1], C[2], C[3], 
          C[4]}]

phi0 = eq0t2 /. consts0[[1]]

Calculate the 0 order pressure

pans01=peq0/.{((Subscript[ϕ, 0]))'[y]->(∂_y) phi0,(Subscript[ϕ, 0])[y] -> phi0,Derivative[3][(Subscript[ϕ, 0])][y]->\(∂\_{y, 3}phi0\)}

We recall that the first term has    (α^(-1))power so we write

pressure0 = pans01/α

Calculate the 0 order shear stress

\(stress0 = 
      \((1/R \((D[phi0, {y, 2}] + ∂\_{y, 2}uave)\))\) /. y -> 0\)

Here we solve the first order problem

Here is the next perturbation term:

eq1

eq1t1 = 
      eq1 /. {(Subscript[ϕ, 0])[y] -> phi0, (Subscript[ϕ, 0])''[y] -> D[phi0, {y, 2}]}

eq1t2 = DSolve[eq1t1 == 0, (Subscript[ϕ, 1])[y], y]

eq1t3 = 
      (Subscript[ϕ, 1])[y] /. eq1t2[[1]]

The first boundary condition is that the top wall normal velocity is 0

bc11 = eq1t3 /. {y -> 1}

The next boundary condition is that the top wall tangential velocity is 0

bc12 = (∂_y) eq1t3 /. y -> 1

At the wavy wall, the normal velocity is 0

bc13 = eq1t3 /. y -> 0

At the wavy wall, the tangential velocity is 0.
At first order we have that ∂(Subscript[ϕ, 1])/∂y = 0.  The average velocity is only in the lowest order α equation.

bc14 = ((∂_y) eq1t3) /. y -> 0

consts1 = 
      Solve[{bc11 == 0, bc12 == 0, bc13 == 0, bc14 == 0}, {C[1], C[2], C[3], 
          C[4]}]

phi1 = eq1t3 /. consts1[[1]]

Calculate the first order pressure at the wall

peq1

pans11=peq1/.{((Subscript[ϕ, 0]))'[y]->(∂_y) phi0,(Subscript[ϕ, 0])[y] -> phi0,Derivative[3][(Subscript[ϕ, 0])][y]->\(∂\_{y, 3}phi0\),
		((Subscript[ϕ, 1]))'[y]->(∂_y) phi1,(Subscript[ϕ, 1])[y] -> phi1,Derivative[3][(Subscript[ϕ, 1])][y]->\(∂\_{y, 3}phi1\)}

pressure1 = (pans11 /. y -> 0)

Thus we can write

pressure0  +  pressure1

Calculate the 0 order shear stress

stress1 = α  1/R (D[phi1, {y, 2}]) /. y -> 0

stress0 +  stress1

Here we solve the second order problem

Here is the next perturbation term:

eq2

\(eq2t1 = 
      eq2 /. {ϕ\_0[y] -> phi0, 
          \(\(ϕ\_0'\)'\)[y] -> D[phi0, {y, 2}], 
		
          ϕ\_1[y] -> phi1, \(\(ϕ\_1'\)'\)[y] -> D[phi1, {y, 2}]}
          \)

eq2t2 = DSolve[eq2t1 == 0, (Subscript[ϕ, 2])[y], y]

eq2t3 = 
      (Subscript[ϕ, 2])[y] /. eq2t2[[1]]

The first boundary condition is that the top wall normal velocity is 0

bc21 = eq2t3 /. {y -> 1}

The next boundary condition is that the top wall tangential velocity is 0

bc22 = (∂_y) eq2t3 /. y -> 1

At the wavy wall, the normal velocity is 0

bc23 = eq2t3 /. y -> 0

At the wavy wall, the tangential velocity is 0.
At second order we have that ∂(Subscript[ϕ, 2])/∂y = 0.

bc24 = ((∂_y) eq2t3 ) /. y -> 0

consts2 = 
      Solve[{bc21 == 0, bc22 == 0, bc23 == 0, bc24 == 0}, {C[1], C[2], C[3], 
          C[4]}]

phi2 = eq2t3 /. consts2[[1]]

Calculate the first order pressure at the wall

peq2

pans21=peq2/.{((Subscript[ϕ, 0]))'[y]->(∂_y) phi0,(Subscript[ϕ, 0])[y] -> phi0,Derivative[3][(Subscript[ϕ, 0])][y]->\(∂\_{y, 3}phi0\),
		((Subscript[ϕ, 1]))'[y]->(∂_y) phi1,(Subscript[ϕ, 1])[y] -> phi1,Derivative[3][(Subscript[ϕ, 1])][y]->\(∂\_{y, 3}phi1\),
			((Subscript[ϕ, 2]))'[y]->(∂_y) phi2,(Subscript[ϕ, 2])[y] -> phi2,Derivative[3][(Subscript[ϕ, 2])][y]->\(∂\_{y, 3}phi2\)}

pressure2 = α (pans21 /. y -> 0)

Thus we can write

pressure0  +  pressure1  +  pressure2

Calculate the 2 order shear stress

stress2 =  α /R (D[phi2, {y, 2}] +  phi0) /. y -> 0

stress0 +  stress1  +  stress2

Here we solve the third order problem

Here is the next perturbation term:

eq3

\(eq3t1 = 
      eq3 /. {ϕ\_0[y] -> phi0, 
          \(\(ϕ\_0'\)'\)[y] -> D[phi0, {y, 2}], 
		
          ϕ\_1[y] -> phi1, \(\(ϕ\_1'\)'\)[y] -> D[phi1, {y, 2}], 

          		ϕ\_2[y] -> phi2, \(\(ϕ\_2'\)'\)[y] -> D[phi2, {y, 2}]}
          \)

eq3t2 = DSolve[eq3t1 == 0, (Subscript[ϕ, 3])[y], y]

eq3t3 = 
      (Subscript[ϕ, 3])[y] /. eq3t2[[1]]

The first boundary condition is that the top wall normal velocity is 0

bc31 = eq3t3 /. {y -> 1}

The next boundary condition is that the top wall tangential velocity is 0

bc32 = (∂_y) eq3t3 /. y -> 1

At the wavy wall, the normal velocity is 0

bc33 = eq3t3 /. y -> 0

At the wavy wall, the tangential velocity is 0.
At second order we have that ∂(Subscript[ϕ, 2])/∂y = 0.

bc34 = ((∂_y) eq3t3 ) /. y -> 0

consts3 = 
      Solve[{bc31 == 0, bc32 == 0, bc33 == 0, bc34 == 0}, {C[1], C[2], C[3], 
          C[4]}]

phi3 = eq3t3 /. consts3[[1]]

Calculate the first order pressure at the wall

peq3

pans31=peq3/.{((Subscript[ϕ, 0]))'[y]->(∂_y) phi0,(Subscript[ϕ, 0])[y] -> phi0,Derivative[3][(Subscript[ϕ, 0])][y]->\(∂\_{y, 3}phi0\),
		((Subscript[ϕ, 1]))'[y]->(∂_y) phi1,(Subscript[ϕ, 1])[y] -> phi1,Derivative[3][(Subscript[ϕ, 1])][y]->\(∂\_{y, 3}phi1\),
			((Subscript[ϕ, 2]))'[y]->(∂_y) phi2,(Subscript[ϕ, 2])[y] -> phi2,Derivative[3][(Subscript[ϕ, 2])][y]->\(∂\_{y, 3}phi2\),
			((Subscript[ϕ, 3]))'[y]->(∂_y) phi3,(Subscript[ϕ, 3])[y] -> phi3,Derivative[3][(Subscript[ϕ, 3])][y]->\(∂\_{y, 3}phi3\)}

pressure3 = (α^2) (pans31 /. y -> 0)

Thus we can write

lampressure = 
      pressure0  +  pressure1  +  pressure2  +  pressure3

Calculate the 3rd order shear stress

stress3 =  
      ((α )^2) 1/R (D[phi3, {y, 2}] +  phi1) /. y -> 0

lamstress = stress0 +  stress1  +  stress2  +  stress3

Pressure and shear variation calculated from the linearized integral momentum equation for a laminar flow.

We start with the integral momentum equation for a single phase (see Lin and Hanratty, IJMF, 12, 79-98., 1986).  We include a perturbation parameter, ϵ, which is giving an idea how large the waves are.  We need ϵ to be small compared to H, the channel height, to make the linearization assumption be valid.

One dimensional mass balance

The mass balance relates changes in channel height to velocity changes for an incompresible fluid.

(∂_x)(u[x] h[x]) == 0

We expand this in terms of the average and fluctuating components.

\(mass1 = 
      ∂\_x\((
          \((U\&_ + ϵ\ u\&^\ \ Exp[I\ k\ x])\)\ 
            \((H\&_ - ϵ\ h\&^\ \ Exp[I\ k\ x])\))\)\)

mass2=Expand[mass1]

Here is the relation between the height change, \(\(h\&^\ \)\) and velocity change, \(\(u\&^\ \)\)for the perturbations

mass3= Expand[Coefficient[mass2,ϵ,1]/Exp[I k x]]

The substitution relation is

\(massans = Solve[mass3 == 0, u\&^]\)

Integral momentum equation

Here is the integral momentum equation in a form that will be expanded to get the average and perturbation equations.

\(momeq = 
      Γ\ 
          ∂\_x\((
              \((H\&_ - ϵ\ h[x])\)\ \((U\&_ + ϵ\ u[x])\)\^2)
              \) + 
	\(\((H\&_ - ϵ\ h[x])\)\ 
            \((dpdx + ϵ\ ∂\_x p[x] + 
                g\ ρ\ ∂\_x\((H\&_ - ϵ\ h[x])\)\ 
                  Sin[θ])\)\)\/ρ + 
        \(τ\_b0 + τ\_w0 + ϵ\ τ\_b[x] + 
            ϵ\ τ\_w[x]\)\/ρ\)

This equation is expanded and the individual powers of ϵ are collected.

ex1=Expand[momeq];

eqs1=Collect[ex1,ϵ];

First solve the basestate equation

basestate1=Coefficient[eqs1,ϵ,0]

basestate2 = basestate1 /. {(Subscript[τ, w0]) -> (Subscript[τ, b0])}

wallshearstress = Solve[basestate2 == 0, (Subscript[τ, b0])]

Calculation of the shear stress

We will need to eliminate the shear fluctuations in terms of average shear and the heights.  The way that this can be done is to assume that the spatial variations in the flow is slow enough that the average relation for the friction factor (or wall shear stress) holds at any location along the wave.  This is implemented by expanding the average relation in terms of a average and a perturbation equation.

Start with the average relation for shear stress and Reynolds number.

taueq1 = τ -  6 μ /ρ/U/H

Now expand.

\(taueq2 = 
      taueq1 /. {τ -> τ\&_ + \ ϵ\ τ\&^, \ 
          H -> H\&_ + \ ϵ\ h\&^, \ 
          U -> U\&_\  + 
					ϵ\ u\&^}\)

taueq3 = Series[taueq2, {ϵ, 0, 2}]

We first recover the average equation.

averagetau = Coefficient[taueq3, ϵ, 0]

We next obtain the first order term

tau10 =  Coefficient[taueq3, ϵ, 1]

tau11 = tau10 /. massans[[1]]

This relation for the perturbation shear stress is the same for both the top and bottom walls.  Note that the perturbation shear stress is twice the average shear stress.

\(tauans = FullSimplify[Solve[tau11 == 0, \ τ\&^]]\)

Now solve the linear perturbation problem

eq1=Coefficient[eqs1,ϵ,1]

We now use a normal mode expansion, that is, assume that the shape of the surface, h[x], can be written as one (or a whole sum of) sinuous modes.  Of course, since the equations are linear, these modes will be independent.   We use complex exponentials since they will cancel out easier.

eq2=eq1/.{\(h[x] -> h\&^\ Exp[I\ k\ x]\),\(τ\_w[x] -> τ\&^\_w\ \ Exp[I\ k\ x]\),\(τ\_b[x] -> τ\&^\_b\ \ Exp[I\ k\ x]\),h'[x]->\(∂\_x\((h\&^\ Exp[I\ k\ x])\)\),p'[x]->\(∂\_x\((p\&^\ Exp[I\ k\ x])\)\),u'[x]->\(∂\_x\((u\&^\ Exp[I\ k\ x])\)\)}

It is worth mentioning at this point, that this is no longer a differential equation.  The restriction of the solution to a periodic domain is enough to remove the first derivatives in x.  The higher order x-derivatives were removed in the "averaging" (i.e., when we decided to just do long waves and when we integrated in the normal direction.

eq3=Expand[eq2/Exp[I k x]]

\(pans = Solve[eq3 == 0, p\&^]\)

\(pressurefluc = p\&^ /. pans[\([1]\)]\)

pf1=Expand[pressurefluc]

pf2=Collect[pf1,k]

We now substitute the result from the mass balance to eliminate \(TraditionalForm\`u\&^\) into the equation for the pressure.

pf3=pf2/.massans[[1]]

Now make this nondimensional

\(pf4 = Expand[pf3\/\(ρ\ U\&_\^2\)]\)

We need to translate wavenumber, k, into the dimensionless wavenumber α.

\(pf5 = pf4 /. k -> α\/H\&_\)

We need to eliminate the stress perturbations in terms of the average pressure drop.

\(pf6 = 
      pf5 /. {τ\&^\_b -> 2\ h\&^\ τ\&_\ /\ H\&_, 
          τ\&^\_w -> 2\ h\&^\ τ\&_\ /\ H\&_}\)

Here we use the basestate solution again.

\(pf7 = 
      Expand[\((pf6\ *H\&_/h\&^\ )\) /. 
          τ\&_\  -> \(-\(\(dpdx\ H\&_\)\/2\)\)]\)

We can recover the Froude number, F.

\(pf8 = pf7 /. g -> \ F\ U\&_\^2/H\&_\)

We can eliminate the pressure drop in terms of the Reynolds number.  First we substitute for the pressure drop.

\(pf9 = pf8 /. dpdx -> \(-\(\(U\&_\ 12\ μ\)\/H\&_\^2\)\)\)

\(pf10 = pf9 /. U\&_ -> \(R\ μ\)\/\(ρ\ H\&_\)\)

We would like to know the value of Γ.

gamma=Integrate[(6*y - 6*y^2)^2,{y,0,1}]

The  final answer for the pressure variation is

pf11 = pf10 /. Γ -> gamma

The  final answer for the stress variation is

\(stressintlam = 
      \(-\((\((τ\&^ /. tauans[\([1]\)])\) H\&_/h\&^)\)\) /. 
        U\&_ -> \(R\ μ\)\/\(ρ\ H\&_\)\)

In terms of the notation for the differential problem we have

pressureintlam = pf10 /. Γ -> gamma

Pressure and shear stress variations from the linearized differential  equation for a polynomial profile (i.e. turbulent) flow.

Base state turbulent channel flow

We choose a polynomial profile. Brodkey and Hershey,Transport Phenomena: A Unified Approach, McGraw-Hill, 1988, p244, for a consistent average profile that fits the boundary conditions.  (Note that the much touted 1/7 power profile does not make sense nor fit the boundary conditions.)  We do not have any Reynolds stresses however.

The form of the turbulent profile translated for a domain of 0 to 1 (from -1 to 1) is  (y-> 2 (y-1/2))

tprof1=1+a1 4 (y-1/2)^2 + a2 2^(2 m) (y-1/2)^(2 m)

Here are the definitions of a1 and a2

tprof2 = 
      tprof1 /. {a1 -> (s - m)/(m - 1), a2 -> (1 - s)/(m - 1)}

We now normalize this so that the average velocity is 1.

tprof3 = 
      tprof2/(
          Integrate[tprof2, {y, 0, 1}, Assumptions -> {  Re[m] > 0}]/
            Integrate[1, {y, 0, 1}])

Here are the expressions for s and m

\(mm = \ Round[\(- .617\) +  .008211\ R^ .786]\ ; 

    ss = \  .585\  +  .003172\ R^ .833; \)

tproftest = (tprof3  /. {s -> ss, m -> mm}) /. R -> 5000

Plot[tproftest, {y, 0, 1}];

Lets look at some profiles,

tproftest2 = 
      ((tprof3 /. {s -> ss,  m -> mm}) /. R -> 4000)

Integrate[tproftest2, {y, 0, 1}]

psi0turb = Integrate[tprof3, y]

One more check for the average velocity being correct,

(∫_0^1) tprof3 \[DifferentialD]y

turbavestress = 1/R D[tprof3, y] /. y -> 0

The pressure drop is

dpturbave1 = 1/R D[tprof3, {y, 2}]

dpturbave2 = Simplify[dpturbave1 /. {m -> mm,  s -> ss}]

We could check to make sure that the "pressure drop" is not a function of y -- effectively it is not-- so set y=1/2

dpturbave3 = dpturbave2 /. y -> 1/2

This concludes the problem for the average flow.

Expansion of the equations for the perturbation problem

Momentum equation to get flow field, from Orr-Sommerfeld

Here is the governing equation, which is the steady Orr-Sommerfeld equation with an eigenvalue of 0.  The Orr-Sommerfeld equation arises for a periodic flow that is governed by the linearzed Navier-Stokes equations.  If this equation is unfamiliar, see the NoteBook: Introduction to linear stability for a derivation of the Orr-Sommerfeld equation.

goveq=((I (α^3) ϕ[y])/R)-(α^2) ϕ[y] u[y]-((2 I α ϕ''[y])/R)+u[y] ϕ''[y]-ϕ[y] u''[y]+((I Derivative[4][ϕ][y])/(α R))

We can expandϕ(y) in powers of α as

ϕexpand = 
      (Subscript[ϕ, 0])[y] + α (Subscript[ϕ, 1])[y] + (α^2) (Subscript[ϕ, 2])[y] + 
        (α^3) (Subscript[ϕ, 3])[y] + (α^4) (Subscript[ϕ, 4])[y] + 
        (α^5) (Subscript[ϕ, 5])[y] + (α^6) (Subscript[ϕ, 6])[y] + 
        (α^7) (Subscript[ϕ, 7])[y]

Now substitute:

expandturb1=goveq/.{Derivative[4][ϕ][y]:>\(∂\_{y, 4}ϕexpand\),ϕ''[y]:>\(∂\_{y, 2}ϕexpand\),ϕ[y] ->  ϕexpand}

expandturb2 = ExpandAll[α expandturb1];

expandturb3 =  Collect[expandturb2, α];

Here are the expanded equations

Here is the 0 order equation

eqturb0 = Coefficient[expandturb3, α, 0]

Here is the 1st order equation

eqturb1 = Coefficient[expandturb3, α, 1]

Here is the 2nd order equation

eqturb2 = Coefficient[expandturb3, α, 2]

Here is the 3rd order equation

eqturb3 = Coefficient[expandturb3, α, 3]

Momentum equation to get pressure

The equation for pressure can be obtained from the x  momentum equation a couple of steps before it is combimed with the y-momentum equation to get the Orr-Sommerfeld equation.  Note that we will then have to set the wave speed , c, to 0.

presseq0=I α press[y]-I α c ϕ'[y]+(((α^2) ϕ'[y])/R)+I α u[y] ϕ'[y]-I α ϕ[y] u'[y]-((Derivative[3][ϕ][y])/R)

presseq1 = presseq0 /. c -> 0

presseq2 = -Expand[(presseq1/(I α))] + press[y]

Now expand this in terms of α as well.

pturbexpand1=presseq2/.{Derivative[4][ϕ][y]:>\(∂\_{y, 4}ϕexpand\),Derivative[3][ϕ][y]:>\(∂\_{y, 3}ϕexpand\),\(\(ϕ'\)[y] :> ∂\_{y, 1}ϕexpand\),ϕ''[y]:>\(∂\_{y, 2}ϕexpand\),ϕ[y] ->  ϕexpand}

pturbexpand2 = ExpandAll[ pturbexpand1];

pturbexpand3 =  Collect[pturbexpand2, α];

Here are the expanded pressure equations

Here is the 0 order equation

pturbeq0 = Coefficient[pturbexpand3, α, -1]

Here is the 1st order equation

pturbeq1 = Coefficient[pturbexpand3, α, 0]

Here is the 2nd order equation

pturbeq2 = Coefficient[pturbexpand3, α, 1]

Here is the 3rd order equation

pturbeq3 = Coefficient[pturbexpand3, α, 2]

We also need the interfacial shear

The shear stress at the wall is given by, 1/R (∂v/∂x+ ∂u/∂y).  We will need to expand this aroung y=0, (domain perturbation) because the wall is not at 0, it is a y=  η and we don't want to have to evaluate the boundary condition there all of the time.  

This should be

τ(0) = 1/R ((α^2) ϕ+∂ϕ/∂y + ∂\(TraditionalForm\`U\&_\)/∂y) @ y=0

Here we solve the 0 order problem

Here is the first perturbation term:

eqturb0

eq0turb1 = DSolve[eqturb0 == 0, (Subscript[ϕ, 0])[y], y]

eq0turb2 = 
      (Subscript[ϕ, 0])[y] /. eq0turb1[[1]]

The first boundary condition is that the top wall normal velocity is 0

bcturb01 = eq0turb2 /. {y -> 1}

The next boundary condition is that the top wall tangential velocity is 0

bcturb02 = (∂_y) eq0turb2 /. y -> 1

At the wavy wall, the normal velocity is 0

bcturb03 = eq0turb2 /. y -> 0

At the wavy wall, the tangential velocity is 0.  Note that the wall is not always at 0 so we need to use domain perturbation to fix this.  It involves expanding the tangential velocity around 0.

The result is that u(η)=0= u'(0) + ∂\(TraditionalForm\`\(\(U\&_\)(0)\)/∂y\)

\(bcturb04a = 
      \((∂\_y eq0turb2 + ∂\_{y, 1}tprof3)\)\)

bcturb04 = bcturb04a /. y -> 0

conststurb0 = 
      Solve[{bcturb01 == 0, bcturb02 == 0, bcturb03 == 0, bcturb04 == 0}, {
          C[1], C[2], C[3], C[4]}]

phiturb0 = eq0turb2 /. conststurb0[[1]]

phi0turbsub = phiturb0 /. {m -> mm, s -> ss};

Short[phi0turbsub]

Calculate the 0 order pressure

pansturb01=FullSimplify[pturbeq0/.{((Subscript[ϕ, 0]))'[y]->(∂_y) phiturb0,(Subscript[ϕ, 0])[y] -> phiturb0,Derivative[3][(Subscript[ϕ, 0])][y]->\(∂\_{y, 3}phiturb0\)}]

We recall that the first term has    (α^(-1))power so we write

pressureturb0 = pansturb01/α

Calculate the 0 order shear stress

\(stressturb0 = 
      \((1/R \((D[phiturb0, {y, 2}] + ∂\_{y, 2}tprof3)\))\) /. 
        y -> 0\)

stressturb01 = stressturb0 /. {m -> mm,  s -> ss};

Short[stressturb01]

Table[stressturb01 /. {m -> mm, s -> ss}, {R, 2000, 10000, 1000}]

<<  Graphics`Graphics`

LogLogPlot[-stressturb01, {R, 2000, 100000}];

Here we solve the first order problem

Here is the next perturbation term:

eqturb1

This next one takes a long time.

\(eq1t1 = 
      Simplify[eqturb1 /. {ϕ\_0[y] -> phiturb0, 
            \(\(ϕ\_0'\)'\)[y] -> D[phiturb0, {y, 2}], 
			
            \(\(u'\)'\)[y] :> \ D[tprof3, {y, 2}], \ u[y] :> \ tprof3}]\)

This DSolve takes a while also.

eq1t2 = DSolve[eq1t1 == 0, (Subscript[ϕ, 1])[y], y];

eq1t3 = 
      (Subscript[ϕ, 1])[y] /. eq1t2[[1]]

The first boundary condition is that the top wall normal velocity is 0

bcturb11 = eq1t3 /. {y -> 1}

The next boundary condition is that the top wall tangential velocity is 0

bcturb12 = (∂_y) eq1t3 /. y -> 1;

Short[bcturb12]

At the wavy wall, the normal velocity is 0

bcturb13 = eq1t3 /. y -> 0

At the wavy wall, the tangential velocity is 0.
At first order we have that ∂(Subscript[ϕ, 1])/∂y = 0.  The average velocity is only in the lowest order α equation.

bcturb14 = ((∂_y) eq1t3) /. y -> 0

conststurb1 = 
      Simplify[Solve[{bcturb11 == 0, bcturb12 == 0, bcturb13 == 0, 
            bcturb14 == 0}, {C[1], C[2], C[3], C[4]}]];

phit1 = eq1t3 /. conststurb1[[1]];

Short[phit1]

Calculate the first order pressure at the wall

pturbeq1

pansturb11=pturbeq1/.{((Subscript[ϕ, 0]))'[y]->(∂_y) phiturb0,(Subscript[ϕ, 0])[y] -> phiturb0,Derivative[3][(Subscript[ϕ, 0])][y]->\(∂\_{y, 3}phiturb0\),
		((Subscript[ϕ, 1]))'[y]->(∂_y) phit1,(Subscript[ϕ, 1])[y] -> phit1,Derivative[3][(Subscript[ϕ, 1])][y]->\(∂\_{y, 3}phit1\),
		u[y] -> tprof3,u'[y] -> D[tprof3, y]};

Short[pansturb11]

pressureturb1 = (pansturb11 /. y -> 0)

Short[pressureturb1]

Thus we can write

totalturbpressure = pressureturb0  +  pressureturb1;

Short[totalturbpressure]

Calculate the 0 order shear stress

stressturb1 = α  1/R (D[phit1, {y, 2}]) /. y -> 0;

Short[stressturb1]

totalturbstress = Simplify[stressturb0 +  stressturb1];

Pressure and shear variation calculated from the linearized integral momentum equation for a turbulent flow.

We start with the integral momentum equation for a single phase (see Lin and Hanratty, IJMF, 12, 79-98., 1986).  We include a perturbation parameter , ϵ, which is giving an idea how large the waves are.  We need ϵ to be small compared to H, the channel height, to make the linearization assumption be valid.

One dimensional mass balance

The mass balance relates changes in channel height to velocity changes for an incompresible fluid.

(∂_x)(u[x] h[x]) == 0

We expand this in terms of the average and fluctuating components.

\(mass1 = 
      ∂\_x\((
          \((U\&_ + ϵ\ u\&^\ \ Exp[I\ k\ x])\)\ 
            \((H\&_ - ϵ\ h\&^\ \ Exp[I\ k\ x])\))\)\)

mass2=Expand[mass1]

Here is the relation between the height change, \(\(h\&^\ \)\) and velocity change, \(\(u\&^\ \)\)for the perturbations

mass3= Expand[Coefficient[mass2,ϵ,1]/Exp[I k x]]

The substitution relation is

\(massans = Solve[mass3 == 0, u\&^]\)

Integral momentum equation

Here is the integral momentum equation in a form that will be expanded to get the average and perturbation equations.

\(momeq = 
      Γ\ 
          ∂\_x\((
              \((H\&_ - ϵ\ h[x])\)\ \((U\&_ + ϵ\ u[x])\)\^2)
              \) + 
	\(\((H\&_ - ϵ\ h[x])\)\ 
            \((dpdx + ϵ\ ∂\_x p[x] + 
                g\ ρ\ ∂\_x\((H\&_ - ϵ\ h[x])\)\ 
                  Sin[θ])\)\)\/ρ + 
        \(τ\_b0 + τ\_w0 + ϵ\ τ\_b[x] + 
            ϵ\ τ\_w[x]\)\/ρ\)

This equation is expanded and the individual powers of ϵ are collected.

ex1=Expand[momeq];

eqs1=Collect[ex1,ϵ];

First solve the basestate equation

basestate1=Coefficient[eqs1,ϵ,0]

basestate2 = basestate1 /. {(Subscript[τ, w0]) -> (Subscript[τ, b0])}

wallshearstress = Solve[basestate2 == 0, (Subscript[τ, b0])]

Calculation of the shear stress

We will need to eliminate the shear fluctuations in terms of average shear and the heights.  The way that this can be done is to assume that the spatial variations in the flow is slow enough that the average relation for the friction factor (or wall shear stress) holds at any location along the wave.  This is implemented by expanding the average relation in terms of a average and a perturbation equation.

Start with the average relation for shear stress and Reynolds number.

\(turbtaueq1 = 
      τ - \  .079\ \((\ \(\(μ\ /ρ\)/U\)/H)\)\^ .25\)

Now expand.

\(turbtaueq2 = 
      turbtaueq1 /. {τ -> τ\&_ + \ ϵ\ τ\&^, \ 
          H -> H\&_ + \ ϵ\ h\&^, \ 
          U -> U\&_\  + 
					ϵ\ u\&^}\)

turbtaueq3 = Series[turbtaueq2, {ϵ, 0, 2}]

We first recover the average equation.

turbaveragetau = Coefficient[turbtaueq3, ϵ, 0]

We next obtain the first order term

turbtau10 =  Coefficient[turbtaueq3, ϵ, 1]

turbtau11 = turbtau10 /. massans[[1]]

This relation for the perturbation shear stress is the same for both the top and bottom walls.  Note that the perturbation shear stress is twice the average shear stress.

\(turbtauans = FullSimplify[Solve[turbtau11 == 0, \ τ\&^]]\)

Now solve the linear perturbation problem

eq1=Coefficient[eqs1,ϵ,1]

We now use a normal mode expansion, that is, assume that the shape of the surface, h[x], can be written as one (or a whole sum of) sinuous modes.  Of course, since the equations are linear, these modes will be independent.   We use complex exponentials since they will cancel out easier.

eq2=eq1/.{\(h[x] -> h\&^\ Exp[I\ k\ x]\),\(τ\_w[x] -> τ\&^\_w\ \ Exp[I\ k\ x]\),\(τ\_b[x] -> τ\&^\_b\ \ Exp[I\ k\ x]\),h'[x]->\(∂\_x\((h\&^\ Exp[I\ k\ x])\)\),p'[x]->\(∂\_x\((p\&^\ Exp[I\ k\ x])\)\),u'[x]->\(∂\_x\((u\&^\ Exp[I\ k\ x])\)\)}

It is worth mentioning at this point, that this is no longer a differential equation.  The restriction of the solution to a periodic domain is enough to remove the first derivatives in x.  The higher order x-derivatives were removed in the "averaging" (i.e., when we decided to just do long waves and when we integrated in the normal direction.

eq3=Expand[eq2/Exp[I k x]]

\(pans = Solve[eq3 == 0, p\&^]\)

\(pressurefluc = p\&^ /. pans[\([1]\)]\)

pf1=Expand[pressurefluc]

pf2=Collect[pf1,k]

We now substitute the result from the mass balance to eliminate \(TraditionalForm\`u\&^\) into the equation for the pressure.

pf3=pf2/.massans[[1]]

Now make this nondimensional

\(pf4 = Expand[pf3\/\(ρ\ U\&_\^2\)]\)

We need to translate wavenumber, k, into the dimensionless wavenumber α.

\(pf5 = pf4 /. k -> α\/H\&_\)

We need to eliminate the stress perturbations in terms of the average pressure drop.

\(pf6 = 
      pf5 /. {τ\&^\_b -> 2\ h\&^\ τ\&_\ /\ H\&_, 
          τ\&^\_w -> 2\ h\&^\ τ\&_\ /\ H\&_}\)

Here we use the basestate solution again.

\(pf7 = 
      Expand[\((pf6\ *H\&_/h\&^\ )\) /. 
          τ\&_\  -> \(-\(\(dpdx\ H\&_\)\/2\)\)]\)

We can recover the Froude number, F.

\(pf8 = pf7 /. g -> \ F\ U\&_\^2/H\&_\)

We can eliminate the pressure drop in terms of the Reynolds number.  First we substitute for the pressure drop for turbulent flow.

\(turbpf9 = 
      pf8 /. dpdx -> 
          \(-\ 0.079\)\ \(\((H\&_\ ρ\ U\&_/μ)\)\^\(- .25\)\) ρ\ 
            U\&_\^2/H\&_\)

\(turbpf10 = turbpf9 /. U\&_ -> \(R\ μ\)\/\(ρ\ H\&_\)\)

We would like to know the value of Γ.

turbgamma=1

The  final answer for the pressure variation is

pressureintturb = turbpf10 /. Γ -> turbgamma

The  final answer for the stress variation is

\(stressintturb = 
      \(-\((\((τ\&^ /. turbtauans[\([1]\)])\) /. 
              U\&_ -> \(R\ μ\)\/\(ρ\ H\&_\))\)\) H\&_/h\&^\)

Results

The analytical expressions

You can just click on this section and load in all of the results without rederiving them above.  Of course, this means that if you change something above, you do not want to click on the same result below because it will get changed to the old result!!

laminar stress

lamavestress=(6/R)

lamstress=-((I (17 (R^3) - 53900 R) 
                    (α^2))/(943250 R)) + 
            ((4 ((R^2) + 1078) α)/(2695 R)) + 
            ((6 I α)/35) + (12/R);

Back to conclusions

integral laminar stress

stressintlam = 12/R;

Back to conclusions.

laminar pressure

lampressure

lampressure=\(\(3\ \((411\ R\^3 - 1331330\ R)\)\ α\^2\)\/\(12262250\ R\) + 
          \((\(-\(\(3\ I\ \((\(-39\)\ R\^2 - 59290)\)\)\/\(13475\ 
                        R\)\)\) - \(6\ I\)\/R)\)\ α - 
	
          54\/35 + \(36\ I\)\/\(R\ α\)\)

Back to conclusions

integral laminar pressure

pressureintlam=F sin(θ) + ((36 I)/(R α)) - 
            (6/5);

Back to conclusions

polynomial stress

turbavestress

turbavestress=-((6 (m - 1) (2 m + 1) 
                  (((4 (((-1))^(2 m - 1)) m (1 - s))/(m - 
                          1)) - ((4 (s - m))/(m - 1))))/(R 
                  (-8 (m^2) - 4 (s - 2) m - 
                      3 (((-1))^(2 m)) + 
                      3 (((-1))^(2 m)) s + s + 3)));

This means nothing to us in this form but...

totalturbstress=(12 (2 m + 1) 
                (-280 
                      (-8 (m^2) - 4 (s - 2) m - 
                          3 (((-1))^(2 m)) + 
                          3 (((-1))^(2 m)) s + s + 3) 
                      (-(((-1))^(2 m)) m + 
                          (((-1))^(2 m)) s m + m - s) + 
                    140 (
                        -8 (m^2) - 4 (s - 2) m - 
                          3 (((-1))^(2 m)) + 
                          3 (((-1))^(2 m)) s + s + 3) 
                      (2 (((-1))^(2 m)) (s - 1) (m^2) + 
                          (-(((-1))^(2 m)) s + 
                                (((-1))^(2 m)) + 1) m - s) + 
                    (I (2 m + 1) R 
                          (32 (
                                  (((-1))^(2 m)) s - 
                                    (((-1))^(2 m)) + 1) (m^6) + 
                              32 s (
                                  6 (((-1))^(2 m)) s - 
                                    6 (((-1))^(2 m)) + 5) (m^5) - 
                              8 (3 (
                                        8 - 56 (((-1))^(2 m)) + 
                                        35 (((-1))^(4 m))) (s^2) + 
                                    (-168 + 
                                        398 (((-1))^(2 m)) - 
                                        210 (((-1))^(4 m))) s + 
                                    5 (
                                        25 - 46 (((-1))^(2 m)) + 
                                        21 (((-1))^(4 m)))) 
                                (m^4) - 
                              8 (420 
                                      (((-1 + (((-1))^(2 m)))
                                        )^2) + 
                                    3 (
                                        56 - 177 (((-1))^(2 m)) + 
                                        140 (((-1))^(4 m))) 
                                      (s^2) + 
                                    (-551 + 
                                        1371 (((-1))^(2 m)) - 
                                        840 (((-1))^(4 m))) s) 
                                (m^3) + 
                              (-3 
                                      (1136 - 
                                        2177 (((-1))^(2 m)) + 
                                        1225 (((-1))^(4 m))) 
                                      (s^2) + 
                                    (6531 - 13753 (((-1))^(2 m)) + 
                                        7350 (((-1))^(4 m))) s + 
                                    7222 (((-1))^(2 m)) - 
                                    3675 (((-1))^(4 m)) - 3547) 
                                (m^2) + 
                              (-840 
                                      (((-1 + (((-1))^(2 m)))
                                        )^2) - 
                                    3 (
                                        1057 - 1425 (((-1))^(2 m)) + 
                                        280 (((-1))^(4 m))) 
                                      (s^2) + 
                                    (4147 - 5955 (((-1))^(2 m)) + 
                                        1680 (((-1))^(4 m))) s)
                                   m + 
                              120 s 
                                ((-5 + 7 (((-1))^(2 m))) 
                                      s - 7 (((-1))^(2 m)) + 7))
                             α)/
                      ((m + 1) (m + 2) (2 m + 3) 
                          (2 m + 5))))/
            (35 R 
                (((-8 (m^2) - 4 (s - 2) m - 
                      3 (((-1))^(2 m)) + 
                      3 (((-1))^(2 m)) s + s + 3))^2));

polynomial pressure

totalturbpressure;

totalturbpressure=-(
                  144 (m - 1) 
                    (-(((-1))^(2 m)) m + 
                        (((-1))^(2 m)) s m + m - s) 
                    ((((((-1))^(2 m)) (1 - s))/(m - 1)) + 
                        ((s - m)/(m - 1)) + 1) (((2 m + 1))^2))/
                ((-8 (m^2) - 4 (s - 2) m - 
                        3 (((-1))^(2 m)) + 
                        3 (((-1))^(2 m)) s + s + 3) 
                    (8 (m^2) + 4 s m - 8 m + 3 (((-1))^(2 m)) - 
                        3 (((-1))^(2 m)) s - s - 3)) - 
            ((144 I 
                ((((-1))^(2 m)) (s - 1) m + m - s) 
                (2 m + 1))/(R 
                (3 (((-1))^(2 m)) (s - 1) + s - 
                    4 m (2 m + s - 2) + 3) α)) - 
            (1/R) 
              (I 
                  (-(
                            144 I (((-1))^(2 m - 3)) m
                               (2 m - 2) (2 m - 1) 
                              (4 m + 7) R (s - 1) 
                              (-(((-1))^(2 m)) m + 
                                  (((-1))^(2 m)) s m + m - s) 
                              (((2 m + 1))^2))/
                          ((m + 1) (m + 2) (2 m + 3) 
                              (2 m + 5) 
                              (((8 (m^2) + 4 s m - 8 m + 
                                    3 (((-1))^(2 m)) - 
                                    3 (((-1))^(2 m)) s - s - 3)
                                  )^2)) - 
                      (216 I (((-1))^(2 m - 2)) m 
                            (2 m - 1) 
                            (16 (m^3) - 64 (m^2) - 190 m - 65) R 
                            (s - 1) 
                            (-(((-1))^(2 m)) m + 
                                (((-1))^(2 m)) s m + m - s) 
                            (((2 m + 1))^2))/
                        ((m + 1) (m + 2) (2 m + 3) 
                            (2 m + 5) (4 m + 2) 
                            (((8 (m^2) + 4 s m - 8 m + 
                                  3 (((-1))^(2 m)) - 
                                  3 (((-1))^(2 m)) s - s - 3))^2))
                           - (
                          9 I R 
                            (768 
                                  ((((-1))^(2 m)) s - 
                                      (((-1))^(2 m)) + 1) (m^7) + 
                                64 (
                                    65 (((-1))^(2 m)) (s^2) - 
                                      52 (((-1))^(2 m)) s + 53 s - 
                                      13 (((-1))^(2 m)) + 13) 
                                  (m^6) + 
                                32 (
                                    65 (
                                        -2 + 15 (((-1))^(2 m)))
                                         (s^2) + 
                                      (949 - 1620 (((-1))^(2 m)))
                                         s + 
                                      645 (-1 + (((-1))^(2 m)))
                                        ) (m^5) + 
                                80 (
                                    (-390 + 
                                        1105 (((-1))^(2 m)) + 
                                        63 (((-1))^(4 m))) (s^2) + 
                                      (1363 - 
                                        2081 (((-1))^(2 m)) - 
                                        126 (((-1))^(4 m))) s + 
                                      976 (((-1))^(2 m)) + 
                                      63 (((-1))^(4 m)) - 1039) 
                                  (m^4) + 
                                16 (
                                    5 (-1105 + 
                                        1347 (((-1))^(2 m)) + 
                                        231 (((-1))^(4 m))) 
                                        (s^2) + 
                                      (12245 - 
                                        13158 (((-1))^(2 m)) - 
                                        2310 (((-1))^(4 m))) s + 
                                      3 (
                                        -2526 + 
                                        2141 (((-1))^(2 m)) + 
                                        385 (((-1))^(4 m)))) 
                                  (m^3) + 
                                2 (
                                    5 (-11280 + 
                                        4205 (((-1))^(2 m)) + 
                                        1827 (((-1))^(4 m))) 
                                        (s^2) + 
                                      (90889 - 
                                        47291 (((-1))^(2 m)) - 
                                        18270 (((-1))^(4 m))) s + 
                                      26266 (((-1))^(2 m)) + 
                                      9135 (((-1))^(4 m)) - 35401)
                                     (m^2) + 
                                (5 (-12106 - 
                                        1647 (((-1))^(2 m)) + 
                                        777 (((-1))^(4 m))) 
                                        (s^2) + 
                                      (80837 - 
                                        795 (((-1))^(2 m)) - 
                                        7770 (((-1))^(4 m))) s + 
                                      105 (
                                        -123 + 
                                        86 (((-1))^(2 m)) + 
                                        37 (((-1))^(4 m)))) m - 
                                15 s 
                                  ((669 + 259 (((-1))^(2 m))) 
                                        s - 7 
                                        (123 + 37 (((-1))^(2 m)))
                                        )) (2 m + 1))/
                        (35 (m + 1) (m + 2) (2 m + 3) 
                            (2 m + 5) 
                            (((-8 (m^2) - 4 (s - 2) m - 
                                  3 (((-1))^(2 m)) + 
                                  3 (((-1))^(2 m)) s + s + 3))^2))
                           + (
                          10368 I (((-1))^(2 m - 1)) m
                             (8 (m^4) - 34 (m^2) - 26 m - 5) R 
                            (s - 1) 
                            (-(((-1))^(2 m)) m + 
                                (((-1))^(2 m)) s m + m - s))/
                        ((m + 2) (2 m + 3) (2 m + 5) 
                            (4 m + 4) 
                            (((8 (m^2) + 4 s m - 8 m + 
                                  3 (((-1))^(2 m)) - 
                                  3 (((-1))^(2 m)) s - s - 3))^2))
                           - (
                          216 I (((-1))^(2 m)) 
                            (104 (m^4) + 84 (m^3) - 202 (m^2) - 189 m - 40)
                               R (s - 1) 
                            (-(((-1))^(2 m)) m + 
                                (((-1))^(2 m)) s m + m - s))/
                        ((m + 1) (m + 2) (2 m + 5) 
                            (4 m + 6) 
                            (((8 (m^2) + 4 s m - 8 m + 
                                  3 (((-1))^(2 m)) - 
                                  3 (((-1))^(2 m)) s - s - 3))^2))
                          ));

integral turbulent stress

stressintturb

stressintturb=0.0395 ((((1/R)))^0.25);

Integral turbulent pressure

pressureintturb

pressureintturb=F sin(θ) + 
            ((0.237 I)/((R^0.25) α)) - 1;

Comparison of laminar differential and integral cases.

The first term of the pressure and shear stress for each of the two cases are identical.  The second term in the pressure fluctuation is about 20% different.

N[(54/35 - 6/5)/(54/35)]

There is no second term in the shear stress for the integral model.  A normalization of the second term of the second term is 6/35 α/(6/R) = fluctation/average stress.  It could be that this is less important that the pressure in causing wave growth, but this difference could be important.

6/35 I α/(6/R)

Back to conclusions

Here are some comparisons of the two differential models

Here are the expressions for s and m

\(mm = \ Round[\(- .617\) +  .008211\ R^ .786]\ ; 

    ss = \  .585\  +  .003172\ R^ .833; \)

<<  Graphics`Graphics`

Shear stress

The real part of the polynomial shear stress normalized with the average shear stress.

turbTRplotRe = 
      LogLogPlot[
        -Re[N[(totalturbstress/turbavestress) /. {s -> ss, m -> mm, 
                  α ->  .01}]], {R, 2000, 20000}];

The imaginary part of the polynomial shear stress

turbTRplotIm = 
      LogLogPlot[
        Im[N[(totalturbstress/turbavestress) /. {s -> ss, m -> mm, 
                α ->  .01}]], {R, 2000, 20000}];

The real part of the laminar shear stress normalized with the average stress.

lamTRplotRe =  
      LogLogPlot[
        Re[N[(lamstress/lamavestress) /. {F -> 0,  α ->  .01}]], {
          R, 2000, 20000}];

The normalized imaginary part of the laminar shear stress.

lamTRplotIm =  
      LogLogPlot[
        -Im[N[(lamstress/lamavestress) /. {F -> 0,  α ->  .01}]]
            , {R, 2000, 20000}];

Stress comparisons

Here is the Real part of the shear stress. The laminar model has a much larger stress variation.

Show[turbTRplotRe, lamTRplotRe, PlotRange -> All];

Back to conclusions

Here is the imaginary part of the shear stress.  The laminar model has a much larger stress variation.

\(\(Show[turbTRplotIm, lamTRplotIm, PlotRange -> All, 

      AxesLabel -> {"Reynolds", "stress variation"}]; \)\)

Pressure

The real part of the polynomial pressure

turbPRplotRe = 
      LogLogPlot[
        -Re[N[( totalturbpressure/2/turbavestress) /. {s -> ss, 
                  m -> mm, α ->  .01}]], {R, 2000, 20000}];

The imaginary part of the polynomial shear stress

turbPRplotIm = 
      LogLogPlot[
        Im[N[(totalturbpressure /2/turbavestress) /. {s -> ss, 
                m -> mm, α ->  .01}]], {R, 2000, 20000}];

The real part of laminar pressure.  Note that it changes signs so we won't worry about the highest Re range

lamPRplotRe =  
      LogLogPlot[
        -Re[N[(lampressure/12 R) /. { α ->  .01}]], {R, 
          2000, 20000}];

There is no imaginary part of the laminar shear stress

lamPRplotIm =  
      LogLogPlot[
        Im[N[(lampressure/12 R) /. { α ->  .01}]], {R, 2000, 
          20000}];

Pressure comparisons

Here is the Real part of the pressure  The agreement is not real good since they head off in opposite directions.

\(\(Show[lamPRplotRe, \ turbPRplotRe, PlotRange -> All, 
      AxesLabel -> 
{"Reynolds", "P variation"}]; \)\)

Back to conclusions

Here is the imaginary part of the pressure.  The quantitative agreement starts off by about 50% then the polynomial model does not change.  This could be fixed by patiently solving for the next term.  Please let me know how it comes out!!

Show[lamPRplotIm,  turbPRplotIm, PlotRange -> All];

Here are some comparisons of the turbulent models

Shear stress

The real part of the polynomial shear stress

turbTRplotRe = 
      LogLogPlot[
        -Re[N[totalturbstress /. {s -> ss, m -> mm, α ->  .01}]], {
          R, 2000, 20000}];

The imaginary part of the polynomial shear stress

turbTRplotIm = 
      LogLogPlot[
        Im[N[totalturbstress /. {s -> ss, m -> mm, α ->  .01}]], {R, 
          2000, 20000}];

The real part of the integral shear stress

intturbTRplotRe =  
      LogLogPlot[
        Re[N[stressintturb /. {F -> 0,  α ->  .01}]], {R, 2000, 
          20000}];

There is no imaginary part of the integral shear stress

Stress comparisons

Here is the Real part of the shear stress, which is the first term (the only term for the integral model).  The agreement is not very good.

Back to conclusions

\(\(Show[turbTRplotRe, intturbTRplotRe, PlotRange -> All, 

      AxesLabel -> {"Reynolds", "stress variation"}]; \)\)

Pressure

The real part of the polynomial pressure

turbPRplotRe = 
      LogLogPlot[
        -Re[N[totalturbpressure  /. {s -> ss, m -> mm, α ->  .01}]]
            , {R, 2000, 20000}];

The imaginary part of the polynomial pressure

turbPRplotIm = 
      LogLogPlot[
        Im[N[totalturbpressure /. {s -> ss, m -> mm, α ->  .01}]], {R, 
          2000, 20000}];

The real part of the integral pressure

intturbPRplotRe =  
      LogLogPlot[
        -Re[N[pressureintturb /. {F -> 0,  α ->  .01}]], {R, 
          2000, 20000}];

The imaginary part of the integral pressure.

intturbPRplotIm = 
      LogLogPlot[
        Im[N[pressureintturb /. {F -> 0,  α ->  .01}]], {R, 2000, 
          20000}];

Pressure comparisons

Here is the out-of-phase (imaginary) pressure comparison, which is the lowest order term.  The agreement is getting somewhat worse at higher Re with about a 20% difference.

\(\(Show[turbPRplotIm, \ intturbPRplotIm, PlotRange -> All, 

      AxesLabel -> {"Reynolds", "pressure variation"}]; \)\)

Here is the in-phase pressure comparison which is the second term.  This never really agrees at all with a factor of 2 at the beginning and then a steady worsening.

Back to conclusions

\(\(Show[intturbPRplotRe, turbPRplotRe, 

      AxesLabel -> {"Reynolds", "pressure variation"}]; \)\)

Conclusions

1.  Comparision of laminar flow differential pressure, different stress, integral pressure and integral stress, shows that the integral model agrees to some extent, but that there can be important differences.  The 20% difference in pressure variation  and the arbitrary differences for the shear component that cause wave growth could be significant. 

2.  The laminar and polynomial differential results have more significant qualitative differences in both the pressure and shear stress.  Note that these comparisons are made after normalization by the average stress-- which is different for laminar versus turbulent flow. 

3.  The polynomial and turbulent integral models give very different answers for the pressure and shear stress.  This is troubling because even though the polynomial model is not necessarily correct, it produces the same average pressure drop as the integral model.  In the absence of a better differential model that reliably produces the correct basestates and has been verified to produce reasonable growth information, we cannot completely discount either of these two.  Thus, at the present time, there is a great deal of uncertainty in the prediction of the stability behavior for turbulent flows -- the only case that matters for process conditions.  

4.  Finally, we have shown above how to solve a difficult ordinary differential equation by a regular perturbation method.  Note that presence of inhomogeneities in the boundary conditions at the lowest order from domain perturbation that allows us to evaluate the interfacial conditions (very much more conveniently) at y=0, instead of the exact position of the interface.