Rotating disk flow and mass transfer

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://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 this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version:  4/12/16

Summary

The notebook examines the flow caused by a rotating disk.

References

The references for this NoteBook are:
F. M. White
Viscous Fluid Flow
McGraw-Hill, 1974
pp 164-170

V. G. Levich
PhysicoChemical Hydrodynamics
Prentice Hall
1962
pp 60-72

Use the full Navier Stokes equations.

The similarity variables given by White and Co. allow the transformation of the pde to an ode.

The interesting result is that the mass transfer coefficient is not a function of position.

Need the NS equations in cylindrical coordinates

The equations are (note the typo in White’s book)

continuityeq = 1/r  D[r vr[r,θ,z],r] + D[vz[r,θ,z],z]

vreq =  vr[r,θ,z] D[vr[r,θ,z],r] + vz[r,θ,z] D[vr[r,θ,z],z] - vθ[r,θ,z]^2/r + 1/ρ D[p[r,θ,z],r]  - ν   (  D[vr[r,θ,z],{r,2}]+1/r  D[vr[r,θ,z],{r,1}]+ D[vr[r,θ,z],{z,2}]-vr[r,θ,z]/r^2)

vθeq =  vr[r,θ,z] D[vθ[r,θ,z],r] + vz[r,θ,z] D[vθ[r,θ,z],z] +vθ[r,θ,z]  vr[r,θ,z]/r   - ν   (  D[vθ[r,θ,z],{r,2}]+1/r  D[vθ[r,θ,z],{r,1}]+ D[vθ[r,θ,z],{z,2}]-vθ[r,θ,z]/r^2)

vzeq =  vr[r,θ,z] D[vz[r,θ,z],r] + vz[r,θ,z] D[vz[r,θ,z],z] +1/ρ D[p[r,θ,z],z]      - ν   (  D[vz[r,θ,z],{r,2}]+1/r  D[vz[r,θ,z],{r,1}]+ D[vz[r,θ,z],{z,2}])

Stream function formulation to satisfy continuity

We will define a stream function from the continuity equation, define a similarity variable involving x and y.  Then transform the equations.

continuityeq/.{
Derivative[a1_,a2_,a3_][vr][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3}))(r ω F[zz])) Sqrt[ω/ν]^a3, Derivative[a1_,a2_,a3_][vθ][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3})) (r ω g[zz]))Sqrt[ω/ν]^a3,
Derivative[a1_,a2_,a3_][vz][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3})) (Sqrt[ω ν] h[zz]))Sqrt[ω/ν]^a3,
Derivative[a1_,a2_,a3_][p][r,θ,z]:> ( (∂_({r,a1},{θ,a2},{zz,a3})) p[r,θ,z])Sqrt[ω/ν]^a3,
vr[r,θ,z]-> r ω F[zz], vθ[r,θ,z]-> r ω g[zz], vz[r,θ,z]-> Sqrt[ω ν] h[zz],p[r,θ,z]->ρ ω ν pp[zz]}

continuityresult= Expand[PowerExpand[%/ω]]

vreq /.{
Derivative[a1_,a2_,a3_][vr][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3}))(r ω F[zz])) Sqrt[ω/ν]^a3, Derivative[a1_,a2_,a3_][vθ][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3})) (r ω g[zz]))Sqrt[ω/ν]^a3,
Derivative[a1_,a2_,a3_][vz][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3})) (Sqrt[ω ν] h[zz]))Sqrt[ω/ν]^a3,
Derivative[a1_,a2_,a3_][p][r,θ,z]:> ( (∂_({r,a1},{θ,a2},{zz,a3})) ρ ω ν pp[zz])Sqrt[ω/ν]^a3,
vr[r,θ,z]-> r ω F[zz], vθ[r,θ,z]-> r ω g[zz], vz[r,θ,z]-> Sqrt[ω ν] h[zz],p[r,θ,z]->ρ ω ν pp[zz]}

PowerExpand[%]

vrresult= Expand[%/r/ω^2]

vθeq

vθeq /.{
Derivative[a1_,a2_,a3_][vr][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3}))(r ω F[zz])) Sqrt[ω/ν]^a3, Derivative[a1_,a2_,a3_][vθ][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3})) (r ω g[zz]))Sqrt[ω/ν]^a3,
Derivative[a1_,a2_,a3_][vz][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3})) (Sqrt[ω ν] h[zz]))Sqrt[ω/ν]^a3,
Derivative[a1_,a2_,a3_][p][r,θ,z]:> ( (∂_({r,a1},{θ,a2},{zz,a3})) ρ ω ν pp[zz])Sqrt[ω/ν]^a3,
vr[r,θ,z]-> r ω F[zz], vθ[r,θ,z]-> r ω g[zz], vz[r,θ,z]-> Sqrt[ω ν] h[zz],p[r,θ,z]->ρ ω ν pp[zz]}

PowerExpand[%]

vθresult= Expand[%/r/ω^2]

vzeq

vzeq /.{
Derivative[a1_,a2_,a3_][vr][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3}))(r ω F[zz])) Sqrt[ω/ν]^a3, Derivative[a1_,a2_,a3_][vθ][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3})) (r ω g[zz]))Sqrt[ω/ν]^a3,
Derivative[a1_,a2_,a3_][vz][r,θ,z]:>  ((∂_({r,a1},{θ,a2},{zz,a3})) (Sqrt[ω ν] h[zz]))Sqrt[ω/ν]^a3,
Derivative[a1_,a2_,a3_][p][r,θ,z]:> ( (∂_({r,a1},{θ,a2},{zz,a3})) (ρ ω ν pp[zz]))Sqrt[ω/ν]^a3,
vr[r,θ,z]-> r ω F[zz], vθ[r,θ,z]-> r ω g[zz], vz[r,θ,z]-> Sqrt[ω ν] h[zz],p[r,θ,z]->ρ ω ν pp[zz]}

PowerExpand[%]

vzresult = Expand[%/Sqrt[ν]/ω^(3/2)]

The equations to be solved, now only a function of the rescaled z are:

continuityresult

vrresult

vθresult

vzresult

So the clever way to create the system of first order equations is to use the continuity results to get rid of the derivatives of h.

vθresult1= vθresult/.Derivative[a1_][h][zz]:>  ((∂_({zz,a1-1}))(-2F[zz]))

vzresult1=vzresult/.Derivative[a1_][h][zz]:>  ((∂_({zz,a1-1}))(-2F[zz]))

Now following White, we can create a system of first order equations using
y1 = h
y2 = F’
y3=F
y4= G’
y5=G
y6 = P

eq1 =y1'[x] == -2 y3[x]

((F[zz])^2)-((g[zz])^2)+h[zz] F'[zz]-F''[zz]

eq2 = y2'[x] == y3[x]^2 - y5[x]^2+ y1[x] y2[x]

eq3 = y3'[x] == y2[x]

2 F[zz] g[zz]+h[zz] g'[zz]-g''[zz]

eq4 = y4'[x] == 2 y3[x] y5[x] + y1[x]y4[x]

eq5 = y5'[x] == y4[x]

-2 F[zz] h[zz]+2 F'[zz]+pp'[zz]

eq6 = y6'[x] == 2 y3[x] y1[x] - 2 y2[x]

We can solve this system with a shooting method.  Here is a pretty good one that uses NDSolve for the integrations.

Here is the numerical solution

Initialize some variables.

finit = {{.7},{-.7}};  (*initial condition vector for the 2 unknown derivatives, F'[0] and G'[0]*) 
youter=6; (* outer value of η *)
deltv = .00001;(* increment on the initial guesses used to generate the Jacobian*)

If you like to know how to pick out a number from a vector, look here.

finit

We need to run through the integration once to get a first values for the error.

The equations are from White’s book..

zz=NDSolve[
{y1'[x]==-2 y3[x],
y2'[x]==-y5[x] y5[x] + y3[x] y3[x] + y2[x] y1[x],
y3'[x]== y2[x],
y4'[x] == 2 y3[x] y5[x] + y1[x] y4[x],
y5'[x] == y4[x],
y6'[x] == 2 y3[x] y1[x] - 2 y2[x],
y1[0]==0,
y2[0]==finit[[1,1]],
y3[0]==0,
y4[0] ==finit[[2,1]],
y5[0] == 1,
y6[0] == 0},
{y1[x],y2[x],y3[x],y4[x],y5[x],y6[x]},{x,0,youter}];eps1=-(zz/.x->youter)[[1,3,2]];
eps2=-(zz/.x->youter)[[1,5,2]];
Print[ "Vr=  ",eps1,"   f''[0]=  ",finit[[1,1]],"  Gerror= ", eps2," G'[0] =   ", finit[[2,1]]];

-(zz/.x->youter)[[1,5,2]]

The first two values for the error are, eps1 (y3==> F, radial velocity) and eps2 (y5==>G, circumferential velocity), which should go to 0 as z->infinity.

eps1

eps2

Here is the main iteration loop.

While[Abs[eps1]+Abs[eps2]>.0000001,
(*first calculate the change due to changing the r velocity gradient guess*)
zz=NDSolve[
{y1'[x]==-2 y3[x],
y2'[x]==-y5[x] y5[x] + y3[x] y3[x] + y2[x] y1[x],
y3'[x]== y2[x],
y4'[x] == 2 y3[x] y5[x] + y1[x] y4[x],
y5'[x] == y4[x],
y6'[x] == 2 y3[x] y1[x] -2 y2[x],
y1[0]==0,
y2[0]==finit[[1,1]]+deltv,
y3[0]==0,
y4[0] ==finit[[2,1]],
y5[0] == 1,
y6[0] == 0},
{y1[x],y2[x],y3[x],y4[x],y5[x],y6[x]},{x,0,youter}];
eps12=-(zz/.x->youter)[[1,3,2]] ;
eps22=-(zz/.x->youter)[[1,5,2]] ;
(*calculate the change due to changing the θ velocity guess*)
zz=NDSolve[
{y1'[x]==-2 y3[x],
y2'[x]==-y5[x] y5[x] + y3[x] y3[x] + y2[x] y1[x],
y3'[x]== y2[x],
y4'[x] == 2 y3[x] y5[x] + y1[x] y4[x],
y5'[x] == y4[x],
y6'[x] == 2 y3[x] y1[x] - 2 y2[x],
y1[0]==0,
y2[0]==finit[[1,1]],
y3[0]==0,
y4[0] ==finit[[2,1]]+deltv,
y5[0] == 1,
y6[0] == 0},
{y1[x],y2[x],y3[x],y4[x],y5[x],y6[x]},{x,0,youter}];
eps13=-(zz/.x->youter)[[1,3,2]] ;
eps23=-(zz/.x->youter)[[1,5,2]] ;
(*calculate the Jacobian*)
alf={{(eps12-eps1)/deltv,(eps13-eps1)/deltv},      
          {(eps22-eps2)/deltv,(eps23-eps2)/deltv}};
(*calculate the new guesses*)
deltvnew = -Inverse[alf] . {{eps1},{eps2}};
finitnew = finit + deltvnew;
(*update everything and save the old numbers*)
finitold= finit;
eps1old = eps1;
eps2old = eps2;
finit=finitnew;
(*run the new numbers to see if the solution has converged*)
zz=NDSolve[
{y1'[x]==-2 y3[x],
y2'[x]==-y5[x] y5[x] + y3[x] y3[x] + y2[x] y1[x],
y3'[x]== y2[x],
y4'[x] == 2 y3[x] y5[x] + y1[x] y4[x],
y5'[x] == y4[x],
y6'[x] == 2 y3[x] y1[x] - 2 y2[x],
y1[0]==0,
y2[0]==finit[[1,1]],
y3[0]==0,
y4[0] ==finit[[2,1]],
y5[0] == 1,
y6[0] == 0},
{y1[x],y2[x],y3[x],y4[x],y5[x],y6[x]},{x,0,youter}];
eps1=-(zz/.x->youter)[[1,3,2]];
eps2=-(zz/.x->youter)[[1,5,2]];
Print[ "Uerror=  ",eps1,"   F'[0]=  ",finit[[1,1]],"  Terror= ", eps2,"  G'[0] =   ", finit[[2,1]]];]

zvelocity = Plot[y1[x]/.zz[[1]],{x,0,youter},PlotRange->All, PlotLegends->{"Z velocity"},PlotStyle->Red]

rvelocity = Plot[y3[x]/.zz[[1]],{x,0,youter},PlotRange->All, PlotLegends->{"r velocity"},PlotStyle->Green]

θvelocity = Plot[y5[x]/.zz[[1]],{x,0,youter},PlotRange->All, PlotLegends->{"θ velocity"},PlotStyle->Blue]

pressure = Plot[y6[x]/.zz[[1]],{x,0,youter},PlotRange->All, PlotLegends->{"pressure"},PlotStyle->Orange]

Show[%,%%,%%%,%%%%]

Now use this for the mass transfer rate

This is easy to solve numerically.

concentrationeq =-Sc H[ζ] Θ'[ζ]+Θ''[ζ]

Now we use the solution from above to get the concentration field.  Since the flow field is not affected by heat transfer, I just use the known values of F’’[0] and G’[0] for the flow field and check that the concentration has converged.

Here are the (converged) initial values for F’ and G’

finit[[1,1]]

finit[[2,1]]

Clear[youter]

youter=70;

Sc=.01;
hppint = .5;
eps=1;

While[Abs[eps]>.0001,
zz=NDSolve[{y1'[x]==-2 y3[x],
y2'[x]==-y5[x] y5[x] + y3[x] y3[x] + y2[x] y1[x],
y3'[x]== y2[x],
y4'[x] == 2 y3[x] y5[x] + y1[x] y4[x],
y5'[x] == y4[x],
y6'[x] == 2 y3[x] y1[x] - 2 y2[x],
y7'[x] == y8[x],
y8'[x] ==  Sc y1[x]y8[x],
y1[0]==0,
y2[0]==finit[[1,1]],
y3[0]==0,
y4[0] ==finit[[2,1]],
y5[0] == 1,
y6[0] == 0,
y7[0] ==1, y8[0] == hppint},{y1[x],y2[x],y3[x],y4[x],y5[x],y6[x],y7[x],y8[x]},{x,0,youter}];xz=(zz/.x->youter)[[1,7,2]];eps=xz;hppint= (hppint/(1-eps));Print["U = ",xz,"    error=  ",eps,"   f''[0]=  ",hppint];]

This is the “z” velocity

Plot[y1[x]/. zz[[1]],{x,0,youter},PlotRange->All,AxesLabel->{"ζ","H(ζ)"}, PlotLegends->{"Sc = .1"},PlotStyle->Dashed]

This is concentration for a specific Sc.

Here are some different Sc values

Show[anssc1,ansscp1,anssc10,anssc100]

Let’s see what we can find out about scaling

y8’[0] == ∂Θ/∂ζ

D ∂Θ/∂ζ Sqrt[ω/ν] == k 

Sh = k Sqrt[ω/ν] /D ~ .6204 Sc^(1/3)

So  = D ∂Θ/∂ζ Sqrt[ω/ν]  Sqrt[ω/ν] /D ~ .62 Sc^1/3

So we see that we need to scale the ζ with Sc^-1/3 to make the concentration profiles collapse

Here is the scaled concentration to see how well it matches the Sc^1/3 prediction for a boundary layer.

anssc100=ParametricPlot[{x  100^(1/3), y7[x]/. zz[[1]]},{x,0,youter/(100)^(1/3)},PlotRange->All,AxesLabel->{"ζ","Θ(ζ)"}, PlotLegends->{"Sc = 100"},PlotStyle->Dashed, AspectRatio->1]

We see that the concentration profiles begin to collapse for Sc > 10 or so, although they look better for >1000

Show[anssc1000,anssc10000,anssc100,anssc10,anssc1,PlotLabel->{ "Profiles showing Sc^1/3 scaling"}]