ENGINEERING. DESIGN, AND PROCESS DEVELOPMENT Literature Cited
(1) Donath, E. E., Trans. Am. Inst. Min. Engrs., 193, 381-5 (1952). (2) Lombardi, A. D., and Nickels, J. E., Mellon Institute and Koppers Co., Pittsburgh, Pa., investigations under terms of cooperative agreement between U. S.Bur. Mines, Demonstration Branch, Louisiana, N o . , and Koppers Co., unpublished data. (3) Orchin, Milton, and Storch, H.H., J . Soc. Chem. Ind. (London), 69, 121-2 (1950).
( 4 ) Petroleum Processing, 7, 1154-60 (1952). (5) Pier, bl., Wissel, K., and Oettinger, W., Erdol u. Kohle, 6, 693, 696 (1953). (6) Woolfolk, E. O., and associates, U. S. Bur. Mines, BUZZ. 487,
3, 1950.
ACCEPTEDJuly 9 , 1954. RECEIVED for review April 23, 1954. Presented as part of the Symposium on Synthetic Liquid Fuels and Chemicals before the Division of Gas and Fuels Chemistry, a t the 125th Meeting, ACS, Kansas City, Mo.
Digital Computer Solution for Heat Transfer to Temperature Probes H. F. KRAEMERI
AND
J. W. WESTWATER
University o f Illinois, Urbana, 111.
W
HENEVER possible, scientists develop exact analytical solu-
tions for problems that arise, but when equations of unusual mathematical complexity are faced, some other attack may be needed. Graphical methods or numerical methods may be used, although the drudgery involved is sometimes discouraging. The digital computer is a powerful tool for attacking problems by numerical techniques. Its value is a result of the high speed with which it can perform. A typical problem that cannot be solved analytically is thnt of a fin or a probe that is undergoing heat transfer by simultaneous convection, conduction, and radiation. The temperature probe is of particular interest. With a probe, a reading of the probe-tip temperature is obtained, but the measurement desired is the temperature of the fluid into which the probe is inserted. This paper is concerned with obtaining exact solutions to the heat transfer equation that can be set up for a simple probe. The exact solutions are compared with approximate solutions presented in a prior paper ( 4 ) .
Equation 1 can be generalized to obtain dimensionless limits of 0 and 1for positions on the probe and dimensionless groups for the other variables. The appropriate substitutions in Equation 1 yield Equation 2.
Equation 2 contains five dimensionless variables. It is a secondorder differential expression and contains the fourth power of the dependent variable, as does Equation 1. Although Equation 1 contains 10 variables, the five variables of Equation 2 are sufficient to describe the problem. I t is possible to reduce the order of the equation by the common technique of letting a new variable be defined, say p = de/dl. A substitution yields a first-order equation which is formidable, because it is not linear. No exact analytical method is available for integrating the new equation, but if a digital computer is available, it can handle the two simultaneous equations de ;fi = P
Digital Computer Handles Differential Equations Relating Fluid and Probe Temperatures
(3)
and
Consider a temperature probe inserted through a duct wall into a fluid. The probe is undergoing heat transfer by conduction, convection, and radiation. The temperature-sensitive element is assumed to give a true reading of the temperature of the tip of the probe. The problem is to calculate the temperature of the fluid. For liquids the problem is often trivial, but for gases the gas temperature may be very different from the probe temperature and the problem becomes important. A short element of the probe, at distance 1: from the tip, may be considered. A heat balance for the length, dx, gives
and thus it is not necessary to deal with one highly complex expression. "eat Transfer at Tip. The heat transfer a t the tip of the probe is important for thick probes and negligible for slender probes. In any case, it is easy to include this factor in the problem to be handled by a computer. A heat balance for the tip area gives
Equation 1 states that the heat received by convection is equal to that lost by conduction and radiation. The equation is general and holds even if the gas is colder than the probe. The fluid is treated as if it were transparent; the base of the probe is assumed to be at the temperature of the durt \Tall, 5°F. The equation cannot be integrated directly. An approximate solution is available ( 4 ) . Bn exact solution can be obtained by using numerical methods with a digital computer. This attack was used.
where ( d T / d x ) z is the temperature gradient a t the tip. Various writers have simplified the tip area problem by assuminga value of zero for ( d T / d x ) ~ .This will not be done in the present treatment. -4much more reasonable assumption is made-namely, that the temperature gradient a t the tip is the same for all positions on the tip. Changing Equation 5 to the generalized, dimensionless variables gives Equation 6, which represents the boundary conditions that at e = BE, p = p ~ .
1
Present address, Ethyl Corp., Baton Rouge, La.
October 1954
(4)
pE
=
Orl*(@E4
- 1) - P I * ( e G - B E )
INDUSTRIAL AND ENGINEERING CHEMISTRY
(6)
2035,
ENGINEERING, DESIGN, AND PROCESS DEVELOPMENT Note that p is the ratio of the tip area to the side area of the probe. Arbitrary values of ,u muet be selected to allow the computer t o handle the additional information. The other boundary conditions necessary for the solution of Equations 3 and 4 are that at the base of the probe the temperature is the same as the duct temperature and that at the tip the temperature is defined as TE. I n the generalized variables
Sets of Numerical Values for Generalized Variables Are Computed
The information represented by Equations 3, 4, 6, 7 , and 8 was fed into the digital computer. A total of 99 sets of numerical values for the generalized variables was selected, and for each case the fluid temperature was calculated by the computer. The range of values selected are listed below. 6~ (ratio of tip temperature to base temperature) from 1.015 to
2.0
01 (radiation dimensionless group) from 0 to 80 p (convection dimensionless group) from 0.36 to 90 p (ratio of tip area to side area) from 0.0 to 0.4
In reality, 01 compares the relative magnitude of heat transfer by radiation to that by conduction, p compares the convection to conduction, and CY/@ compares radiation to convection. The range of 01/b was from 0 to 100. The calculated absolute temperature of the fluid waa found to have a range of from 1.0307 to 7.15 times the absolute temperature of the duct walls-that is, 00 varied between these limits.
hc hC
Figure 1.
8
DIMENSIONLESS
+hR
Range of Variubles
Since it is axkward to visualize the meaning of the dimensionless ranges, Figure 1 illustrates the values in familiar terms. All possible combinations of the wall temperature, tip temperature, and fluid temperature can be represented by ordinate values from 0 to 1. All possible combinations of convection and radiation heat transfer rates can be represented by abscissa values from 0 to 1. All real cases must lie in Figure 1, although the particular values of probe length and tip area cannot be conveniently shown on the graph. The cases solved on the computer are indicated in the figure. They represent a planned survey of the infinite number of possible cases. In Figure 1 a value of 1 for the abscissa corresponds to cases with no radiation occurring. An abscissa of zero corresponds to no convection. Values of the 45 O line can be shown to correspond to no conduction. Values inside the triangle correspond to cases that include all three modes of heat transfer.
2036
Runge-Kutta Method. The numerical solutions of the heat transfer equations were made with a large scale electronic computer, the ILLIAC, located a t the University of Illinois. The two equations, Equations 3 and 4, were integrated simultaneously using as boundary conditions Equations ti, 7 , and 8. Since boundary conditions are specified at each end of the probe, an iterative process must be used to obtain the eigenvalue, SQ, that satisfies all conditions. For a specified set of thermocouple parameters, 01, 0,,u, and a thermocouple temperature, BE, the computer selects limiting values of the gas temperature, Bc. A temperature profile of 0 versue I is generated with the boundary conditions including the slope, p~ (from Equation ti), and with the starting point being the end temperature, 68. Forty intervals were taken along the length of the probe, and each step of the temperature profile was calculated by the method of Runge-Kutta (1, 6). In general, the temperature a t the wall, ew, obtained a t the end of the integration does not equal the boundary value of 1. The computer then iteratively selects closer estimates of BO and generates corrected temperature profiles until the boundary condition a t the wall is satisfied. Convergence is assumed to occur when high and low The computaestimates of either 6~ or 6w deviate less than tions were carried out to eight significant figures. The calculation of 00 is completed by the computer in about 25 seconds and involves approximately 10 iterations. The roundoff and truncation error6 of the numerical method can be shown to no greater thnn 0.04% of the absolute temperature. Computer Solutions Are Used To Test Assumptions of Other Investigators
Inasmuch as problems solved by a digital computer are solved by the manipulation of actual numbers, the results consist of numerical answers to a finite number of cases. Difficulty arises in reaching generalizations about the solutions. Since the final tables of numbers for the present problem are space consuming, they are not given in this paper. However, the tables do lead to eome clear conclusions. -4n approximate solution to the probe problem was derived previously by West and Restwater (4). The solution was based on assumptions that the probe tip area was negligible and that the heat transfer coefficient for radiation, h ~was , constant along the probe length. Neither of these two assumptions was necessary for the digital computer; therefore, the reliability of the assumptions could be tested. It was found that the assumption of constant hB leads to values of the calculated fluid temperature that are too low in all cases. However, the error in the absolute gas temperature was not objectionable, being no more than 0.2% under the most extreme conditions. The other assumption, that of negligible tip area, leads to calculated gas temperatures that are also too low. The discrepancy amounts to several per cent (on the absolute temperature scale) for the extreme cases. This error is enough to be undesirable. Prior writers have suggested that the probe tip (or fin tip) problem might be simplified for ease of mathematical attack by the introduction of two compensating suppositions. The tip area is treated as if it were zero and this is balanced by treating the probe as if it were longer than its actual size. Harper and Brown ( 2 ) state that the additional length should be sufficient to increase the side area by an amount equal to the actual tip area. Jakob ( 3 )also treats the problem, using the method of finite differences. In each case, hon-ever, radiation is not included in thc mathematical treatment. The equation of West and 1Test.rvater can be niodified t o treat the tip area by the scheme above of Harper and Brown. The equation then becomes
INDUSTRIAL AND ENG INEERIWG CHEMISTRY
VOl. 46, No. 10
ENGINEERING, DESIGN, AND PROCESS DEVELOPMENT eE=i.t 13 = a
ee=i.03
2-
/3
d
d
n1.44 =0.072
0.2-
d '80
USING b C ll ttlZ# # U l)
0.1
-
USINCb(lt/4)
1
- -- - - -- - -
-DIGITAL 7 - -COMPUTER ----
a
DIGITAL COMP./
-I
5
-
-
USING b
___(*D
-
-0.1
\WCW,
K
;
-2
USING L E N C T W b
Y
-2
-2
-0.2-
c I
-3 0
0.1 /A,TIP AREA
-3 0
0.3
0.2
+ SIDE A R E A
0.2
0.1
-3 -3
0.3
o
ai
0.2
P
Figure 2.
-0.3 -0.3
0
0.1
0.2
0.3
0.4
A
C
B
A
0.3
P
D
Comparison of Methods for Handling Probe Tip Areas W & W = West and Westwater ( 4 )
+
b' is the fictitious probe length defined as b' = b( 1 p), where b is the actual probe length and p is the ratio of the tip area to side area. Values of the gas temperature were calculated by Equation 9 and were compared t o the digital computer values. The improvement in accuracy is striking. In not one of the 99 cases solved did an error as large as 1% occur. In some cases the tip area was as great as 40% of the side area. Figure 2, A-D, shows the reliability of Equation 9, compared to the digital computer solution of the problem. Four sets of conditions are considered, and in each set the probe tip area is allowed to vary from 0 to 32% (or more) of the side area. The original equation of West and Westwater ( 4 ) yields negative errors (gas temperature too low) in all cases. The improved equation, Equation 9, decreases the error greatly. The error becomes zero a t some value of p and then becomes slightly positive. The advantage of using the fictitious length is obvious. A few calculations were made using a revised length equal to b(l 2 p ) t o see if this might prove of value. The graphs show that this gives an undesirable overcorrection. Figure 2, D, which represents extreme conditions of large radiant heat transfer, shows that the method of correcting for the tip area is not crucial for this case; in fact, no correction is needed. It is recommended that Equation 9 be used for probe calculap). For a tions. The probe length should be treated as b ( l circular rod this is equal to the actual length plus one half of the radius.
+
+
m would usually occur. The expression fits only the curious case of a single fin surrounded by enclosing walls that are at the same temperature as the base of the fin.
Nomenclature
a = transverse cross-sectional area, sq. ft. A E = area a t tip, sq. ft. b = length of probe or fin, ft. b' = b ( 1 H ) , ft. F A = radiation angle factor, dimensionless FE = radiation emissivity factor, dimensionless It, = convection heat transfer coefficient, B.t.u./(hr.)(sq.ft.)
+
hR
k
Z
m n
p p~ P
Q
T TE TG Tw
5 01
P
A fin having a rod shape can be treated by Equation 9 to calculate the tip, base, or fluid temperature. The temperature profile along such a fin can be calculated from Equation 10.
T G- T w -= T - Tw
("___h R )
cosh nb' (cosh nb' - cosh nx
=
(h,P/n)(Tc - T w ) tanh nb'
=
= = = = = = = = =
6G P
=
U
or fin, dimensionless
= perimeter of transverse cross section, ft.
= =
heat transfer rate, B.t.u./hr. temDerature of Drobe a t location x. ' R. temperature a t tip, ' R. temperature of fluid, ' R. temperature of enclosing duct or walls, O R. distance along probe from tip, ft. uFAF.PbZTwS/(ka), dimensionless hoPb2/(ka),dimensionless T / T v , dimensionless T z / T w , dimensionless T G / T w ,dimensionless ratio of tip area to side area = An/(Pb),dimensionless Stefan-Boltzman constant = 0.173 X 10-8 B.t.u./(hr.)(sq.ft.)(' R.4)
(1) Gill, S., Proc. Cambridge Phil.SOC.,47, 96 (1951). (2) Harper, D. R., and Brown, W. B., Natl. Advisory Comm. Aeronaut., Rept. 158, 1922. (3) Jakob, &Phil. I.,Mag., 28, 571-8 (1939). (4) West, W. E., Jr., and Westwater, J. W., IND. EXG.CHEM.,45,
215243 (1953).
( 5 ) Wilkes, M. V., Wheeler, D. J., and Gill, S., "Preparation of Programs for an Electronic Digital Computer," pp. 32-3, 56-7, 86-7, 132-4, Addison-Wesley Press, Cambridge, 1951.
(11)
This equation is the same as the ordinary one for fins experiencing no radiation, except for the appearance of the term, n, where
October 1954
+
= d(hc hR)P/(ka),ft.-l = d@/dl,dimensionless = value of dO/dZ a t tip of probe
Literature Cited
The digital computer results in this paper suggest that the equation should be satisfactory. The heat transferred through a fin is of interest. The portion flowing by conduction at the base of the fin is obtained by differentiating the Equation 10 to find the temperature gradient at the base.
Q (at base)
thermal conductivity of probe or fin, B.t.u./(hr.)(ft.)(' F.) = x / b , dimensionless = d h , P / ( k a ) , ft.-l =
BE
6 Method M a y Be Applied to Fins
c p F.)
= radiationheat transfer coefficient,B.t.u./(hr.)(sq.ft.)(' F.)
RECEIVED for review February 19, 1954. ACCEPTED May 14,1954. Presented before the Division of Industrial and Engineering Chemistry a t the 126th Meeting of the h f E R 1 C A N CHEMICAL SOCIETY, New York, N. Y.
INDUSTRIAL AND ENGINEERING CHEMISTRY
2037
