Infinite Rotating Disk

>

Consider the steady fluid flow generated when the infinite plane z = 0 , immersed in a Newtonian viscous fluid, rotates about the axis r = 0 with a constant angular velocity omega . The viscous drag of the rotating surface creates a swirling flow toward the disk. The motion is characterized in terms of the pressure, p , and the three components of the velocity, v^r , v^theta , v^z , in cylindrical coordinates. The radial symmetry of this problem eliminates theta as a independent variable. Thus, with rho and nu denoting the density and kinematic viscosity of the fluid and writing partial derivatives as subscripts, the equations of continuity and conservation of momentum reduce to [6]:

1/r (r*v^r)[r] + (v^z)[z] = 0 ,

v^r*(v^r)[r]+v[z](v^r)[z] - 1/r (v^theta)^2 = -1/rho p[r] + nu ( (v^r)[rr] + 1/r (v^r)[r] + (v^r)[zz] - v^r/(r^2) ) (9)

v^r*(v^theta)[r]+v^z*(v^theta)[z] + 1/r v^r*v^theta = nu ( (v^theta)[rr] + 1/r (v^theta)[r] + (v^theta)[zz] - v^theta/(r^2) )

v^r*(v^z)[r]+v^z*(v^z)[z] = -1/rho p[z] + nu ( (v^z)[rr] + 1/r (v^z)[r] + (v^z)[zz] )

The boundary conditions are chosen (for all r >0) to enforce no slippage at the interface with the disk:

(v^r)(r,0) = (v^z)(r,0) = 0,

(v^theta)(r,0) = r*omega , (10)

p(r,0) = 0

and no non-axial viscous effect in the far-field:

limit((v^r)(r,z),z = infinity) = 0 ,

limit((v^theta)(r,z),z = infinity) = 0 , (11)

limit((v^z)[z](r,z),z = infinity) = 0 .

>

Similarity solutions to this equation were found by Karman [9], who noted that each of v^r/r , v^theta/r , v^z/r , and p depends only on the distance from the disk, z . The system of PDEs reduces to a system of ODEs with the introduction of z *= z*sqrt(omega/nu) as the dimensionless independent variable for the dimensionless functions F , G , H , and P defined according to

v^r = r*omega*F(z*`*`) ,

v^theta = r*omega*G(z*`*`) , (12)

v^z = sqrt(omega*nu)*H(z*`*`) ,

p = rho*omega*nu*P(z*`*`) ,

The substitution of (12) into (9), (10), and (11) can be completed in a manner analogous to that used in the Blasius problem. These steps are omitted here in the interest of space; the resulting BVP is:

H*`'` = -2*F ,

F*`''` = F^2-G^2+F*`'`*H ,

G*`''` = 2*F*G+H*G*`'`

P*`'` = -H*H*`'`+H*`''` , (13)

F(0) = H(0) = P(0) = 0 ,

G(0) = 1 ,

limit(F(z*`*`),z*`*` = infinity) = limit(G(z*`*`),z*`*` = infinity) = 0 .

Note that (13) can be solved as a two-point BVP for F , G , H . Once this solution is known, the pressure equation can be integrated to yield P =- 1/2 H^2 + H '. In addition, note that any physically realistic solution must have F , G , P > 0 and H < 0 for all z * [6, p. 164].

>

Assuming the boundary condition at z* = infinity can be applied at a finite distance from the disk, shoot can be used to obtain an approximate solution to this system. First, introduce the first derivatives of F and G as new dependent variables:

> restart;

> with(Shoot):

> with(plots):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the assigned name shoot now has a global binding

Warning, the name changecoords has been redefined

> FNS:={ F(Z), G(Z), H(Z), Fp(Z), Gp(Z) }:

and reformulate the system as a first-order system of ODEs:

> ODE:={diff(H(Z),Z) = -2*F(Z),
diff(F(Z),Z) = Fp(Z),
diff(Fp(Z),Z) = -G(Z)^2+F(Z)^2+Fp(Z)*H(Z),
diff(G(Z),Z) = Gp(Z),
diff(Gp(Z),Z) = 2*F(Z)*G(Z)+H(Z)*Gp(Z)}:

with initial conditions:

> IC:={ F(0)=0, G(0)=1, H(0)=0,

> Fp(0)=alpha, Gp(0)=beta }:

and boundary conditions (again, z *=10 can be shown to be in the far-field):

> BC:={ F(10)=0, G(10)=0 }:

Note that, as in the first example, the boundary condition at infinity has been moved to a finite position. The value z *=10 is chosen, as before, on the basis that this is already in the far field. (In fact, truncating the computations at z *=7$ yields essentially the same solution.) Note also that since there are two boundary conditions at the second boundary point, there are two parameters to be determined in the shooting method ( i.e. , m[2] = 2 ). This can be handled by shoot without modification:

> infolevel[shoot]:=1:

> S:=shoot( ODE, IC, BC, FNS,

> [alpha=0.51, beta=-0.62] ):

newton: Step # 1

newton: Parameter values : alpha = .51, beta = -.62

newton: Step # 2

newton: Parameter values : alpha = .5100033193, beta = -.6161487475

newton: Step # 3

newton: Parameter values : alpha = .5102117353, beta = -.6159070528

newton: Step # 4

newton: Parameter values : alpha = .5102135912, beta = -.6159096495

>

> P:=Z->-H(Z)^2/2-2*F(Z):

> C:=odeplot(S,[ [Z,F(Z)], [Z,G(Z)],

> [Z,H(Z)], [Z,-P(Z)] ], 0..10 ):

> seq( op(0,lhs(f))=rhs(f), f=select(has,S(3/4),{F,G,H}) ):

> [%,P=subs(%, H^2+2*F)]:

> L:=textplot( [ seq([3/4,rhs(f)+1/45,lhs(f)],f=%) ],

> font=[TIMES,ITALIC,18], align=ABOVE ):

> display( [C,L], labels=[`z*`,``],

> title=`Infinite Disk: `

> ||`alpha=0.51, beta=-0.62` );

[Maple Plot]

This plot confirms that the solution satisfies the original boundary conditions at z *= infinity . (In fact, the solution obtained when the computational domain is truncated at z *=7.) Compare the solution in the preceding plot with the solution obtained when the starting values of the control parameters are alpha = .50 and beta = -.61 -- each just 0.01 less than the first attempt:

> S2:=shoot( ODE, IC, BC, FNS,

> [alpha=0.5, beta=-0.61] ):

newton: Step # 1

newton: Parameter values : alpha = .5, beta = -.61

newton: Step # 2

newton: Parameter values : alpha = .5181506110, beta = -.5834922896

newton: Step # 3

newton: Parameter values : alpha = .5178653205, beta = -.5900697917

newton: Step # 4

newton: Parameter values : alpha = .5019889184, beta = -.5940268798

newton: Step # 5

newton: Parameter values : alpha = .5040798104, beta = -.5927957600

newton: Step # 6

newton: Parameter values : alpha = .5042719728, beta = -.5927086094

newton: Step # 7

newton: Parameter values : alpha = .5042737798, beta = -.5927079136

The different values of alpha and beta suggest that the shooting method has returned a different solution. Is this also a solution to the problem?

> C:=odeplot(S2,[ [Z,F(Z)], [Z,G(Z)],

> [Z,H(Z)], [Z,-P(Z)] ], 0..12 ):

> seq( op(0,lhs(f))=rhs(f), f=select(has,S2(3/4),{F,G,H}) ):

> [%,P=subs(%, H^2+2*F)]:

> L:=textplot( [ seq([3/4,rhs(f)+1/45,lhs(f)],f=%) ],

> font=[TIMES,ITALIC,18], align=ABOVE ):

> display( {C,L}, labels=[`z*`,``],

> title=`Infinite Disk: `

> ||`alpha=0.50, beta=-0.61` );

#read `/k1/faculty/meade/www/maple/worksheets/shoot`; [Maple Plot]

Note that this solution does not satisfy the aforementioned sign constraints for F , G , H , and P for a physical solution to (13). Furthermore, while the boundary conditions at z *=10 are satisfied, it is extremely unlikely that this solution will satisfy the boundary conditions at z *= infinity . This is a spurious solution [6, p. 167]. Automating the detection of spurious solutions is extremely difficult. In some instances, this difficulty can be avoided by the use of one of the other solution techniques for boundary value problems mentioned previously.

>