by user

Category: Documents





Journal of Applied Mathematics and Physics (ZAMP)
Vol. 29, 1978
Birkh~iuser Verlag Basel
Free-Convection Boundary Layer on
an Isothermal Horizontal Cylinder
By D. B. Ingham, Dept. of Applied Mathematical Studies, The University of Leeds,
Leeds, England
1. Introduction
Solutions of the equations for the free-convection boundary layers have mainly
been confined to those body shapes for which the partial differential equations can
be reduced to ordinary differential equations by a suitable similarity transformation.
These solutions are reviewed by Ostrach [1].
The free-convection boundary layer on a two-dimensional, or axisymmetric body,
whose temperature T' is suddenly increased from that of the surrounding fluid at
time t' = 0 has been studied by Elliott [2]. He obtained a series solution in terms of
powers of the square of the time and this theory was applied to a horizontal circular
cylinder. The procedure for extrapolating to infinite time, as devised by Shanks [3]
and which has successfully been applied to several problems in fluid dynamics by
Van Dyke [4], was employed.
Merkin [5] used an accurate numerical step-by-step method to solve the steadystate boundary-layer equations on a heated horizontal circular cylinder. These results
showed that the heat transfer and skin friction as extrapolated by Elliott were in good
agreement with the exact numerical values especially at small values of 0, where 0 is
the angle measured from the downward vertical through the axis of the cylinder.
Hermann [6] obtained experimental results to this problem. The numerical results of
Merkin in the range 0 ~ < 0 < 150 ~ for the heat transfer and skin friction on the plate
agree well with these experimental results.
One serious deficiency in the numerical results of Merkin is that the solution violates the boundary condition of symmetry about the plane 0 = 180 ~ Obviously the use
of the boundary-layer equations does not allow us to be able to enforce a boundary
condition at two stations of 0 and Merkin solved the boundary-layer equation starting
with the analytical solution at 0 -= 0 ~ The small-time solution of Elliott does, however, satisfy the symmetry condition at 0 = 180 ~
Further the numerical solution of Merkin shows that the flow reaches 0 = 180 ~
without separating and that the boundary layer has finite thickness at this point. This
suggests that the two boundary layers which form on each side of the cylinder collide
at 0 = 180 ~ and are deflected through an angle of 90 ~ to form a buoyant plume.
The problem of the colliding boundary layers is similar to that encountered on a
D.B. Ingham
sphere which is rotating about a diameter, in an otherwise undisturbed fluid. This
problem has been investigated by many authors, for example Howarth [7], Nigam [8],
Banks [9, I0], Stewartson [11], Manohar [12], lllingworth [13], Benton [14] and Barrett
[15]. In this case fluid is sucked into the boundary layer from near the poles and ejected
at the equator. The collision of the boundary layers, which grow from the poles, at the
equator has been the subject of many postulates. Nigam proposed a solution in which
the boundary layers were empty of fluid at the equator but Banks and Stewartson
showed that this solution was not possible. Banks also attempted to use an inviscid
theory as proposed by Stewartson but concluded that the interaction at the equator
was a viscous problem. Smith and Duck [16] postulated that there is a large (compared
with the boundary layer thickness) recirculating region near the equator. This flow
has not been observed experimentally.
As in the rotating sphere problem the nature of the flow in the free convection
problem near 0 = 180 ~ has not yet been resolved. The main aim of this paper is to
develop a numerical method which will deal with the collision of the boundary layers
which occurs near 0 = 180 ~ and to investigate the generation of the buoyant plume in
this vicinity. A numerical technique is developed which enforces the symmetry
condition about 0 = 180 ~
2. Equations
Consider a circular cylinder of radius a which is at a constant temperature To. It
is placed such that its axis is horizontal in an infinite expanse of homogeneous viscous
fluid which is initially at rest and at a constant temperature To. The fluid is assumed
to be a Boussinesq fluid and the gravity force, g, is taken to act in the vertical direction.
At time t' = 0 the temperature T' of the cylinder is impulsively increased to Tw,
a constant. On the assumption that Tw - To = AT << To and the Grashof number
Gr >> 1 {where Gr = (gfiATaa)/v ~ and/? is the volumetric coefficient of expansion and
v is the kinematic viscosity of the fluid}, then the governing equations are the usual
incompressible free-convection time-dependent boundary-layer equations.
We take cylindrical polar coordinates (r', 0, z') with the axis of the cylinder at
r ' -- 0, 0 = 0 to be the downward vertical and assume all quantities are independent
ofz'. The motion may then be described by the dimensionless velocities (u, v, 0) which
are obtained by dividing the actual components by (t~gaAT) ll2. A dimensionless
vorticity ~ ( = [dimensional vorticity]/[~gAT/a]ll2), dimensionless temperature T
(-= [T' - To]/AT) and dimensionless time t ( = t'[~gAT/a] zl~) are also introduced.
Using the transformation ~: = in (r), where r = r'/a, then the velocity components
can be expressed as
u = e-~ a~
v = - e - e ~-~,
where ~b is the non-dimensional stream function ( = [dimensional stream function]/
Vol. 29, 1978
Free-Convection Boundary Layer
[{3gaSAT]~/~). The full Navier-Stokes equations governing the motion can now be
-~ + e- ~[--~-~
-~ "~l - GP'2 ~Of2 + ~02] e-2r
or) ,
0f - ~ J - Pr Gr t'2 \ e e ~ + -g0-if] e-~r
(o r
+ e- ~ sin 0-~- + cos 0 ~
_ Jo or
a~- + e
[gO 0f
o or]
~ 2 + ~O2~b = - g e 2r
where Pr is the Prandtl number which is the ratio of the kinematic viscosity to the
thermal diffusivity of the fluid.
Equations (2), (3) and (4) have to be solved subject to the boundary conditions
T= 1 onf=0,
~, T--~O
as f - + ~ ,
t _> 0,
O--if=0 o n 0 = 0 ~ and 180 ~, f >__0, t >__0.
In the early development of the flow a boundary layer of thickness proportional
to (t/Grl/2) lj2 forms on the cylinder and therefore the coordinate f normal to the
cylinder is transformed by
x = ~/k,
k = 2~/t/GP/4.
This leads to the following scalings for the dependent variables
~b = 88 ~/2 k3~F,
~ = 88 ~/2 kw,
T = ~.
Substituting expressions (9) and (10) into Eqns. (2), (3) and (4) gives
2 + 2x-b-~x - 2w
. 2[a~F &o
4 t -~- + 4t ~-~--ff~xx
0iF ~w) e_2~x
- 4 e-k~ (sin 00~2-b--xx+ kcos 0~-(0U2)- k2 ~ff0%~e_2~x,
D.B. Ingham
e -~k~ {a2f~
Pr kax 2 + 2 x - ~ x
at2 +
Of2 + O~F Of2] e_2k x
4t2k-g-ff ax
ax aO]
k2 a2"q e-2kx,
Pr a02
ax---T + k z a0 2 = -oJ e 2kx.
The power-series solution for small-time for the boundary-layer equations, as
obtained by Elliott, can be obtained by setting k = 0 and looking for a power-series
solution in t 2 of the Eqns. (11), (12) and (13).
The method of series truncation, as described by Collins and Dennis [17], has
been employed in this problem. We assume that, for all values of k,
T = ~ g,(x, t) sin nO,
o~ = ~ hn(x, t) sin nO,
f~ = fl(x, t) + ~ f~+l(x, t) cos nO.
These forms of solution have been taken so as to satisfy the boundary conditions (8).
Substituting expressions (14) into the boundary-layer equations, which are obtained
by setting k = 0 in Eqns. (11)7 (12) and (13), gives
x Oh~
Ox2 + 2 ~ -
2 h ~ - 4 t g---i~
= 2t 2
In - mlgl._ml Ox
kra = l,raCn
mh~ -~ glm-nl sgn (n - m) -
1 ~2fl +
Pr ax 2
Pr 0x
= 2t 2 2ngn ~
-- nlglm_n
fa = 1 r o m a n
mfm+l- ~ x g l m - n l
+ -~x Jm+z) J'
(m + n)gm+~ Ofm+~
m =1
4t afn+l
~ mhm Ogm+.]
n > 1,
= 2 t 2 | ~ ' km g,.
1 a2f.+a
2x 0f~+l
Pr 7x ~- + Pr ax
m=1 (m + n)gm+.
sgn (m - n) ,
n > 1
Vol. 29, 1978
Free-Convection Boundary Layer
Ox 2 =
n > 1
where sgn (m - n) denotes the sign of (m - n), with sgn (0) = 0 and 3i,s = 0 if
i :# j and 8~,j = 1 if i = j. In practice the terms in expressions (14) must be truncated
by setting identically zero all terms with subscript n which is greater than a prescribed
integer no. This defines a truncation of order no. Thus in Eqs. (15), (16) and (17) all
functions with a subscript greater than no are set zero.
The boundary conditions (5), (6) and (7) reduce to (for n > 1),
g" = Ox
f~ = 8~,1
on x = O, t >_ O,
f., h. -+ 0
as x -+ 0% t _> 0,
~g" > 0
as x --~ 0% t > 0,
Collins and Dennis found that the use of an integral condition rather than condition
(21) was easier to implement and they were equivalent. Thus integrating equation (17)
from x = 0 and x = m and using conditions (18) and (20) gives
h. dx
The 3 no Eqns., (15), (16) and (17), have now to be solved subject to the
boundary conditions (18), (19), (20) and (22) to give an approximate solution.
For simplicity we will henceforth restrict our attention to Pr = 1.
3. Numerical Method
To describe the method of approximating eq'laations (15) and (16) it is convenient
to suppress the subscripts and write a typical equation in the form
4t ~ t = A-5-~x2 + B-~x + C k + D,
where k, A , B, C and D are identified from Eqns. (15) and (16). Central differences in
space are used with a constant mesh size h. A Crank-Nicolson type approximation is
used by integrating over the time interval t - H to t, where H is the constant mesh
size in the t direction. Using these approximations Eqn. (23) can be approximated by
t)} + k ( x t)
t)} - 7 k ( x - h, t ) { A ( x , t) - 89
= 89
t) + k ( x , t - a){4(t - 8 9 - 2 7 A ( x , t - H ) - 8 9
t - n)}
+ 7 k ( x + h, t - H ) { A ( x , t - H ) + 89
t - H)} + 7 k ( x - h, t - H )
• { A ( x , t - H ) - 89
t - H)} + 8 9
t - H)
where 7 = H/(2h2).
The numerical solution of the Eqns. (24) for each function f , and h, is required
+ h, t ) { A ( x , t) + 89
x {4(t - 89
+ 2 7 A ( x , t) - 8 9
D.B. Ingham
subject to the boundary conditions (18), (19), (20) and (22). A direct method could be
used to solve the tri-diagonal matrices of the form (24) but as the equations are
non-linear and therefore require an iterative method of solution, an exact inversion
of the matrices is probably not the most efficient approach.
The solution at the start of the motion is found by putting t = 0 in Eqns. (11),
(12) and (13). The analytical solution to these equations, which are given by Elliott,
can be used to start the integration. Rather than use these analytical expressions the
governing equations at t = 0 were solved numerically using exactly the same techniques as those used when solving the time-dependent equations. This forms a check
on the numerical method. Using even a small number of mesh points in the x direction,
N, of 40 and h = 0.1 gave results which agreed to four decimal places with the
analytical results of Elliott.
In the initial stages of the time integrations only a few terms of the series (14)
are required to describe the motion but more terms become necessary as the timedependent integrations proceed. This is to be expected since the small-time expansion
in powers of t 2 introduces higher periodic terms in the expansions (14). The advantage
of the numerical method over the power-series method is that the numerical method
can be continued long after the power-series method fails to be valid.
The numerical integration is started with no = 2, H = 10 -3, h = 0.05 and N = 80.
The right-hand sides of the equations of the type (24) were evaluated for each o f f ,
and h, and one sweep through these equations was performed, for each value of n,
using the boundary conditions (18), (19), (20) and (22). As is usual in this type of
problem, see Collins and Dennis, an under relaxation factor, ~, has to be used. In this
case the minimum value of o~ used was 0.5 which contrasts with a value of 10 -2 as
used by Collins and Dennis. Simpson's Rule was used to satisfy the integral condition
(22), evaluate dg,/dx and thence to obtain g, for each value of n. In order to work to
high order accuracy throughout a special starting procedure as described by Milne [18]
was used when using Simpson's Rule. It was found that the value ofh,(0) was the last
quantity to converge so the test for convergence was that, for each value of n, the
difference between two successive iterations of h,(0) was less than some prescribed
value. This was taken to be 10 -5 in the results presented here. I f this condition is not
satisfied then the iterative process is repeated.
Errors arising from using finite differences in the t direction were kept small by
covering the step from t to t - H in first one and then two steps and insisting that
the difference between the two solutions be less than a prescribed value. This value
was taken to be 5 x 10-5 in the results presented in this paper. I f this accuracy is not
achieved then the step length H is halved and the process repeated. During the early
stages of the calculations a smaller time-step is required than at later times. Thus if
the solution converges very rapidly in going from time t to t + H, say in less than four
iterations, then a test is performed to see if the time-step can be doubled. This is done
by covering the step from t to t + 2 H in first one and then two steps and testing if the
difference between the two solutions is within the tolerance.
Vol. 29, 1978
Free-Convection Boundary Layer
In order to keep the calculations down to a minimum the integration was started
with no = 2 but at the completion of each time-step integration the values of each of
the functionsf~o(X), g~o(X) and h~o(X), for each value of x, were tested to see if any of
their values were greater than 10 -6 . If they were then no was increased by one. When
the calculations were terminated no was eighty.
The results presented in this paper correspond to h = 0.05. A complete calculation was performed with h = 0.1 and these results do not differ substantially from those
with h = 0.05. At small values of t the maximum value of x, xm say, may be taken as 4
but at larger values of t this value has to be increased. When the calculations were
terminated Xm was 12.
4. Results and Conclusions
The numerical integrations were performed up to t -~ 3.5 for Pr = 1. At this
stage there were 240 mesh points across the boundary layer and the order of the
truncation was 80. The integrations were terminated not because of any numerical
difficulties, as was the case experienced by Collins and Dennis, but because of the
excessive amount of computer time and storage required.
At t = 0 the numerical results were in excellent agreement with the analytical
results of Elliott. At t _-_ 0.5 the three term power expansion of Elliott is correct
to three decimal places. The agreement between the small-time analytical solution
and the numerical solution indicates that the numerical method is an appropriate
Figure 1 shows the temperature at 0 -- 30 ~ 90 ~ and 180 ~ at various values of t.
Also shown are the steady-state profiles as obtained by Merkin. These figures indicate
that at 0 = 30 ~ and 90 ~ the unsteady calculations approach the steady-state results
at large values of time. This approach towards the steady-state occurs for all values of
0 except near 180 ~ The discrepancy near 0 = 180 ~ is to be expected since the formulation of the unsteady and steady problems are different. In this work the correct
symmetry boundary conditions (8) have been used whereas in the steady-state solution
these boundary conditions have not been enforced.
Figure 2 shows the angular component of velocity at 0 = 30 ~ 90 ~ and 150 ~ for
several values of t and the steady-state results of Merkin. These again show the profiles
tending towards the steady-state results. In fact the smaller the value of 0 the quicker
are the steady-state results being approached as one would expect since the equations
are parabolic in the 0 direction. It is the temperature difference which drives the flow
and hence one would expect the temperature profiles to approach their steady-state
values quicker than the velocity. This is clearly shown in Figures 1 and 2.
At time t = 0 the analytical solution of Elliott illustrates that the radial velocity
at the outer edge of the boundary layer is inwards for 0 ~ <_ 0 < 90 ~ and outwards for
90 ~ < 0 _< 180 ~ In the steady-state solution of Merkin the radial velocity is always
inwards but, physically, we expect the fluid to be sucked into the boundary layer
D , B. I n g h a m
r'- a~ Gr o 25
Figure l(a)
The t e m p e r a t u r e at 8 = 30 ~ at v a r i o u s t i m e
i n t e r v a l s as a fu nction o f r a d i a l distance.
a J
F i g u r e l(b)
The t e m p e r a t u r e a t 0 = 90 ~ a t v a r i o u s t i m e intervals as a f u n c t i o n of r a d i a l distance.
( r 'a- a ~/ G r ....
F i g u r e l(c)
T h e t e m p e r a t u r e at 8 = 180 ~ at v a r i o u s t i m e intervals as a func t i on of r a d i a l di s t a nc e .
(r~-a~ G r ....
Figure 2(a)
~ a /
T h e angular c o m p o n e n t o f velocity at 0 = 30 ~ at v a r i o u s t i m e intervals as a f u n c t i o n o f radial
~t= 3"5
F i g u r e 2(5)
T h e a n g u l a r c o m p o n e n t o f velocity at O = 90 ~ at v a r i o u s t i m e intervals as a function o f radial
D. B. l n g h a m
t =3,0
(~)Gr ~
Figure 2(c)
The angular c o m p o n e n t of velocity at 0 = 150 ~ at various time intervals as a function of radial distance.
for all values of 0 except near 0 = 180 ~ Near 0 = 180 ~ we expect a thin buoyant
plume to form above the cylinder. Figure 3 shows the value of 0, say 0~, at which the
radial velocity at the outer edge of the boundary layer changes from being radially
inwards to being outwards as a function of time. This shows a narrowing of the region
of outflow. Also shown in Figure 3 is the maximum radial velocity at 0 = 180 ~ Urn,
as a function of time. This shows that as time increases this velocity increases rapidly.
-wo o
.130 ~
9o ~
Figure 3
The m a x i m u m radial velocity at the outer edge of the boundary layer at 0 = 180 ~ and the value of
0 at which there is no radial flow at the edge of the boundary layer as functions of time.
Vol. 29, 1978
Free-Convection Boundary Layer
~r r=l
9 Hermann's results
Figure 4
The variation of the heat transfer as a function of 0 at various time intervals.
These two results support the idea of two boundary layers originating at 0 = 0 ~
colliding at 0 = 180 ~ and forming a thin buoyant plume above the cylinder at 0 = 180 ~
Figures 4 and 5 show the variation of the non-dimensional heat transfer,
) ( ,=1
a G r - ' / 4 [~5;~
= 0.5t -112[~fln~163
-~- N" UJn+l
g/0] x=0 }
/~,=, Ox cos
and the non-dimensional skin friction on the plate,
for several values of time as functions of 0. The steady state values as obtained by Merkin
and the experimental results of Hermann are also shown. The heat transfer settles
down to the steady-state values earlier than the skin friction for all values of 0,
except near 0 = 180 ~ as one would have expected. At 0 = 180 ~ the limit as t ~ oo
obviously will not tend to the steady-state values as obtained by Merkin because in the
work presented here the boundary conditions (8) have been imposed. As a consequence
the skin friction at 0 = 180 ~ cannot have the non-zero value which is predicted by the
steady-state solution. Experimentally Hermann found, as expected, that the skin
D.B. Ingham
9 Het'moot~
~ 4 " f
Figure 5
The variation of the skin friction at the cylinder as a function of 0 at various time intervals.
friction at 0 = 180 ~ to be zero. The heat transfer at 0 = 180 ~ as predicted by the present
numerical calculation is seen to overshoot the steady-state value as predicted by
Merkin and keeps on decreasing with time. It is interesting to observe that the experimental results o f H e r m a n n predict that the heat transfer is zero at 0 = 180 ~ Although
the numerical integrations have not been performed for sufficiently large times both
Figures 5 and l(c) show definite signs of a rapid decrease in heat transfer at 0 = 180 ~
with time. A matching approach with a region near 0 = 180 ~ in which the temperature
is a constant might well be the appropriate mathematical method in dealing with this
U p to the time at which the integration was terminated no sign o f a recirculating
flow was observed. In fact, the skin friction near 0 -- 180 ~ increases as time increases
thus indicating no tendency for flow reversal.
The numerical method can easily be extended to include finite G r a s h o f numbers
but, because o f the large a m o u n t o f computer time required to obtain the results for
the boundary-layer equations, these calculations have not been performed here.
Part o f this work was done while the author was a visitor at the University o f
Western Ontario, London, Canada, and the author would like to thank Professor
S. C. R. Dennis for arranging this visit and making computing time available. This
part of the work was carried out with the assistance of a grant from the National
Research Council of Canada.
Vol. 29, 1978
Free-Convection Boundary Layer
[1] S. OSTRACH,'Laminar Flows with Body Forces' (Section F of Theory of Laminar Flows,
ed. F. K. Moore), Princeton University Press (1964).
[2] L. ELLIOTT, Free Convection on a Two-Dimensional or Axisymmetric Body, Quart. J. Mech.
Appl. Math. 23, 153 (1970).
[3] D. SHANKS, Non-Linear Transformation of Divergent and Slowly Convergent Sequences, J.
Math. Phys. 34, 1 (1955).
[4] M. VAN DYKE, Perturbation Methods in Fluid Dynamics, Academic Press, New York (1964).
[5] J. H. MERKIN, Free-Convection Boundary Layer on an Isothermal Horizontal Cylinder,
A.S.M.E.-A.I.Ch.E. Heat Transfer Conference, St. Louis, Mo. (1976).
[6] R. HERMANN, Wiirmeiibertragung bei freier Striimung am waagerechten Zylinder in zweiatomigen Gasen, Z. Ver. dt. Ing. 379, 1 (1936) {also N.A.C.A. Tech. Mere. 1366 (1954)}.
[7] L. HOWARTI-I,Note on the Boundary Layer on a Rotating Sphere, Phil. Mag. 42, 1308 (1951).
[8] S. D. NIGAM, Note on the Boundary Layer on a Rotating Sphere, Z. angew Math. Phys. 5, 151
[9] W. H. H. BANKS, The Boundary Layer on a Rotating Sphere, Quart. J. Mech. Appl. Math. 18,
443 (1965).
[10] W. H. H. BANKS, The Laminar Boundary Layer on a Rotating Sphere, Acta Mechanica 24, 273
[11] KSTEWARTSON,On Rotating Laminar Boundary Layers, Boundary Layer Research, I.U.T.A.M.
Symposium, p. 59 (1957).
[12] R. MANOHAR,The Boundary Layer on a Rotating Sphere, Z. angew Math. Phys. 18, 320 (1967).
[13] C. R. ILLINGWORTH,Boundary Layer Growth on a Spinning Body, Phil. Mag. 45, 1 (1953).
[14] E. R. BENTON, Laminar Boundary Layer on an Impulsively Started Rotating Sphere, J. Fluid
Mech. 23, 611 (1965).
[15] K. E. BARRETT, Or/the Impulsively Started Rotating Sphere, J. Fluid Mech. 27, 779 (1967).
[16] F. T. SMITH and P. W. DUCK, Separation of Jets or Thermal Boundary Layers from a Wall,
Quart. J. Mech. Appl. Math. 30, 145 (1977).
[17] W. M. COLLINS and S. C. R. DENNIS, Flow Past an Impulsively Started Circular Cylinder,
J. Fluid Mech. 60, 105 (1973).
[18] W. E. MILNE, Numerical Solution of Differential Equations, p. 48 (1953).
The free convection on a horizontal circular cylinder whose temperature is suddenly increased
is studied at large Grashof number. An accurate numerical method is described for determining
the solution of the time-dependent boundary-layer equations. The development of the various
physical properties of the flow are calculated and compared with their values at small and large
times as obtained from the previously obtained analytical solutions.
Cette 6tude consid~re la convection libre &nombre de Grashof 61ev6 au dessus d'un cylindre
horizontal de section circulaire lorsqu'on augmente brusquement sa temp6rature. On pr6sente
une m6thode num6rique qui permet de d6terminer avec pr6cision la solution des 6quations
instationnaires de la couche limite. Le d6veloppement des caract6ristiques physiques de l'6coulement est calcul6 et les r6sultats sont compar6s avec ceux qui ont 6t6 d6j~t obtenus analytiquement
pour les petites et grandes valeurs de la variable du temps.
(Received: February 27, 1978; revised: September 11, 1978)
Fly UP