WADC TECHNICAL REPORT 54-255 April, 1954 HEAT TRANSFER TO LAMINAR FLOW IN A ROUND TUBE OR FLAT CONDUIT THE GRAETZ PROBLEM EXTENDED JOHN SELLARS MYRON TRIBUS JOHN KLEIN WRIGHT AIR DEVELOPMENT CENTER, U. S. AIR FORCE CONTRACT NO. AF 18(600)-51, E. 0. NO. 462 Br-1

SUMMARY The complete set of eigenvalues and eigenfunctions for the classical Graetz problem is presented and the solution is extended to cover arbitrary wall temperature or heat-flux variations. WADC TECHNICAL REPORT 54-255 iii

SYMBOLS AND NOMENCLATURE A = Coefficient occurring in Equation (10). B = Coefficient occurring in Equation (10). b = Half width of flat duct, ft. Cn = Coefficient in Equation (3). Cp = Unit heat capacity at constant pressure, Btu/lb-~F. D = Coefficient occurring in Equation (18). E = Coefficient occurring in Equation (18). F = F(s, r+), laplace transform of Graetz solution. g = g(x+, r+) Integrating kernel for heat-flux problems, see Equation (44). G = G(s, r+) Laplace transform of g. h = h(x+) Integrating kernel for heat-flux problems, see Equation (34). H = H(s) Laplace transform of h. i = 4-1. J = Bessel function of first kind, zero order. J1/3 = Bessel function of first kind, 1/3 order. J-1/3 = Bessel function of first kind, -1/3 order. k = Thermal conductivity of fluid, Btu/sec-ft2 (~F/ft). Kn = Coefficient occurring in Equation (A-3). Pr = Prandtl modulus, dimensionless, ((ICp/k)(3600 gc). q = q(x) heat flux per unit wall area, Btu/hr-ft2. WADC TECHNICAL REPORT 54-255 iv

Q = Laplace transform of (kq/ro). r = Radius, ft. o = Tube radius, ft. r+ = (r/ro). R = R(r+). Re = Reynolds modulus, dimensionless, (2Um rOp/4igc) or (4Um bp/pgc ). S = Transform variable. t = t(x+, r+) temperature, ~F. T = T(s, r+), Laplace transform of t. u = Velocity of fluid, ft/sec. Um = Average fluid velocity in tube, ft/sec. x = Distance along tube, ft. x = (x/ro)(RePr)-1 or (x/b)(RePr)1. y = Distance from duct wall, ft. =(y/b). z = Distance from tube wall, ft. += (Z/rO). 7 = Zero of H(s). X = Eigenvalue. = Viscosity of fluid, lb-sec/ft2. p = Fluid density, lb/ft3. = Dumny variable. = Dummy Variable. = Gamma function. WADC TECHNICAL REPORT 54-255 v

HEAT TRANSFER TO LAMINAR FLOW IN A ROUND TUBE OR FLAT CONDUIT THE GRAETZ PROBLEM EXTENDED INTRODUCTION The problem considered here is posed by a system in which a fluid of constant properties flows in steady laminar motion in a round tube or flat duct. The velocity profile is fully established and parabolic. Up to a point (x = 0) the fluid is isothermal. After this point a prescribed heat flux or temperature is given at the wall of the conduit and the problem is to find the temperature distribution, as well as the connection between heat flux and wall temperature. The application of this solution to practical problems of heat exchange has already been so well established that further comment is unnecessary. The problem has been considered in detail by a number of workers and an excellent review is contained in the book "Heat Transfer" by M. Jakob.l The problem readily reduces to the finding of eigenvalues and prior to this paper only the first three eigenfunctions and the first four eigenvalues have been known. A recent paper2 has brought out the importance of obtaining more eigenvalues, and by using the complete set of eigenvalues and the methods of reference 2 the classical "Graetz Problem"3 is extended to more complicated boundary conditions. The problem can be stated in mathematical terms as follows: Given tt t(%,A^) za 4. \ A (1) X = ECIARE(AR 5 and for t=to X<o WADC TECHNICAL REPORT 54-255 1

with either t(z, S,) A FI)) X >o or jt-(S, J^ ) = %( find t(x,r) and the relation between q(x) and tw(x). The nondimensional form of the equations is (2) The boundary conditions are "t tt. and either t (t, I) tw() or t (%+ =X In view of the linearity of Equation (2), it is necessary to have only the fundamental solution, known as the Graetz solution, to construct all other needed solutions. Therefore, the initial step is the completion of the Graetz solution. THE GRAETZ SOLUTION The problem considered by Graetz and most other workers is Equation (2) with boundary conditions. t -( < x<o tt 0-o % >0 Led ~ be a solution of Equation (2), then 0 c0 R ) (3) where the Xn are the eigenvalues required to make the solution to the following differential equation WADC TECHNICAL REPORT 54-255 2

V R~ + R o, ~~. (-^;- ) n,o (4) satisfying the boundary conditions Rn(1) = 0, Rn(O) = 1. The coefficients Cn are determined from the relation CI' (,_77,2 d* X (5) The eigenfunctions and eigenvalues have been given only for n = 1, 2, 3. The higher modes of Equation (4) are very difficult to calculate for large values of X. Therefore, to obtain kn and Cn for n > 3, a solution is sought which will be valid as n -+ o. It will be found that the resulting formulae will provide good answers even when Xn is small. First, look for a solution in the form Rse and find that g(r+) satisfies I + (4.t - o (6) Now an asymptotic solution is sought in the form a=' = a * +h>:, (7) Substitution in Equation (6) and equating powers of X gives Io = ~ i T-t (8) 3. = -; w (9) Since X is large, the remaining terms in Equation (7) are neglected. Substitution of Equations (8) and (9) in Equation (7) gives for R ACe 3 e R = — e ~ "',j, —_ (10) Equation (10) is the so-called WKB approximation and is valid for o<r+ < 1 for sufficiently large X. Now the coefficients A and B must be determined so that WADC TECHNICAL REPORT 54-255 3

Equation (10) will correspond to the regular solution of Equation (4), where r+ is small. For small r+ Equation (10) is As +se_ e (11) Tw ( I - nt ) Inspection of Equation (4) shows that when r+ is small enough so that \2 (1 - r+2) +> 2, the classical solution behaves as Jo (Xr+), since Equation (4) then becomes a Bessel equation. For large xr+, even if r+ is small, the asymptotic expression for Jo(Xr+) is TJ(xOt) J- J^ (2P ( X_ 5-^) (12) and thus, it is seen that to make Equations (11) and (12) equal for r+ small, it is required that A =J e = / et (13) and for 0 < r+ < 1 (.^t) = vTTE _7V4 (14) Equation (14) is not a good approximation to the solution as r+ + 1, since it has a singularity there. Because a boundary condition is to be imposed at r+ = 1, the development of an alternate solution, valid near r+ = 1, is considered. By patching it on to Equation (14) the solution over the range O' r+ 1 is obtained. The following change of variable is made 3 -,-+ and Equation (4) becomes A - -- 0+t( (15) Now consider 0 < z+ < 1 and define a new variable WAD = TEHIA R Rt a(16) Substitution of Equation (16) into Equation (14) yields for large x WADC TECHNICAL REPORT 54-255 4

IL (17) which has the solution o 3,(h^ a )t (~Y tA(1) (18) The constants D and E are to be so chosen that for small z+ Equations (18) and (14) are equivalent. Change the variable from r+ to z+ in Equation (14) and perform the integration +t f Si- AI = i j +- -a + t- (19) For small z+ Equation (19) yields + 3/ So At 13 =^-c a/ ^3 (20) so that Equation (14) for small z+ is Co <><a (j V\xa+'_i (X-S') X (21) For large \z+, even if z+ is small, Equation (18) becomes R(3) = D (7I'L S=; _ (4) E C(st 2) Expanding the cosines of differences of angle occurring in Equations (21) and (22) yields the simultaneous equation iD an are lt ~r Eut (18) 4 WADCeN T~ ECHNICt = ~AL REP~ORT~ 5-(25) from which D and E are evaluated. Therefore, Equation (18) is WADC TECHNICAL REPORT 54-255 5

3 -1/3 volving J becomes constant. Therefore, the coefficient of J must be zero if R = O at z+ = O. The values of Xn must therefore be given by A ~ Vla, ok (25) The equations for Rn are therefore for small r+ (center of pipe) R., (A') " o (Am.At) (26) for medium r+ R, ^r7).. V T I - R+ 0-~).tL (27) (,_,t.),4 and for small z+ = 1 - r+ (near the wall) R3(3<) -. C (-1) ^ T () (28) Equations (24) to (28) contain all the information essential to the problem solution. The coefficients Cn in Equation (3) are found from Equation (24) in accordance with Equation (5). Thus it is found that (/>)A.\_ = no (0 ok 1, a(29) and therefore Go e, (- Fib) >3 ~' ad, ",, (30) The derivative of R at the wall (z+= O) which is ^ " 3 -@ = r 3.",, - ~... (31) ~'(~I): - =.... ~:*~(o r (%) 3~/L will be required later. WADC TECHNICAL REPORT 54-255 6

Table I shows the first ten eigenvalues and the important constants for the case of flow in a round tube. Table II gives the same data for a flat duct with opposite walls at the same temperature. The development of the flatduct system is similar to the round duct and the equations are given in the appendix, numbered to correspond with the text. The previously known eigenvalues given by Jakob are shown in Table III for comparison. Since the solution presented here is valid for large \X, and in view of the agreement even at moderate values of Xn, it has been concluded that all the eigenvalues and functions are now sufficiently accurately known. The heat flux at the wall is computed from the equation ( =i(-t - -)-t _ _ X-/t) (32) Equation (32) is presented in the above form to bring it into agreement with Jakob.1 ARBITRARY WALL-TEMPERATURE VARIATIONS If the wall-temperature variation is given by tw(x), then, as shown by Tribus and Klein,2 the principle of superposition may be applied and the solution may be written in a Fourier-type Stieltjes integral t-o =To E' - [ e( - ).] dUW (5) (33) t=0 where 9 is the solution to Equation (2) defined by Equation (3). The temperature of the wall and fluid for x+ < 0 is to. The Stieljes integral in Equation (33) is evaluated by substituting (dtw/dS) dS for dtw wherever tw is continuous and substituting [1 - 0 (x+ - %1, r+)][t (j1+) - t (Ei-)] as the contribution of the integral wherever tw (x+) has a discontinuity. (See Tribus and Klein2 for a more detailed discussion.) The heat flux is computed from (34) W (, C — EC i CA (R E-P ) 5d-w 25( WADC TECHNICAL REPORT 54-255 7

TABLE I FIRST TEN EIGENVALUES AND THE IMPORTANT CONSTANTS FOR THE CASE OF FLOW IN A ROUND TUBE n Cn /n2 Cn -2 Cn Rn (1) 0 2 2/3 7.1129 +1.47989 0.7303 1 6 2/3 44.489 -0.80345 0.553810 2 10 2/3 113.785 +0.58732 0.460074 3 14 2/3 215.121 -0.474993 0.413743 4 18 2/3 348.457 +0.404448 0.381785 5 22 2/3 513.793 -0.355345 0.357853 6 26 2/3 711.129 +0.318858 0.338988 7 30 2/3 940.465 -0.290488 0.323555 8 34 2/3 1201.8 +0.267691 0.310596 9 38 2/3 1495.1 -0.248895 0.29950. (-,) -srv) -l - C - CC-. v (2/3) 3 x p - R- r(3) 3 = 3,7\ 3 am s i +^ 2s', -., *, W, 6 _ 2 C, g(nu) e.2: ^ W- TZE H N e- (4f5) WADC TECHNICAL REPORT 54-255 8

TABLE II FIRST TEN EIGENVALUES AND THE IMPORTANT CONSTANTS FOR THE CASE OF FLOW IN A FLAT DUCT WITH OPPOSITE WALLS AT THE SAME TEMPERATURE n n2 Kn -Kn Yn' (1) 0 1.667 2.779 +0.50.685 1 5.667 32.11 -0.121.454 2 9.667 95.45 +0.0648.380 3 13.67 186.9 -0.0431.338 4 17.67 312.2 +0.0319.311 5 21.67 469.6 -0.0253.291 6 25.67 658.9 +0.0207.274 7 29.67 880.3 -0.0174.262 8 33.67 1134 +0.0150.251 9 37.67 1419 -0.0131.242 3 383 r() 3yX - (- o 3 X - -7Trr ~L \^ _ + + t/3'~! ~ 0 ^ ~ I, * * 6 = = /< Y\ r )) e WADC TECHNICAL REPORT 54-255 9 WADC TECHNICAL REPORT 54-255 9

H TABLE III o COMPARISON WITH PREVIOUSLY KNOWN EIGENVALUES Results Obtained AV1 Sellars, Tribus, Klein JakobAnalogue Computer n C Rn (1) cn. -n Rn' (1) -nCn Rn (1) n An Cn 2...... 2n Cn 2 2 2 0 2.667 +1.47989 0.7303 2.705 +1.477 0.749 2.71 1.46 0.735 0 1 6..667 -0.80345 0.5381 6.66 -0.810 0.539 6.69 -0.809 0.533 2 10.667 +0.58732 0.4601 10.3 +0.385 0.179 10.62 +0.592 0.444 3 14.667 -0.47499 0.4137 14.67* -0.479* 14.58 -0.51 0.598 *Attributed to Lee, Nelson, Cherry and Boelter.

HEAT FLUX AT THE WALL GIVEN The inverse problem; namely, "Given the heat flux at the wall, what is the temperature?", may be solved with the aid of the Laplace transform theory Define the following transforms T(s,+t) -5e &Sz-(tt t) 6SPA") e- t - to) je (35) T (S) = T(, 1) (36) 00 f(-_ht) = i F,- 8t^ ^)] e -5% (37) -s. (38) H(s) T^+ (s, I) = - I B^ ( l C d (58) 00 Q(S) =f & e0) (39) Applying the Faltung theorem to Equations (33) and (34) yields; T(s, A+) = F(, A+) S Tw (s) (40) (s) = H(s)s r(s) (41) If the heat flux is finite, tw(x+) will be continuous. Eliminating tw(s) from the above equations, T(s, A/) aF(s, AL ) q (42) WAIDC TECHNICAL REPORT 54-255 11 WADC TECHNICAL REPORT 54-255 11

Now define G(s,+) = (S)) (43) and let g(x+, r+) be the inverse transform of G(s, r+). Then, for arbitrary heat flux at the wall, the temperature is given by = So 3( it-5, )5()1t (44) Thus, the problem is reduced to finding g(x+, r+), which is given by C+ joo,(.) = a e (YS) AS (45) o' an-t J.H(s) C- 6 (S Returning to Equations (37), (38), and (3) it is found that F S~, C r".,, (A+) rS.. S+X, (46) H(S ) -(., ( %o 5+_ a S+ (47) Because F and H have poles at s = -X2n, the quotient F/H has no poles except at S = 0 and the zeroes of H(s) and the zeroes of H(s) must be found numerically. Because H'(s) is monotonic, it is found that the zeroes of H(s) occur between the -X2n. Letting y2m be the values satisfying H(-Y2m) = O, from the theory of residues _., AD _ X ) _: e -_ z'c^, Re^% E (48) ii /0) Table IV gives the values of y2ml H'(-72m) for the first three values of m. The term H(O) has been shown by others to be given by H(o) - + 1 (49) hence, the wall temperature may be easily calculated with the aid of -e X 5v7ADI,- 42 H'EPO) (50) WADC TECHNICAL REPORT 54-255 12

TABLE IV THE VALUES OF 2m, H'(-r2m) FOR THE FIRST THREE VALUES OF m Roots of H(s) = O, Values of H'(-72m) - N z canw I (1) Hs) = - ct I (I) _. _ 1, i 1 25.639 8.854 x 1o-3 4.405 2 84.624 2.062 x 1-3 5.7508 3 176.4 9.45 x O1-4 6.oo84 WADC TECHNICAL PEPORT 54-255 15 WA~DC TECHNICAL REPORT 54-255 13

TABLE V VALUES OF V2, H'(-72) FOR THE FIRST THREE VALUES OF m FOR A FLAT DUCT Roots of H(s) = O, Values of H'(-72m) ('s) = -Z kOY'(I) K, Y,/(') (S + x3' Is) = ~ Y'(,). 4, X, /, 2, m 72 -H' (-72m) 2 (In m} 72m H(1 49.345 7.45 x 10-4 27.2 2 185.94 1.67 x 10-4 32.1 3 409.45 6.89 x 10-5 35.4 WADC TECHNICAL REPORT 54-255 14

A SAMPLE CALCULATION FOR CONSTANT WALL HEAT FLUX By way of illustration consider the computation of the asymptotic value of the Nusselt modulus for the case of constant heat flux at the wall. Combining Equations (44) and (48) with q(S) = q = constant, the following is obtained. = e;nv~ (%^S4)+) x' (%+J +) ht^ ^ -^ -^ (c X7 ~t =' 8) S(_ - X CR ~ X VW (51) Letting Px = x+, where P = ixk/2WCp, and integrating Equation (51) gives tl+ A/t) -to - t A~,, v; { _ i e W^ - H'(- Y < (52) -z~.,. which may be rewritten as t(+ at+) -to = J X - E, t-' E 5i- e_ x (53) E+q ( s H' fr dc t priv Equation (55) shows that far down the pipe (x+ co) the derivative of t with respect to x is independent of x or r+; i.e., 2it _h _ao r _t r (54) Substituting this quantity into Equation (2) leads to WADC TECHNICAL REPORT 54-255 15

' / JL- I ) (55) R =(A(I-A;) ^+ l ^+/ which may be integrated directly to give t(tA+) -t( = / ( ) (56) Now the mixed mean temperature along the pipe is given by -e>5, -t e (57) but the mixed mean temperature is also defined by I t (^+) = u c p t(? X A) 2T < J A! _ L~jD~p_ ii~t~)2~,jt~j~t(58) W C. Substituting Equation (56) into Equation (58) and integrating results in, ) - -( o) -= 4 — (59) Combining Equations (59) and (56) - (A+ = > __ _ A. ":,, - Z 4i7 (60) and substituting Equation (57) into Equation (60) A;)(. AD-___ EA + -*' 7 (61) at r+ = 1. This expression reduces to t (X ) - o = + (62) but from Equation (52), since Rn(1) = 0, t If ) t-o _ L-t, =-' I( (63) Hence, WADC TECHNICAL REPORT 54-255 16

^- ~1' —----- = i0. AfS Y H'( - _ + (64) (Note that the first three terms sum to approximately -0.27.) Now, the Nusselt modulus is given by Nu = = (65) X ( tW "tMI) From Equation (60) t ( -t I),,, ='If * (66) which when substituted into Equation (65) gives Nu = % 1.36 (67) Substitution of (64), (58) and (53) into (65) yields the local value of the Nusselt Modulus for the case q(x+) = constant. v___1~~~~ / ~~(68) Nu =: -J ~ N~. CALCULATION FOR LINEARLY VARYING WALL TEMPERATURES In similar fashion the use of the boundary condition Tw(x*) - To A x* where A = any constant, gives: a^x+) = AA s JA RI (') a I (69) -,, ('- = n — _ _ (71) Jr' n,, -?, - + aj ~ —- 7 — s e WADC TECHNICAL REPORT 54-255 17

APPROXIMATIONS FOR SMALL x+ Whenever x+ is small, a large number of the terms in the series, Equation (3) must be taken. The Leveque solution is a good approximation for such cases. As shown by Tribus and Klein,2 the wall temperature and heat flux for such a case are related by en, ((5" (X- 5) t( (72) 3 Cf4) v ip7kJI o,0 J 0tw(s) and 31r(t1) ( w (73) For flat ducts (. ) = 3 / (74) for round ducts -,tu 44/,(75) Substitution of Equation (75) into Equations (72) or (73) (and noting that the mixed mean temperature of the fluid is essentially equal to its inlet value at small values of x*) gives for the three cases under consideration: For constant wall temperature:. /3 3 +/ _ -a _ -,1.3S-S' X + < o. ol (76) J' r(/3) For constant heat flux " s"3 r(3) +-'~ -3 hlua S - "3 _ — 1= - 1.6373 X - (77) For linearly varying wall temperature 3. X=fS (78) WC T ICA REPO = R &"/ r'({%) WADC TECHNICAL REPORT 54-255 18

Figure 1 shows a graph of the functions Ro, R1, and R2 compared with solutions given by Jakob. Figure 2 shows the variations in Nusselt modulus for three cases (1) wall temperature constant, (2) heat flux constant, and (3) wall temperature increasing linearly along the pipe wall. The Nusselt modulus is defined by the equation NL. = O (79) i, (X) - t_'. The mixed-mean temperature, tmm, is determined by integrating the heat flux from the origin (X+ = O) to the position where q(x) is known. CONCLUSIONS The methods used in this paper have a wide applicability. For example, the liquid metals systems analyzed by Poppendiek4 could be treated by the methods used here. The unsymmetrical boundary conditions treated by Yih and Cermak5 can also be readily treated by these methods. The authors are somewhat surprised at the fact that whereas the asymptotic formulae are all supposed to be valid only for very large X, in actuality values of n as small as 4 seem to give excellent results. The reasons for the good results are not now clear. APPENDIX A The equations for a flat duct system with walls at y = + b (A-l) Defining Re = 4Um pb/[, x+ = (x/b)(RePr)-1, y+ = (y/b) the equation to be solved is 3 vt = I t(A-2) which has a solution WADC TECHNICAL REPORT 54-255 19

Q 1. Q R0 FIG. i p'.8 0 -.86V _- ^ EIGEN FUNCTIONS FOR GRAETZ SOLUTION -.6 -I~ -I JAKOB -I. 0 _ _ _ _ _ _ _ _ 0.1.2.3.4.5.6.7.8.9 1.0 PRESEN PAPE

0Q 1.00 S I ot - 3 17 - i. TVX+)-To=Ax+ x3 I x 2. q(x+) = constant E E 3Fig.2. Tw(xaminr -Tlow of Constant CJ ui 10 -3.0001.1:3.01. HI-u ~ ~ ~ ~ ~ F 092 d Fig.2. Laminar Flow of a Constant Property Fluid in a Round Tube

- o= Z K 1 Ue (A-3) satisfying ~ = 1 at x+ = 0, 9 + O, x + O, if y(y+) satisfies Y "+ X - t)Y=o (A-4) with Y'(O) = Y(l) = O, Y(O) = 1. Xn is the value of X to permit Yn(1) = O. The coefficients Kn are given by -, 1< OV ) (A-5) By the methods in the text the WKB approximation is found to be Y(t) = c (I....). ( (A-14) for 0' y+ < 1. Defining z = 1 - y+, the solution of A(A-4) for z << 1 is found to be y 3, ( ( f * (' -((^ z- 3( - S T ( z,) The eigenvalues are \X = Lt+ /3 (A-25) ( X t ( / IT, (A-29) K, (-3) 3 - _ _r' - (A-30) 3, r (13) WADC ^1TE CHNICA =REPORT (54 )+(A-51) 3 (< -l= i It t 2 Y - 8 (A-32) WADC TECHNICAL PEPORT 54-255 22

To obtain the fluid temperature for a given heat flux use e -zto = ) S? -a'd)'6 )'7 (A-44) The integrating kernel, g, is given by 2 - e _ z K Y^(*)s at y^ (A-48) e- e where the -y2m are the zeroes of H 15) - - Y - ad (\ ) _(A-47) H (s) = - (_ _. REFERENCES 1. Jakob, Mas, Heat Transfer, vol. 1, Wiley, 1949. 2. Tribus, M., and Klein, J., "Forced Convection from Nonisothermal Surfaces", Heat Transfer: A Symposium Held at the University of Michigan during the Summer of 1952, University of Michigan, Engineering Research Institute, 1953, p. 211-35. 3. Graetz, L., "Uber die Warmeleitungsfahigkeit von Flissigkeiten", Ann. d. Phys. Chem., 25, 337-357 (1885). 4. Poppendiek, H. F., "Forced Convection Heat Transfer in Thermal Entrance Regions, Part I", Oak Ridge National Laboratory, Tenn., ORNL-913, Series A, Physics, March, 1951. 5 Yih. S., and Cermak, J. E., "Laminar Heat Convection in Pipes and Ducts" Civil Eng. Dept., Colorado Agriculatural and Mechanical College, Fort Collins, Col., September, 1951. (ONR Contract No. N90 nr 82401, NR063-071/119-49). WADC TECHNICAL REPORT 54-255 23

UNIVERSITY OF MICHIGAN 3 9015 03524 4410