Stability of rotating flow

This notebook has been written in Mathematica by 

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

Mark.J.McCready.1@nd.edu
http://chemeprof.com

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

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

Version:  5/2/18

Analysis of the Taylor-Couette instability following the discussion in Leal: Advanced Transport Phenomena, Cambridge, pp 829. and Drazin and Reid:  Hydrodynamic Stability (also Cambridge).

First consider the steady state solution.

For an inner cylinder R1 and rotation rate, Ω1 and an outer cylinder R2 and rotation rate, Ω2, the steady solution for rotating flow is

vtheta0 = (Ω2 R2^2 - Ω1 R1^2)/(R2^2-R1^2) r + 1/r (Ω1-Ω2) R1^2 R2^2/(R2^2-R1)

Plot[vtheta0/.{R1->1,R2->2,Ω1->-1,Ω2->1},{r,1,2}]

This centrifugal force is balanced by a pressure gradient as can be seen in the “v r” equation.

dpdr = ρ vtheta0^2/r

Integrate[ρ ( c1^2 r + 2 c1 c2/r + c2^2/r^3),r,Assumptions->{R1>0, R2>0, Im[Ω1] ->0, Im[Ω2] -> 0, r>0}]

Integrate[ρ ( c1^2 r + 2 c1 c2/r + c2^2/r^3),{r,R1,r},Assumptions->{R1>0, R2>0, Im[Ω1] ->0, Im[Ω2] -> 0, r>0}]

Start of stability analysis

vθ0 = C1 r + C2/r

This section leaves the equations in dimensional form

Here is the r equation

rmom = ρ ( D[vr[r,z,t],t]  - vθ[r,z,t]^2/r  )+ D[p[r,z,t],r]- 
μ(D[1/r D[r vr[r,z,t],r],r] + D[vr[r,z,t],{z,2}])

We first linearize

rtemp1=rmom/.{ Derivative[a1_,a2_,a3_][vr][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3}))(  ϵ vr1[r,z,t])) , Derivative[a1_,a2_,a3_][vθ][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (vθ0+ ϵ vθ1[r,z,t] )),
Derivative[a1_,a2_,a3_][vz][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (vz0+ ϵ vz1[r,z,t] )),
Derivative[a1_,a2_,a3_][p][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (p0[r]+ ϵ p1[r,z,t] )),vθ[r,z,t]-> (vθ0 + ϵ vθ1[r,z,t]),vr[r,z,t]-> ( ϵ vr1[r,z,t]),vz[r,z,t]-> (  ϵ vz1[r,z,t])}

rtemp2=Expand[rtemp1]

rtemp3 = Collect[rtemp2,ϵ]

At lowest order, the pressure balances the centripetal force

zeroorder = Coefficient[rtemp3,ϵ,0]

The first order disturbance equation is

vreq1 =Coefficient[rtemp3,ϵ,1]

Do the same for the θ direction

θmom = ρ ( D[vθ[r,z,t],t] + vr[r,z,t] D[ vθ[r,z,t],r]  +vr[r,z,t] vθ[r,z,t]/r  + vz[r,z,t] D[ vθ[r,z,t],z])+ D[p[r,z,t],Θ]- 
μ(D[1/r D[r vθ[r,z,t],r],r] + D[vθ[r,z,t],{z,2}])

θtemp1 =θmom/.{ Derivative[a1_,a2_,a3_][vr][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3}))(  ϵ vr1[r,z,t])) , Derivative[a1_,a2_,a3_][vθ][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (vθ0+ ϵ vθ1[r,z,t] )),
Derivative[a1_,a2_,a3_][vz][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (vz0+ ϵ vz1[r,z,t] )),
Derivative[a1_,a2_,a3_][p][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (p0[r]+ ϵ p1[r,z,t] )),vθ[r,z,t]-> (vθ0 + ϵ vθ1[r,z,t]),vr[r,z,t]-> ( ϵ vr1[r,z,t]),vz[r,z,t]-> (  ϵ vz1[r,z,t])}

θtemp2 =Expand[θtemp1]

θtemp3 = Collect[θtemp2,ϵ]

The rotating flow is an exact solution to the θ equation at 0 order

Coefficient[θtemp3,ϵ,0]

The linearized disturbance equation is

vθeq1 =Coefficient[θtemp3,ϵ,1]

zmom = ρ ( D[vz[r,z,t],t] + vr[r,z,t] D[ vz[r,z,t],r]   + vz[r,z,t] D[ vz[r,z,t],z])+ D[p[r,z,t],z]- 
μ(1/r D[r D[ vz[r,z,t],r],r] + D[vz[r,z,t],{z,2}])

Now do z

ztemp1 =zmom/.{ Derivative[a1_,a2_,a3_][vr][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3}))(  ϵ vr1[r,z,t])) , Derivative[a1_,a2_,a3_][vθ][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (vθ0+ ϵ vθ1[r,z,t] )),
Derivative[a1_,a2_,a3_][vz][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (vz0+ ϵ vz1[r,z,t] )),
Derivative[a1_,a2_,a3_][p][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (p0[r]+ ϵ p1[r,z,t] )),vθ[r,z,t]-> (vθ0 + ϵ vθ1[r,z,t]),vr[r,z,t]-> ( ϵ vr1[r,z,t]),vz[r,z,t]-> (  ϵ vz1[r,z,t])}

ztemp2 =Expand[ztemp1]

ztemp3= Collect[ztemp2,ϵ]

Coefficient[ztemp3,ϵ,0]

The disturbance equation is

vzeq1=Coefficient[ztemp2,ϵ,1]

finally the continuity equation

continuity= 1/r   D[ r vr[r,z,t],r]    + D[ vz[r,z,t],z]

contemp1 = continuity/.{ Derivative[a1_,a2_,a3_][vr][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3}))(  ϵ vr1[r,z,t])) , Derivative[a1_,a2_,a3_][vθ][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (vθ0+ ϵ vθ1[r,z,t] )),
Derivative[a1_,a2_,a3_][vz][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (vz0+ ϵ vz1[r,z,t] )),
Derivative[a1_,a2_,a3_][p][r,z,t]:>  ((∂_({r,a1},{z,a2},{t,a3})) (p0[r]+ ϵ p1[r,z,t] )),vθ[r,z,t]-> (vθ0 + ϵ vθ1[r,z,t]),vr[r,z,t]-> ( ϵ vr1[r,z,t]),vz[r,z,t]-> (  ϵ vz1[r,z,t])}

contemp2 =Expand[contemp1]

contemp3 =Collect[contemp2,ϵ]

Coefficient[contemp3,ϵ,0]

conteq1 = Coefficient[contemp3,ϵ,1]

Here we nondimensionalize the equations.

There is nothing profound to be gained by doing this.  However, if you wish to examine the inviscid problem, then it is convenient to get the Reynolds number into the equations and then take the limit of Re-> ∞.

conteq2 =conteq1/.{Derivative[a1_,a2_,a3_][vr1][r,z,t]:> (∂_({s,a1},{ζ,a2},{τ,a3}))ur1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),  
                                        Derivative[a1_,a2_,a3_][vθ1][r,z,t]:> (∂_({s,a1},{ζ,a2},{τ,a3}))uθ1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),
                                       Derivative[a1_,a2_,a3_][vz1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))uz1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),  
                                                              vr1[r,z,t]->(Subscript[R, 1])Ω1 ur1[s,ζ,τ],         
                                                              vθ1[r,z,t]->(Subscript[R, 1])Ω1 uθ1[s,ζ,τ],                       
                                                              vz1[r,z,t]->(Subscript[R, 1])Ω1 uz1[s,ζ,τ],
Derivative[a1_,a2_,a3_][p1][r,z,τ]:>(∂_({s,a1},{ζ,a2},{τ,a3}))P[s,ζ,τ] (R^(-a1)) ρ ((Subscript[R, 1])Ω1)^2 ,r->(Subscript[R, 1]) s}

conteq3 = Expand[s conteq2/Ω1]

vreq2 = vreq1/.{Derivative[a1_,a2_,a3_][vr1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))ur1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),         
                               Derivative[a1_,a2_,a3_][vθ1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))uθ1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),
                               Derivative[a1_,a2_,a3_][vz1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))uz1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),  
                                                     vr1[r,z,t]->(Subscript[R, 1])Ω1 ur1[s,ζ,τ],
                                                     vθ1[r,z,t]->(Subscript[R, 1])Ω1 uθ1[s,ζ,τ],
                                                     vz1[r,z,t]->(Subscript[R, 1])Ω1 uz1[s,ζ,τ],
Derivative[a1_,a2_,a3_][p1][r,z,t]:>(∂_({s,a1},{ζ,a2},{τ,a3}))P[s,ζ,τ] ((Subscript[R, 1])^(-a1)) ρ ((Subscript[R, 1])Ω1)^2 ,r->(Subscript[R, 1]) s}

vreq3 =Expand[vreq2/ρ/((Ω1)^2) /(Subscript[R, 1])]

vreq4 = Expand[vreq3/.{μ-> 1/re ρ Ω1 (Subscript[R, 1]^2), C1->C1p Ω1, C2->C2p Ω1 (Subscript[R, 1]^2)}]

vθeq2 = vθeq1/.{               Derivative[a1_,a2_,a3_][vr1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))ur1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),         
                               Derivative[a1_,a2_,a3_][vθ1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))uθ1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),
                               Derivative[a1_,a2_,a3_][vz1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))uz1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),  
                                                     vr1[r,z,t]->(Subscript[R, 1])Ω1 ur1[s,ζ,τ],
                                                     vθ1[r,z,t]->(Subscript[R, 1])Ω1 uθ1[s,ζ,τ],
                                                     vz1[r,z,t]->(Subscript[R, 1])Ω1 uz1[s,ζ,τ],
Derivative[a1_,a2_,a3_][p1][r,z,t]:>(∂_({s,a1},{ζ,a2},{τ,a3}))P[s,ζ,τ] ((Subscript[R, 1])^(-a1)) ρ ((Subscript[R, 1])Ω1)^2 ,r->(Subscript[R, 1]) s}

vθeq3 =Expand[vθeq2/ρ/((Ω1)^2) /(Subscript[R, 1])]

vθeq4 = Expand[vθeq3/.{μ-> 1/re ρ Ω1 (Subscript[R, 1]^2), C1->C1p Ω1, C2->C2p Ω1 (Subscript[R, 1]^2)} ]

vzeq2 = vzeq1/.{Derivative[a1_,a2_,a3_][vr1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))ur1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),         
                               Derivative[a1_,a2_,a3_][vθ1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))uθ1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),
                               Derivative[a1_,a2_,a3_][vz1][r,z,t]:>  (∂_({s,a1},{ζ,a2},{τ,a3}))uz1[s,ζ,τ](Subscript[R, 1])Ω1 ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3),  
                                                     vr1[r,z,t]->(Subscript[R, 1])Ω1 ur1[s,ζ,τ],
                                                     vθ1[r,z,t]->(Subscript[R, 1])Ω1 uθ1[s,ζ,τ],
                                                     vz1[r,z,t]->(Subscript[R, 1])Ω1 uz1[s,ζ,τ],
Derivative[a1_,a2_,a3_][p1][r,z,t]:>(∂_({s,a1},{ζ,a2},{τ,a3}))P[s,ζ,τ] ((Subscript[R, 1])^(-a1)) ((Subscript[R, 1])^(-a2)) ((Ω1)^a3) ρ ((Subscript[R, 1])Ω1)^2 ,r->(Subscript[R, 1]) s}

vzeq3 =Expand[vzeq2/ρ/((Ω1)^2) /(Subscript[R, 1])]

vzeq4 = Expand[vzeq3/.μ-> 1/re ρ Ω1 (Subscript[R, 1]^2)]

Now we start the normal mode analysis

conteq3 /.{ Derivative[a1_,a2_,a3_][ur1][s,ζ,τ]:>((∂_({r,a1},{z,a2},{t,a3}))( ur1h[r] Cos[α z] Exp[σ t])) ,              
                      Derivative[a1_,a2_,a3_][uθ1][s,ζ,τ]:>((∂_({r,a1},{z,a2},{t,a3})) (uθ1h[r] Cos[α z] Exp[σ t] )),
                      Derivative[a1_,a2_,a3_][uz1][s,ζ,τ]:> ((∂_({r,a1},{z,a2},{t,a3})) (uz1h[r] Sin[α z] Exp[σ t] )),
                          Derivative[a1_,a2_,a3_][P][s,ζ,τ]:> ((∂_({r,a1},{z,a2},{t,a3})) (p0[r]+ ϵ p1h[r,z,t] )),     
                                            uθ1[s,ζ,τ]-> (                                       uθ1h[r] Cos[α z] Exp[σ t] ),
                                            ur1[s,ζ,τ]-> (                                       ur1h[r] Cos[α z] Exp[σ t] ),
                                            uz1[s,ζ,τ]-> (                                       uz1h[r] Sin[α z] Exp[σ t] ),s->r}

This let’s us eliminate uz

uzans=Solve[%==0,uz1h[r]]

Here we have uz in terms of ur for later elimination.

uzhsubs=uz1h[r]/.uzans[[1]]

take the r and z equations next

crossr =vreq4/.{ Derivative[a1_,a2_,a3_][ur1][s,ζ,τ]:>((∂_({r,a1},{z,a2},{t,a3}))( ur1h[r] Cos[α z] Exp[σ t])) ,              
                      Derivative[a1_,a2_,a3_][uθ1][s,ζ,τ]:>((∂_({r,a1},{z,a2},{t,a3})) (uθ1h[r] Cos[α z] Exp[σ t] )),
                      Derivative[a1_,a2_,a3_][uz1][s,ζ,τ]:> ((∂_({r,a1},{z,a2},{t,a3})) (uz1h[r] Sin[α z] Exp[σ t] )),
                          Derivative[a1_,a2_,a3_][P][s,ζ,τ]:> ((∂_({r,a1},{z,a2},{t,a3})) (p0[r]+ ϵ p1h[r,z,t] )),     
                                            uθ1[s,ζ,τ]-> (                                       uθ1h[r] Cos[α z] Exp[σ t] ),
                                            ur1[s,ζ,τ]-> (                                       ur1h[r] Cos[α z] Exp[σ t] ),
                                            uz1[s,ζ,τ]-> (                                       uz1h[r] Sin[α z] Exp[σ t] ),s->r}

crossz = vzeq4/.{ Derivative[a1_,a2_,a3_][ur1][s,ζ,τ]:>((∂_({r,a1},{z,a2},{t,a3}))( ur1h[r] Cos[α z] Exp[σ t])) ,              
                      Derivative[a1_,a2_,a3_][uθ1][s,ζ,τ]:>((∂_({r,a1},{z,a2},{t,a3})) (uθ1h[r] Cos[α z] Exp[σ t] )),
                      Derivative[a1_,a2_,a3_][uz1][s,ζ,τ]:> ((∂_({r,a1},{z,a2},{t,a3})) (uz1h[r] Sin[α z] Exp[σ t] )),
                          Derivative[a1_,a2_,a3_][P][s,ζ,τ]:> ((∂_({r,a1},{z,a2},{t,a3})) (p0[r]+ ϵ p1h[r,z,t] )),     
                                            uθ1[s,ζ,τ]-> (                                       uθ1h[r] Cos[α z] Exp[σ t] ),
                                            ur1[s,ζ,τ]-> (                                       ur1h[r] Cos[α z] Exp[σ t] ),
                                            uz1[s,ζ,τ]-> (                                       uz1h[r] Sin[α z] Exp[σ t] ),s->r}

Cross differentiate to remove pressure

cross1 =D[crossr,z]- D[crossz,r]

cross2 = Expand[cross1/Exp[σ t]/Sin[α z]]

This could have been done more elegantly...

cross3 = cross2/.{Derivative[a1_][uz1h][r]:>D[uzhsubs,{r,a1}] ,uz1h[r]->  uzhsubs }

We get what is presumably Leal’s eq. 12-125

rzeqsubs = Expand[-cross3 re α]

Coefficient[rzeqsubs, re]

Now do the same transformations on the θ equation.

vθtemp1=vθeq4/.{ Derivative[a1_,a2_,a3_][ur1][s,ζ,τ]:>((∂_({r,a1},{z,a2},{t,a3}))( ur1h[r] Cos[α z] Exp[σ t])) ,              
                      Derivative[a1_,a2_,a3_][uθ1][s,ζ,τ]:>((∂_({r,a1},{z,a2},{t,a3})) (uθ1h[r] Cos[α z] Exp[σ t] )),
                      Derivative[a1_,a2_,a3_][uz1][s,ζ,τ]:> ((∂_({r,a1},{z,a2},{t,a3})) (uz1h[r] Sin[α z] Exp[σ t] )),
                          Derivative[a1_,a2_,a3_][P][s,ζ,τ]:> ((∂_({r,a1},{z,a2},{t,a3})) (p0[r]+ ϵ p1h[r,z,t] )),     
                                            uθ1[s,ζ,τ]-> (                                       uθ1h[r] Cos[α z] Exp[σ t] ),
                                            ur1[s,ζ,τ]-> (                                       ur1h[r] Cos[α z] Exp[σ t] ),
                                            uz1[s,ζ,τ]-> (                                       uz1h[r] Sin[α z] Exp[σ t] ),s->r}

vθtemp2 = Expand[vθtemp1/Exp[σ t]/Cos[α z]]

Here is Leal’s eq 12-126.  (It matches.)

vθsubs = Expand[- re vθtemp2]

Some further transformations would be needed but let’s just skip ahead

Here is the thin gap approximation from Leal

leal147a = D[(D[ur[s],{s,2}]- α^2 ur[s]),{s,2}]- α^2 (D[ur[s],{s,2}]- α^2 ur[s]) -σ (D[ur[s],{s,2}]- α^2 ur[s]) -  (1+β s)uθ[s]

leal147b = D[uθ[s],{s,2}]- α^2 uθ[s] -σ uθ[s] == -tt α^2 ur[s]

We want to get a single ode in uθ to solve numerically.   First get an equation for ur.

Solve[leal147b,ur[s]]

ur1 = ur[s]/.%[[1]]

leal147a/.{ur[s]->ur1, ur'[s]-> D[ur1,s], ur''[s]->D[ur1,{s,2}], ur'''[s]->D[ur1,{s,3}],ur''''[s]->D[ur1,{s,4}]}

ode=Expand[- % tt α^2]

Make a notation change that will allow use of Chebyshev polynomials in the solution.

goveeq=ode/.{  Derivative[a1_][uθ][s]:>  Derivative[a1][vθ][ξ,i] 1^a1,uθ[s]:>vθ[ξ,i],s-> ξ}

Now we start the numerical spectral analysis.

Amazingly, 12 modes is enough.   Of course the original solution by Chandrasekar ultimately is the {1,1} element as an approximation to the secular determinant!

n1 = 12;

order 0 coefficient

coeff0 =Coefficient[goveeq, vθ[ξ,i]]- ξ /vθ[ξ,i]Coefficient[goveeq,ξ]

For vθ[ξ,i]].  Note that the flow domain is (0,1) so we adjust the argument of the Chebyshev function to account for this.  We also need to adjust the weighting function and coefficient.

exp0 = Sum[a[i] ChebyshevT[i,2ξ-1],{i,0,n1}]

We know what this should give, just the a[i]’s so this would be just a check of the “4” and the domain change to (0,1).

z0d=Table[a0d[i]=((4 (∫_0^1)exp0 ChebyshevT[i,2 ξ-1] (1/(Sqrt[1-(((2 ξ-1))^2)]))\[DifferentialD]ξ)/(π c[i])),{i,0,n1,1}]/.c[0]->2/.Table[c[i]->1,{i,1,n1,1}]

It works!

z0d=Table[a0d[i]= a[i],{i,0,n1,1}]

order 0 ξ coefficient

coeff0ξ=Coefficient[goveeq,ξ]/vθ[ξ,i]

exp0ξ = Sum[a[i] ChebyshevT[i,2 ξ-1],{i,0,n1}]ξ

z0dξ=Table[a0dξ[i]=((4 (∫_0^1)exp0ξ ChebyshevT[i,2 ξ-1] (1/(Sqrt[1-(((2 ξ-1))^2)]))\[DifferentialD]ξ)/(π c[i])),{i,0,n1,1}]

z0dξ1 = z0dξ/.c[0]->2/.Table[c[i]->1,{i,1,n1,1}]

order 2 derivative

coeff2 =Coefficient[goveeq, vθ''[ξ,i]]

exp2 =D[exp0,{ξ,2}]

z2d=Table[a2d[i]=((4 (∫_0^1)exp2 ChebyshevT[i,2ξ-1] (1/(Sqrt[1-(((2 ξ-1))^2)]))\[DifferentialD]ξ)/(π c[i])),{i,0,n1}]

z2d1=Expand[z2d/.c[0]->2/.Table[c[i]->1,{i,1,n1,1}]]

order 4 derivative

coeff4 =Coefficient[goveeq, vθ''''[ξ,i]]

coeff4 =-3 (α^2)-2 σ

exp4 =D[exp0,{ξ,4}]

z4d=Table[a4d[i]=((4 (∫_0^1)exp4 ChebyshevT[i,2 ξ-1] (1/(Sqrt[1-(((2 ξ-1))^2)]))\[DifferentialD]ξ)/(π c[i])),{i,0,n1,1}]

z4d1=Expand[z4d/.c[0]->2/.Table[c[i]->1,{i,1,n1,1}]]

order 6 derivative

coeff6 =Coefficient[goveeq, vθ''''''[ξ,i]]

exp6 =D[exp0,{ξ,6}]

Here we also check the timing!

Timing[z6d=Table[a6d[i]=((4 (∫_0^1)exp6 ChebyshevT[i,2 ξ-1] (1/(Sqrt[1-(((2 ξ-1))^2)]))\[DifferentialD]ξ)/(π c[i])),{i,0,n1,1}]]

z6d1=z6d/.c[0]->2/.Table[c[i]->1,{i,1,n1,1}]

This is the array of n1 + 1 equations from the complete substituted ODE that have been “orthogonalized” with each of the Chebyshev modes.

dcoefs=z6d1+z4d1*coeff4+z2d1*coeff2 + z0d*coeff0   + z0dξ1 coeff0ξ

Here are the boundary conditions.

vexpand=(∑_({{i=0}, {Δ i=1}})^n1)a[i] ChebyshevT[i,2ξ-1]

bc1=vexpand/.{ξ->0}

bc2=vexpand/.{ξ->1}

bc3=Expand[(D[vexpand,{ξ,2}]  )/.{ξ->0}]

bc4=Expand[(D[vexpand,{ξ,2}]  )/.{ξ->1}]

bc5=Expand[ (D[vexpand,{ξ,3}] - (α^2 +σ)D[vexpand ,ξ])/.{ξ->0}]

bc6=Expand[ (D[vexpand,{ξ,3}] -  (α^2 +σ)D[vexpand ,ξ])/.{ξ->1}]

We solve the first 4 BCs for the top 4 coefficients.  We use this result to eliminate the top 4 “a[i]’s” for all subsequent calculations.  We save BC5 and BC6 and replace them in the coefficient matrix from the ODE expansion.  This is done because α and the eigenvalue is in the last two boundary conditions which prevents solving linearly for any lower “a[i]’”s,

bcsolve=Solve[{bc1==0,bc2==0,bc3==0, bc4==0},{a[n1],a[n1-1],a[n1-2],a[n1-3]}]

To allow enforcement of the 6 boundary conditions, we off the six highest equations and then paste in bc5 and bc6.  The other 4 boundary conditions provide a substitution list for the 4 highest “a[i]’s”.

coefs1=(Table[dcoefs[[i]],{i,1,n1-5}])

Here we attach the two boundary conditions.

coefs2 = Join[coefs1,{bc5,bc6}]

Here we substitution for the 4 highest a’s.

coefs=Expand[coefs2 /.bcsolve[[1]]]

There should be n1 + 1 - 4 equations....

Length[coefs]

Here we construct L (λ) * x =0,  So we have to pick out all of the L coefficients.  Note that the “x” is the vector of “a[i]’s.]

AB=Table[Table[Coefficient[coefs[[j]],a[i]],{i,0,n1-4}],{j,1,Length[coefs]}]

AB={{(19428147200/3)-((40832000 (α^2))/3)+tt (α^2)-((17600 (α^4))/9)-(α^6)+(1/2) tt (α^2) β-((81664000 σ)/9)-((70400 (α^2) σ)/27)-2 (α^4) σ-((17600 (σ^2))/27)-(α^2) (σ^2),(1/4) tt (α^2) β,(640221675520/99)-((1344716800 (α^2))/99)-((568000 (α^4))/297)-((2689433600 σ)/297)-((2272000 (α^2) σ)/891)-((568000 (σ^2))/891),0,(622985216000/99)-((1293584384 (α^2))/99)-((495872 (α^4))/297)-((2587168768 σ)/297)-((1983488 (α^2) σ)/891)-((495872 (σ^2))/891),0,(60668968960/11)-((120193024 (α^2))/11)-((38752 (α^4))/33)-((240386048 σ)/33)-((155008 (α^2) σ)/99)-((38752 (σ^2))/99),0,(114440929280/33)-((203187200 (α^2))/33)-((52544 (α^4))/99)-((406374400 σ)/99)-((210176 (α^2) σ)/297)-((52544 (σ^2))/297)},{(1/2) tt (α^2) β,(275904921600/67)-((912384000 (α^2))/67)+tt (α^2)-((190080 (α^4))/67)-(α^6)+(1/2) tt (α^2) β-((608256000 σ)/67)-((253440 (α^2) σ)/67)-2 (α^4) σ-((63360 (σ^2))/67)-(α^2) (σ^2),(1/4) tt (α^2) β,(272424370176/67)-((896071680 (α^2))/67)-((175104 (α^4))/67)-((597381120 σ)/67)-((233472 (α^2) σ)/67)-((58368 (σ^2))/67),0,(246900326400/67)-((782622720 (α^2))/67)-((129600 (α^4))/67)-((521748480 σ)/67)-((172800 (α^2) σ)/67)-((43200 (σ^2))/67),0,(163590045696/67)-((465960960 (α^2))/67)-((61056 (α^4))/67)-((310640640 σ)/67)-((81408 (α^2) σ)/67)-((20352 (σ^2))/67),0},{(33810022400/3)-((78848000 (α^2))/3)-3520 (α^4)-((157696000 σ)/9)-((14080 (α^2) σ)/3)-((3520 (σ^2))/3),(1/4) tt (α^2) β,(1114180157440/99)-((2596894720 (α^2))/99)+tt (α^2)-((116480 (α^4))/33)-(α^6)+(1/2) tt (α^2) β-((5193789440 σ)/297)-((465920 (α^2) σ)/99)-2 (α^4) σ-((116480 (σ^2))/99)-(α^2) (σ^2),(1/4) tt (α^2) β,(1084719104000/99)-((2500198400 (α^2))/99)-((103552 (α^4))/33)-((5000396800 σ)/297)-((414208 (α^2) σ)/99)-((103552 (σ^2))/99),0,(105879961600/11)-((232775680 (α^2))/11)-((24576 (α^4))/11)-((465551360 σ)/33)-((32768 (α^2) σ)/11)-((8192 (σ^2))/11),0,(200970076160/33)-((394649600 (α^2))/33)-((11200 (α^4))/11)-((789299200 σ)/99)-((44800 (α^2) σ)/33)-((11200 (σ^2))/33)},{0,(190749081600/67)-((830269440 (α^2))/67)-((142560 (α^4))/67)-((553512960 σ)/67)-((190080 (α^2) σ)/67)-((47520 (σ^2))/67),(1/4) tt (α^2) β,(188445229056/67)-((816721920 (α^2))/67)+tt (α^2)-((733824 (α^4))/335)-(α^6)+(1/2) tt (α^2) β-((544481280 σ)/67)-((978432 (α^2) σ)/335)-2 (α^4) σ-((244608 (σ^2))/335)-(α^2) (σ^2),(1/4) tt (α^2) β,(171550310400/67)-((717373440 (α^2))/67)-((113280 (α^4))/67)-((478248960 σ)/67)-((151040 (α^2) σ)/67)-((37760 (σ^2))/67),0,(115489898496/67)-((430940160 (α^2))/67)-((273984 (α^4))/335)-((287293440 σ)/67)-((365312 (α^2) σ)/335)-((91328 (σ^2))/335),0},{7064780800-((69386240 (α^2))/3)-((7040 (α^4))/3)-((138772480 σ)/9)-((28160 (α^2) σ)/9)-((7040 (σ^2))/9),0,(232830402560/33)-((2285731840 (α^2))/99)-((233248 (α^4))/99)-((4571463680 σ)/297)-((932992 (α^2) σ)/297)-((233248 (σ^2))/297),(1/4) tt (α^2) β,(226990489600/33)-((2209464320 (α^2))/99)+tt (α^2)-((250880 (α^4))/99)-(α^6)+(1/2) tt (α^2) β-((4418928640 σ)/297)-((1003520 (α^2) σ)/297)-2 (α^4) σ-((250880 (σ^2))/297)-(α^2) (σ^2),(1/4) tt (α^2) β,(66954854400/11)-((207585280 (α^2))/11)-((20800 (α^4))/11)-((415170560 σ)/33)-((83200 (α^2) σ)/33)-((20800 (σ^2))/33),0,(43287838720/11)-((356249600 (α^2))/33)-((29312 (α^4))/33)-((712499200 σ)/99)-((117248 (α^2) σ)/99)-((29312 (σ^2))/99)},{0,(73574645760/67)-((638668800 (α^2))/67)-((47520 (α^4))/67)-((425779200 σ)/67)-((63360 (α^2) σ)/67)-((15840 (σ^2))/67),0,(72757149696/67)-((3149217792 (α^2))/335)-((257472 (α^4))/335)-((2099478528 σ)/335)-((343296 (α^2) σ)/335)-((85824 (σ^2))/335),(1/4) tt (α^2) β,(66762178560/67)-((565125120 (α^2))/67)+tt (α^2)-((80640 (α^4))/67)-(α^6)+(1/2) tt (α^2) β-((376750080 σ)/67)-((107520 (α^2) σ)/67)-2 (α^4) σ-((26880 (σ^2))/67)-(α^2) (σ^2),(1/4) tt (α^2) β,(46869774336/67)-((1751887872 (α^2))/335)-((211392 (α^4))/335)-((1167925248 σ)/335)-((281856 (α^2) σ)/335)-((70464 (σ^2))/335),0},{2523136000-17121280 (α^2)-((3520 (α^4))/9)-((34242560 σ)/3)-((14080 (α^2) σ)/27)-((3520 (σ^2))/27),0,(7560232960/3)-((564162560 (α^2))/33)-((118784 (α^4))/297)-((1128325120 σ)/99)-((475136 (α^2) σ)/891)-((118784 (σ^2))/891),0,(7385907200/3)-((548208640 (α^2))/33)-((168640 (α^4))/297)-((1096417280 σ)/99)-((674560 (α^2) σ)/891)-((168640 (σ^2))/891),(1/4) tt (α^2) β,2202009600-((158945280 (α^2))/11)+tt (α^2)-((43520 (α^4))/33)-(α^6)+(1/2) tt (α^2) β-((105963520 σ)/11)-((174080 (α^2) σ)/99)-2 (α^4) σ-((43520 (σ^2))/99)-(α^2) (σ^2),(1/4) tt (α^2) β,1495531520-((94289920 (α^2))/11)-((66496 (α^4))/99)-((188579840 σ)/33)-((265984 (α^2) σ)/297)-((66496 (σ^2))/297)},{-(4026880/9)-((3200 (α^2))/27)-((3200 σ)/27),(16727040/67)+((6400 (α^2))/67)+((6400 σ)/67),-(12049408/27)-((8960 (α^2))/81)-((8960 σ)/81),(81779712/335)+((5376 (α^2))/67)+((5376 σ)/67),-(11497472/27)-((7168 (α^2))/81)-((7168 σ)/81),(13977600/67)+((3584 (α^2))/67)+((3584 σ)/67),-(1042432/3)-((512 (α^2))/9)-((512 σ)/9),(39831552/335)+((1536 (α^2))/67)+((1536 σ)/67),-(1690112/9)-((640 (α^2))/27)-((640 σ)/27)},{(4026880/9)+((3200 (α^2))/27)+((3200 σ)/27),(16727040/67)+((6400 (α^2))/67)+((6400 σ)/67),(12049408/27)+((8960 (α^2))/81)+((8960 σ)/81),(81779712/335)+((5376 (α^2))/67)+((5376 σ)/67),(11497472/27)+((7168 (α^2))/81)+((7168 σ)/81),(13977600/67)+((3584 (α^2))/67)+((3584 σ)/67),(1042432/3)+((512 (α^2))/9)+((512 σ)/9),(39831552/335)+((1536 (α^2))/67)+((1536 σ)/67),(1690112/9)+((640 (α^2))/27)+((640 σ)/27)}}

n1=12

Even though we will go further below, you can get the neutral stability conditions by setting σ = 0 and then looking a Det[AB] == 0.

For the case of the inside cylinder rotating and the outside fixed, β=-1.

This is apparently right on the mark with the “exact” answers in Drazin and Reid at about 2*1694.5.

Table[{tt/2.,Det[AB/.{α->3.127,β->-1,σ->0}]},{tt,3380,3400,1}]

Likewise if β=0, both cylinders rotating, the Taylor number matches Drazin and Reid.

Table[{tt,Det[AB/.{α->3.117,β->0,σ->0}]},{tt,1700,1720,1}]

Now construct the quadratic eigenvalue problem in a form that can be converted into a generalized eigenvalue problem.  Wikipedia had a simple procedure for this.

First (σ^0)

A0=AB/.σ->0;

(σ^1)

A1 = Coefficient[AB, σ, 1];

(σ^2)

A2= Coefficient[AB, σ, 2];

Dimensions[A2]

MatrixForm[A2]

We need a square matrix of 0’s to attach to the left of -A2.

zeromatriz = Table[0,{i,1,n1-3},{j,1,n1-3}];

Attach the square matrix of 0' s to attach to the left of - A2.

zza1 =Transpose[Join[zeromatriz,Transpose[-A2]]];

It worked!

MatrixForm[zza1]

Now we need a square identity matrix on the left of a square matrix of 0’s.

eys = Transpose[Join[IdentityMatrix[n1-3],Table[0,{i,1,n1-3},{j,1,n1-3}]]];

We now attach the identity and 0 matrix to the bottom of the 0 ,-A2 matrix.  zza is the coefficient matrix for “λ” or in our case σ in the newly created generalized eigenvalue problem.

zza =Join[zza1,eys];

It looks like this .  Note that we have 2 rows of 0’s in the middle because the boundary conditions that we included in the coefficient matrix have σ to only the first power.  Fortunately, the Mathematica eigenvalue solver can handle this singularity.

MatrixForm[zza]

For the main matrix, we attach A0 to the left side of A1.

zz =Transpose[Join[Transpose[A0],Transpose[A1]]];

MatrixForm[zz]

We need a matrix of 0’s and an identity matrix to attach below the A0-A1 matrix.

eys2 =  Transpose[Join[Table[0,{i,1,n1-3},{j,1,n1-3}],IdentityMatrix[n1-3]]];

Here we attach them

zzx =Join[zz,eys2];

Now solve the eigenvalue problem numerically.

We choose conditions where we know (from the literature or above) close to the neutral point.

Here is just the inside cylinder rotating.

nzzx=(zzx/.{tt->2 1694.95, β->-1,α->3.127});

nzza=N[(zza/.{tt-> 2 1694.95,β->-1,α->3.127})];

ans =Eigenvalues[{nzzx,nzza}]

ListPlot[Transpose[{Im[ans],Re[ans]}],PlotStyle->Blue]

Here is a growth curve of the disturbance as a function of wavenumber.

subs2={tt-> 2 1694.95, β->-1}

Table[{α,Re[Eigenvalues[{zzx/.subs2,zza/.subs2}][[18]]]},{α,2,5,.05}]

ListPlot[%]

Let’s see if we can find another mode at higher rotation speed?

subs22={tt-> 4 2 1694.95, β->-1,α->4.1}

Eigenvalues[{zzx/.subs22,zza/.subs22}]

Eigenvalues[N[{zzx/.subs2,zza/.subs2}]]

subs2={tt-> 4 2 1694.95, β->-1}

Table[{α,Re[Eigenvalues[{zzx/.subs2,zza/.subs2}][[17]]]},{α,2,5,.05}]

ListPlot[%,PlotRange->All]

Here we have both cylinders rotating.

nzzx=(zzx/.{tt->1708, β->0,α->3.117});

nzza=N[(zza/.{tt-> 1708,β->0,α->3.117})];

ans =Eigenvalues[{nzzx,nzza}]

Here are all of the eigenvalues.

ListPlot[Transpose[{Im[ans],Re[ans]}],PlotStyle->Red]

Here is the transition in terms of Taylor number

subs={α->3.13, β->0}

Table[Eigenvalues[{zzx/.subs,zza/.subs}][[18]],{tt,1650,1740,10}]

subs3 ={ tt->1708,β->0}

Table[{α,Re[Eigenvalues[{zzx/.subs3,zza/.subs3}][[18]]]},{α,2,5,.05}]

The growth curve for β=0.

ListPlot[%]

We can calculate the eigenfunction for the most unstable mode.

We chose conditions close to neutral stability and substitute for all of the variables, including the eigenvalue.  This is coefficient matrix of L(λ) which with λ possessing the special value, has a solution for all of the “a’s” in terms of one of the a’s.  We want to pick a[0].

ABsubs = AB/.{tt-> 2 1694.95,β->-1,α->3.127,σ->-0.000050079276562195155`}

Now we construct a set of equations from this matrix that we wish to solve, L[λ=...]. x=0.

ABsubs.Table[a[i],{i,0,8}]

eqsforas = Table[%[[i]]==0,{i,1,9}]

It is so nice that Mathematica gives the result in terms of a[0]!

ansforas = Solve[eqsforas,Table[a[i],{i,0,n1}]]

exp0 = Sum[a[i] ChebyshevT[i,2ξ-1],{i,0,n1-4}];

eigenfunction =(exp0/.ansforas [[1]])/.{a[10]->0, a[0]->1, a[11]->0, a[12]->0,a[9]->0}

We get the function!   Well maybe more like  yawn...

Plot[eigenfunction,{ξ,0,1}]

We can calculate the eigenfunction for the most unstable mode.

We chose conditions close to neutral stability and substitute for all of the variables, including the eigenvalue.  This is coefficient matrix of L(λ) which with λ possessing the special value, has a solution for all of the “a’s” in terms of one of the a’s.  We want to pick a[0].

ABsubs = AB/.{tt->1708, β->0, α->3.1,σ-> 0.0012902216419651525`}

Now we construct a set of equations from this matrix that we wish to solve, L[λ=...]. x=0.

ABsubs.Table[a[i],{i,0,8}]

eqsforas = Table[%[[i]]==0,{i,1,9}]

It is so nice that Mathematica gives the result in terms of a[0]!

ansforas = Solve[eqsforas,Table[a[i],{i,0,n1}]]

exp0 = Sum[a[i] ChebyshevT[i,2ξ-1],{i,0,n1-4}];

eigenfunction =(exp0/.ansforas [[1]])/.{a[10]->0, a[0]->1, a[11]->0, a[12]->0,a[9]->0}

We get the function!   Well maybe more like  yawn...

Plot[eigenfunction,{ξ,0,1}]