INVERSE PROBLEMS RELATED To” TRANSIENT CONVECTION :INIHE‘ ’- ; THERMAL ENTRANCE REGION’BEIWEE‘N; 15-52}.:§:;.:’[":_3.1L3 ' PAMLLLL PLATES * ._ . I ' "es“ t‘°’ L“L'Dee‘gree of Ph, D." g ’ ' “ ' .MLCH-IGANSTAIE u-mvmsm- a RAYMOND 78,;00LLADAY; . . :1969 ' ' LIBRAP" THESIS Michigan State"; Universmy This is to certifg that the thesis entitled INVERSE PROBLEMS RELATED TO TRANSIENT CONVECTION IN THE THERMAL ENTRANCE REGION BETWEEN PARALLEL PLATES presented by RAYMOND S. COLLADAY has been accepted towards fulfillment of the requirements for Doctora] _ Mech. Engr. _..___ degree m —_ Major profes or Dam September 10, 1969 0-169 ABSTRACT INVERSE PROBLEMS RELATED TO TRANSIENT CONVECTION IN THE THERMAL ENTRANCE REGION BETWEEN PARALLEL PLATES BY Raymond S. Colladay Transient heat transfer in the thermal entrance region between parallel plates with a fully developed laminar velocity profile is studied for low Peclet number flows. Axial heat conduction is in- cluded and its effect on the temperature profiles upstream of the heated region as well as in the entrance region itself is noted. The energy equation is solved numerically using an alternating direction implicit finite difference method. Two boundary condition cases are presented, (1) a uniform wall heat flux and (2) a uniform wall temperature. The plate boundary upstream of the heated region is insulated in both cases. Steady state and transient results are presented for a range of Peclet numbers between 1 and 50. These results: have not been previously obtained. Using this solution as a basic building block in superposition, a number of inverse problems are formulated and presented in terms of integral equations. Particular emphasis is given to the inverse problem of estimating the mean velocity from temperature measurements at the wall and the solution of the energy equation. The optimum wall heat flux profile and Optimum axial location of wall thermocouples which give the best estimate of this velocity parameter are studied. An example of the nonlinear estimation procedure used in calculating Raymond S. Colladay the mean velocity is also presented. This work is the first to formulate these inverse problems related to convection. INVERSE PROBLEMS RELATED TO_TRANSIENT CONVECTION IN THE THERMAL ENTRANCE REGION BETWEEN PARALLEL PLATES By 4 l L' |' Raymond Si(Colladay A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 1969 Coco/330 L; — L. 7 0 ACKNOWLEDGEMENT The author wishes to express his sincere appreciation to his major professor, Dr. James V. Beck, for his valued guidance and encouragement during the period of research and graduate study. He also wishes to thank Dr. Charles R. St. Clair for his helpful advice throughout the years at Michigan State University. Thankful acknowledgement is extended to the other members of the Guidance Committee: Dr. John F. Foss, Dr. Merle C. Potter, Dr. David H. Yen, and Dr. George E. Mase. The author is grateful for financial support received from a NASA Traineeship and supplementing funds from the University. To his wife Judy, the author dedicates this dissertation for her understanding and cooperation during graduate study and research. ii TABLE OF CONTENTS Page Chapter I . DESCRIPTION OF THE PROBIEM ................... 1 II. MATHEMATICAL DESCRIPTION . . . ............. . . . . . 6 III. NUMERICALILIETHOD...” ....... ..... 13 IV. DISCUSSION OF RESULTS FOR CASES I AND 11 ..... 31 V . OPTIMUM EXPERIMENT FOR DETERMINING THE MEAN VEIOCITY FROM THERMAL MEASUREMENTS ...... . . . . . 52 VI . SOME INTEGRAL INVERSE PROBLEMS ............... 78 VII . CONCLUSIONS ................ . ................. 9O BIBLIOGRAPHY...................... ..... 93 Appendix A. TRUNCATION ERROR 98 B. STABILITY ANALYSIS FOR ADI NUMERICAL METHOD . . 101 C. DERIVATION OF FULLY DEVEIOPED TEMPERATURE PROFILES AND SENSITIVITY COEFFICIENTS . . . . . . . . 104 iii LIST OF FIGURES Figure Page 2.1 Plate Geometry ..... ................................ 7 3.1 Finite Difference Spatial Grid for Parallel Plate Geometry 0.0.0.0... 000000 O0.0.0....0.00.00.00.00...O 14 3.2 Y-Direction Calculation ..................... ....... 22 3.3 x-Direction calculation 00.00.000.00... .......... 0.. 23 3.4 Energy Balance on A Differential Element ........... 25 3.5 Development Length - Case I ...... ..... .. ....... .... 29 4.1 Steady State Bulk Temperature in the Entrance Region 000......O0....OOOOOOOCOOOCOOOOOOOO... ..... .0 32 4.2 Steady State Temperature Profiles at x = O - Case I 34 4.3 Steady State Nusselt Number in the Entrance Region - caSeI OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO0. 35 4.4 Functional Dependence of Nusselt Number on Peclet Nuwer .00...0.00.0...OOOOOOOOOOOOOOOOOOOO0.0.0.0000 36 4.5 Steady State Temperature Profiles in the Entrance Region-casel 0.0.0000000000000000000000000000...0 38 4.6 Transient Nusselt Number in the Entrance Region - caseI OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO0.0.0...0 40 4.7 Nusselt Number in the Entrance Region at Various Times-C8881 oeoo00000000000000.0000.-oeoooo ...... 42 4.8 Steady State Bulk Temperature in the Entrance Region-C3861]: 0.000000000000000. ccccc oooooo-ooooo 46 4.9 Steady State Temperature Profiles at x = O - case II .00...OOOOOOOOOOOOOOOOOOOOOOOOO0.. .......... 47 4.10 Steady State Temperature Profiles in the Entrance Region for Pe = 6 - Case II ..... ........ .......... 48 4.11 Steady State Temperature Profiles in the Entrance Region for Pa = 50 - Case 11 ................. ..... 49 4.12 Steady State Nusselt Number in the Entrance Region - case II 000......0.0000000000000000000000000....O... 50 4.13 Fully Developed Transient Nusselt Number - Case II . 51 4.14 Hydrodynamic Entrance Region for a System of Parallel Plates O..00...OOOCOOCOOOOOOOOOOOOO00....O. 44 iv 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.1 6.2 6.3 C.l C.2 Illustration of Sum of Squares Function ............ Relation of the Objective Function to Heating Length and Thermocouple Location ................... Surface Temperature for a Finite Heating Length .... Steady State Sensitivity Coefficient for Case I Boundary condition OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO. Transient Sensitivity Coefficient for Case I Boundary condition 0..OOOOIOOOOOOOOOOOOO00.0.0000... Illustration of Discrete Optimum Heat Flux Profile . Transient Wall Heat Flux Which Maintains a MBximum Surface Temperature Rise of ATmax ................. Objective Function in the Entrance Region for Optimum.Transient Heat Flux ........................ Transient Objective Function for Optimum Transient Heat Flux .0.00....0.0000000000000000000000000...... Flow Between Heated Parallel Plates ................ Differential Element in a Conducting Plate ......... Illustration of Radiation Between Parallel Plates .. Energy Balance on a Fluid Element ...... .......... .. Entrance Region Control Volume ....... ......... ..... 57 59 62 64 65 68 7O 71 72 79 83 87 105 106 LIST OF TABLES Table Page 5.1 60 vi Re Pr Pe h CX NOMENCLATURE_ Temperature Velocity components Physical coordinates Time Half-plate separation distance Absolute viscosity Specific heat Thermal conductivity Density Thermal diffusivity = k/pCp Mean velocity Reynolds number = pua/u Prandtl number = uCp/k Peclet number = RePr = ua/a = Local heat transfer coefficient Nu = Local Nusselt number = hCXa/k T s Dimensionless physical coordinate x/aPe Dimensionless physical coordinate y/a Dimensionless time = ta/a2 Dimensionless velocity = u/u Wall heat flux Uniform upstream temperature and initial temperature Wall temperature + II T = T-TO/fli2 = Dimensionless temperature for uniform heat flux boundary condition T+ = T-To/Ts-To = Dimensionless temperature for uniform wall temperature boundary condition Tb Bulk temperature T = Duration of experiment x = Location of wall thermocouples K II L Heating length Tmax = Maximum temperature constraint “Tmax - Max1mum temperature rise constraint - Tmax-To § = Objective function - basic criterion to be extremized Pe aT/aPe = Sensitivity Coefficient with respect to the Peclet number. 1. DESCRIPTION OF THE PROBLEM 1.1 Introduction The problem of extracting information from a set of measured data obtained under optimum experimental conditions has recently gained wide interest. In the area of convection, where it may be difficult to measure the velocity with a probe, parameters associated with the velocity profile can be estimated from thermal measurements at the boundary - an application which has considerable appeal. An optimum experiment insures that the parameters are determined more accurately than by any other similar experiment with the same errors in the temperature measurements. The solution of the partial differential equation describing a process is often labeled as an inverse type problem.when conditions are overspecified at a given boundary. The energy equation, with both the heat flux and temperature at the boundary known and the velocity an unknown coefficient to be determined falls into this category. In this work, the basic objective is to investigate inverse problems related to the thermal entrance region between parallel plates under transient conditions. Consequently, an accurate transient solution to the thermal entrance region problem is needed, particularly in the liquid metal-small Peclet number range. Many extensions to the original Graetz problem (1) of steady laminar heat transfer to a constant prOperty fluid in the thermal 1 entrance region of a constant wall temperature tube with a parabolic velocity profile have appeared in the literature. The corresponding flat plate geometry has been considered by various investigators (2, 3, 4); their work.was followed by the treatment of such additional effects as transient conditions (5, 6), uniform heat flux boundary condition (7-9), hydrodynamically undeveloped flow (7-15), nonuniform prOperties (16, 17), wall surface resistance (18, 19), flow in an annulus (20), viscous dissipation (ll, 21), power velocity profiles (22, 23), and nonuniform boundary conditions (24). In all of these studies, the one fundamental assumption made is that the fluid conduction in the axial direction is negligible. In many cases this assumption is justified. However, when the Pe number is small (say less than 100), 93 when results in the entrance region in the immediate vicinity of, and including x = 0 (location where heating begins) are of interest, these solutions are not applicable. Another case where axial conduction should not be neglected is wherever the boundary conditions at the plate vary substantially over small regions. In many modern heat exchanger devices envolving large heat fluxes such as nuclear reactors, the use of liquid metals has become more popular. Due to the good thermal conductivity of these metals, solutions that include axial conduction are required. In reviewing the literature, it is apparent that very few publications treat this problem. References (25-31) include the axial conduction term in the energy equation, but apply a uniform temperature boundary condition at x = O. This forces the solution at this location to be the same as for the case of negligible axial conduction. However, it is at the start of the heating section, more than anywhere else in the entrance region, that the effect of axial conduction is most significant. Some authors (32-36) have considered the "complete" problem; allowing axial conduction upstream of the x = 0 location by assuming a uniform temperature profile far upstream of the heated region. A uniform velocity profile is assumed in (32-34). Only Agrawal (35) and Hennecke (36) treat this case for a parabolic velocity profile; both under steady state conditions. The former assumes a flat plate geometry and the latter treats a tube geometry. Using analytical methods, Agrawal treats a constant temperature boundary condition at the wall; that is, the regions x > O and x.< O at the wall are maintained at constant temperatures Ts and To respectively. In addition to this boundary condition, Hennecke treats the wall con- dition of constant heat flux for x 2 0 and insulated upstream. An analytical solution of the energy equation assuming a fully developed parabolic velocity profile leads to a differential equation whose eigenfunctions are not orthogonal. Representing the eigenfunctions by infinite Fourier sine series, this difficulty is reflected in the fact that to date only the first five eigenvalues have been reported, (35). These few values are not sufficient for accurate evaluation of the temperature distribution in the neighborhood of X = 0 since at this location, the convergence of the series is very slow. Hence, Agrawal's solution in this region is not reliable. The methods of both (35) and (36) require numerically matching the temperature solutions from the upstream and downstream regions at the interface x = 0. In transient calculations this could require a substantial amount of computer time. To the author's knowledge, no results for the complete problem of transient laminar flow between parallel plates have appeared in the open literature. Once the wall temperatures are available from the solution of the entrance region problem, an optimum experiment can be designed to estimate the velocity (or some parameter of the velocity profile) from temperature measurements at the wall. The experiment is optimized with respect to the heat flux profile and location of thermocouples at the wall for a given duration of the experiment and maximum temperature rise. Optimum conditions permit the estimation of the parameters more accurately than any other similar experiment with the same constraints. Nonlinear estimation - the procedure for calculating unknown parameters in a differential equation describing a physical process - has its first known reference by Gauss (47). The method, as applied to Optimum experimental conditions, has been deveIOped from a statistical approach by Box and coworkers (48-50) and recently extended to applications involving partial differential equations by Beck (51-53). 1.2 Problem Description In this dissertation, the thermal entrance region between parallel plates is studied. The energy equation, with axial con- duction and transient effects included, is solved numerically using an alternating direction implicit (ADI) method. Assuming a hydro- dynamic development section upstream of the heated region, the velocity profile is parabolic. Two sets of boundary conditions at the wall are treated, (I) a uniform heat flux is applied to both plates in the region x 2 0 (where x is the spatial coordinate in the axial direction) at time t = O, and insulated in the region x < 0; and (2) the plates are maintained at the same uniform temperature in the region x 2 0 beginning at t = O and insulated upstream. Axial conduction is allowed upstream by assuming a uniform temperature profile at x = «n. The ADI numerical method does not require matching the solutions to the upstream and downstream regions at x = 0. Both the regions are treated together as a single domain. Since the thermal conductivity and the velocity are assumed independent of temperature, the energy equation is linear in x and t, making superposition of solutions possible. A number of inverse problems are formulated using Duhamel's superposition integral with the uniform heat flux boundary condition solution as a basic building block. Particular emphasis is given to the inverse problem of estimating the mean velocity from temperature measurements at the wall. A basic criterion is developed which gives a measure of the effectiveness of an experiment for determining the mean velocity. This criterion is then maximized with respect to the heat flux profile (a function of x and t) and wall thermocouple location to establish an optimum experiment. The effect of errors in the temperature measurements is in- vestigated. Also, an example using the nonlinear estimation procedure under optimum conditions is presented. II. MATHEMATICAL DESCRIPTION 2.1 Energy Equation Transient convection for fully developed laminar flow of an incompressible viscous fluid with constant properties is considered in the region between parallel plates. Two cases are investigated, Case I: a constant wall heat flux for x 2 0 (see Figure 2.1) and insulated upstream, x < 0; and Case II: a constant surface temper— ature for x 2 O and insulated upstream. For negligible viscous dissipation and flow work and constant thermal conductivity, the equation expressing energy conservation is, T a: a1- 2T bi pcpg: + u ax + v 3y) — 14:: + ayz) (2.1) For fully developed flow, u and v are given by, u = % G[1 — ($21 (2.2) v = o (2.3) a - '_l where the mean veloc1ty, u - a $0 udy (2.4) Equation (2.1) is linear in temperature, T since k, p, Cp, and G are assumed to be independent of T. Since the thermal entrance region is of particular interest, the axial conduction term azT/ax2 must be retained. Solutions to Equations (2.1), (2.2), and (2.3) are obtained for the two sets of boundary conditions on T(x,y,t) given below. ‘ /////A l L l i L ’L X f/f/W T 11 T f UNIFORM HEAT FLUX - CASE I ///4/1 S y L. X 7///// T S UNIFORM WALL TEMPERATURE - CASE II 2a FIGURE 2.1: PLATE GEOMETRY CASE I: UNIFORM HEAT FLUX T(X.y.0) = T0 T(-°°.y.t) = TO 2 a T(”1_Y2t) = 0 5x2 LTKX goitl = 0 BY kaT(X,a,t) = 0 for X < 0 By q" for x 2 0 CASE II: UNIFORM WALL TEMPERATURE T(X.y.0) = TO T('°°9Y:t) = TO AT («.550 = 0 ax BT(X20:t) = 0 BY T(x,a,t) = TS for x 2 O aT(:;a,t)= O for x < 0 It is convenient to express the energy equation and its (2.5a (2.5b (2.5e (2.5d (2.5e (2.6a (2.6b (2.6c (2.6d (2.6e (2.6f boundary conditions in dimensionless form. Let the variables be non- dimensionalized as follows, using the half plate separation distance, (1) a, as the characteristic length. l ( ) Note that some references for the characteristic length use the hydraulic diameter, defined D = 4x(Flow Area/Perimeter). the parallel plate case, D = 4a. hen comparison is made with results based on DH’ a factor of 4 will be introduced. for Case I for Case II o for x+ 0 (2.10e) + + + 5T (X $1,: ) = o for x*'< o (2.101) BY 2.2 Other Useful Relations The local Nusselt number, Nu, used in expressing the local heat transfer coefficient, hx’ in dimensionless form is of particular . . . + interest, espec1ally 1n the entrance region near X = O. The defining equation for hx is, q" = hx(Ts - T ) (2.11) b where Tb is the bulk temperature as defined by Equation (2.16). Also, the heat flux is given by, .. = 3:1; 2.12 q k(ay)y=a ( ) The Nu number is then formed as follows, g"a as; h a k a(a ) =a (2.13) k T -Tb (TS-Tb) Expressing Equation (2.13) in terms of the appropriate dimen- sionless variables of Equation (2.7), we have: 11 for Case I: Nu = l = l T{To Tb-To T+-T: 3"a ' g"a S k k or based on DH = 4a, 4 Nu = T+'T+ (2.14) s b for Case II: T-T a 0 + ( _ ) 5T 5(1) Ts To 1=1 ( 4) +_ Nu = a a = M T -T T -T + s o) _ b o) l - Tb (T -T (T -T s o s o and based on D H + 4 d1. + +_ Nu = 5y 1 " (2.15) + l - Tb The bulk temperature, as used in the above equations, is defined as T ‘lr' dA b(x’t) AGIAUT c c c where Ac is the cross sectional area normal to the direction of flow. For the given geometry, =L aquy (2.16) Tb ad 0 For the particular boundary conditions of this problem, Tb in dimensionless form becomes, using Equation (2.7), Case I: T 12 II Since it [a udy = l and g_§_= constant. Therefore, u o k T-T +'_ l. a .__Jl b afi j‘o ( g"a)dy k + _ 1 +1+ + or Tb - f0 u T dy (2.17) Case II: T+ = T-To = 1_ Ya T _ 1_ Ya To y b Ts-TO a6 0 u T -T ad 0 T -To and again, + 1+++ =1 Tb o u T dy (2.17) The Nusselt number as defined in Equation (2.13) is a function of x and t. One can consider an average Nu over an axial length L, or time T, or both. §E(c) =-%.f: Nu dx (2.18) fi3(x) = %-j; Nu dt (2.19) --_ l_ L T Nu.-T1‘IOIO Nu dt dx (2.20) For both cases, the Nusselt number is zero for x < O. In the region x 2 O, a numerical method of solution is used and is described in the next chapter. III. NUMERICAL METHOD 3.1 Finite Difference Expressions Expanding a function f(x,y,t) about a point x in a Taylor Series gives for (x + Ax,y,t), 2 2 f(x+Ax,y,t) = f(x,y,t) + Ax W + (£21?) a f(Xiy,t) . 8X 3 3 + mx) 5 mew) + (Axf‘ aafémuti + (3 1) 3! 3 4! a 0.. 0 ax BX and for (x-Ax,y,t), 2 2 f(X"AX,y,t) = f(x,y,t) _ AX lag—M + %{L Lf®21Y1t) 6x 3 3 4 4 (AX) a f(x,y,t) (AX) a_f(x.y.£2 3! 3 + 4! 4 + ... (3.2) 5x 3x Subtracting we have, 2 3 afgxzyfi) = f(X+AXQY2t) '7 f(x-Ax,y,t) - (AX) a f(x’y’t)‘ + ... ax 2Ax 6 M3 (3.3) and adding, 62f(x2}'1t) = fCX+AX:Y:t) ' 2f(X:Y:t) + fg’AX335t) 2 2 5x (AX) (1502 4f( t) - a. X’y’ + ... (3.4) 12 4 5x A rectangular grid with a mesh size Ax by Ay and At in time is imposed on the region of interest (see Figure 3.1). If the 13 l4 E ‘— H— H— I‘— ‘— I‘— 0) e. w— +. *_ I I T I I I i I I : I I Mdr—+-—-—F —————— , —————————— fi—-4 I I I I I I I I I I I j I I I ' _L A a r—f——*‘T ————— fl —————————— (“fl I I I l I I IN" : I I I 1 FM -7 -—--, ——————— r —————————— +—--+ I AAY' I I I o J I L:; L (h) I 11_ X 'Y 0 '1 MK 1 NX-l. NX MX = Number of intervals into which the upstream length, L is divided. NX = Number of intervals into which the entire x-domain, L is divided. NY = Number of intervals into which the half-plate separation distance a is divided. x = O = MX(Ax) = Location where heating begins. FIGURE 3.1: FINITE DIFFERENCE SPATIAL GRID FOR PARALLEL PLATE GEOMETRY 15 indices i,j, and n correspond to nodal points in x,y, and t respectively, then the temperature at a typical grid point in time and space is, T(xi,yj ,tn) = TEiAX’jAYSHAt] Using Space indices as subscripts and the time index as a super- script we have, n Ti,j - T(xi’yj’tn) Equations (3.3) and (3.4) then become in particular, r ------- 1 n n n 3 n - I I 5%,]- , Ti+1,j Ti-m . ML 1211' . - I (3-5) H ______ 5i-4 2 n n n n F""‘-"‘5'fi"7 - I ' 5 TL; ,___ 1+1,j ”m + Ti-1,1u_ $92 a Tm: (3 6) 2 2 : 12 4 I . ax (AX) . 5x ' .__ _________ J Similarly for the y derivative 2n n n n "'"’"'”ZYF“ . . . . - 2 . . + . . ' 2 2 2 2 I 12 4 I ' ay (Ay) . By J L...___._-_ ...... From Equation (3.1) the time derivative can be expressed as, —-----—--—- n n+1 n = n : EEEEJ 2 Ti:3 - Tia ' di- $A£l.gL ‘ (3.8) at At . 2 at2 : t... _________ J The boxed in terms in the above equations are retained only to give an order of the truncation error. See Appendix A. Adopting the following notation for the first and second order central difference operators we have, 16 n n n ~ n n = - ~ }’ - n c . i-(Tisj) Ti+%aj Ti'%:j 2(Ti+19j Tl'laJ) (3 9a) 2 n n n n = - + 3. sicri,j) Ti+l,j 2Ti,j Ti-l,j ( 9b) 2 n _ n n n Siam) " Ti,j+1 ”m + Tia-1 (3'9” 3.2 Alternating Direction Implicit Method The numerical solution of finite difference approximations to parabolic equations can be accomplished in a number of ways depending on the time step at which the difference operators in Equation (3.9) are evaluated. A number of these approximations are described by Douglas (42). The method used in this work is an alternating direction implicit (ADI) method similar to that of Douglas (42). This latter method is a modification of the method first prOposed by Peaceman and Rachford (43), and later modified by Douglas and Rachford (44). To the author's knowledge, this numerical method has only been used in solving diffusion-type equations. Pearson and Serovy (41) used the method in solving the Navier-Stokes equations but did not include the first order derivatives in the alternating direction scheme. Substituting Equation (3.9) into (2.8), (drOpping the super- script "+" for simplicity) and alternately evaluating the x and y derivatives at the unknown or future time level, we have for the intermediate time step n + %, n+% n T - T iii i,j .l_ n =.___l____ 2 n At/2 + Ax uj 61“) (”021,82 an) +--———1 52,(T“+!5) (3-10) am 2 J 17 and for the (n+1) time step, Tn+1 _ Tn+3§ i ' i ' 1 n+1 l 2 n+1 +—-U.6.(T )=T5.(T ) At/Z Ax J 1 (Ax) Fe 1 iv + 1 2 57:(T“+2) (3.11) (Ay) J Formulation in this way allows for the separation of unknowns into a y-direction calculation with NY unknowns and an x direction calculation with NX unknowns. Thus, the maximum number of simul- taneous equations to solve is Max(NX,NY) rather than (NX) X (NY) for a fully implicit method. Combining Equations (3.10) and (3.11) to eliminate the inter- mediate time step results in the following equation for the total time step, Tn+1 _ T i ' 1 1 n n+1 _ l 2 n n+1 At +2(Ax)uj6i(T +T )———-——2 251(1' +T ) 2(Ax) Fe + 1 2 agar“ + In”) +————§t 2 2 a: agcr“ - Tn+1) 2(Ay) 4(Ax) (Ay) Fe 3 2 + - —AtT u. 5. sicr“ - Tn 1) (3.12) 4(Ay) Ax J J The formulation as just described is that due to Peaceman and Rachford. If, however, the modified method of Douglas is used, the end reSult for a total time step is the same, but the particular x and y direction equations are much easier to solve. Rather than . . . n+% * . defining a half time step temperature T , let Ti j be an Inter- , mediate temperature defined by, 18 1* n 1,1 ' i,j 1 n 1 2 n + _ = _ At Ax uj 51(T ) (6x)2Pez 51(T ) -+-———l-§ agcr“ +-T*) (3.13) 2 (by) + and let T: ; be defined by, 9 n+1 n . . - T. . 1,3 1,1 + 1 u, 6,(Tn+1 + Tn) = l 62(Tn+l + Tn) At 2(Ax) J 1 2(AX)2Pe2 1 * + ——l-—§ 6?(T + Tn) (3.14) 2(Ay) Subtracting Equation (3.14) from (3.13) gives n n+1 At uj 61(T - T ) 4--————-§—-§ 2(Ax) Pe * n+1 At T. . = T. . - 1:] 1).] 2(AX) Tn+1 5:(T“ - ) (3.15) Substituting Equation (3.15) into (3.14) eliminating T* results in exactly the same equation as (3.12) for the Peaceman-Rachford formu- lation, so the two methods are equivalent for the total time step. The truncation error, which gives an indication of the error resulting from the replacement of the differential equation by its finite difference approximation, is dealt with in Appendix A. It is seen that the error associated with Equation (3.12) as given by Equation (A-8) is of order O[(Ax)2,(Ay)2,(At)2]. 2 2 +L ERROR = 1A%1_ [T 3 + u T 2 "‘l2'T 2 21n +'Lé%l— Tnaz + t txy Pe tx y Y (szz 1 n+% 6 E 2 T 4 ‘ u T 33 Fe X X A subscript notation has been used to denote derivatives, e.g. Tn3 = a3Tn/at3. t 19 Referring to Appendix B, one can see that a stability analysis predicts that the ADI method as formulated by Equation (3.12) is un- conditionally stable for any size Ax, Ay, and At. However, for a mesh that is too coarse, the solution of the finite difference equation, though stable, differs significantly from the solution of the differ- ential equation. Define, and Y At Then Equations (3.13) and (3.15) can be put in the form, X-Direction: + (a +Axu.6 -—l—5‘:‘)T‘i‘1— J Fe 3.] x l 2 n axTi,j + (Ax ujéi - 2 61)Ti,j (3-16) Pe Y-Direction: 2 * a - 6. T. . = ( Y J) 1:] 0’ 2 0’ 2 2 n (a - 2Ax—1u,5,+——-—16.+6.)T. . (3.17) y 01X J1 Pezozx 1 J 1.3 Written in matrix form, Equations (3.16) and (3.17) result in a tridiagonal system of the form, ._ ._ +1 ._ B +A F} =Q j = 0,1,...,NY (3.16a) x x J x. j J 3 +37”? =6y i = 1,2,...,NX (3.17a) Y y 1 1 20 Bx and E; are boundary condition vectors for the x and y directions reSpectively, ”6.1 o 0 Ex f 3 ET = 0 o th’L EN; 6; and a; are known vectors for the x and y direction equations respectively, t * qu F0 q1 q1 a=. a=' X. y. J qi 1 qj [508 for each j LfNZJ for each i where a typical element of 6; is j * n 2 n n = _ +—_ _ qi “xTi.j RJTi-1,j P82 Ti,j SjTi+1,j and a typical element of 6; is i n n n n n . = AX U. T. . ‘ T. . + E T. . + a - 2E T. . + E T. . q] B J( 1'13] 19]) 1'13] ( y ) l’J 1+laJ _ l Ax where R. -'——§ +'—— U. (3.18a) 2 J Pe _ l Ax S, _ m_§ ..E_ u, (3.18b) J Pe J a 3 =1 (3.18c) ax E _ 2 ““§ B (3.18d) 21 AK and A; are tridiagonal coefficient matrices of order NX-l and J 1 NY reSpectively,( ) r A -S 0 ‘1 X J -R A -S 0 J X J A = where A = a +-—g§ x. x x J Pe -R. A 'S O J X J O -R. A L J PA 0 0 O -l A -1 Y A = where A = a + 2 and y Y y -l A -1 c and d depend on O y 0 d A the boundary con- -n+l -* T and T are the unknown temperature vectors for the x and y directions respectively, n+11 1,J' Tn+1 2,J — + T”,1 1 = J = 0,1, ,NY J +1 T? . 1,] n+1 T NX:J (1) A; is of order NX-l since the boundary temperatures TO j j , j = 0,1,...,NY are known for all times and therefore do not have to be calculated. Also, Afl will be of order NY-l for x 2 O in Case II since Ti NY 1 = MX,...,NX is the specified surface temperature 9 equal to l for all time. 22 a”; T. L 1,113 3.3 Method of Solution Assume the temperatures over the first n time steps are known. To calculate the temperatures for the (n+1) time step, Equation (3.17a) is solved for each i = l,...,NX for the temperatures * T , i.e. a tridiagonal system of NY equations and unknowns is solved NX-l times. See Figure (3.2) m T G A G v J=NY 1 o O 090 o 00‘ 0 --—0J —e—-—_ g j=0 O 1 1 NX FIGURE 3.2: Y DIRECTION CALCULATION * Then with T known, Equation (3.16a) is solved for each j = 0,1,...,NY + for the temperature Tn 1. This requires the solution of a tri- diagonal system of NX-l equations and unknowns NY times. See Figure (3.3). 23 T—eg: A- e vHE—J—NY 8 o to .. o 01—; «Leg ,i; e ”v 9:2,J=o o 1 i Nx FIGURE 3.3: X DIRECTION CALCULATION The tridiagonal matrices are easily upper triangularized to solve for Tn+1 or T* by Gauss elimination. The solution time required to run the calculation to steady state with a typical mesh size of NX = 60 and NY = 10 on a CDC 3600 computer is from 30 seconds to 1 minute. Pe = 1 required NX = 108 with a run time of 1.5 minutes. It should be mentioned that at X+i= 0 where Ax must be very small for accurate calculation of the Nu number, the surface and bulk temperatures have been extrapolated to their respective T-intercept values as (Ax)2 * 0. Since the truncation error is of order (Ax)2, (Ay)2 and (At)2 (see Appendix A), for a given (by)2 and (At)2 the temperature is linear in (Ax)2 provided Ax is small enough. Two calculations for the same Ay and At and different Ax establishes the clope C1 of the straight line, T=C(Ax)2+C(A)2+C(At)2+T 1 2 y 3 0 C2 and 03 are nearly zero in the range of Ay and At used, so at X+ = O we can express TS and Tb as 24 2 Ts — Cls(Ax) +Tso and _ 2 Tb Clb(Ax) +Tb0 For the larger Pe numbers, Tso and Tbo are used in Equation (2.14). For smaller Fe and positions downstream from X+i= 0, C1 is negligible in the Ax range used. 3.4 Finite Difference Boundary Conditions It is more instructive to formulate the boundary conditions from a conservation of energy principle on a finite element centered on a boundary node than from the Taylors series approach used for an interior node. In this section, let the difference Operators at a boundary be denoted by equating the subscript to the appropriate boundary 2 2 2 6. , 62 are the second order index. For example, 6i=0’ 6i=NX’ J=0 j=NY difference operators at the respective boundaries. CASE I: Uniform Heat Flux The upstream boundary condition on X+ (Equation (2.9b)) is applied at X+'= -L where L is large enough that the fluid temperature is not effected by the heated plates, (see Figure 3.1). This distance must be larger for the smaller Pe numbers. 4. At X -L = ~MX(Ax): - O,1,2,...,NY 0,1,... :3 (.... I To,j=0 At y = O 25 6j=o(T ) = 0 2 _ 1 - 1,2, ,NX 6j=o(T ) - 2(T ,l - T1,o) - 0,1,. . . + . Cons1der the plate surface at y = l, J = NY. ll Cl1 II qz l ‘1 l c l l l .. Y —: q = NY d -——-—O‘ ——————>- . -‘3.°9-_9.-_-.-__g_ AV _..- fol‘d J = NY - 3:; conv -——u- r——-—->qconv qcond. o o 0 NY - 1 o o 0 j + l .1J<———£DK-—+4 O M O O J o O O J " 1 i-l i i+l FIGURE 3.4: ENERGY BALANCE ON A DIFFERENTIAL ELEMENT For an interior node, an energy balance applied to the element centered on the i,j node leads to the same finite difference equation derived before, Equation (3.12). 26 However, when considering a boundary node, the resulting "half" element or cell associated with this node should have fluxes centered on the half cell at a location Ay/4 from the plate surface (row of temperatures designated by 9 connected by a dotted line in Figure 3.4). The net rate of energy into the element is given as follows: 1 C d t' : ° - ' ‘ ( ) on uc ion q q i,NY-% cond cond i-%,J - qcond i+%,J k A1 Ti-1,J ' Ti,J _ (91 T1,.) ' Ti+1,J AZ (2 ) Ax 2 ) Ax Ti,NY ' Ti,NY-1:l - Ax Ay where J designates the position yJ = a - g2 . (2) Convection dconv.|i;%,J ’ é‘i+%,J = T + + Ax [1g Ti-1,J Ti,J Ti+1,J] 2 puJCpAz 2 2 - . . . = Ali? Rate of energy accumulation Within cell pCpAzAx 2 6t i,J Completing the energy balance, T. - T LT + 1+1,J i-l,J = pCp at i,J pCpuJ 2Ax T - T + - k i+1,J 2 i,J Ti-ILJ Ti,NY TilNY’l n l— + n .1— 3 1 2 '2 2 +q1Ay quy('9) (AX) (M) 3 1 = -— + _ Let TJ 4 TNY 4 TNY-l or in general = + TJ ClTNY CZTNY-l 27 where 3/4 0 ll 1/4 0 II The nodes along the boundary y = a fall into three categories depending on the heat flux at that point. + (1) For those boundary nodes which fall in the region X < 0, (1 s i S MX-l), q; = q3 = O. + (2) For those in the region X > O (MX < i s NX), _ n: u (3) For the boundary node located at the discontinuity in + heat flux at X = 0, q; = q" and q; = O. Expressing Equation (3.19) in the dimensionless ADI form we have, for the X-Direction l 2 n+1 * l 2 n [ax + Ax ujai - 2 611%,J - OtXTi,J + [Ax uJéi - 2 511T,“J (3.20) Pe Pe In terms of C1 and C2, 1 2 n+1 _ 1 2 n+1 CIEGX + Ax uNYéi - 2 6iJTi,NY CZEaX +'Ax uNY-lsi - 2 i]Ti,NY-1 Be Pe * l 2 n + - ....— + C1{°’xTi,NY [AX uNYéi Pe2 6i1Ti,NY} * l 2 n + Cz{°’xTi,NY-1 + ”X ”NY-151 ' P82 6iJTi,NY-'1} (3'21” Note the ease with which this boundary condition can be implemented using the ADI numerical method. The right hand side of (3.21a) is known since the TP+1 1,NY-1 1 = 1’ ° 0 - ,NX rOW Of temperatures are calculated before the NY row. 28 . _ 2 _ For 1 - NX, 51 - o and 51=Nx (T) becomes %(TNX,J _ TNX-1,J) in (3.21a). ll q" If 62 (T) = ZCT - T ) + Ay .El +-—g then, j=NY i,NY-1 i,NY q" q" using (3.18c) and (3.18d) we have for the Y-Direction * 2 x 2 n n - = - +-E + c T + ayTi,J 5j=NY Ti,NY Edy AXE 61 5i](C1Ti,NY 2 i,NY-l) 2 n 6j=NY Ti,NY (3.21b) 2 For i = NX, 6i O and 61=NX(T) = 35(TNX,J - T NX-1,J) 1“ (3.21b). For this case, the downstream x-boundary condition, 52T+(L-L,y+,t+)/ax+2 = O can be applied before the development length, xgév is reached without effecting the resultsin the region near 3 + X = 0. Therefore, when a smaller mesh size for greater accuracy + + at X = O is required, (LsL) may be less than XDev° The de- velopment length, defined as being that length from the entrance at which the local steady state Nu number is within 5% of its final value, is shown verses Pe number in Figure 3.5. This shows that X 9 O , is important the Pe number dependence, other than in X+I= 35; for Pe < 20. The fully developed value of xgév = 0.186 agrees with Cess and Schaffer (46). CASE 11: Uniform Tempreature The boundary conditions for this case are the same as for case n . I except for the constant wall temperature, Ti NY = 1 1 = MX,...,NX, , and the downstream x-boundary condition. 29 OJ.-’ 0’ i l, L_:Lv1 5 2| 1 32 7:0 20 40 Fe FIGURE 3.5 DEVELOPMENT LENGTH - CASE I 30 Applying the boundary condition (2.10c) at a location upstream of Xgév is more critical to the results in the entrance region than in Case I. However, with the following treatment of this boundary, the region about X+I= 0 appears to be relatively insensitive to the placement of the downstream condition. 2 Let 6i=NX = 6i=NX = O and Tn = Tn NXJ Nx‘laj This treatment of the downstream boundary on x eliminates the back-propagation of boundary anomalies into the entrance region resulting from not applying the condition far enough downstream. IV. DISCUSSION OF RESULTS FOR CASES I AND II In this chapter we shall investigate the effects of allowing axial conduction upstream of the heated region. Of particular + interest are the results at the entrance X = 0, since it is here thatconsiderable deviation from previous studies exists. 4.1 Case I If axial conduction is neglected, the steady state bulk temperature increases linearly with X+ from zero. As can be seen from Figure 4.1, axial conduction has the effect of increasing the bulk temperatures. Since the region X+I< 0 is insulated, the heat conducted upstream is convected downstream, raising the fully developed, dimensionless bulk temperature, T+, by 1/Pe2 as compared to the case with no axial conduction. The steady state, fully developed temperature profile, T+, and T+ are derived in Appendix b B from an energy balance on a control volume which includes the . + + . . entrance and upstream regions. For X >XDev the dimenSionless + steady state temperature T is, 2 4 + _ + l 3 + .l + 39 T‘x+ 2+4y '8y ‘280 (4‘2) Fe and the bulk temperature is, + + T = X +--l¥- (4.3) b 2 Fe 31 32 FIGURE 4.1 STEADY STATE BULK TEMPERATURE IN THE ENTRANCE REGION - CASE I 33 At X+ = 0 and the region immediately downstream, the bulk temperatures are slightly greater than their extrapolated asymptotic or fully developed values and decay exponentially for X+-< 0. Note that the develoPment length for T:' is much shorter than that used to define the entrance region, (i.e. xgév ). The temperature profiles at X+I= 0 for various Pe numbers are given in Figure 4.2. Notice how axial conduction from the heated region alters the temperature profile of the incoming fluid. This clearly shows the error in those papers which make a point of in- cluding axial conduction and then, through the boundary condition, force a uniform temperature profile at the entrance. This forces the temperature here to be the same as for the case where axial conduction effects are neglected. However, it is at X+ = 0, more then anywhere else in the entrance region, that these effects are greatest. Also, note that since T: and T+ at X+ = 0 are both b zero only as Fe a'm, the Nu as defined by Equation (2.14) x+-o is finite for finite Fe. The Nu number in the entrance region for Fe = l,2,6,10, 20,50, and m is given in Figure 4.3. For X+-< 0.02, the Nu number decreases with decreasing Pe, while the Opposite is true for X+i> 0.02. This can also be seen from Figure 4.4 which gives the Nu number as a function of Fe for X+'= O and several positions downstream of the entrance. The fully developed value of 8.225, as shown in Figure 4.3 is in agreement with Kays (45). It is apparent from these figures that axial conduction can have a considerable effect on the Nu number. The greater the Pe + number the smaller the region about X = 0 where the Nu number 34 LC) 0.8—- m -0 .0 o b/ 2 _- —~ Fe 006 — ...... I! pn- l 6.4 - / -— C>2.- __ .. 0,2 -— __ «’CL4-r— \\ .. 49.6 —- , .... ....{.O l I l I l - I . I 1 l l FIGURE 4.2 STEADY STATE TEMPERATURE PROFILES AT X=0 - CASE I 35 N H mmo 50 this deviation could be significant when dealing with a wall heat flux varying with X+. Since the energy equation is linear, superposition of the basic kernel for a step change in q" can be used to generate the reaponse to any variable heat flux q"(X+). In this way, the results in a neighborhood of, and including X+ = 0 would be distributed throughout the region of interest. If the axial cone duction is not included the resulting surface temperature and Nu number for such a case could be in substantial error even for large Pe number. For a further discussion on the superposition of solutions, see Chapter 6. The temperature profiles for Pe = 6 in the entrance region are shown in Figure 4.5. The profiles should converge to the fully developed profile envelope as given by Equation (4.2). As can be seen, when X+ = 0.45, the numerical results lie on this envelope. The transient heat conduction solution for an infinite slab of thickness 2a in the y- direction with a constant wall heat flux, q" is (54), T=_935 m6 maze: .8 mg v.0 MO dd HO 0 H6 I 38 1a.? i so- ..IL «.0! ll .1 N.O I 3.3 3638» m3 :1an 1 ... do 3 ... x o 8.0 w u .. ofx ‘\.. L 3 J 1.3.0 L 39 For flow in the entrance region downstream of the heat flux dis- continuity at X+I= 0, convection has no influence on the temperature response during early times since aT+75X+' is zero then and Equation (2.8) reduces to the heat conduction equation. Therefore, if one were to use the temperatures given by Equations (4.5) and (4.6) in the usual expressions for the bulk temperature (Equation (2.17)) and Nu number (Equation (2.14)), the slab transient heat conduction problem should predict the early transient behavior some distance (dependent on the time considered) downstream of X+'= O. This is indeed the case as can be seen from Figure 4.6, which gives the transient results before steady state is reached. Recall the expression for the bulk temperature, 2 + 1 ++ + 3 1 + + + - = -2- J‘Oa-y )T dy (4.7) - d Tb O u T y Substituting Equation (4.5) into (4.7) and integrating gives, 2 2 + + + - T=t -—1 “fit b 15 (4’8) [‘18 6 1 +7 Te U n n 1 Equations (2.14), (4.6), and (4.8) are combined to give an effective slab Nu number of, _ 2 Nu — .. 2 2 + (4.9) .1; .. 1— E (L + 2- L)e-n TT t 5 2 2 2 4 N n=1 n n n For steady state, Equation (4.9) becomes Nu = 10 (4.10) Notice that the velocity profile enters into the slab Nu number only as the weighting function in the bulk temperature H mm xDev.' We have seen that the steady state, fully developed Nu number is 8.225. However, if the distance downstream is sufficiently large (X+ > 0.3 for Fe = 6) the Nu number reaches a secondary steady state dwell at a value of 10 predicted by the transient heat conduction problem, Equation (4.10). Under these conditions, the length of time for the heat conduction response (Equation (4.9)) to reach steady state is shorter than the time required for a slug of fluid experiencing the entrance effects to be convected downstream to this location. Hence, the Nu number remains at this platau of 10 until the convection effects act to reduce its value to the final or primary steady state value of 8.225. Figure 4.6 shows the numerical results following exactly the curve of Equation (4.9) until convection effects are felt. For the region in the immediate downstream vicinity of X+D= 0, (e.g. the curves shown for X+-< 0.06) convection is significant even for very early times, resulting in transient Nu number curves which differ widely from Equation (4.9). 42 H mm H< onwmm muz935 wé mxszm 63 uwm . ...to L J l I l I l l I 47 L 8| “‘3: T+ 06 FIGURE 4.9 STEADY STATE TEMPERATURE PROFILES AT X=0 - CASE II 48 .1 fl .1 LO-l‘lnnrl. .|.14 0.2. 0.4 T+ 0.6 0.8 1.0 FIGURE 4.10 STEADY STATE TEMPERATURE PROFILES IN THE ENTRANCE REGION FOR Pe = 6 49 I -CL2.k ~CIOOQ '0MOOH I : 0 -1'0 7- I J n L L L l L l j l 0.2. 0.4 TI” 06 0.8 FIGURE 4.1] STEADY STATE TEMPERATURE PROFILES IN THE ENTRANCE REGION FOR Pa = 50 l l l i 5.. 50 HH mm355 NHé mag: ZOZUDOZOU .272 oz .. 6$ bad ..x «.0 so 5 was was So moo NE 5d 0 i_ 4: fl _ ..H._.H..i.__._._. r J [j J / 8 l I i o. i 1 «use ...U IJ J Baum i o ' --” b-i' 1'1“?" " I; 3 dz 51 HH mm = z [Axon (t) - 9(t)].dt (5.1) . o i s 1 i=1 is to be minimized with respect to Fe. The temperature Ts(t)|i a function of Pe, is the calculated numerical solution to Equation . . + + . . . (2.8) at pOSition Xi’ y = l, and time t, for a spec1fied heat flux, q"(x,t), (derived from the Case I solution by superposition). The corresponding experimental temperature is given by ei(t). The quantity Ai’ subsequently taken to be unity, is a weighting factor included here for generality. If the statistics of the errors in the temperature measurements are known, Ai is Often taken to be the inverse matrix of the variance-covariance matrix of the errors. A number of procedures have been suggested for minimizing F, some of which are the method of steepest descent and modifications to the Taylor series approach (47, 48). Under Optimum conditions 54 the Taylor series approach yields very rapid convergence (see the example in Section 6.4) and will be used in this work. Though the calculated temperature is a nonlinear function of the Pe number, the Taylor series method is an iterative procedure which assumes that at each step the temperature is a linear function of Pe. If PeO is an initial estimate Of the Pe number, then M (rec) S TS(PeO) + APe aPe (5-2) Ill TS(Pe) where APe = Pe - PeO is the first correction in the Pe number. x . . * . . . . 5F . Noting that at the pOint Pe at which F is a minimum, aPe IS zero, Equations (5.1) and (5.2) can be solved to give the first correction in Pe. n aT (Pe ) aT .(Pe ) 117—: T ____9_ _ ng 0 aPe 2 1:1 IOETS(PeO) + APe aPe e1i aPe dt * At the minimum value of F (i.e. F ) we have, n aTS i(Peso) d =1; I08” (Pee) 931 are t n 5T 82(Peo) + g fot——§P———]: APe dt i-l Solving for APe, TS i(Peo ) iii‘fo [e - TS (Fe 0)]: aPe dt APe - n 6T (Pee ) (5.3) i21IT\-2____o_]23 dt aPe i For one thermocouple making NT discrete measurements in time T, Equation (5.3) becomes, 55 NT aT (Pe) s,k 0 £21[9 - T3(Peo)]k aPe APe = NT BT(PeO) 2 (5.4) z -——---1 k=l aPe k T = NT(At) An improved value of Fe is given by Pe Pe + APe l o and the iteration procedure is continued until some convergence criterion is met, say £9< 0.0001 Pe 5.2 Objective Function for Determining Optimum Experiment In order to determine an optimum experiment it is necessary to find some basic criterion which when extremized yields the Optimum experiment. The criterion is frequently called an objective function. An experiment can be an optimum one in various aSpects and thus there may be different optimums depending upon the desired Objectives and constraints. For this problem, the best heat flux boundary condition is determined to insure that, for a random distribution of errors in the wall temperature measurement, the resulting error in the mean velocity is a minimum. In maximizing this Objective function the following conditions are satisfied. 1. The maximum surface temperature rise is to be ATmax This can be chosen small enough to make the assumption of k, p, Cp, and u 2. 56 independnet of temperature realistic. The number of thermocouples, n, along the plate is fixed. Specifically, optimum conditions for one thermocouple are found. 3. NT discrete measurements in Assume that the value F(Pe) * sum of squares function, F , There is a fixed duration of the experiment, T. x is known and is designated Pe . For time, T = NT(At). of the Pe number at the minimum of The minimum value of the need not be zero but it is assumed that the errors in the temperature measurements are small. Let us examine F(Pe) * x in the local region near F = F(Pe ). The sum of squares function is, n _ T 2 F(Pe) - 33. Tons“) - e(t)]idt (5.1) 1—1 * d * 2 d2 * F(Pe + APe) ’5 APe d—g-g + (A381 F2 (5.4) dPe * aT* n dF __ = T * _ __s dPe 0 2 E To“, 9&5... d (5.5) i-l 2 * 52 * n T d F * 2=2 2 flag-eh Zdt dPe i=1 aPe n aT T 2 + 2 z fo(§%)i dt (5.6) i=1 If the error in the temperature measurements is small, then the first term on the right hand compared to the second term. into (5.4) we have, side of Equation (5.6) is negligible Substituting Equations (5.5) and (5.6) 57 * T * * 2 n T a 2 F(Pe + APe) = F + (APe) 2 IO(;§§91 i=1 aT* * APe 2 n T * s 2 =F+(——;) z (Pa—.dt Pe i=1 Io aPe 1 x n T * aTS Let ¢- E fowe find: (5.7) 1-1 * * Then F(Pe + APe) - F = (ABE->295 (5.8) Pa F I I ,/.Large ¢ I I I I I I I I I i I I Small ¢ l I * I F __________ ._ I I 1* rib Pe Pe FIGURE 5.1: ILLUSTRATION OF SUM OF SQUARES FUNCTION Referring to Figure 5.1, we see that for a well defined minimum in F(Pe), ¢ should be a maximum. Considerable information can be gained from the sensitivity . . . 5T . . . . coefficient defined as Pe aPe (51). It is the prinCiple quantity of interest in the objective function since it gives in a relative way a measure of how sensitive the temperature is to a change in Pe number. Obviously, if we are to calculate the Pe number from temperature measurements we would like this sensitivity coefficient as large as possible. Also, it can be eSpecially instructive when 58 considering two or more parameters since it establishes to which parameter the temperature is most sensitive. The objective function, ¢, must be modified to include the constraints of maximum temperature rise, AI = T - T , and max max 0 given duration of experiment, T. This can be accomplished by dividing . . . 2 . the right hand Side of Equation (5.7) by T ATmax' The resulting Objective function T, given by Equation (5.9), is a dimensionless criterion which is to be maximized with respect to the boundary heat flux, q"(X+,t+), and thermocouple location. n 5T 1 2 0 = -——2 .2: [Ewe ~34), dt (5.9) 1' AT i=1 5P9 1 max For one thermocouple and NT discrete measurements over time T, Equation (5.9) is, NT 3T 0 =——-1-2-— z (Pe—P-E-lidt (5.10) T AT k=1 5 max ++ Before investigating the optimum heat flux, q"(X ,t ), which ++ makes é a maximum, consider a simple but useful q"(X ,t ) profile. Suppose that at any instant the heat flux is uniform over a heating length x and that the plates are insulated elsewhere. Also, over L a given duration of the experiment, let the heat flux have a constant value such that the maximum temperature rise is less than or equal to ATmax' The following question is then answered: for a given experiment duration, T, what is the length of heating, xL, and placement of one thermocouple, x , which will make T a maximum 0 for this constant q"? Once xL is determined, then the magnitude of q" which satisfies the AT constraint for a given T can max be determined. This can be done for each experiment duration, 59 giving a full time spectrum Of optimum transient experiments for a constant q". These results are given in Table 1. Both xL and x increase continuously as T increases, but due to the nature of a finite difference method, the smallest increment detected is Ax. (For Table l, AX+ = 0.015 and Fe = 6, hence Ax/a = Pe AX+ = 0.09) This explains why xL and x9 in the table appear to be constant for some time and increase in increments of 0.09. For example, a particular entry in Table l is illustrated in more detail in Figure 5.2. This figure shows that if the duration of the experiment is T+ = 0.12, the Optimum heating length is xL/a = 0.54 and the Optimum placement of the thermocouple is xe/a = 0.09. 252 --04 -02 °’ 00 M: 0.2 I I I I I I II T ' Pe:6 : - 0.002— I*=O.54 5 _. i> 02.1w 5 ‘ 0.001 F- I .... IENGTH: XLIe' X961-op'TIM0M - I PLAcEMI-Zm OF one - __ I THERMOCOUPLE WHEN : XL/xL by superposition. The surface temperature does not reach a maximum at the end of the heating length as would be the case if axial conduction upstream‘were not allowed. As can be seen from Figure 5.3, AT = Tmax - To’ occurs some distance upstream of + XL. The smaller the Pe number, the greater this distance. Notice that the position of maximum temperature shifts downstream‘with increasing time. Also, recall from Section 4.1 that the further one goes downstream, the longer is the time required to reach a steady state condition. As a consequence of these two facts, the optimum heating length becomes larger as the duration of the experiment increases, (see Table 5.1). In terms Of the dimensionless variables of Equation (2.7) for Case I, Equation (5.10) becomes, + NT 3T8 l "a 2 @= + 2 (9k) 21(Pe 7:2) At+ AT k=l 5e T max NT=2,3,... '* Ihmm.- To Since ATmax - q"a k , we have + NT aTs M+ =-———:2—'k Z 1(Pe --): (5.11) NT = 2,3,... 62 TT :52“: wzH._.3rw V0.0u+P «.0 mud 3:35:23; 3.0 63 BTS The sensitivity coefficient, Fe 33;, is evaluated by solving Case I for two different Pe numbers differing by a small value, such as APe = 0.005. Then, 5T n Tn .(Pe + APe) - Tn .(Pe) Pe ——S '3 Pe 5’1 $41 ape i APe i = 0,1,...,NX The coefficient with reSpect to u is actually the quantity of interest. However, ; M = Pe $08:an au x aPe Ia,a,x . . . + Ax . . . 1f the dimen31on1ess Ax = 353' in Case I is changed pr0port10nate1y when the Pe number is changed so that T(Pe +-APe) and T(Pe) are evaluated at the same physical location. The steady state sensitivity coefficient for a step change in surface temperature is given in Figure 5.4. The entrance results converge to the fully deve10ped value of -(X+'+--£§) derived in Appendix C. Fe The transient results for the same boundary condition are given in Figure 5.5. Note that for early times a maximum occurs near X+'= 0, and shifts downstream with increasing time. This causes the optimum placement of the thermocouple to move from the vicinity of X+u= 0 for early times, downstream as the duration of the experiment increases. The sensitivity coefficient for an arbitrary variation in q"(X+,t+) is obtained by superposition of the step change case in + Peal”— is linear in aPe the same way as with the temperature, since T T T I I U I T I ‘ I I ' j— I T j 4 CH3 - - 4 ' Pe =6 '7 0.64— a L -Pe 55: J are P I- q o4 - 4 . , J ’l T“ ’1 * 2 .— , FULLY DEVELOPED: x + 5;, F x, J (12»— q ~ + y... .. o l 1 l l L J l I i 1 E l l j 1 l -04 o 0.7. 0.4 0.6 08 FIGURE 5.4 STEADY STATE SENSITIVITY COEFFICIENT FOR CASE I BOUNDARY CONDITION 65 'l'er'Ijl'T'l'TjTT Y' FTB=€5 0.19.— +. . tf-OAO CADE— + ..k T P o 0 0.081- 0.06:»— o.” 0.0!“— r C11C) 0.02}— - 0.05 L o . ... l A _L 1 L 1 ._ l I 1 l L -0.1 O 0.1 x+ 0.2. 0.4 FIGURE 5.5 TRANSIENT SENSITIVITY COEFFICIENT FOR CASE I BOUNDARY CONDITION 66 X+ and t+ for a given Pe number. Let us investigate what happens to the objective function as the duration of the experiment increases. Since Q is a time average of the Square of the sensitivity coefficient, it takes much longer to approach its final value of '12 F 313' Pe-——§ aPe x+ 9 AT+ L. max _J than for the temperature response to reach its steady state value. + In fact, §'~ 1 as T a m. To gain some insight into what happens + . . . as T d'm, suppose there 18 no axial conduction upstream; the + . . . , and p081t10n of max1mum 9 + temperature occur at the end of the heating length, XL. The above optimum location of the thermocouple, X expression of the summand in Equation (5.11) can be expressed in aT terms of the fully deve10ped values of Fe 332' and T:' given in + Appendix C, sinse XL is large for large T+1 Namely, + BTS + 2 Pa aPe X+ XL + """2' e -—> _ Fe + + l 17 ATmax xL + Pez + 35 The optimum heating length maximizing this quantity as T+ increases is x::~’w. Allowing axial conduction upstream will not alter these limiting results. Obviously, an infinite heating length is not tolerable, physically. Since no optimum time duration exists for this case as formulated, the duration of the experiment is left to the choice of the investigator. Once it is chosen, however, xL, xe, and 67 q" a constant which make é a maximum can be found using Table 1. We would like T to be small in order that 3' may be estimated even though it is time varying. Also, from a physical consideration, it would be desirable to have the heating length short, and since the optimum x increases with the duration of the experiment this L means that T should be made as small as possible. Incidently, the nonlinear estimation procedure requires repeated transient calculations of the Case I solution over the duration of the experiment. Hence, a smaller choice of T requires less computer time. We now wish to extend the optimum conditions to include a transient heat flux profile. At a particular instant, let the heat flux be constant over the heating length, XI. However, in addition to maximizing Q with resPect to the heating length and placement of one thermocouple, allow q" to vary with time in such a way as to make Q a maximum for a given experiment duration while satisfying the constraint. AT max 8T Consider Equation (5.10). The quantity, Fe 33%, for a given heat flux variation in time can be calculated, using as a basic building block in superposition, the sensitivity coefficient result- ing from a unit pulse in the heat flux from t+ = 0 to t+ = At. let Wk be this basic building block. Then if Pe g;% is the sensitivity coefficient due to a unit step change in q" at t+ = 0 (T: = T+YX+,l,t+) is the solution to Case I at y+ = l), aT+ aT+ wk = (Fe 35% k - (Pe Sifik-l + represents the sensitivity coefficient evaluated at t = k At 68 + due to a unit pulse in q" for the duration At beginning at t+r= 0. The coefficient resulting from a similar pulse in heat flux of magnitude q" is then q; Wk. The superposition principle gives the following expression for Fe SEE at time k At+ due to a heat flux varying with time as illustrated in Figure 5.6. 5T (Fe -—8) + q" W = " +...+ " 5.13 aPek ql Wk 2 k-l q W ( ) k 1 H n q']: l H V + k=0 1 2 3 4 t Akfl FIGURE 5.6: ILLUSTRATION 0F DISCRETE OPTIMUM HEAT FLUX PROFILE aTS 2 k k —— = II II (Fe aPe k E {3 clj qp wk-j+1 wk-p+1 (5'14) —1 p-l Hence, 1 NT k k Q = _""— Z 2 2 q". q" W _. W _ (5.15) ‘1' AT2 k=1 j=1 p=l J P 1‘ 3+1 1‘ PH max If W has the same sign for all k, the objective function, k Q, will be a maximum*when each term in the sum over k is a 69 maximum, (i.e. when (Pe aTS/aPe): in Equation (5.14) is a maximum for any k). The coefficient Wk is always negative since aT:75Pe monotonically decreases in time from its zero initial value. The optimum transient heat flux is one which varies with time in such a way that the surface is maintained throughout the experiment at its maximum permissible value. This insures that the heat flux is as large as possible at every instant. If the q" illustrated conceptually in Figure 5.6 is this heat flux then its value at any time, q", can not be increased or the ATmax constraint will be violated. If qfi for some k were smaller than the maximum permissible, then noting that W is independent of q", the terms in the sum in Equation (5.14) would not be maximized. For example, the optimum transient heat flux for a heating + length of x = 1.02 a, (XL = 0.17; Fe = 6) is given in Figure 5.7 L per degree The position of the maximum temperature varies ATmax' as before, moving downstream with increasing time until a steady state position of x = 0.72 a is reached. (Keep in mind that the heat flux at any given instant is constant over the region: + + Osx sx). L + The heating length of XL = 0.17 is an optimum for a choice + in experiment duration of T = 0.8. From Figure 5.8, with 6 given by Equation (5.15), we see that the optimum location of one thermocouple for this example is x;'= 0.06 or since Pe = 6, = 0.36 . X6 8 The choice of the experiment duration is somewhat arbitrary as mentioned before. It can be seen from Figure 5.9, which gives the transient behavior of the objective function, that once the 70 to $9 x<2h< mo mmua mascamMQZme ma <1ng (6.1) T 0567 + for —m < X < w + t > 0 The terms in brackets in the above equation account for the non-zero initial value in q" as well as its non-zero value at X+ = 0 and the end of the heating length, x:. If no discontinuities in q" exist, (i.e. if q"(X+,0) = q"(0,t+) = q"(X-£,t+) = 0) then Equation (6.1) reduces to XL 2 T<+—§,y+,t+-T>d§d~r (6.4) 81 Equations (6.1) and (6.4) differ from the result for no axial con- duction in that the limit of integration on g is over the entire length of heating rather than from 0 to X+. This is due to the fact that results upstream are effected by downstream conditions. To the author's knowledge, Equation (6.4) does not appear anywhere in the open literature. As indicated in chapter 4, the temperature resulting from a space variable heat flux could be in significant error unless axial conduction is taken into account. The results at X+ = 0, where the greatest error exists in neglecting axial conduction, are distributed throughout the X+-domain. Equation (6.4) becomes an inverse problem if the surface temperature is a known function of position and time and the heat flux is to be determined. The resulting integral equation, with q” in the integrand unknown, is of the first kind; a Fredholm type on X+ and a Volterra type on t+ (60). Parenthetically, it is noted that with this formulation, the constant heat flux results - (Case I) can be used to determine the response to a unit step + change in surface temperature for X 2 0 (Case II), a/k 1 = TS-TO X + 0L J“: q"(i,T)K1(X+-§,1,t+-T)d§d'r (6.5) There are a number of ways to solve integral equations of this type, some of which are mentioned in references (57, 60). To obtain a solution numerically, suppose the heat flux in Equation (6.5) has been determined for the first (n-l) time steps. At the nth time step an algebraic difference equation can be written for each node in the x direction over the heated region (say m 82 nodes), then summed over the first n time steps. Each such equa— tion in turn contains a sum over the m nodes in the heated region - a 'l + + 1eav1ng (m+1) equations and (m+l) unknowns. Once q (X ,t ) is determined from Equation (6.5), it can be used in Equation (6.4) to calculate the temperature response throughout the region between the plates. No reported work has been done on this inverse problem. When using experimental data, severe problems with stability might arise (59). 6.2 Conducting Parallel Plates Case I assumes that the heat flux applied to the outer sur- face is conducted through the plate to the fluid without axial conduction losses. It is of considerable practical interest to investigate the case where the plates themselves conduct heat in the axial direction. Let ks be the thermal conductivity of the plates and let q; be the heat applied to the outer plate Surface over a length XE. WI 1 1 l" 1‘ |-——-x———~l k—dx T it: . L —k5LT_.. k s; _, sax x T M M FIGURE 6.2: DIFFERENTIAL ELEMENT IN A CONDUCTING PLATE The following assumptions are made, 1. Steady state 2. k8 = constant 3. The Biot number, h6/ks, is small, say less than 0.1. Assumption (1) is made in order to direct attention to the conduct- ing plates. The results of section 6.1 could be used to consider transient effects. An energy balance on the element shown in Figure 6.2 gives, 2 T qux - q"dx + ksé :;3 dx = 0 84 By assumption (3), q; + kss a 23 = q" (6.6) X where TS is the fluid-plate interface temperature. Note that if the transient case is considered the heat capacity of the plate must be included. This adds an additional term to the left hand side of Equation (6.6), namely, -(pCp)S aTs/at. + In dimensionless X+ and y form, Equation (6.6) is, 2 k 6 a T q: as. = m M a Pe 5X T =T(x+.1) Note that Ts has the dimension of temperature. Equation (6.4) for steady state conditions reduces to, T0<+,y+)-To = if: q"<§>K2d§ (6.8) on where now the kernel is KZCX+:§,y+) = aT+(X+-§,y+)/ax+ with T+(X+,y+) being the known steady state solution to Case I. The function q"(X+) is the variable heat flux applied to the fluid (not the plate). Note that the limits of integration on g must be over the entire X+ domain from «m to m, since q"(X+) has a non-zero contribution over this range. Introducing Equation (6.7) into (6.8) gives, X+ + + _ a L n +_ + m .y >—To - RIO qo K2d§ k 6 2 s 1 a T(§,1) +_ + + k a —Pe2 j: 5:2 K2d§ Let q: = constant and note that 85 aT+(x+-§,y+2 = , aT+£X+-§,y+) + as = K2 (X+'§ 237+) - ax Then we have the following integral equation of the second kind. 'I _qoa + + + T(X+,y+)'To - _k— [T (X ,y ) ~ T+(X+-X:,y+)] ks6 l QZngzl) + + “LEFI: F32 K2(X -§.y )dé (6.9) e a The temperature T+(X+,y+), its derivative K£X+,y+), To’ and the heat flux, qg, are knOWn while T(X+,y+), appearing both on the left hand side and in the integrand (for y+ = l), is to be determined. One possible method for solving this equation would be to write the finite difference equation for each node at the wall (y+ = 1) over the enitre x domain (say M nodes). The equation for a given node contains all M unknown wall temperatures. One from the left hand side of Equation (6.9) for that node plus the M temperatures resulting from the summation of the second derivative. Hence, there are M equations and M unknowns. Once the wall temperatures are known, Equation (6.9) reduces to an ordinary integral with the integrand known and the fluid temperature throughout the region between the plates can easily be calculated. The problem of conducting plates has not been considered in this manner previously to the author's knowledge. 6.3 Combined Convection and Radiation Between Parallel Plates. Radiation has become of general concern in many applications recently. This section considers the situation where the plates 86 are of sufficiently high temperature that radiation as well as con- vection is an important mode of heat transfer. If the fluid is a non-interacting medium, then radiation between the plates is a boundary phenomenon dealing with the entire x-domain from -m to +m. Hence, the solution to Case I serves as a natural basic building block in formulating the integral relation for this situation. The following assumptions are made: (1) 'The plates are gray (2) The fluid between the plates is non—absorbing and non-emitting (3) Steady State (4) Heat flux, qg, is applied for X+n> 0 From Love (57), the total heat flux both emitted and re- flected leaving the surface is given by, R = so T:(X+) + p [Emmmimdn (6.10) + where R(X ) = the radiosity o = Stefan-Boltzmann constant e = emissivity of the gray surface p = reflectivity T = absolute temperature in oR TS = steady state plate temperature = T(X+,1) n = dummy X+I variable The kernel F, a function of the geometry, is related to the shape factor. It is the fraction of radiation leaving an elemental strip at n on one plate which strikes an elemental strip at X+ on the other plate per unit length (in n direction) of emitting surface. See Figure 6.3. 87 Elemental strip ///////////1 Latl l l/l l l .k—dx* |—— x* d1\——-| ”WW? 1% I 1 T r 1 1 FIGURE 6.3: ILLUSTRATION OF RADIATION BETWEEN PARALLEL PLATES + As given in Ref. (57), the kernel F(X ,n) for this geometry is, 2 2 + _ 1 (2a 2 2a F(X an) — __E 1 = 2 r3 C(X+-fi)2 +4a233/2 Hence, from Equation (6.10), the radiosity for either plate is, + 4 + 2 R(md‘fl R(X ) = 60 T (x ) + Zoe [1 (6.11) s [(X+-fi)2 + 43213/2 The total heat flux absorbed by a differential strip dX+ + at X by radiation from the other plate is, + 2 m R(fl)dn + (1"(X ) = 203- " R(X) (6-12) 1‘: I-oo [$414.02 + 4a213/2 where a = absorptivity, (a = e for a gray surfact). Finally, the total heat flux added to the fluid is, 88 q"(X+) = q32 + 4a2] with R(X+) in turn given by the integral Equation (6.11). Note the limits of integration on Equation (6.14) must run from «n to 00. Equations (6.11), (6.14), and (6.15) can be solved for R, q”, and T by an iteration procedure. At y+ = 1 these equations are of the inverse type. However, once q" and R have been determined the temperatures throughout the region between the plates . can be calculated by a straightforward integration. Let the x domain be divided into M nodes. Then for y+ = l we have 3M equations and 3M unknowns. If Equation (6.14) is written in finite difference form for each of the M nodes, there results M equations each containing M unknown q"'s and one TS = T(Xj,l). Hence, there are 2M unknowns and M equations. Similarly, Equation (6.11) yields M equations each containing M unknown R's and one T:, or altogether 2M unknowns with M equations. Finally, Equation (6.15) gives M equations each containing M unknown R's and one q" for a total of 2M unknowns. The three equations simultaneously supply the 3M equations needed to solve for 3M unknowns in q", R, 89 and TS. One poSsible iteration procedure would be to pick TSCX+) + + and solve for R(X ) in Equation (6.11). Then q"(X ) can be calculated from Equation (6.15) which in turn allows new surface temperatures to be calculated using Equation (6.14). VII. CONCLUSIONS In summarizing, the following conclusions can be drawn: 1. Axial conduction has a substantial effect on the temperature distribution and resulting Nusselt number in the entrance region for Fe < 50 in both boundary condition cases. 2. In the neighborhood of, and including the location where heat- ing begins, the results are influenced significantly by axial conduction even for Fe numbers much greater than 50. The larger the Pe number, the smaller this region of influence becomes. 3. If axial conduction is important, then the temperature profile in the region upstream of the heated section is significantly altered. Hence, one can not assume a uniform temperature profile at the entrance cross-section in this case. 4. Axial conduction has the effect of increasing the development length. As the Pe number decreases, the development length increases. 5. The early transient response following the application of the wall heat flux can be described by the heat conduction equation if the location of interest is sufficiently far downstream from the entrance cross-section. This location may still be well within the entrance region, however. 6. At axial positions located in the fully developed thermal region, the transient Nusselt number decreases in time until reaching a 90 91 "secondary" dwell at a constant value predicted by considering only conduction between the plates. It remains at this intermediate level until the entrance effects are convected downstream to this location. Thereafter, the Nusselt number decreases to its final steady state, fully developed value. 7. An objective function is defined which, when maximized and satisfying given constraints, establishes an optimum experiment to estimate the mean velocity from temperature measurements at the wall. 8. The Optimum transient wall heat flux profile is one which de- creases with time in such a way as to maintain the maximum wall temperature at its constrained value throughout the duration of the experiment. At any instant the heat flux is constant over a finite heating length in the axial direction. 9. The Optimum heating length for a given experiment duration and Fe number is determined. Also, the optimum location of one thermocouple at the wall for this case is found. 10. No optimum experiment duration exists when the objective function is expressed in the given form. This design factor is left to the choice of the investigator. However, there are several important considerations effecting this choice. 11. The optimum heating length and location of one thermocouple are also presented for a simpler wall heat flux which is not only constant over the heating length, but also constant throughout the duration of the experiment. The magnitude of this constant heat flux decreases with increasing experiment duration so as not to exceed the maximum wall temperature constraint. 92 12. A number of inverse problems are formulated using as a basic building block in the superposition principle, the uniform heat flux case. 10. ll. BIBLIOGRAPHY Graetz, L., "fiber die‘warmeleitfahigkeit von Flfissigkeiten", Ann. Phys., 18, 79, 1883, 25, 337, 1885. Prins, J.A., Mulder, J., and Schenk, J., "Heat Transfer in Laminar Flow Between Parallel Plates", Appl. Sci. Res., A2, 341-348, 1951. Schenk, J. and Beckers, H.L., "Heat Transfer in Laminar Flow Between Parallel Plates", Appl. Sci. Res., A4, 405, 1954. Mercer, A. McD., "The Growth of the Thermal B.L. in Laminar Flow Between Parallel Flat Plates", Appl. Sci. Res., Hague A8, No. 5, 357, 1959. Siegal, R. and Sparrow, E.M., "Transient Heat Transfer for Laminar Forced Convection in the Thermal Entrance Region of Flat Ducts", Trans. ASME, Journal of Heat Transfer, Vol. 81, 29-36, 1959. Siegel, R., "Heat Transfer for Laminar Flow in Ducts with Arbitrary Time Variations in.Wall Temperature", Journal of Applied Mechanics, E 27, No. 2, 241, 1960. Kays, W.M., "Numerical Solutions for Laminar-Flow Heat Transfer in Circular Tubes", Trans. ASME, Vol. 77, 1265-1274, 1955. Siegel, R. and Sparrow, E.M., ”Simultaneous Development Of Velocity and Temperature Distributions in a Flat Duct with uniformTWall Heating", AIChE Journal, Vol. 5, No. 1, 73-75, 1959. Ulrichson, D.L. and Schmitz, R.A., "Laminar-Flow Heat Transfer in the Entrance Region of Circular Tubes", International Journal of Heat and Mass Transfer, Vol. 8, 253, 1965. Tien, Chi and Pawelek, R.A., "Laminar Flow Heat Transfer in the Entrance Region of Circular Tubes", Appl. Sci. Res., A13, 317. Hornbeck, R.W., "An All-Numerical Method for Heat Transfer in the Inlet of A Tube", ASME Publication, 65-WA/HT-36, 1965. 93 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 94 Mercer, W.E., Pearce, W.M., and Hitchcock, J.E., "Laminar Forced-Convection in the Entrance Region Between Parallel Flat Plates”, Trans. ASME, Journal of Heat Transfer, August 1967. Han, L.S., "Simultaneous Developments of Temperature and Velocity Profiles in Flat Ducts", Interantional Developments in Heat Transfer, International Heat Transfer Conference, Boulder, Colorado, 591-597, 1961. Sparrow, E.M., Lin, S.H., and Lundgren, T.S., "Flow Develop- ment in the Hydrodynamic Entrance Region of Tubes and Ducts", Phys. Fluids, Vol. 7, No. 3, 338-347, March 1964. Sparrow, E.M., "Analysis of Laminar Forced-Convection Heat Transfer in Entrance Region of Flat Rectangular Ducts", NACA, TN3331, 1955. Diessler, R.G., "Analytical Investigation of Fully Developed Laminar Flow in Tubes with Heat Transfer with Fluid Properties Variable Along the Radius", NACA TN2410, 1951. Deissler, R.G. and Presler, A.F., "Analysis of Developing Laminar Flow and Heat Transfer in A Tube for A Gas with Variable Properties", NASA TM X-52171. van der Does de Bye, J.A.W..and Schenk, J., "Heat Transfer in Laminar Flow Between Parallel Plates", Appl. Sci. Res., A3, 308-316, 1953. Schenk, J., Dumoré, J.M., "Heat Transfer in Laminar Flow Through Cylindrical Tubes", Appl. Sci. Res., A4, 39-51, 1954. Lundberg, R.E., Reynolds, W.C., and Kays, W.M., NASA TN D-1972, Washington, D.C., 1963. Tyagi, V.P., "Laminar Forced Convection of A Dissipative Fluid in A Channel", Journal of Heat Transfer, 88, No. 2, 161, 1966. Wolf, H., "Heating and Cooling Air and CO in the Thermal Entrance Region of A Circular Duct with Large Gas to Wall Temperature Differences", Journal of Heat Transfer 81, No. 4, 267, 1959. POppendick, H.F., "Liquid-Metal Heat Transfer", Heat Transfer Symposium, University of Michigan Press, 1953. Sellers, J.R., Tribus, M., and Klein, J.S., "Heat Transfer to Laminar Flow in a Round Tube or Flat Conduit - The Graetz Problem Extended", Trans. ASME, Vol. 78, 441-448, 1956. Schneider, P.J., "Effect of Axial Fluid Conduction on Heat Transfer in the Entrance Regions of Parallel Plates and Tubes", Trans. ASME, Vol. 79, 765, 1957. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 95 Hwang, C.L. and Fan, L.T., "Finite Difference Analysis of Forced Convection Heat Transfer in Entrance Region of A Flat Rectangular Duct", Appl. Sci. Res., A13, 401-422, 1964. MOMordie, R.K. and Emery, A.F., "A.Numerical Solution for Laminar-Flow Heat Transfer in Circular Tubes with Axial Con- duction and Developing Thermal and Velocity Fields", ASME Paper No. 66-WA/HT-7. Hsu, C.J., "An Exact Mathematical Solution for Entrance-Region Laminar Heat Transfer with Axial Conduction”, Appl. Sci. Res., A17, NO. 4-5, 1967. Singh, S.N., "The Determination of Eigenfunctions of a Certain Sturm-Liouville Equation and its Application to Problems of Heat Transfer", Appl. Sci. Res., A7, 237-250, 1957. Singh, S.N., "Heat Transfer by Laminar Flow in a Cylindrical Tube", Appl. Sci. Res., A7, 325-340, 1957. Eraslan, A.H. and Eraslan, H.F., "Heat Transfer in Magneto- hydrodynamic Channel Flow", The Physics of Fluids, Vol. 12, No. 1, 120-128, 1969. Harrison, W.B., "Forced Convection Heat Transfer in Thermal Entrance Regions: Part III, Heat Transfer to Liquid Metals", ORNL-915, 1954. ' Schneider, P.J., "Effect of Axial Fluid Conduction on Heat Transfer in the Entrance Regions of Parallel Plates and Tubes", Trans. ASME 79, 766-773, 1957. Burchill, W.B., Jones, B.G., and Stein, R.P., "Influence of Axial Heat Diffusion in Liquid Metal-Cooled Ducts with Specified Heat Flux", Trans. ASME, Journal of Heat Transfer, 90C, No. 3, 1968. Agrawal, H., "Heat Transfer in Laminar Flow Between Parallel Plates at Small Peclet Numbers", Appl. Sci. Res., A9, 177-196, 1960. Hennecke, D.K., "Heat Transfer by Hagen-Poiseuille Flow in the Thermal Development Region with Axial Conduction", University of Minnesota Heat Transfer Lab., 1968. Hishida, M., "Temperature Distribution in the Entrance Region of Fluid Flow in A Pipe", Trans. Japan Society of Mechanical Engineers 31, 222, 1965. Larkin, B.K., "Some Finite Difference Methods for Problems in Transient Heat Flow", Chem Engrg. Progress Symposium, Series v61, No. 59, 1-11, 1965. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 96 Lundburg, D.D. and Miller, J.A., "Laminar Convective Heat Transfer in Entrance Region Between Parallel Flat Plates", U.S. Naval Postgraduate School, Research Paper No. 54, 58, 1965. Wang, Y.L. and Longwell, P.A., "Laminar Flow in the Inlet Section of Parallel Plates", AIChE Journal (10), 323-329, 1964. Pearson, J.D. and Serovy, G.K., "Numerical Solution of the Navier—Stokes Equations for the Entrance Region of Suddenly Accelerated Parallel Plates", Iowa State University Publication, Ames, Iowa. Douglas, J., Jr., "A Survey of Numerical Methods for Parabolic Differential Equations?, Advances in Computers 2, 1-54, 19 Peaceman, D.W. and Rachford, H.H., Jr., "The Numerical Solution Of Parabolic and Elliptic Differential Equations", Jour. o the Soc. for Ind. and Appl. Math., 3, 27-41, 1955. Douglas, J., Jr. and Rachford, H.H., "0n the Numerical Solution of Heat Conduction Problems in Two and Three Dimensions", Trans. American Mathematical Society 82, 421-439, 1956. Kays, W.M., Convection Heat and Mass Transfer, McGraw-Hill Book Co., Inc., New York, 1966. 61. f Cess, R.D. and Shaffer, E.C., "Heat Transfer to Laminar Flow Between Parallel Plates with a Prescribed Wall Heat Flux", Appl. Sci. Res., A8, 339, 1959. Gauss, K.F., "Theoria Motus Corporum Coelestium", Werke, 8 (1809). Box, G.E.P. and Coutie, G.A., ”Application of Digital Computers to the Exploration of Functional Relationships", The Inst. Elec. Eng., Vol. 103, 100, 1956. 5 of Box, G.E.P. and Lucas, H.L., "Design of Experiments in Non- linear Situations", Biometrika, Vol. 46, 77-90, 1959. Box, G.E.P., "Fitting Empirical Data", Annals N.Y. Academy Sciences, Vol. 86, 792, 1960. Beck, J.V., "The Optimum Analytical Design of Transient of Experiments for Simultaneous Determinations of Thermal Con- ductivity and Specific Heat", Ph.D. Thesis, Michigan State University, East Lansing, Mcihigan, 1964. Beck, J.V., "Analytical Determination of Optimum, Transient Experiments for Measurement of Thermal Properties", Proc. the Third Int. Heat Transfer Conf., Vol. IV, 74-80, 1966. of 53. 54. 55. 56. 57. 58. 59. 60. 61. 97 Beck, J.V., "Transient Determination of Thermal Properties", Nuclear Engineering and Design, Vol. 3, 373-381, 1966. Carslaw, H.S. and Jaeger, J.C., Conduction of Heat in Solids, Second Edition, Oxford University Press, London, 1959. Forsythe, G.E. and Wasow, W.R., Finite-Difference Methods for Partial Differential Equations, John Wiley and Sons, New York, 1960. Chapman, A.J., Heat Transfer, Second Edition, The Macmillan Company, New York, 1967. Love, T.J., Radiative Heat Transfer, Merrill Publishing Company, Columbus, Ohio, 1968. Arpaci, V.S., Conduction Heat Transfer, Addison-Wesley Publish- ing Company, Reading, Mass., 1966. Beck, J.V., "Surface Heat Flux Determination Using an Integral Method", Nuclear Engineering and Design 7, 170-178, 1968. Hildebrand, F.B., Methods of Applied Mathematics, Prentice- Hall, Inc., Englewood Cliffs, N.J., 1956. Hald, A., Statistical Theory with Engineering Applications, John Wiley and Sons, Inc., 1952. APPENDI CES APPENDIX A Truncation Error in Approximating Equation (2.8) by the ADI Finite Difference Equation (3.12). Let the difference, D, between (2.8) evaluated at the point [iAx,jAy,(n+%)At] and (3.12) be the truncation error at a typical node. Omitting the superscript "+" from the dimensionless variables for simplicity, 2 2 n+% Tn+1 - n 51., nT____1 £11-11 123' 12J+ uj,’i‘6('r +Tn+1) at 5x Fe? axz By? 1,j At 2(AX) 1 1,1 1,1 1 2 + ' €>i(n :11 - 1 (Tn+T 1) 2(AX) Pe ’3 ’3 2 t 662 2 + t +1 - 2A 2 2 616j(TI j - 2 ;) +--—-A———— uj6j5i(TI j - T: ji] = D 4(AX) (Ay) Pe ’ ’ 4(Ay) Ax 1: , (A-l) + 1 - where Tn % = 3(Tn + Tn+1) Using Equations (3.5) through (3.8) 2 —-------l 2 . .1: . +T12+¥> eta: «r 1 . +T .125 2(AX) ,J .1 X ,1 X X ,1 X .J 4T where the variable subscripts denote differentiation, e.g. T 4 = h—Zu x ax Similarly, 2 + ,. + +1 —1-—2-6? :41": i) =§ (A-3) 2(A)’) J 9.] 3.] y 9.] y L] y a] y L] l n+1 ~ 1 n n+1 (? 222 n+1 —-— + A-4 2(AX) 6191.1 Mm) 2‘ x‘iu Tx‘i j) + (T X31: j +TX3‘1,J‘) ( ) 99 Tn+1 _ i,j i,j g n At n _ 1 n n+1 . _ T +.__ _._ At t i,j 2 th i,j 2(Tt‘i,j t‘i,j) +1 +2-£ (A-S) t 13.] t ,J , 1 . 2 2 n+1 N n+1 52‘ 2 6i6j(TI j ' Ti j) = (T 2 2‘? j ' 2 21i j) (Ax) (AV) ’ ’ x y ’ x y ’ M n n+1 mi n n+1 + 12 (T 4 211,, ' 4 211,1) + 12 (T 2 411,3- ' T 2 411,3) (A'6) y X y X y X y (where the higher order terms have been omitted) 2 + +1 2 éiéj I j ' T: ;) = (T 2‘: j ’ T 2‘: j) (AV) (AX) ’ ’ xy ’ xy ’ + MCI ‘11 T ‘n-H') + MCI! ‘n ln-HL) (A 7) 6 321,j 3Zi,j 12 4i,j 4i,j x y x y xy xy Introducing into (A-S) through (A-7) the fact that n n+1 n (At)2 n (T11.j ' T‘i.j) = 'At Tt‘i.j ' 2 TtZ‘i,j and neglecting terms higher than second order we have, Tiff;- - ngj '5 lCII ‘n ‘n+1) L—LAt 2 T ‘n (A-5a) At 2 t i,j t i,j 4 t3 i,j 1 2 2 n n+1 n At 2 n 2 2 616. T1 . - i .) ‘ -At T 2 2 11 . -‘£—§l—'T 2 2 2 (A-6a) (AX) (BY) J :J 3.] X y t 2.] X y t 2 2 + “’1‘2—“5i5jarilj ‘Ti1 3}) =‘AtT 21: ' “LIX—2LT 2212' “'73) (Ay) M , ’ Xy t ’J Xy t 3.] Substituting (A-2), (A-3), (A-4), (A-Sa), (A-6a), and (A-7a) into (A-l) the truncation error becomes, 100 2 At n n 1 n D='£_i)_[T3+UT 2'—"2T22] t txy Pe txy i,j 2 + Ax _1_ Tn+35 _ u TM“?! + munfik] “(A-8) 6 2 4 3 , . 6 4 . . Pe X X :J y 13.] APPENDIX B Stability Analysis for ADI Numerical Method Let E(x,y,t) be the error in the temperature, T, at a point (x,y) and time t due to the accumulation of round off errors in the numerical calculations. * E(x,y,t) = T (x,y,t) - T(x,y.t) * where T is the calculated temperature with a small error and T is the corresponding exact temperature. Changing the usual notation Tn to Tn in order to avoid confusion with the imaginary i, i,j p,q the error expressed in difference notation is, n n * n E = E x t = T - T -1 PA ( P’yq’ n) p.q paq (B ) Using the Fourier series method deve10ped by von Neumann, we ' assume that the distribution of errors in the x,y plane at t = O can be expanded in a Fourier series. iBmx ivy E(x,y) =22Amne e m 11 Since the energy equation is linear, and therefore solutions are additive, we need consider only the propagation of the error due , i x i . . . to a Single term, say e B e YY (the coeff1c1ent 18 constant and may be neglected). To investigate the propagation of this error as t increases, a solution must be found which reduces to elfixele as t a 0. Assume 101 102 _ in ivy fit E(x,y,t) — e e e The solution is said to be stable if the initial error does not grow with time. This condition is met if ‘entl s 1 (B-2) Substituting (B-l) into (3.12) we see that the error function satisfies the same difference equation. En+1 En ' 2 PigAt qu +_2(Ax) u 6 (En + En+1) ='____l§__§ 6 (En + En+1) q p 2(Ax) Pa 1 + 2 2 + + 2 520::n +En 1) + gt 2 2 5 5 (En - E" 1) 2(Ay) q 4(Ax) (Ay) Pe P ‘1 _ At2 u 62(En _ En+l) 4(Ay) Ax q q 2 2 Defining ax = Zififil— and ay =.§L%%l_ as before and substituting E: q = elepelyyqefltn we have after cancelling the exponential ’ ' factors common to each term, [enAt _ 1] = _ 5&5 u [enAc + 1][eiBAx _ e-iBAx] + ‘-l-§[enAt + 1][eIBAX - 2 + e-IBAX] «LienAt + 1][e1YAy - 2 + e'lyAy] axPe ay 1 t o - o o - o - -—--§[enA - 1][elBAx - 2 + e 1BAX][elYAy -2 + e 1yAy] axayPe u At At i AX -i Ax i A -i A +;;2Ax‘;n “1119-8 'eBJEeyy-ZH Yy] (3-3) Making use of the tigonometric relation, 103 iBAx -iBAx 8 Te =SinBAx=2$ianosgAl 21 2 2 we have (elBAx - 2 + e-IBAX) = - 4 Sin2 2%5 Using these relations in (B-3) and grouping the real and imaginary parts we obtain, “At _ A +-iB _ e -C+iD (BA) with A,B,C, and D given below. C = 1 + 4 Sin2 EAE-+-£-‘— Sin2 XAX‘+ 16 Sin2 EAE Sin2 KAI A 2 -'a 2 2 0 Fe y aid Pe D - 2At u Sin £95 Cos QAE’EA Sin2 XAX'+ l] B 2 2 2 '— \Cl 2 W M 2 IBI Therefore, lA + iB| s ‘C + iD‘ ‘ and, ‘enAt| S 1 Thus, the numerical method is always stable. APPENDIX C. Derivation of Fully Developed Temperature Profiles and Sensitivity Coefficients. For the fully developed region, the steady state energy equation reduces to 2 = 3% (C_]_) BY elc 3P4 Kays (45) defines the region of fully developed temperature profiles by the criterion,, L<§t¥> =0 ex 3 b Differentiating, dT T - T dT T - T dT fl=__5___§____3_+_i____b. (C-2) ax dx T - T dx T - T dx 3 b s b Also, since Nu is constant, h is constant and for x>xDev ' the constant heat rate, Case I, q" = h(TS - T = constant. b) Ts - Tb = constant de dTb 3‘21?” (0'3) Substituting (C-3) into (C-2) d d 5:. = .32 = i (G-A.,) ax dx dx ' 105 Equation (C-l) becomes, 3- y_2 u=-2'u(l-az) fi=12:(1_fi)fl 2 a 2 o2 dx BY a Integrating twice, applying the boundary conditions BY we have, -— d 2 2 T = T + §2._32(2_.- Zfi__.- 2E.) 3 2a dx 2 12a2 12 Nondimensionalizing according to Equation (2.7) + 2 4 dT + + T T5+2dx+(2 '12 12 (CS) Consider an energy balance on an element in the fully deve10ped region, &q"Ax I I I I ._ : : __ dTb pCp(2a)u Tb ~H AX lr-v- pCp(2a)u(Tb + E;— AX) *qlle FIGURE C.1: ENERGY BALANCE ON A FLUID ELEMENT 106 '_ dTb —.—= " pCpua dx q dT n or 8 Fe a;§-= 3E3 dT+ —__b_=1 + dX Hence, (C-S) becomes 2 + + 1 + = +— T TS 8(6y + _ l +>+ + Tb-Y uTfl (C-6) (C-7) Substituting the temperature profile given in (C-6) and the parabolic velocity profile and integrating, the bulk temperature becomes, + Tb + =T .- S + To evaluate the unknown TS a control volume which includes the 11_ 35 (C-8) consider an energy balance on entrance and upstream regions. ..Q a .l__La_.aaa3gL_.a...........;&..___1_¢?__O _ .. X-..” 0::0 a T f--- .52» lx>xmw FIGURE C.2: The energy balance gives, q"x+f:kg%dy=j|: ENTRANCE REGION CONTROL VOLUME - d pCpu(T TS) y 107 In dimensionless form, using (C-4) we have where + Integrating and solving for T8 :1ng 1 + + + x+-—-j'-—dy32:fT(1-y)dy Pe dx+ dT; —:=1 from (C-6), dX From (C-9) T =‘x +—2-+— (C-9) 2 4 ..+ _1_. 2+ 1+ 12. - ‘X + 2+4), "sy “280 (cm) Pe =x++—1—2- (c-11) Pe we can obtain the expression for the fully de- veloped steady state sensitivity coefficient, or + T .a._.s. - (a... X+) L.(_.]:._ 8P8 aPe x a ape 2 Pe + is 2 _ as .1. _ _2_ aPe a Pe2 Pe3 + Pe 23; = -(X+'+-£-) (C-12) aPe 2 M. IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII mu«1|mmmunxuwmwwwum