An assessment of the accuracy of numerical solutions to the graetz problem

INT. COMM. HR~TMASS TRANSFER 0735-1933/85 $3.00 + .00 Vol. 12, pp. 209-218, 1985 @Pergamon Press Ltd. Printed in the United States

AN ASSESSMENT OF THE ACCURACY OF NUMERICAL SOLUTIONS TO THE GRAETZ PROBLEM Nancy Conley, Adeniyi Lawal and Arun S. Mujumdar Department of Chemical Engineering McGill University Montreal, Quebec Canada (Ctamalnicated by J.P. Hartnett and W.J. Minkowycz) ABSTRACT I t has been w e l l documented i n p r i o r s t u d i e s t h a t n u m e r i c a l s o l u t i o n s t o t h e c l a s s i c a l Graetz problem i n t h e e n t r a n c e r e g i o n o f p i p e s produce e s t i m a t e s o f N u s s e l t number t h a t are greatly in error. I n t h i s r e p o r t , we examine t h i s problem and a t t e m p t t o p r o v i d e p l a u s i b l e e x p l a n a t i o n f o r t h e d i s c r e pancy.

Introduction In spite of the gross idealizations involved in its formulation,

the

classical Graetz problem surprisingly, provides valuable guidance in the design and development of compact heat exchangers.

Of greater significance

is the fact that this problem has served for decades as a model for testing various solution techniques since its solution has been adequately verified and widely disseminated [I]. Since the historical attempt by Graetz [2] in 1883 considerable effort has been expended on the solution of this problem or one of its many variations; the discussion here will he limited to only the most relevant ones. Hmploying the separation of variables technique, a closed-form solution is obtainable in the form of an infinite series, the eigenvalues and eigenfunctions of which satisfy the Sturm-Louiville system. were evaluated by Graetz.

The first two terms

Sellars et al. [3] subsequently extended this

work by providing the first ten eigenvalues and constants.

The major draw-

back of this solution approach is the poor convergence behaviour of the series near the entrance to the duet where even with 121 terms the series does not converge

[13.

The difficulty was partially resolved by Leveque ~4~ 209

210

N. Conley, A. Lawal and A.S. Mujumdar

Vol. 12, No. 2

who employed the 'flat plate' solution near the singularity at the entrance. Transport problems commonly encountered in industrial practice are almost always more involved than the Graetz problem.

Furthermore, the re-

sulting mathematical models do not generally admit of closed-form solutions. Usually recourse must be made to numerical solution techniques.

Inherently

numerical methods suffer from a major weakness -- the absence in most situations of a means of internally verifying the accuracy of the solutions which they generate.

This becomes even more pertinent since for every new

combination of parameters, their repetitive nature demands that the computer codes be validated.

The Graetz problem and other such standard prob-

lems that possess exact or closed-form solutions have served this useful purpose. The extended Graetz problem which included axial heat conduction was solved numerically by Schmidt and Zeldin [53.

For very high Peclet numbers

this problem reduces to the classical Graetz problem.

Thus it is expected

that the numerical solutions at very high Peclet numbers should closely

approach the closed-form solutions.

However, as noted in

[13, close

to the

duct entrance a deviation in local Nusselt number of up to 25% was observed making the numerical results in this region (for any Pe number) highly suspect.

Another numerical effort in the solution of the Graetz problem was

reported by Grigull and Tratz [6~; here again the values of the mean Nusselt number were about 16% lower than those evaluated from the closed-form solutions.

The implications of this general observation for other more comp-

lex problems that do not possess verifiable analytical solutions are greatly disturbing. Intuitively, the discrepancy between numerical and exact solutions can be attributed to the singularity at the entrance and the associated method of treatment by the numerical

algorithm[13. There

in the literature in favour of this argument.

is no compelling evidence

Another possible explanation

is rooted in the observation made by Collins [73 that I percent general errors (even greater near entrance) are common in the energy balance when the wall thermal condition is explicitly imposed numerically.

He therefore

suggested replacing the discretization equation at the wall node but one by an integral energy balance reminiscent of the overall continuity equation for mass conservation.

This hypothesis is tested in this report.

Because

of the large axial gradients in the vicinity of the entrance and also the occurrence of rapid radial changes in the thermal profile one might want to appropriately concentrate the grids in these regions.

The obvious

VOl. 12, No. 2

NIIMERICAL S(X;/rlCNS TO THE GRAETZ ~

211

advantage of t h i s procedure i s t h a t i t affords the placement of a v a i l a b l e computing power i n a r e a s where i t i s d e s i r e d most. f u n c t i o n s i n t h e r a d i a l d i r e c t i o n were t e s t e d ,

Two c o o r d i n a t e s p a c i n g

one o f which had been p r e -

v i o u s l y employed by Schmidt and Z e l d i n f o r a s i m i l a r problem. The a u t h o r s a r e unaware o f any e f f o r t s t o n u m e r i c a l l y r e p l i c a t e t h e e x a c t r e s u l t s o f t h e Graetz problem ( e s p e c i a l l y i n terms o f t h e l o c a l Nuss e l t number) i n t h e e n t r a n c e r e g i o n . on t h e v a l i d i t y

Such an e x e r c i s e w i l l shed some l i g h t

of various numerical procedures.

Herein l i e s t h e prime m o t i -

vation for this study. Analysis Of i n t e r e s t

i s the entrance region laminar heat t r a n s f e r to viscous i n -

c o m p r e s s i b l e c o n s t a n t p r o p e r t y Newtonian f l u i d f l o w i n g i n a c i r c u l a r t u b e w i t h a f u l l y developed e n t r a n c e v e l o c i t y p r o f i l e .

The f l u i d e n t e r s t h e

duct w i t h a u n i f o r m t e m p e r a t u r e T and i t i s t h e n h e a t e d or cooled t h r o u g h e t h e duct w a l l whose t e m p e r a t u r e i s m a i n t a i n e d u n i f o r m a t Tw ( d i f f e r e n t from the entrance value).

The n e g l e c t o f such e f f e c t s as v i s c o u s d i s s i p a t i o n ,

i n t e r n a l h e a t g e n e r a t i o n and a x i a l h e a t c o n d u c t i o n completes t h e d e f i n i t i o n of the c l a s s i c a l

Graetz problem,

For t h e model i d e a l i z e d above t h e t h e r m a l energy e q u a t i o n r e d u c e s t o : 1 3

( r 3T

= u 3T

(la)

w h i l e t h e t h e r m a l boundary c o n d i t i o n s a r e g i v e n by: T ( r , 0 ) = Te; T(ro,X ) = Tw and -3T ~ (0,x) = 0

(lb)

With t h e i n t r o d u c t i o n o f t h e f o l l o w i n g d i m e n s i o n l e s s v a r i a b l e s ; R = r/r ° 0 = T-Tw_ T -T e w

X = (x/2ro)/P e t h e e n e r g y e q u a t i o n a l o n g w i t h t h e boundary c o n d i t i o n s become: 2 ~ (R 30 30 3--R " -~) = (l-R2) * ~X 30 O(R,O) = 1 . 0 ; O(1,X ) = 0 and ~-~

(2a)

(0,X*) = 0

(2b)

where t h e v e l o c i t y p r o f i l e has assumed i t s u s u a l p a r a b o l i c form. finitions

The de-

o f o t h e r v a r i a b l e s have been r e l e g a t e d t o t h e n o m e n c l a t u r e .

212

N. Conley, A. Lawal and A.S. Mujumdar

Vol. 12, No. 2

Method of Solution: To exercise control on the grid near the wall we introduce a new variable such that variable spacing in the R direction corresponds with uniform spacing in ~. R

=

If ~ and R are related as below:

cA

(3)

with A suitably chosen we can generate a mesh that is fine near the wall and progressively coarser near the center.

We refer to this as COORDI.

If

the necessary substitution is made, Eq. (2a) can be written as: __ ~2-2A 3 2 0 + A2 3~2

gl-2A (~A-~---) 3 0 _ i (i_ 2A) 30 3~ 2 3X*

(4)

In order to confine grids near the wall, Schmidt and Zeldin employed a logarithmic relationship (herein called COORD2) which can be generalized as: R

=

1

In (I + (ea-l)~)

(5)

The parameter a serves the same purpose as previously stated, and the energy equation rewritten in the new coordinate becomes:

a2 (1+ E ( e a - 1 ) ) 2

__320

(ea-1) 2

a 2 (1 + E ( e a - 1 ) )

3~2 +

(1 +

(ea-1)

1

30

In (1 + ( e a - 1 ) ~ ) ) - ~ ' =

2 (1 - ( l n ( l + ( e a - 1 ) ~ ) ) 2) 830 -~

(6)

A s t a n d a r d g r i d network w i t h t h e dependent v a r i a b l e

s t o r e d a t each g r i d

p o i n t i s l a i d out on t h e c o m p u t a t i o n a l domain h a v i n g a u n i f o r m r a d i a l o f Ag.

Second o r d e r c e n t r a l

all partial

derivatives

the first-order

difference

in the radial

derivative

spacing

scheme i s employed i n a p p r o x i m a t i n g

direction

and t h e upwind scheme f o r

in the axial direction.

The r e s u l t i n g

discretiza-

t i o n e q u a t i o n s a r e s o l v e d u s i n g t h e TDMA ( T r i d i a g o n a l M a t r i x A l g o r i t h m ) . As remarked e a r l i e r

i n s t e a d o f forming a d i s c r e t i z a t i o n

w a l l node b u t one u s i n g Eq. ( l a ) , which r e l a t e s the vicinity

we employ t h e i n t e g r a l

equation at the

energy balance

t h e a x i a l change i n b u l k t e m p e r a t u r e t o t h e t e m p e r a t u r e s i n of the wall.

point using regular radial

The f i n a l

form o f t h e e q u a t i o n f o r t h i s node

s p a c i n g in R i s g i v e n b y :

Ob(X* ) _ 0b(X*_AX* )

=

4~X ARR (30(M) - 4 0 ( M - l )

+ 0(M-2))

(7) In t h i s c a s e , however, t h e d i s c r e t i z a t i o n

e q u a t i o n s a r e n o t amenable

t o s o l u t i o n s by TDMA. T h e r e f o r e , t h e s u c c e s s i v e o v e r - r e l a x a t i o n

(SOR)

Vol. 12, NO. 2

NL~ERICAL~GNS

p r o c e d u r e was u t i l i z e d . two q u a n t i t i e s

TO THEGRAETZ P ~

Upon o b t a i n i n g t h e t e m p e r a t u r e p r o f i l e we c a l c u l a t e

of s p e c i a l i n t e r e s t ,

t h e l o c a l N u s s e l t number, t h e l a t t e r

t h e d i m e n s i o n l e s s b u l k t e m p e r a t u r e and o f which i s d e f i n e d a s :

~_.__RRIR=I -2~e

Nu . = X

213

(8)

eb

The bulk temperature was computed using Simpson's rule while the thermal gradient at the wall was numerically evaluated from the 3-point rule except otherwise implied.

The other choices for the approximation of the latter

quantity include the 5-point rule and the log relationship presented in detail i n [5]. Results and Discussion In Fig. 1 is depicted the effect of the number of axial steps on Nu

at X*=I0 -6 for varying degrees of radial grid refinement. Clearly with X* 21 regular radial grid points the Nu which is substantially in error, X*'

15£

M= 101

lOG

. . . . . . . . . . . . . . . . EXACT . . . . . . . . . . . SOLN. ...............................

M= 301

Nu x, M=21

5(:;

0

I

I

5

No oFIOAxIAL

I

15 STEPS N

FIG. 1 Effect of axial step size for regular radial grids *

with X =10

-6

.

20

214

N. Conley, A. Lawal and A.S. Mujumdar

shows no noticeable dependence on the axial grid distribution.

Vol. 12, No. With M=301

the numerical solution is almost indistinguishable from the exact solution except for N
Keeping the number of axial steps fixed and varying M

produced the effects illustrated in Fig. 2 from which we infer that no matter what value M attains an optimal number of axial steps is required for

. . . . . . . ~NA~Y-gb-Cg-- -

1OO

~

~

=

~

5O

o

~

I~o

,~o No

OF

~o

2~o

36o

RADIAL POINTS M

FIG. 2 Effect of radial grid refinement for three values of *

N with X =I0

reasonable accuracy.

-6

.

The profiles presented on these figures provide con-

vincing evidence that the conventional test of numerical convergence of finite difference solutions may not be conclusive especially with regard to derived properties such as NUx,.

For illustration, in Fig. 2 there is

no a p p r e c i a b l e are greatly

c h a n g e i n Nu . when M l i e s b e t w e e n 75 a n d 125 y e t t h e v a l u e s X in error. The s a m e a p p l i e s t o F i g . 1 w h e r e t h i s b e h a v i o u r i s

strikingly clear for M=21. The radial grid points can be conveniently concentrated near the wall with the parameter A
Figure 3 examines its effect on Nu .. X Evidently in contrast to the regular grid an appropriate choice of A can assuredly replicate the exact solution even for a coarse grid of M=21 but the @b (not shown) begins to deviate quite significantly from the exact solution for A<0.6 and M
The behaviour of Nu , for different values X of parameter a is not dissimilar from that for A as can be observed from Fig. 4.

In fact when a~0.0 and A=I.0 (regular grid) the NUx, values are

Vol. 12, No. 2

~ C A L

SOLUTIONS X D T H E I

t

I

•5

of parameter * -6 M with X =i0 . 1

.

.

.

.

.

.

.

M with X =I0

I

.7

!

!

/

I

~B

A

.9

I

1J

values

-I

of

"

. . . .

FIG. 4 a on Nu . f o r X

.

........

.

of parameter -6

I

.6

215

i

FIG. 3 A o n Nu . f o r t h r e e X

Effect

*

1

-~--~:

10(

Effect

I

GRAETZ P N 0 ~

three

values

of

216

N. Conley, A. Lawal and A.S. Mujumdar

Vol. 12, No. 2

* -6 Ten uniform axial steps were taken to reach X =I0 .

the same.

Since the intent of this work is to examine Nu ~nd @b near the entrance, further results are presented in Table 1 for the t~ree types of grid employed in the analysis.

The initial step size AX

was 10 -7 , and for each order of

TABLE 1 Ob and NUX. distributions using different radial grid distributions

Regular Grid M

21

X

0b

,000001 .000004 .000008 .000020 .000060 .000100 X

.000001 .000004 101 .000008 .000020 .000060 .000100 X

.99979 .99919 .99842 .99634 .99092 .98685

Nu . X 59.345 57.452 55.075 49.002 35.443 27.951

@b

NUx*

.99937 130.561 .99842 71.557 .99748 54.674 .99542 40.142 .99048 27.066 .98665 22,558 @b

.000001 .000004 .000008 301 .000020 .000060 .000100

99937 99841 99747 99540 99047 98664

magnitude change in

Nu . X 109.811 68.375 53,492 39.757 26.985 22.519

X*

COORD1 (A=. 6) 0b

(a=l. O)

.98656 .98565 .98462 .98239 .97745 .97365

Nu . X 96.607 84.536 72.004 51.218 29.879 24.014

.99968 .99880 .99779 .99552 .99055 .98677

Ob

NUx*

Ob

.99747 113.856 .99652 69.103 .99558 53.813 .99351 39.903 .98857 27.050 ,98475 22.568 @b

Shah & London

COORD2

NUx*

.99886 109.252 .99790 68,313 .99696 53.483 .99496 40.122 .98996 26.996 .98613 22.529

Ob

.99957 .99841 .99748 .99541 .99047 .98664 @b .99937 .99841 .99747 ,99540 .99047 .98664

Nu . X 91.195 81.418 70.849 51.833 30.075 23.940 Nu . X 114.463 69.098 53.756 39.841 27.007 22.525 Nu . X 109.225 68.283 53.457 39.746 26.982 22.517

Ob .99936 .99839 .99745 .99534 .99039 .98657 Ob .99936 .99839 .99745 .99534 .99039 .98657 @b .99936 .99839 .99745 .99534 .99039 .98657

a corresponding change in AX* was made.

Nu . X 106.538 66.731 52.763 38.637 26.544 22.275 Nu . X 106.538 66.731 52.763 38.637 26.544 22.275 Nu . X 106.538 66.731 52.763 38.637 26.$44 22.275

With M=21 and

X*=I0 -4 the Nu . values show a maximum deviation of about 26% from Shah and X London's results which were obtained from analytical solution of Leveque's problem.

However, when M=I01 the numerical results regardless of the radial *

grid employed can be accepted as accurate to within 3.6% for X >8XI0

-5

.

The results of Table 2 indicate that for coarse radial grids Nu . is X dependent on the expression employed in approximating the wall thermal gradient, the arbitrariness so introduced further attests to the unsuitability of such *

grids even though this becomes less pronounced for X =i0

-4

.

The integral

energy balance was incorporated into the code and the Nu . evaluated at X

Vol. 12, No. 2

Nu . a t X

~

N~CAL

selected

values

SOLUTIGNS TO THE GRAETZ P R ~

, TABLE 2 of X using three thermal

gradient

different

methods for wail

evaluation

M

3-Point

X*=10 -6 5-Point

LOG

21

59.345

82.017

38.674

27.951

25.304

22.641

ioi 130.561 ~ 301 109.811

116.340

109.474

22.558

22.514

22.517

108.761

109.026

22.519

22.517

22.517

,-,

21

96.607

131.316

63.548

24.014

20.904

22.802

~

101

113.856

103.208

109.203

22.568

22.560

22.560

"-~ 301

109.308

109.071

109.082

22.529

22.529

22.529

21

91.195

124.395

59.870

23.940

19.963

22.643

101

114.463

101.572

109.130

22.525

22.517

22.547

301

109.225

109.012

109.074

22.517

22.517

22.527

0°5

~ O II or.~ , .m w

X*=10 -6 for the coarse grid of M=21.

217

X*=10 -4 3-Point 5-Point

LOG

There is no noticeable improvement

(Nu ,=59.244) when compared with the corresponding value in Table I, an outX come similar to that obtained when the control volume method was utilized in discretizing the energy equation E83.

It is appropriate to comment that all

the results require about the same number of CPU seconds so there are no cost differences as far as computing is concerned. Based on the above, it may be concluded that very close to the entrance an optimal grid size, both in the radial and axial directions, is required and a simple redistribution of grid points may not eliminate the inaccuracy of the results.

As no exact solutions exist in general for the more compli-

cated entrance problems such as the d e v e l o p i n g

flow and heat transfer

in ducts, the findings cannot strictly be generalized, except in so far as they provide usable guidelines.

Prior results obtained with few grid-

lines may not be sufficiently accurate very near the entrance. Acknowledgements The authors are indebted to Professor W. Thorpe for providing extensive computer time on the Amdahl V7A of the McGill Computing Centre. Nomenclature a

- Coordinate spacing parameter

A

- Coordinate spacing parameter

218

N. Conley, A. Lawal and A.S. Y~/jumdar

M

Number of radial grid points

N

- Number of axial steps

Nu .

Local Nusselt number

Pe X

Peclet number = 2roPCpUm/k

r r

Vol. 12, No. 2

- Radial coordinate o

R

- Tube radius - Radial coordinate (dimensionless)

T

- Temperature

Te

- Entrance temperature

Tw

- Wall temperature

um

- Mean velocity

x

- Axial coordinate

X*

- Axial coordinate (dimensionless)

Greek Symbols: - Thermal diffusivity (pCp/k) - Radial coordinate @

- Temperature (dimensionless)

@b

- Bulk temperature (dimensionless) References

i.

R.K. Shah and A.L. London, Laminar Forced Convection in Ducts, Advances in Heat Transfer, New York; Academic Press.

2.

L. Graetz, On the Thermal Conductivity of Liquids, Part I; Ann. Phys. Chem. 18, 79 (1883).

3.

J.R. Sellers, M. Tribus and J.S. Klein, Heat Transfer to Laminar Flow in a Round Tube or Flat Conduit -- the Graetz Problem Extended, Trans. ASME 78, 441 (1956).

4.

M.A. Leveque, Les lois de la transmission de chaleur par convection, Ann. Mines, Mem., Ser. 1 2 1 3 , 201 (1928).

5.

F.W. Schmidt and B. Zeldin, Laminar Heat Transfer in the Entrance Region of Ducts, Applied Scientific Research 2_33, 73 (1970).

6.

U. Grigull and H. Tratz, Thermischer einlauf in ausgebildeter laminarer Rohrstromung, Int. J. Heat Mass Transfer 8, 669 (1965).

7.

M.W. Collins, Finite Difference Analysis for Developing Laminar Flow in Circular Tubes Applied to Forced and Combined Convection, Int. Journal for Numerical Methods in Engineering, 15, 381 (1980).

8.

N.A. Poirier, Numerical Solution of the Graetz Problem in the Entrance Region, Internal Report Chemical Eng. Dept. McGill University (1984).