Similarity Solution for Axisymmetric Flow in a Tube with an Accelerating Surface Velocity
​Author: Housam Binous

The following BVP describes the similarity solution for axisymmetric flow in a tube with an accelerating surface velocity ​​
1
r

r
r

r
f'
r
-β=ℛe
2
f'
r
-
f
r
f"
r
-
f'
2
r

​​with​​
f(0)-f'(0)=0,f(1)=0,f'(1)=1
​​The constant β is given by​​
β=f'(1)-ℛe-f"(1)+f'''(1)
​​The function f(r) is related to the velocity field:
v=v
e
r
+u
e
z
​​​
u=zf'(r)/r,v=-f(r)/r
​​Brady and Acrivos(1981) were the first to show that the solution space exhibits a turning point when β is plotted versus ℛe, near ℛe=10.25. That is, for ℛe>10.25, there is no steady axisymmetric solution until ℛe >186.​Note there is an analytically solution for ℛe=0. The solution to​​
1
r

r
r

r
f'
r
-β=0
​​subject to BC: f(1)=0 gives​​
f
0
(r)=
1
2
(
4
r
-
2
r
),β=8
​​We use chebyshev orthogonal collocation in order to discretize the nonlinear BVP. With the above scheme and arc-length continuation, we were able to reproduce Fig. 7 of Ref. 1 or Figure 7 of Ref. 2 (i.e., a plot of β vs. ℛe). The turning point occurs at
ℛe≈10.25
. Also along the lower branch
β=-3
6.6 at
ℛe=6
, which is in good agreement with the value found by Brady and Acrivos [1] as well as Higgins and Binous [2].​​References:​​[1] J. F. Brady and A. Acrivos, Steady flow in a channel or tube with an accelerating surface velocity. An exact solution to the Navier-Stokes equations with reverse flow, J. Fluid Mech., 112, 127-150 (1981). DOI: http://dx.doi.org/10.1017/S0022112081000323.​[2] Higgins, B. G., & Binous, H. (2010). A Simple Method for Tracking Turning Points in Parameter Space. JOURNAL OF CHEMICAL ENGINEERING OF JAPAN, 43(12), 1035–1042. DOI: https://doi.org/10.1252/jcej.10we122.​
In[]:=
Np=20;
In[]:=
y[i_]:=Cos[iPi/Np]//N;
In[]:=
d[0,0]=(2Np^2+1)/6//N;​​d[Np,Np]=-(2Np^2+1)/6//N;​​d[i_,i_]:=-y[i]/(2(1-y[i]^2))//N;​​d[j_,k_]:=If[j!=k&&k!=0&&k!=Np&&j!=0&&j!=Np,(-1)^(j+k)/(y[j]-y[k])]//N;​​d[Np,0]=-1/2(-1)^Np;​​d[0,Np]=1/2(-1)^Np;​​d[j_,0]:=If[j!=Np&&j!=0,-1/2(-1)^j/(1-y[j])];​​d[j_,Np]:=If[j!=Np&&j!=0,1/2(-1)^(j+Np)/(1+y[j])];​​d[0,j_]:=If[j!=Np&&j!=0,2(-1)^j/(1-y[j])];​​d[Np,j_]:=If[j!=Np&&j!=0,-2(-1)^(j+Np)/(1+y[j])];​​D1=Table[d[i,j]//N,{i,0,Np},{j,0,Np}];D2=D1.D1;​​D3=D2.D1;
In[]:=
r[i_]:=(y[i-1]+1)/2
In[]:=
Table[y[i],{i,0,Np}]
Out[]=
{1.,0.987688,0.951057,0.891007,0.809017,0.707107,0.587785,0.45399,0.309017,0.156434,0.,-0.156434,-0.309017,-0.45399,-0.587785,-0.707107,-0.809017,-0.891007,-0.951057,-0.987688,-1.}
In[]:=
Table[r[i],{i,1,Np+1}]
Out[]=
{1.,0.993844,0.975528,0.945503,0.904508,0.853553,0.793893,0.726995,0.654508,0.578217,0.5,0.421783,0.345492,0.273005,0.206107,0.146447,0.0954915,0.0544967,0.0244717,0.00615583,0.}
In[]:=
1/rD[rD[f'[r]/r,r],r]//Simplify
Out[]=
′
f
[r]-r
′′
f
[r]+
2
r
(3)
f
[r]
3
r
In[]:=
eq[i_]:=1/r[i]^3(2(D1.Table[f[j],{j,1,Np+1}])[[i]]-r[i]4(D2.Table[f[j],{j,1,Np+1}])[[i]]+r[i]^28(D3.Table[f[j],{j,1,Np+1}])[[i]])==beta+Rey((2(D1.Table[f[j],{j,1,Np+1}])[[i]]/r[i])^2-f[i]/r[i](4(D2.Table[f[j],{j,1,Np+1}])[[i]]/r[i]-2(D1.Table[f[j],{j,1,Np+1}])[[i]]/r[i]^2))
In[]:=
BC={f[Np+1]==0,2(D1.Table[f[i],{i,1,Np+1}])[[Np+1]]==0,f[1]==0,2(D1.Table[f[i],{i,1,Np+1}])[[1]]==1};
In[]:=
eqsup=beta==2(D1.Table[f[i],{i,1,Np+1}])[[1]]-Rey-4(D2.Table[f[i],{i,1,Np+1}])[[1]]+8(D3.Table[f[i],{i,1,Np+1}])[[1]];
In[]:=
system=Flatten@Join[{Table[eq[i],{i,3,Np-1}],eqsup,BC}];
In[]:=
Length[system]
Out[]=
22
(*solutionatℛe=Reynoldsnumber=0*)
In[]:=
Rey=0;
In[]:=
unknown=Flatten[Join[{Table[{f[i],0.5},{i,1,Np+1}],{{beta,-20.5}}}],1]
Out[]=
{{f[1],0.5},{f[2],0.5},{f[3],0.5},{f[4],0.5},{f[5],0.5},{f[6],0.5},{f[7],0.5},{f[8],0.5},{f[9],0.5},{f[10],0.5},{f[11],0.5},{f[12],0.5},{f[13],0.5},{f[14],0.5},{f[15],0.5},{f[16],0.5},{f[17],0.5},{f[18],0.5},{f[19],0.5},{f[20],0.5},{f[21],0.5},{beta,-20.5}}
In[]:=
Length@unknown
Out[]=
22
In[]:=
FindRoot[system,unknown]
Out[]=
{f[1]5.42029×
-25
10
,f[2]-0.00606156,f[3]-0.0230037,f[4]-0.0473913,f[5]-0.0743949,f[6]-0.0988817,f[7]-0.116515,f[8]-0.124593,f[9]-0.122435,f[10]-0.111278,f[11]-0.09375,f[12]-0.073126,f[13]-0.0525583,f[14]-0.0344883,f[15]-0.0203378,f[16]-0.0104933,f[17]-0.00451774,f[18]-0.00148054,f[19]-0.000299254,f[20]-0.0000189464,f[21]0.,beta8.}
In[]:=
Table[1/2(r[i]^4-r[i]^2),{i,1,Np+1}]
Out[]=
{0.,-0.00606156,-0.0230037,-0.0473913,-0.0743949,-0.0988817,-0.116515,-0.124593,-0.122435,-0.111278,-0.09375,-0.073126,-0.0525583,-0.0344883,-0.0203378,-0.0104933,-0.00451774,-0.00148054,-0.000299254,-0.0000189464,0.}
(*arc-lengthcontinuationallowsonetoplotβvs.ℛe*)
In[]:=
Rey=.
In[]:=
system1=Flatten@Join[{Table[eq[i],{i,3,Np-1}],eqsup,BC}]/.beta->beta[s]/.Rey->Rey[s]/.f[1]f[1][s]/.f[2]f[2][s]/.f[3]f[3][s]/.f[4]f[4][s]/.f[5]f[5][s]/.f[6]f[6][s]/.f[7]f[7][s]/.f[8]f[8][s]/.f[9]f[9][s]/.f[10]f[10][s]/.f[11]f[11][s]/.f[12]f[12][s]/.f[13]f[13][s]/.f[14]f[14][s]/.f[15]f[15][s]/.f[16]f[16][s]/.f[17]f[17][s]/.f[18]f[18][s]/.f[19]f[19][s]/.f[20]f[20][s]/.f[21]f[21][s];
Auxeq=beta'[s]^2+Rey'[s]^2+Sum[f[i]'[s]^2,{i,1,Np+1}]==1;(*auxiliaryequation*)
In[]:=
Auxeq
Out[]=
2
′
beta
[s]
+
2
′
Rey
[s]
+
2
′
f[1]
[s]
+
2
′
f[2]
[s]
+
2
′
f[3]
[s]
+
2
′
f[4]
[s]
+
2
′
f[5]
[s]
+
2
′
f[6]
[s]
+
2
′
f[7]
[s]
+
2
′
f[8]
[s]
+
2
′
f[9]
[s]
+
2
′
f[10]
[s]
+
2
′
f[11]
[s]
+
2
′
f[12]
[s]
+
2
′
f[13]
[s]
+
2
′
f[14]
[s]
+
2
′
f[15]
[s]
+
2
′
f[16]
[s]
+
2
′
f[17]
[s]
+
2
′
f[18]
[s]
+
2
′
f[19]
[s]
+
2
′
f[20]
[s]
+
2
′
f[21]
[s]
1
In[]:=
init={Rey[0]==0,f[1][0]==5.420294834209781`*^-25,f[2][0]==-0.006061559927098614`,f[3][0]==-0.023003707754474047`,f[4][0]==-0.04739129078141423`,f[5][0]==-0.07439486317129335`,f[6][0]==-0.09888167382396222`,f[7][0]==-0.11651544951362046`,f[8][0]==-0.12459324509542181`,f[9][0]==-0.1224353861247772`,f[10][0]==-0.11127758171556887`,f[11][0]==-0.09374999999905077`,f[12][0]==-0.07312602141672996`,f[13][0]==-0.05255826182772709`,f[14][0]==-0.034488317574284395`,f[15][0]==-0.020337838976681143`,f[16][0]==-0.010493326175474276`,f[17][0]==-0.004517738874315583`,f[18][0]==-0.0014805370846029263`,f[19][0]==-0.00029925375429063676`,f[20][0]==-0.000018946401674099253`,f[21][0]==0.`,beta[0]==8.000000000095842`}
Out[]=
{Rey[0]0,f[1][0]5.42029×
-25
10
,f[2][0]-0.00606156,f[3][0]-0.0230037,f[4][0]-0.0473913,f[5][0]-0.0743949,f[6][0]-0.0988817,f[7][0]-0.116515,f[8][0]-0.124593,f[9][0]-0.122435,f[10][0]-0.111278,f[11][0]-0.09375,f[12][0]-0.073126,f[13][0]-0.0525583,f[14][0]-0.0344883,f[15][0]-0.0203378,f[16][0]-0.0104933,f[17][0]-0.00451774,f[18][0]-0.00148054,f[19][0]-0.000299254,f[20][0]-0.0000189464,f[21][0]0.,beta[0]8.}
In[]:=
inidev=Join[Table[f[i]'[0]==0.0,{i,1,Np+1}],{Rey'[0]==0.01,beta'[0]==0.01}]
Out[]=
{
′
f[1]
[0]0.,
′
f[2]
[0]0.,
′
f[3]
[0]0.,
′
f[4]
[0]0.,
′
f[5]
[0]0.,
′
f[6]
[0]0.,
′
f[7]
[0]0.,
′
f[8]
[0]0.,
′
f[9]
[0]0.,
′
f[10]
[0]0.,
′
f[11]
[0]0.,
′
f[12]
[0]0.,
′
f[13]
[0]0.,
′
f[14]
[0]0.,
′
f[15]
[0]0.,
′
f[16]
[0]0.,
′
f[17]
[0]0.,
′
f[18]
[0]0.,
′
f[19]
[0]0.,
′
f[20]
[0]0.,
′
f[21]
[0]0.,
′
Rey
[0]0.01,
′
beta
[0]0.01}
In[]:=
SYS=Join[system1,Flatten@{Auxeq,inidev,inidev}];
In[]:=
Length[SYS]
Out[]=
69
In[]:=
unknown1=Join[Table[f[i][s],{i,1,Np+1}],{Rey[s],beta[s]}];