Linear Stability of Pressure Driven Channel Flow

This notebook has been written in Mathematica by 

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

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

Mathematica aside

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}] -
		(Subscript[ρ, 1])(-D[(Subscript[ϕ, 1])[x,y,t],t]-U1 D[(Subscript[ϕ, 1])[x,y,t],x]-g η[x,t])+
			(Subscript[ρ, 2])(-D[(Subscript[ϕ, 2])[x,y,t],t]-
					U2 D[(Subscript[ϕ, 2])[x,y,t],x]-g η[x,t])

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

dc=-(γ (∂_({x,2}))η[x,t])-(Subscript[ρ, 1]) (-(∂_t)(Subscript[ϕ, 1])[x,y,t]-U1 (∂_x)(Subscript[ϕ, 1])[x,y,t]-g η[x,t])+(Subscript[ρ, 2]) (-(∂_t)(Subscript[ϕ, 2])[x,y,t]-U2 (∂_x)(Subscript[ϕ, 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.

This notebook gives the Orszag-tau spectral linear stability calculation for channel flow.

Reference:  S. A. Orszag (1971) "Accurate solution of the Orr-Sommerfeld stability equation", Journal of Fluid Mechanics, 50 pp 689-703.

The coefficients of the algebraic equations are computed directly from the orthogonality property of the Chebyshev polynomials using integration. 

Reference:  R. Miesen and B. J. Boersma (1995) "Hydrodynamic stability of a sheared liquid film", Journal of Fluid Mechanics, 301 pp 175-202.

In the tau method, the solution is approximated with Orthogonal polynomials that are valid over the entire domain, as in the Galerkin methods.  However, in many problems because of the complexity of the boundary conditions, it is not possible to make a convenient set of functions exactly fit the boundary conditions.  To solve this problem, the tau method enforces the boundary conditions as additional algebraic equations along with those that originate with the differential equation.  To obtain the extra degrees of freedom that this requires, a requisite number of the highest "frequency"  (i.e. highest n) modes generated from the differential equations are not enforced.  As long as the number of modes can be made large (>10-20) this method will probably work fine.  For the long wave region, numerical resolution limits the number of modes that can be used in this formulation and so there could be some problem getting the tau method to give reliable results.

This section expands the terms of the Orr-Sommerfeld equation in terms of the Chebyshev polynomials and gets the coefficients using the orthogonality.

We need use only even modes because the boundary conditions are symmetric.  There are then two boundary conditions, no slip and no flow through the wall, thus we will drop off the two highest equations produced from substitution of the Chebyshev's in the ODE.

Here is the start of the calculations.  
In this example we use nn/2 +1 modes.

nn=8;

nx=(nn/2)+1;

Here is the second derivative term of the Orr-Sommerfeld equation

TD=-2 (α^2)D[vv[y,i],{y,2}]

td1=TD/.{vv[y,i]->a[i] ChebyshevT[i,y],(∂_({y,2}))vv[y,i]->(∂_({y,2}))(a[i] ChebyshevT[i,y]),(∂_({y,4}))vv[y,i]->(∂_({y,4}))(a[i] ChebyshevT[i,y])}

temp2=Sum[td1,{i,0,nn,2}]

temp2a=Simplify[temp2]

Here is where we use orthogonality

z2d=Table[a2d[i]=((2 (∫_(-1)^1)temp2a ChebyshevT[i,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y)/(π c[i])),{i,0,nn,2}]

Here is the fourth derivative term of the Orr-Sommerfeld equation

FD=D[vv[y,i],{y,4}]

fd1=FD/.{vv[y,i]->a[i] ChebyshevT[i,y],(∂_({y,2}))vv[y,i]->(∂_({y,2}))(a[i] ChebyshevT[i,y]),(∂_({y,4}))vv[y,i]->(∂_({y,4}))(a[i] ChebyshevT[i,y])}

temp4=(∑_({{i=0}, {Δ i=2}})^nn)fd1;

temp4a=Simplify[temp4]

z4d=Table[a4d[i]=((2 (∫_(-1)^1)temp4a ChebyshevT[i,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y)/(π c[i])),{i,0,nn,2}]

Here is the zero derivative term of the Orr-Sommerfeld equation

ZD=(α^4) vv[y,i]

z00d=Table[a00d[i]=(α^4) a[i],{i,0,nn,2}]

Here is the rest of the Orr-Sommerfeld equation

XD=-I α rr(((Overscript[U, ─])[y]-λ) ((∂_({y,2}))vv[y,i]-(α^2) vv[y,i])-(∂_({y,2}))(Overscript[U, ─])[y] vv[y,i])

We can substitute the average profile and the (Overscript[U, ─])''(y)

xd1=XD/.{(Overscript[U, ─])[y]->1-(y^2),(Overscript[U, ─])''[y]->(∂_({y,2}))(1-(y^2))}

xd2=xd1/.{vv[y,i]->a[i] ChebyshevT[i,y],(∂_({y,2}))vv[y,i]->(∂_({y,2}))(a[i] ChebyshevT[i,y]),(∂_({y,4}))vv[y,i]->(∂_({y,4}))(a[i] ChebyshevT[i,y])}

tempxd=(∑_({{i=0}, {Δ i=2}})^nn)xd2;

tempxda=Simplify[tempxd]

z0d=Table[a0d[i]=((2 (∫_(-1)^1)tempxda ChebyshevT[i,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y)/(π c[i])),{i,0,nn,2}]

Here are the boundary conditions.

vexpand=(∑_({{i=0}, {Δ i=2}})^nn)a[i] ChebyshevT[i,y]

bc1=vexpand/.{y->-1}

bc2=(∂_y)vexpand/.{y->-1}

Tau solution method

We solve for the two highest coefficients because the two highest modes are dropped from the matrix of the differential equation mode coefficients (the Tau method).  The boundary conditions replace these modes assure that the BC's are satisfied.  However, the BC equations do not contain the eigen value and thus must be eliminated algebraically before solving the eigenvalue problem.

bcsolve=Solve[{bc1==0,bc2==0},{a[nn],a[nn-2]}]

Here is the table of coefficient equations from the O-S equation

dcoefs=z4d+z2d+z00d+z0d

Here is where we use the boundary conditions to replace the two highest coefficients.

dcoefsub=dcoefs/.bcsolve[[1]]

Here we pick out the number of valid equations

coefs=Table[dcoefsub[[i]],{i,1,nx-2}]

temp1=Expand[(coefs/.c[0]->2)/.Table[c[i]->1,{i,2,nn,2}]]

Here we construct the matrix that comprises the generalized eigenvalue problem 
AB = (AA  -   B  λ)

AB=Table[Table[Coefficient[temp1[[j]],a[i]],{i,0,nn-4,2}],{j,1,nx-2}]

The left side matrix will be AA

AA=AB/.λ->0;

Now we get BB

B=((AB-AA)/λ)

At this point we start doing things numerically.
First choose α and Reynolds number

alfx=1.02056;
rrx=5772.22;

Here is the calculation of Binverse to solve the GEP.

Binv=Inverse[-B];

Here is the matrix of the now regular eigenvalue problem

Binv.AA = AZ

AZ=Binv.AA;

Now make it numerical for sure

N[AZ/.{α->alfx,rr->rrx}]

Get the eigenvalues.

Eigenvalues[N[AZ/.{α->alfx,rr->rrx}]]

These are not very accurate!!  (See below for the good ones)

Here is the Mathematica Solver for the generalized eigenvalue problem

Eigenvalues[{AA,-B}/.{α->alfx,rr->rrx}]

Try a different Reynolds number and wavenumber.

alfx=1;
rrx=10000;

Here is the calculation of Binverse to solve the GEP.

Binv=Inverse[-B];

Here is the matrix of the now regular eigenvalue problem

Binv.AA = AZ

AZ=Binv.AA;

Now make it numerical for sure

N[AZ/.{α->alfx,rr->rrx}]

Get the eigenvalues.

Eigenvalues[N[AZ/.{α->alfx,rr->rrx}]]

Eigenvalues[N[{AA,-B}/.{α->alfx,rr->rrx}]]

Here is the same calculation but made to run faster numerically by using recursion formulas for the coefficients.

nn=72;

nx=(nn/2)+1;

Here is the second derivative term of the Orr-Sommerfeld equation

TDcoef=-2 (α^2)-I α rr+I α λ rr

The recursion formulas to construct the Matrix eigenvalue problem are available in Orszag's paper. There are several more that are useful when solving the two-layer interfacial eigenvalue problem are given in:

Primary and Secondary Interfacial Disturbances in Horizontal Cocurrent Flows
Ph. D. Thesis, 
William C. Kuru
Department of Chemical Engineering
University of Notre Dame
1995.

Note how these match the ones obtained by integration.

z2d=Table[TDcoef (∑_({{p=i+2}, {Δ p=2}})^nn)((p ((p^2)-(i^2)) a[p])/(c[i])),{i,0,nn,2}];

Here is the fourth derivative term of the Orr-Sommerfeld equation

z4d=Table[(∑_({{p=i+4}, {Δ p=2}})^nn)((p ((p^2) ((((p^2)-4))^2)-3 (i^2) (p^4)+3 (i^4) (p^2)-(i^2) ((((i^2)-4))^2)) a[p])/(24 c[i])),{i,0,nn,2}];

Here is the zero derivative term of the Orr-Sommerfeld equation.  I now have all of the zero derivative terms here.

z00d=Table[a00d[i]=((α^4)-2 I α rr+I (α^3) rr-I (α^3) λ rr) a[i],{i,0,nn,2}];

Here is the y^2 vv[y] term

zxt=Table[(c[i-2] a[i-2]+2 a[i]+a[i+2]) (1/1),{i,0,nn-2,2}];

z2y=-(1/4) I (α^3) rr Append[zxt,a[nn-2]+a[nn] 2];

Here is the y^2 v''[y] term

z0d=I α rr Table[((i (i-1) a[i]+(∑_({{j=i+2}, {Δ j=2}})^nn)j ((j^2)-(i^2)-2) a[j])/(c[i])),{i,0,nn,2}];

dcoefs=z4d+z2d+z00d+z0d+z2y;

Here are the boundary conditions

vexpand=(∑_({{i=0}, {Δ i=2}})^nn)a[i] ChebyshevT[i,y];

bc1=vexpand/.{y->-1};

bc2=(∂_y)vexpand/.{y->-1};

bcsolve=Solve[{bc1==0,bc2==0},{a[nn],a[nn-2]}];

dcoefsub=dcoefs/.bcsolve[[1]];

coefs=Table[dcoefsub[[i]],{i,1,nx-2}];

Here we replace all of the coefficients.

temp1=Expand[ (coefs/.c[0]->2)/.Table[c[i]->1,{i,2,nn,2}]];

Here we are back at the matrix

AB=Table[Table[Coefficient[temp1[[j]],a[i]],{i,0,nn-4,2}],{j,1,nx-2}];

Table[c[i]=0,{i,-3,nn+1,2}];

AA=AB/.λ->0;

B=((AB-AA)/λ);

alfx=1;
rrx=10000;

BN=B/.{α->alfx,rr->rrx};

AAN=N[AA/.{α->alfx,rr->rrx}];

We are exceedingly happy that Mathematica finds these all and sorts them.  The last one is the one that we want.

eigs=Eigenvalues[{AAN,-BN}]

alfx=1.02056;
rrx=5772.22*1.5;

BN=B/.{α->alfx,rr->rrx};

AAN=N[AA/.{α->alfx,rr->rrx}];

Eigenvalues[{AAN, -BN}]

Here is the growth curve in Im[λ] vs. α space

Initialize functions to use in the calculations.  If this is done, all that is needed is that alfx and rrx be assigned values, and then "lambda" gives the eigenvalues

Clear[BN,Binv,AAN,AZ,zeigs,Neig,lambda,rr,rrx]

BN:=N[B/.{α->alfx,rr->rrx}];
AAN:=N[AA/.{α->alfx,rr->rrx}];
zeigs:=Eigenvalues[{AAN, -BN}]
Neig:=Length[zeigs];
lambda:=zeigs[[Neig]];

rrx=10000;

Do[{alfx=0.5+i .05;alfz[i]=alfx;imlamb[i]=Im[lambda];relamb[i]=Re[lambda];},{i,0,25}]

imtt=Table[{alfz[i],imlamb[i]},{i,0,25}]

ListPlot[imtt,Joined->True]

The break occurs when a different mode becomes the last eigenvalue.  (Thus we should be more careful about the sorting of eigenvalues than Mathematica is doing).

rett=Table[{alfz[i],relamb[i]},{i,0,25}];

ListPlot[rett,Joined->True]

rrx=5772.22,  right at neutral stability!

rrx=5772.22;

Do[{alfx=0.5+i .05;alfz[i]=alfx;imlamb[i]=Im[lambda];relamb[i]=Re[lambda];},{i,0,25}]

imtt=Table[{alfz[i],imlamb[i]},{i,0,25}]

ListPlot[imtt,Joined->True]

rett=Table[{alfz[i],relamb[i]},{i,0,25}]

ListPlot[rett,Joined->True]

Here is the stability boundary in α-Re space

This module employs the secant method (adapted from Press, Teukolsky, Vetterling and Flannery  Numerical Recipes in Fortran, 2nd Ed. Cambridge Press, 1992, p 351.)  (Since there is no Numerical Recipes in Mathematica.

Do[{alfx=0.9+i .01;dx=10;xacc=.01;rrx=6800;x1=rrx;fl=Im[lambda];rrx=6500;f=Im[lambda];x2=rrx;If[Abs[fl]<Abs[f],{rtsec=x1;x1=x2;swap=fl;fl=f;f=swap},{xl=x1;rtsec=x2}];While[Abs[dx]>xacc,{dx=(((xl-rtsec) f)/(f-fl));xl=rtsec;fl=f;rtsec=rtsec+dx;rrx=rtsec;f=Im[lambda];Print["rrx= ",rrx," f= ",f]}];rr[i]=rrx;alfz[i]=alfx;},{i,0,19}];

ans=Table[{rr[i],alfz[i]},{i,0,19}]

ListPlot[ans,Joined->True,AxesLabel->{"Re","α"},PlotLabel->"Stability boundary"]