UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
TIMEDEPENDENT ENHANCED HEAT TRANSFER IN OSCILLATING PIPE FLOW By GUOJIE ZHANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1988 To my beloved motherland ACKNOWLEDGEMENTS The author wishes to express his deep gratitude to the chairman of his committee, Dr. Ulrich H. Kurzweg, for the valuable assistance and advice in guiding this research work to its completion and for thoroughly reviewing the entire manuscript leading to the realization of this dissertation. Also he wishes to express his appreciation to Drs. E. Rune Lindgren, Lawrence E. Malvern, Arun K. Varma and David W. Mikolaitis for the many helpful discussions in the formulation of the problem and constructive suggestions for overcoming many difficulties in the solution process. Thanks are expressed here also to Dean Eugene R. Chenette, to the department chairman, Dr. Martin A. Eisenberg, and to Dr. Charles E. Taylor for their support in allowing the author to pursue his educational goals within this lovely country. Part of the work presented here was funded by a grant from the National Science Foundation, under contract number CBT8611254. This support is gratefully acknowledged. iii TABLE OF CONTENTS Page ACKNOWLEDGEMENTS . LIST OF FIGURES . LIST OF TABLES . KEY TO SYMBOLS . ABSTRACT . CHAPTERS I INTRODUCTION II FORMULATION OF THE PROBLEM . Governing Equations . Boundary Conditions . Model 1 . . Model 2 . . Model 3 . . Initial Conditions . Calculation of Tidal Displacement Effective Heat Flux . III NUMERICAL TECHNIQUES EMPLOYED . . Transformation . . CrankNicolson Method for Momentum Equation . ADI Method for Axisymmetric Heat Equations. . Convergence Criteria . . Grid Generation . . IV NUMERICAL RESULTS AND DISCUSSION . Part 1. Oscillatory Pipe Flow Features . Velocity Profiles . . Lagrangian Displacements . . Tidal Displacements . . Phase Lags . . Part 2. Enhanced Heat Transfer Investigation Periodic Temperature Buildup in Thermal Pumping Process . . . . iii . . . vi . . . ix . . . x ix X xii 18 18 . 24 24 28 30 . 33 35 . 37 r r r r rr Temperature Distribution in Model 1 . Temperature Distribution in Model 2 . Temperature Distribution in Model 3 . Heat Flux versus Tidal Displacement (Model 2) . . Influence of Thermodynamic Properties . Heat Flux versus Tidal Displacement (Model 1) . . Influence of Wall Thickness . Influence of Pipe Diameter . Variation of Axial Temperature Gradient In Model 3 . . Comparison of Enhanced Oscillatory Heat Transfer and Heat Conduction . Enhanced Feat Flux as a Function of Wormersley Number . . Tuning Curves . . V CONCLUDING REMARKS . . APPENDIX: ETP COMPUTER CODE . . REFERENCES . . . BIOGRAPHICAL SKETCH . . 95 . 100 . 108 . 111 . 115 . 119 . 124 . 126 . 130 . 133 . 135 * 137 . 143 * 149 . 157 . 161 LIST OF FIGURES Ficrure Page 11 F(a) Curve . . 7 21 Thermal Pumping Device . .... 19 22 Model 1Fixed End Temperature Model. 25 23 Model 2Periodic Heat and Cold Sources on Insulated Wall . 29 24 Model 3Pipe with Extended Conducting Sections 31 31 Grid System Used in the Numerical Simulations. 39 32 Coordinate Transformation ... 42 41 1D Velocity Profiles in Oscillating Flow for Wormersley Number a = 1, 10, 100, and 1000 65 42 Velocity Profiles (a = 5, after Uchida) 66 43 Magnified View of Velocity Profile near Wall 67 44 Lagrangian Displacement for a = 0.1, 1.0 70 45 Lagrangian Displacement at a = 10 71 46 Dimensionless CrossSection Averaged Displacement versus Time (a = 0.1 1.0) 73 47 Dimensionless CrossSection Averaged Displacement versus Time (a = 2 50). 73 48 Relationship Between Dimensionless Tidal Displacement AX and Wormersley Number a 77 49 Relationship Between Tidal Displacement Ax and Exciting Pressure Gradient in Water 78 410 Phase Variation Along Radius for Different Wormersley Numbers . 84 411 Temperature Buildup Process in Oscillating Flow (Model 1, a = 1, Ax = 2cm). .. 89 412 Temperature Buildup Process in Oscillating Flow (Model 1, a = 1, Ax = 5cm) ... 90 413 Temperature Buildup Process in Oscillating Flow (Model 1, a = 1, Ax = 10cm) 91 414 Temperature Buildup Process in Steady Flow (Model 2, Ugve = 1.5 cm/sec) . 93 415 Buildup Time versus Tidal Displacement (Model 2, a = 1). . 94 416 Temperature Distribution in Oscillating Pipe Flow (Model 1, a = 1, Ax = 1cm) 96 417 Temperature Distribution in Oscillating Pipe Flow (Model 1, a = 1, Ax = 2cm) 97 418 Temperature Distribution in Oscillating Pipe Flow (Model 1, a = 1, Ax = 5cm) 98 419 Temperature Distribution in Oscillating Pipe Flow (Model 2, a = 1, Ax = 1cm) 101 420 Temperature Distribution in Oscillating Pipe Flow (Model 2, a = 1, Ax = 5cm). 102 421 Temperature Distribution in Oscillating Pipe Flow (Model 2, a = 1, Ax = 10cm) 103 422 Temperature Distribution in Oscillating Pipe Flow (Model 2, a = 1, Ax = 20cm) 104 423 Temperature Distribution in Oscillating Pipe Flow (Model 1, a = 1, Ax = 30cm) .. 105 424 Temperature Distribution in Steady Flow (Model 2, Uave = 0.5 7.5 cm/sec) 106 425 Temperature Distribution in Oscillating Pipe Flow (Model 3, a = 1, Ax = 10cm) 109 426 Magnified View of Temperature in the Central Pipe Section (Model 3, a = 1, Ax = 10cm) 110 427 Heat Flux in Oscillating Flow and Steady Flow (Model 2, a = 1, Water as Working Fluid). 112 vii 428 Influence of Thermodynamic Properties of H20 on the Enhanced Heat Flux (Model 2, a = 1, Ax = 10 cm). . . 117 429 Heat Flux versus Tidal Displacement (Model 1, a = 1, Pr = 7.03). 120 430 Heat Flux versus Tidal Displacement (Model 1, a = 3, Pr = 7.03). 121 431 Influence of Wall Thickness on Axial Heat Flux (Model 1, WaterGlass, a = 1, Ax = 5cm). 126 432 Influence of Pipe Diameter on Heat Flux for Fixed Frequency (Model 3, Waterglass, a = 3, Ax = 10cm) . . 128 433 Typical IsoTemperature Contour in Oscillating Pipe Flow (Model 3, WaterGlass, a = 3, Ax = 10cm). . 129 434 Variation of Temperature T1 and T2 versus Ax (Model 3, WaterGlass, a = 3) ... 131 435 Comparison of Enhanced Heat Transfer and Heat Conduction in Oscillating Pipe Flow (Model 3, WaterGlass, a = 3) . 134 436 Variation of Axial Heat Flux versus Wormersley Number (Model 3, H20Glass, HgSteel, Ax = 10cm) . . 136 437 Computed Tuning Curves (Model 3, H20Steel and HgSteel, Ax = 10cm) . 138 438 Tuning Curve versus Wormersley Number (after Kurzweg). . 140 439 Ratio of Heat Flux due to Conduction to Enhanced Heat Flux versus Wormersley Number (Model 3, H20Steel, HgSteel, Ax = 10cm) 141 viii LIST OF TABLES Tables Pages 41 Dimensionless Tidal Displacement at Different Wormersley Numbers . .... 72 42 Phase Lags Along Radius (Working Medium: H20, Ax = 10cm, a 0.1 20) . 80 43 Comparison of Phase Lags with Different Working Mediums (Ax = 10cm, a 1, 5). 83 44 Comparison of Enhanced Heat Flux Using Numerical Velocity with Heat Flux Using Analytic Velocity (Model 3, H20Glass, R1 = 0.1cm, R2 = 0.15cm, Pr = 7.03, a = 3) 86 45 Enhanced Axial Heat Flux via Tidal Displacement. 113 46 Enhanced Axial Heat Flux in Steady Flow 114 47 The Influence of Properties of Water on the Enhanced Axial Heat Flux . 116 48 Variation of the Axial Temperature Gradient versus Wormersley Numbers (WaterGlass, Ax = 10 cm) . . 132 KEY TO SYMBOLS x x coordinate r radial coordinate t time 5 coordinate (x) in transformed plane 1 coordinate (r) in transformed plane T transformed time L pipe length R1 pipe inner radius R2 pipe outer radius w oscillating frequency 6 boundary layer thickness P pressure p density c specific heat v Kinematic viscosity p dynamic viscosity K thermal diffusivity k thermal conductivity Ke coefficient of enhanced heat diffusivity Pr Prandtl number A = l/p ap/axI a measure of the maximum axial pressure gradient (cm/sec2) a Wormersley number a = Ju/ov T Temperature I 7 = aT/ax timeaveraged axial temperature gradient S = r/R1 dimensionless radial distance g radial temperature distribution function U velocity Uo representative velocity f velocity shape function X Lagrangian displacement DX dimensionless tidal displacement Ax dimensional tidal displacement Qtotal time averaged total enhanced axial heat flow over pipe crosssection axial heat flux Subscript f fluid w wall h hot c cold th thermal eq equivalent adj adjacent Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy TIMEDEPENDENT ENHANCED HEAT TRANSFER IN OSCILLATING PIPE FLOW By Guojie Zhang April 1988 Chairman: Ulrich H. Kurzweg Major Department: Aerospace Engineering, Mechanics, and Engineering Science The problem of timedependent enhanced heat transfer in an incompressible viscous laminar fluid subjected to sinusoidal oscillations in circular pipes which are connected to a hot reservoir at one end and a cold reservoir at the other end has been examined numerically in detail. Three models were designed for the investigation of such an enhanced thermal pumping process and a computer code (ETP) was developed to implement all the numerical calculations. To increase the understanding of the mechanism of thermal pumping, the periodic velocity profiles and Lagrangian displacements as well as tidal displacements at various Wormersley numbers (from a = 0.1 to 1000) were studied. Some transient problems of enhanced axial heat transfer in oscillating pipe flow such as the periodic final temperature buildup process in oscillating pipe flow were also examined. The timedependent temperature distribution xii in the different models was numerically studied in detail. The enhanced axial heat flux magnitude versus different tidal displacements with water and mercury as the working fluids bounded by pipe walls of different material were observed and the quadratic coefficients found. The influence of the variation of water properties on the enhanced axial heat flux was numerically examined and the results show that the enhanced axial heat flux can vary about 150 percent even within the temperature range from 0C to 1000C. The effects of wall thickness and pipe diameter in enhanced thermal pumping were also studied and the optimum wall thickness was found to be about 20 percent of the pipe radius in the waterglass combination. The tuning effect in the watersteel and the mercurysteel cases was examined and the results show good agreement with analytic predictions. A comparison of the enhanced axial heat flux with the axial heat flow due to heat conduction at various tidal displacement and Wormersley numbers shows that the latter is quite small and negligible provided the tuning condition is satisfied. This study has shown that the enhanced thermal pump is indeed a very effective tool for those problems where large amounts of heat must be transported without an accompanying convective mass exchange. The investigation also indicates that turbulent flow in the reservoirs is preferable to laminar conditions and should receive more attention in future studies. xiii CHAPTER I INTRODUCTION Enhanced heat transport in a viscous laminar fluid subjected to sinusoidal oscillations in a very long pipe which connects a hot fluid reservoir at one end and a cold fluid reservoir at the other end (Fig. 21) has been recognized and studied recently by Kurzweg [15, 16, 17, 22]. The results obtained show that with this oscillatory pipe flow the heat transferred axially from the hot end to the cold end can be orders of magnitude larger than that obtained by pure molecular conduction in the absence of oscillations. In addition, the more important thing of interest is that this heat transfer process involves no net convective mass transport. Major assumptions made in the above cited studies on enhanced heat diffusion are that a constant timeaveraged nonzero axial temperature gradient is always present in the oscillating flow and that the axial molecular conduction along the wall and in the oscillating fluid is negligible. Discovery of this enhanced heat transport phenomenon was made possible by earlier studies on axial dispersion of contaminants within steady laminar flows through capillary tubes by Taylor in 1953 [32], and Aris in 1956 [4]. These 2 earlier studies show that when a small quantity of a contaminant is introduced into a circular pipe, the dispersion of the resultant contaminant cloud is greatly enhanced by the flow of the fluid. Bowden's 1965 results show that similar dispersion effects occur in oscillatory flow [8]. This enhanced axial dispersion of contaminants in the presence of laminar oscillatory flow within capillary tubes was studied in 1975 by Chatwin who suggested that the assumption of constant timeaveraged axial contaminant gradient can be made [10]. Recently, Bohn et al. extended this work to the study of gas component transfer in binary gas mixtures when these are confined to single tubes and a sinusoidal pressure variation is applied [7]. Further studies in 1983 by Watson [38] show that the effective diffusion of contaminants is proportional to the square of the tidal displacement. This has been experimentally verified by Joshi et al. [13] and by Jaeger [12], both in 1983, and most recently by Kurzweg and Jaeger [19], in 1987. All these results show that the contaminant would spread axially in both steady and oscillatory laminar pipe flow at rates as much as five orders of magnitude higher than in the absence of fluid motion. The first significant research work extending these enhanced axially diffusion studies to the heat transfer problem in oscillatory flow within very slender pipes or flat plate channels is due to Kurzweg [15, 16, 22]. In 3 early 1983 Kurzweg suggested that a similar dispersion process should occur in the heat transfer area because of the similarity in both the governing diffusion and heat transfer equations [16], and the first preliminary theory was formulated in 1985 [15, 16], in which, referring to Chatwin's idea, a time averaged constant axial temperature gradient assumption was used. The instantaneous temperature distribution was taken to be of the form [15, 22]. T = y[x + R1 g(n)e ] (11) where y = aT/ax is the timeaveraged axial temperature gradient, R1 is the tube radius, x is the axial distance along the capillaries (with x = L/2 and x = L/2 denoting the ends), L is the pipe length under consideration, w is the oscillating frequency of fluid, g(S) is a radial temperature distribution, and = r/R1 is the dimensionless radial distance. The theoretical analysis shows that under certain conditions such enhanced axial heat diffusivity can indeed be significantly larger than the axial molecular conduction [15, 16] and this has been verified by some experimental measurements by Kurzweg and Zhao [22]. In order to better understand the physical mechanisms of this interesting and potentially useful heat transfer technique, we shall examine in greater detail the thermal pumping model shown in Fig. 21. It is assumed that a bundle of very thin and long tubes is connected to a 4 reservoir which supplies unlimited hot liquid at one end to a second reservoir which supplies unlimited cold liquid at the other end. The liquid in the pipes is oscillating axially with an amplitude such that none of the liquid which is originally in the middle portion of pipes ever runs into either reservoir. That is, there is no net convective mass exchange between two reservoirs. The largest axial fluid dimensionless displacement (when crosssection averaged) is referred to as the nondimensional tidal displacement and is denoted by "AX" (it should not be confused with the dimensional tidal displacement "Ax" frequently used in the present study). At time t = 0, the fluid within the pipes is set into axial oscillations at angular frequency w and tidal displacement AX. After a short transient, this oscillatory motion will lead to very large axial heat flows which can be readily made to exceed those possible with heat pipes. Before exploring the mechanisms of this enhanced heat transport, it is necessary to introduce some new concepts which are commonly used in the study of this type of oscillatory motion. As is well known [36], for high frequency viscous laminar axial oscillations within fluid flow along rigid pipes, the nonslip boundary condition creates a very thin Stokes' viscous boundary layer of thickness 6 = 72 v/ (12) 5 where v is the fluid kinematic viscosity. For room temperature water at a frequency of 10 Hz, this viscous boundary layer is approximately 1.7.102 cm. The corresponding thermal boundary layer thickness is about 6th = 6//Pr (13) where Pr is the Prandtl number. Note that both 6 and 6th decrease in thickness as the oscillating frequency w increases. In the theoretical analysis, it is always assumed that a fully developed velocity profile of the oscillating flow exists within the pipes. At high frequency w, this flow consists essentially of a slug flow over most of the fluid core bounded by a thin boundary layer of width S as discussed by Uchida [36]. Neglecting end effects, the fully developed oscillating laminar velocity profile in pipes, due to a periodic axial pressure gradient, is found to be [19] U(S,t) = Uof(C)eiW (14) where Uo is a representative velocity, r = r/R1 is again the dimensionless radial distance, f() the velocity shape function, and w the angular velocity of the oscillatory flow. The explicit form of f(r) is af() = 1 (1iaC ) (15) 6 where A = R12 8p/a8x/pvUo is the nondimensional pressure gradient maximum acting along the capillaries, a = R1J/wlv is the Wormersley number measuring the ratio of inertia to viscous forces, v is the fluid kinematic viscosity of the carrier liquid, and p is the fluid density. This velocity profile will reduce to the familiar Poiseuille parabolic shape as the angular frequency w becomes small, while at moderate frequency, f() has the shape demonstrated by Uchida [36]. Another new term commonly used when dealing with oscillating flow is the crossstream averaged dimensionless tidal displacement AX which can be mathematically defined as 4Uo 1i AX =I4 J~ff(s)dr (16) and on integration, yields [18] ap/ax 2 AX = 1/2 1 + F() (17) where the complex function F(a) has the form F(a) =i J( ia) with the prime denoting differentiation. Using the definition of Kelvin function: Jo(,/ia) = ber a + i bei a ( u OD 00000 + ++ + 4. co r 0d Y 1 , I iCOM 'cT nt) o d LL.  30 d1 8 [1], the complex function F(a) (Fig. 11) can be further written as: F(a) = Fr(a) + i Fi(a) = i ber a + i bei'a (18) ber a + i bei a This dimensionless tidal displacement is related to the maximum of the periodic pressure gradient via [18] 2 a << 1 AX 1/2 p2 J2 (19) 1  a >> 1 Apparently, for any fixed tidal displacement Ax and oscillating frequency w, the axial pressure gradient required is proportional to the inverse square of the Wormersley number when a is small; however, it is independent of a when a is very large. This also implies that the exciting axial pressure gradient ap/axI is approximately proportional to the fluid kinematic viscosity v and inversely proportional to the square of the pipe radius R1 when Wormersley number is very small (it happens only at low w, small R1, and large v, for example, oil), while it is almost independent of the fluid kinematic viscosity and the pipe radius when the Wormersley is very large (it happens only at very large R1, high w, small v, 9 for instance, a liquid metal). Note, if the tidal displacement AX is fixed, while allowing the oscillating frequency to change, the axial pressure gradient ap/ax can become very large when the oscillating frequency becomes large. This is mainly due to inertial effects and not so much due to viscous drag forces which dominate the oscillatory flow at small Wormersley number. With the definition of the above quantities, we are now in a position to explore the details of the enhanced axial heat transfer in oscillatory flow within pipes. It is assumed that a constant temperature gradient exists along the pipe in the axial direction and that a very large time dependent radial temperature gradient variation is superimposed. When the fluid moves towards the cooler side (we term this the positive stroke), the hotter fluid within the pipe core which is initially brought into the pipe from the hot reservoir produces a large radial heat flow via conduction to the cooler portions of the fluid within the Stokes' boundary layer and to the cooler pipe wall; while during the negative (or reverse) stroke, i.e., when the fluid moves towards the hotter side, the higher temperature in the boundary layer and the pipe wall conducts the heat back into the cooler fluid core. This coupled radial heat conduction with an axial convective transport leads to an enhanced axial heat flux along the entire length of the pipes. 10 Further, from the system point of view, the heated fluid near the cold reservoir will eventually be ejected into the cold reservoir and mixed there with the lower temperature liquid. Contrarily, near the hot reservoir side the fluid within the pipes which has been cooled during each positive stroke is pushed out into the hot reservoir and mixes with the higher temperature liquid. This process of thermal pumping is what leads to a timeaveraged heat flow from the hot reservoir into the cold reservoir. It differs essentially from the working principle of a normal pump. For a normal pump one can draw the analogy with transport of a oneway vehicle which transports passengers as well as the carrier from one point to another. For the thermal pump one can draw the analogy with a twoway busline which periodically loads and unloads the passengers (heat) from the hot reservoir to the cold reservoir, heat can be continually transferred, and the carrier, in the time averaged sense, does not move. This property is particularly important for those systems where a large amount of heat transfer is needed while the working fluid is required to remain in the system (as in nuclear reactors). Note that the axial heat conduction, in general, is assumed to be very small in this thermal pumping process compared with the enhanced axial heat flow [16]. Apparently, the heat transport rate in thermal pumping is governed by both the thermal properties of the 11 working medium and pipe wall and the characteristics of the oscillatory motion of the fluid. The enhanced axial heat flow does increase with increasing oscillating frequency as this thins out the boundary layer and leads to an accompanying increase in the radial temperature gradient. This observation holds only as long as the thermal properties of the liquid and the wall are compatible. If the molecular conduction of the fluid in the radial direction is very small, then even high frequency oscillatory motion will not produce a large increase in the rate of enhanced axial heat flow. This is because such a system fails to supply the "passengers" enough time to be loaded onto the "bus" and to be unloaded from the "bus". Obviously, the system just wastes energy. On the other hand, if the molecular conduction of the fluid in the radial direction is very large, but the frequency of the oscillatory motion is low, once again one can not expect that there will be an efficient enhanced axial heat transfer between two reservoirs because the "bus" is now moving too slowly. The above observation can be confirmed by an analysis of the performance of a waterglass combination (i.e., water is the working fluid medium within a glass tube) and of a mercurysteel combination. For the former, it is necessary to employ rather small diameter tubes and low frequency oscillations with large tidal displacement, for it has 12 relatively poor heat diffusive properties as compared with the mercurysteel case. This small tube diameterlower oscillation frequency setup is necessary in order to ensure that there is sufficient time to transfer the excessive heat content of the bulk water core to the tube wall during the positive stroke of each period and to permit the transfer of excess heat from wall to the cooler core fluid during the negative stroke. Otherwise, the water would just carry a portion of its heat content back and forth in the pipe and the condition for achieving optimum enhanced axial heat transfer could not be met. For the mercurysteel case, one can chose relatively large pipe diameters and higher oscillation frequencies with smaller tidal displacement because of the higher thermal conductivity. This assures that only very short times are needed to exchange the heat between the core of the fluid and the wall. It should be pointed out that the suggestion of using smaller tidal displacement is purely due to the mechanical considerations and that one always tries to keep Ax as large as possible in order to produce large axial heat flows. From the above discussion, it can be concluded that the process of enhanced heat transfer via oscillatory pumping requires a precise tuning of parameters governing the enhanced heat transport. Indeed there is an expected "tuning effect" as discussed in references [16, 20, 21]. 13 The tuning effect is a very important concept in the study of the presently considered heat transfer process. It shows that there will be an optimum combination between thermal properties of the working medium and wall and the characteristics of the oscillatory motion. The qualitative aspects of the tuning effect have been observed earlier for both the case of a flat channel and that of the cylindrical pipe [15, 16]. From Fig. 438, one can see that an optimum for axial heat transfer occurs only at or near the tuning point which depends on the oscillating frequency and the thermodynamic properties of the fluid and wall. As has been pointed out by Kurzweg [20], in order to obtain the optimum enhanced heat transfer one has to carefully select suitable values for the pipe size, the material of pipe and the working medium as well as the manner of oscillatory motion. The nondimensionalized enhanced heat diffusivity is defined as Ke A = w~X2 (110) where Ke = 0/ypc is the coefficient of enhanced heat diffusivity, 0 is the axial heat flux, 7 is the time averaged axial temperature gradient, p is the density of the fluid, and c is the specific heat. One can show that this nondimensional enhanced heat diffusivity is a function of both the Wormersley number and the Prandtl number [21] and hence that the dimensional axial heat diffusivity Ke is a 14 function of the tube radius R1, the oscillating frequency w, the kinematic viscosity v, and the square of the dimensionless tidal displacement AX. This can be explained from the fact that the radial heat flow is proportional to the product of the representative radial temperature gradient YAX/Sth and the surface area per unit depth of rAX available for crossstream heat transport. The use of large tidal displacement is always beneficial in the enhanced axial heat transfer within oscillating pipe flow. However, in order to avoid the direct convective net mass exchange between the two reservoirs, the tidal displacement must be limited to less than about one half of the pipe length. As has already been predicted by theoretical studies and will be confirmed by the present numerical simulations, the axial heat transfer will be further enhanced if the rigid surface (part of the rigid wall with finite thickness) has a nonzero thermal diffusivity and hence heat storage capability. Note that the existing considerations are restricted to laminar flow. Turbulent flow conditions can occur in oscillating pipe flow at higher values of wAx2/,v [26, 27] and apparently would destroy the assumptions of the current analytic model of the thermal pumping process. Fortunately, the condition for optimum enhanced heat transfer in such oscillating pipe flow obtained at the tuning point requires 15 very slender pipes, such that the Reynolds number is usually small enough so that the oscillating motion falls within the laminar range [22]. The theoretical aspects of the oscillatory enhanced axial heat transfer process have been developed much further than its experimental and numerical counterparts. The theoretical predictions are quite limited and consider only cases under certain simplifying assumptions [8]. Numerical work is necessary in order to not only to examine the correctness of the theoretical analysis but also to further the development of advanced enhanced thermal pumping devices. Numerical studies are not only fast, economical and accurate, but they also offer a handy way to access complex geometries which can not be handled analytically. It is the purpose of this study to extend the analytic work on thermal pumping by a detailed numerical study. We intend first to examine some transient problems of axial heat transfer in oscillatory pipe flow, such as the development of the velocity profile at various Wormersley numbers in contrast with those of reference [36], where only several special cases with intermediate Wormersley number a were discussed, and to examine numerically the relationship between the tidal displacement and the required corresponding strength of the periodic pressure gradient as a function of Wormersley number a and of tidal displacement Ax. 16 Next, we examine the buildup process of the temperature distribution in a pipe which connects a hot reservoir at one end to a cold reservoir at the other end and see whether there actually exists a constant time averaged temperature gradient along the pipe axis when the final periodic state is eventually reached. Note that a timeaveraged linear temperature distribution along the axial direction is an essential assumption in the existing theoretical studies. The third part of this investigation which forms the main effort, is a computeraided numerical simulation of the thermal pumping technique, including an investigation of the variation of the enhanced axial heat flux versus the tidal displacement, the variation of enhanced axial heat flux versus different Wormersley numbers, and a study of the variation of heat flux versus different Prandtl numbers. It also includes a study of the influence of wall thickness and pipe diameter as well as the change of the fluid properties on such an enhanced axial heat flux and an examination of the tuning effect in the conducting wall case. Further, we compare axial heat transfer in oscillating flow with that in the steady flow, and compare quantitatively this enhanced heat flux with that induced by the pure axial molecular conduction under various Wormersley numbers and different fluidsolid wall combinations. Finally, the author examines the phase lag phenomenon in oscillatory flow in order to 17 further reveal the mechanism whereby this new heat transfer technique functions. The computational work presented here was completed on the Vax11/750 electronic computer in the Department of Aerospace Engineering, Mechanics, and Engineering Science, and on Vax11/780 in the Center for Instructional and Research Computing Activities, University of Florida, Gainesville, Florida. The numerical approach used, for solving the presently considered heat transfer problem, employed a secondorder CrankNicolson scheme and an Alternating Direction Implicit method (ADI) with Thomas algorithm. A Fortran computer code named ETP (Enhanced Thermal Pumping) was developed to implement all the calculations. CHAPTER II FORMULATION OF THE PROBLEM A single tube with inner radius R1, outer radius R2 and length L (such that L >> R1 and R2) connects a large reservoir of hot fluid at temperature T = Th at one end, to another large reservoir of the same fluid at cold temperature T = Tc at the other end (Figs. 21 and 22). The tidal displacement Ax and the oscillating frequency w are adjusted to assure that the turbulent flow does not occur. The tube is oriented in such a manner so that the effect of gravity on the oscillatory motion of the fluid in the pipe is negligible. Water and mercury are employed as the working mediums. They are understood to be Newtonian and incompressible fluids. The variations of thermal properties of the fluid with temperature during the heat transport process are assumed to be negligible. With these assumptions the problem of enhanced timedependent heat transfer induced by a simple harmonic oscillatory laminar fluid motion in a very long circular pipe can then readily be formulated. Governing Equations The use of a very slender capillary circular tube with constant crosssectional area and neglecting end 0 J LL LU 3 I0 W ~ 0 I Z )8 I CL~ I0 Qz L 0 0 o 0U u, o J5 O 0 C.. Q (LV a \. \0 1r >0 C) a I o cn 0r w 1 4 EC 01 41 r10 20 effects insures a laminar axisymmetrical onedimensional timedependent motion. It is convenient to employ cylindrical coordinates for this problem and we denote the coordinate in the axial direction of the tube by x, and the radial direction by r. The axial velocity can be taken to be independent of x and, in order to satisfy the requirements of continuity the other velocity components must vanish. We shall further assume that the pressure gradient induced by moving the piston (Fig. 21) is harmonic and has the form [29] 1 ap S= A cos wt (21) where A = l/p ap/ax is a constant which measures the maximum pressure gradient existing along the xaxis. Clearly this equation implies that we are now dealing with a timedependent sinusoidal pressure gradient which is constant over the pipe crosssection at any instant and the pressure varies linearly along the xaxis. The simplified NavierStokes equation for this problem is [36] aU v 8 aU t~ Acos(wt) + B [r(a)J 0 < r < R1 (22) where U(r,t) is the time and radially dependent axial velocity component. 21 The corresponding temperature T(x,r,t) of the fluid within the pipe is governed by the heat conduction equation [15] aT OT 1 a aT a2T at U 8x + Kr r (r + a x2 0 < r < R1 (23) where R1 is the inner radius of the tube and xf the fluid thermal diffusivity which is related the thermal conductivity kf by kf = pc Here p is the fluid density, and c is the specific heat. Note that the viscous heating term has been neglected in equation (23) since it is very small for most experimental conditions; this is justified provided one does not deal with very high Prandtl number fluids such as oils. The temperature in the wall can be determined from the solution of 8T 1 aa T a2T at w[I r arr8ar) + ~ax\ R1 < r < R2 (24) where Kw is the thermal diffusivity of the conducting wall in R1 < r < R2 and is defined by 22 kw PwcW where kw, py, and cw, are the wall conductivity, density, and the specific heat, respectively. By introducing the following nondimensional terms, equation (22), (23) and (24) can be treated more easily. t x t / x 6 r U r* U*  6 A/w T T = Th Tc where 6 = J12i/ is again the fluid viscous boundary layer thickness. The dimensionless governing system of equations can then be written as aU* 1 a2U* 1 au* = cos(t*) + 2 [ r*2 + w 0 < r < R1 (25) aT* OT* 1 82T* 1 aT* 82T* at = CU* x + 2Pr + + a x 0 < r < R1 (26) and T* 1 a2T* 1 T* a2T* a*= 2S + + ax r * R1 < r < R2 (27) where Cf Pr 'UCw S  Cw and A C= ~^ 076 This nondimensionalization has some advantages in the computing process to be carried out below. The dimensionless velocity and its distribution over the cross section found from the momentum equation (25) are expected to be universal for any Wormersley number a = R1Jw p/ and any associated quantities, such as the tidal displacement and Lagrangian displacement. Its final periodic form is of the form given by (14). The dimensionless temperature in the pipe is only related to the Prantdl number and the dimensionless velocity, while that in the wall is related to the ratio of wall heat diffusivity to the kinematic viscosity, as seen from equation (26) and (27). The governing dimensionless equations (25), (26), and (27) are a set of secondorder parabolic type of partial differential equations expressed in cylindrical coordinates. To solve this set of simultaneous equations, a corresponding set of boundary conditions and initial conditions are required. Since the velocity is assumed to be a function of r and t only, just two boundary conditions 24 are needed for solving the momentum equation, while for the temperature T(x,r,t) the heat conduction equations require two boundary conditions in both the r and x directions as well as compatibility conditions along the interface between the fluid and the solid wall. It should be pointed out that the initial conditions are less important than the boundary conditions if one seeks only a final periodic state. Boundary Conditions The boundary conditions for the velocity are the usual ones for viscous flow, namely, that the velocity is zero at the inner surface of the wall (r* = Rl*). Also by symmetry, the normal derivative of velocity at the axis is zero. That is, U*(R1 *, t*) = 0 (28) and 8U*(0, t*) ar = 0 (29) The boundary conditions for the temperature depend on the particular model investigated. Model 1 In this model, it is assumed that the temperature at each end of the pipe is equal to that of the connecting fluid filling the reservoirs when the fluid moves into the capillary tube (i.e., during one half of each cycle). The r( 0 E E XO 0 3Wr 0 k II a, X * 26 end boundary temperature at x = 0 and x = L will take the values of temperature of the adjacent pipe fluid elements. With this model, we tried to simulate the real experimental observations where it can be clearly seen that a fluid jet exits the pipe during the outstroke, while the elements of the fluid enter the tube in a funnel pattern during the instroke. This observation shows that there is enough time during each oscillating cycle to allow the fluid elements within the exiting jet to fully mix with the fluid particles which originally are in the reservoir before the next instroke, provided the oscillating frequency is not too high. The end temperature boundary conditions are here taken as T*h (when fluid enters pipe) T*(0, r*, t*) = T adj (when fluid leaves pipe) (210) [T*c (when fluid enters pipe) T (L r t ) = T*adj (when fluid leaves pipe) (211) where T*c, T*h are the nondimensional temperatures of cold and hot reservoirs, respectively. While T*adj is the 27 temperature of the fluid elements which are adjacent to the corresponding ends at a particular instant. The thermal boundary conditions along the outer surface r* = R2 of the pipe wall are taken to satisfy the insulating wall condition, and along the axis of the tube the temperature is assumed to meet the symmetric boundary requirement, i.e., the radial derivatives of temperature along axis are equal to zero. We thus have 9T* r I R = 0 (212) and 8T* S* = 0 (213) r =0 The compatibility conditions along the interface between the fluid and the solid wall are that the radial heat flux and temperature are constant across the interface. That is a8T* aT* kf r = kw r* at r = R1 (214) fluid wall and T* = T* at r = R1 (215) fluid wall Since numerically the same nodes are chosen along the interface, condition (215) will be automatically satisfied. Model 2 Here it is assumed that a heat source rim of width 2b is directly mounted on the interface between the solid wall and the fluid at x = L/2, while two small cold rim sources of width b each are mounted at x = 0 and x = L (see Fig. 2 3). The wall thickness is assumed to be zero. This model is intended to simulate the enhanced heat transfer process of oscillatory fluid in an infinitely long pipe which is heated and cooled by the alternative evenly distributed heat and cold sources along the outer surface of a very thin wall which possesses infinite heat conductivity. Apparently, this geometry can be well simulated by the present model with periodic boundary conditions. Mathematically the boundary conditions are T*(x*, Rl*, t*) = T*h X2* < X* > X3* (214) 0 < x* X1* T (x*, R1*, t*) = Tc or (215) X4 : x* x< L* aT* X1* < x* < X2* = 0 or (216) X3* < X* < X4* where Xl*, X2* X3*, X4*, X5* and X6* are nondimensional coordinates of points which correspond to X = 0, b, L/2b, L/2+b, Lb, and L, respectively in Fig. 23. The periodic end boundary conditions are given as U) 0 V H 0 U ^ :3 o H O H II V 0 0) IWU a 0C Oar' Pfr 30 T*(0, r*, t*) = T*(L*, r*, t*) (217) and 8T* aT* ax X* ax *(218) x 0 x =L Model 3 In model 1, the boundary conditions for the heat equation at X = 0 and L are somewhat artificial. The temperature distribution in the vicinity of pipe ends may not precisely match the constant temperature end conditions. This is particularly true in the large tidal displacement and/or very high oscillating frequency cases since a temperature drop occurs at the hot end and a rise at the cold end as fluid particles exit the pipe. This leads to a discontinuity of temperature which will lead to unexpected dispersion errors if an even order numerical method is used [3, 24, 30]. To improve the end boundary conditions presented in model 1, the following extended model has been considered. It is assumed that an extension pipe which is some 5 times the length of the central pipe is attached to the original pipe used in model 1. The heat and cold source are assumed to be located on the outer surface of the wall in the extended sections of the pipe as well as at both ends of the tube as shown in Fig. 24. This model is used to 0 IS 0 0 .I CH o U II cr 0 u >N 0 Ir 32 investigate the situation where the heat contained in the jet is exchanged by pure heat conduction with the surrounding fluid elements in each reservoir without any convective mixing. In this case, the boundary conditions along the outer surface of the wall are T* = T*h T = T*c O s X* < X1* X2* X3* 8T* ir^ 0 Xl X* < X2* (221) where Xl*, X2*, and X3* are nondimensional coordinates of points which correspond to X = 5L, 6L, and 11L, respectively, in Fig. 24. At both ends, we have T *(0, r*, t*) = T*h (222) T* (L*, r*, t*) = T*c (223) As in the other models, the symmetry condition along the axis of pipe requires that aT* ar r=O= 0 (224) and (219) (220) and 33 and the radial flux and temperature continuity condition along the interface are aT*  T* kf far li = kw w (225) fluid wall and T* = T* (226) T fluid wall Here again kf and kw are the same as described in the previous models. Initial Conditions The initial condition chosen depends on the problem under consideration. If one's goal is to investigate the periodic quasisteady solution only, the initial conditions chosen should be as close to the quasisteady state as possible so that a rapid convergence rate at low CPU time cost is achievable. If one intends to study the transient process, then various initial conditions should be supplied according to purpose of the investigating cases selected. For the velocity initial condition we choose for all our studies U*(x*, r*, 0) = 0 (227) For the temperature, if the purpose of investigation is to examine the buildup process of the periodic quasisteady 34 state in the thermal field, the initial condition should be taken as T*(x*, r*, 0) = 0 (228) However, for the other cases, in order to gain faster convergence, the initial temperature can be assumed to have a linear distribution in the axial direction, and thus have the form (T** T*h)xx* T*(x*, r*, 0) = T*h + LL* (229) where LL* is the dimensionless length within which the axial linear temperature distribution is assumed to hold, and xx* is the dimensionless distance measured from the "origin" which is chosen only for the purpose of this linear temperature initialization. Both the LL* and the origin selected depend on the model considered. For Model 1 and Model 3, LL* is equal to L*, and the origin is taken at the left end (Model 1), or at the left intersection between the central pipe and the left extension pipe (Model 3). However, for Model 2, equation (229) is valid only for the right half portion of the pipe as LL* is taken to be L*/2, and the origin is chosen at the middle section of the pipe. The initial temperature of the left half pipe can then be found from the plane symmetric condition about the origin cross section. 35 Calculation of Tidal Displacement An important quantity encountered in the study of the enhanced heat transfer process in oscillatory pipe flow is the tidal displacement, which is usually required to be smaller than one half of the total pipe length in order to avoid any convective mass exchange occurring between the two reservoirs. It has already been defined in the introduction chapter (16) and in the present computation the dimensional tidal displacement takes the form 2 R1 AX = R12 f x(r,w/w)rdr (230) where R1 is the inner radius of the pipe and x(r,w/w) is the Lagrangian displacement of the fluid elements located along a radius within the capillary tube at t = w/w. It is assumed these elements are initially lined up at axial position x = L/2 (Model 1 and Model 2) or x = 5.5L (Model 3), half way between the tube ends. Numerically, the dimensional Lagrangian displacement at time t is computed via the equation t x(r,t) = U(r,t)dt (231) 0 It is obvious that the dimensional Lagrangian displacement x(r,t) is a function of both time t and the radial position r. Note that, since the existence of phase lags between the 36 stimulating pressure gradient and the displacements vary for different Wormersley numbers, in actual calculations, the tidal displacement is equal to the sum of the absolute maximum crosssection averaged Lagrangian displacement and the absolute minimum crosssection averaged Lagrangian displacement within a period. With the same non dimensionalizaton procedure used in the previous section for the governing equations, the dimensionless Lagrangian and tidal displacement can be written as R* AX = 2 f X*(r*,r/w)r*dr* (232) The nondimensional Lagrangian displacement within the period can then be computed by X* (r*,t*) = U*(r*,t*)dt* (233) *0 where = X (A/w2) and AX AX = It is pointed out that this dimensionless tidal displacement Ax* differs from AX defined by equation (17) by a constant 1/2. 37 Effective Heat Flux The time averaged total effective axial heat flow within the pipe has the form Qtotal = cpw dt I UTrdr (234) where c is the specific heat, p is the density and again, w is the oscillatory frequency. The per unit area effective heat flow (termed heat flux) can then be written as = (235) rR1 Equation (234) follows from the fact that pcUT is the convective heat flux. The dimensionless total effective heat flow can be written as 2n R1* Q*total = cpw dt* I U*T*r*dr* (236) and the dimensionless heat flux by Q total S= *2 (237) CHAPTER III NUMERICAL TECHNIQUES EMPLOYED Equations (25), (26) and (27) which were derived in Chapter II can not generally be solved analytically in any simple manner. Therefore, it is necessary to seek a numerical solution approach to the problem. As is well known, in order to numerically solve a set of simultaneous governing differential equations more accurately and efficiently, an optimized grid system plays a very important role. For the present purpose the grid net works were generated in the following way: within the tube, a nonuniform grid (see Fig. 31) is clustered along the radial direction in the vicinity of the inner surface of the pipe where the larger gradients of velocity and temperature are expected to be present, while along the xaxis the grid is distributed uniformly except for those calculations dealing with model 3. The solution process was carried out in the computational plane rather than directly in the physical plane. Thus, a transformation which converts the governing equations as well as the boundary condition plays an essential role. 39 r 1). For Model 1 1). For Model I wall 2). for Model 3 Fig 31 Grid Systems Used in the Numerical Simulations 40 The secondorder implicit unconditionally stable CrankNicolson scheme and the ADI method (2, 31] are employed to break up the transformed equations into finite difference form. A computer code designated ETP (for Enhanced Thermal Pumping) has been developed for obtaining the desired results. Details of this procedure will be described in this chapter. Transformation The numerical solution of a system of partial differential equations can be greatly simplified by a well constructed grid. It is also true that an improper choice of grid point locations can lead to an apparent instability or lack of convergence. Early work using finite difference methods was restricted to some simple problems which can be numerically solved in the physical domain. As experience was gained, general mappings were used to transform the physical plane into the computational domain as well as the governing differential equations [31]. Such a grid transformation technique brought to the numerical simulation numerous advantages, and the computational work was no longer restricted to a few simple geometries. A bodyfitted nonuniform grid in the physical domain could be used, which retains the uniform spaced grid system features in the computational domain [33, 34, 35]. However, several requirements must be placed on the transformation: first, the mapping must be one to one, second, the grid lines must 41 be smooth and have small skewness in the physical domain, third, grid node point spacing should be small where large numerical errors are expected (i.e., large solution gradient regions) in the physical domain (Fig. 32). In present study, the grid system in the physical space is numerically generated by solving an algebraic equation which clusters the grid lines in the region where large gradients of the physical quantities are expected so as to gain higher resolution of these physical quantities. Because of the simple pipe geometry, the requirements of grid smoothness and orthogonality are not a serious problem. For the present purpose, the transformation is a simple one which transforms a nonuniform grid network in the physical plane onto a uniform grid system in the computational plane. x* = x*(e) r* = r*(s) (31) t* = t*(r) The inverse transformation can be found as = ( (x*) S= P7 (r*) (32) r = r (t*) where x*, r* are the dimensionless coordinates and t* the dimensionless time in the physical domain, ,. 1 are the transformed coordinates and r is the transformed time in the 3 C i' 0 e 0 0 H o 0 4V 4. 41 0 r. u 0. 1I 'I computational domain. With this transformation, the transformed governing equations have the form: au a2U* au* f(r) + a()) g + b(,q) a (33) aT* BT* a2T* 8T* a2T* r TP1 + P2a + P3 + P4 (34) and OT* aT* 82T* 8T* 82T* r w + w2 2+ w3 + w4 (35) where at* ar f(r) = Or cos (t*) (36) 1 at* 1 a() = 2 Or (ar*/a,)2 (37) 1 at* 1 1 a2r* b() 2 r*(ar*/,) (ar*/a,)3 (38) and at CU* 1 a2x* P1 r tx*/a~ 2Pr(9x*7/a)3 + 2 + (39) at* 1 P2 = a7 2Pr(ax*/aO2) (310) at* 1 1 1 a2r* P3 r 2Pr r(r (r*/8Oj) (ar*/)3 q2 (311) at* 1 P4= r 2Pr (ar /a)2 (312) 44 at* 1 1 a2x* W1 ar 2PrW (xl 0 8 (313) 72Pr;(ax*/aa) a3 " at* 1 1 W2 = 87 2Prw (ax*/a)2 (314) at* 1 1 1 a2r* W3= r 2PrW (r*(ar*/a8,) (ar*/a?)3 7J (315) at* 1 1 W4 ~9 2Pr; (ar*/a8)2 (316) As mentioned in chapter II, equations (33), (34) and (35) are a set of secondorder parabolic partial difference equations in cylindrical coordinates. In addition, since the oscillating flow is considered as incompressible in the present studies, the momentum equation (33) can then be independently solved at each time step. As a result, the timedependent update velocity found can then be substituted into the heat conduction equation (34) as a known quantity at the same time step level. Eventually, the coupled heat conduction equations (34) and (35) are solved simultaneously to obtain the temperature distribution both in the fluid and in the wall at any time. To best solve this set of equations in terms of accuracy and efficiency, the proper choice of numerical technique and grid network is dictated by an understanding of the physical aspects of the problem. 45 The same transformation should be also applied on all the boundary conditions proposed in the three different models. For the temperature compatibility conditions along the interface between the fluid and the wall (225) one has the following transformed forms: aT* 1 aT* 1 a17 (Or a) fluid a8q (ar /l) 'wall or 3T* aT* Slow= Ka wall (318) where kw ar*/8a pipe Kg = r*/a wall (319) kf ar*/8 1,ll To make the subsequent form of the corresponding finite difference governing equations less cumbersome, the superscript will be dropped from the variables T*, r*, t*, and x*, and in addition, U* will be replaced by V. In the process of deriving the finite difference governing system equations, the secondorder central differences in the domain and forward or backward differences along the boundaries or interface of the fluid and wall have been employed at each nodal point. 46 CrankNicolson Method for Momentum Equation The secondorder accurate CrankNicolson method is quite well known and widely used. It is an unconditional stable, implicit scheme for solving the parabolic types of partial difference equations. When the CrankNicolson method is applied to equation (33), the finite difference algorithm at a typical node k in the radial direction and at the time step n assumes the simpler form: (Bkk1 + DV + A Vk+ n+ = C SkOkl kk kk+lj k k = 1, 2, .....km (320) where C= f n+l+ fn Bk1 + EkVk Akk+ ]n (321) 1 Ak 4 ( 2 a + b )k (322) 1 Bk = 4 (2 a + b )k (323) Dk = 1 + ak (324) Ek = 1 ak (325) 47 The boundary conditions along the interface (i.e., the inner surface of the wall, k = kmid) and the axis (k = 1) then become V k = kmid= 0 (326) 8V 8ri k = 1= 0 The initial conditions of the velocity at all nodal points is taken to be zero. k in = 1 k = 1, 2, 3, .....kmid (327) Equation (320) associated with the boundary conditions (3 26) can then be written in the matrix form. This yields a set of linear system algebraic equations which can be solved in terms of the nodal values of velocity in the capillary tube by using either an iterative method or a Thomas algorithm at each time step. The explicit form is A* n n+l n D A V C 1 1 1 1 B2 D2 A2 V2 C2 B3 D3 A3 V3 C3 Bkd1Dkd1 Akd1 kd1 kd1 Bkd Dkd Vd kd (328) where A = A + B kd = kmid1 ADI Method for Axisymmetric Heat Equations The governing PDEs (33), (34) and (35) are all of the secondorder parabolic type. Thus it might be suggested that the CrankNicolson scheme used in solving the momentum equation (33) can also be applied to the axisymmetric heat equations (34) and (35), and one can then take advantage of the tridiagonal matrix form while using this unconditionally stable technique. However, when attempting to use such a formulation, one immediately finds that the resulting system of linear algebraic equations is no longer of the tridiagonal type (320) but rather a nontridiagonal matrix system requiring substantial CPU time to solve. This difficulty can be avoided by applying the unconditionally stable Alternating Direction Implicit method (ADI), 49 developed by Peaceman and Rachford and Douglas in 1955. According to this scheme, the entire solution process at each time step is "split" into two portions, i.e., the first half of solution processes for kth column (radial direction), while the other half processes for the jth row (axial direction). With the ADI scheme, secondorder central differences are used to approximate the values of derivatives at each nodal point in equations (34), (35). The finite difference algorithm for those equations during the first half of each time step for the jth column are then IBP T + DP T + AP T n+1/2 n BPjjk jjkl + DPjjk jjk + APjjk jk+1 n+/2 PXjj,k k = 1, 2, ....kmid (329) and BWjjk jj,kl + DWjjkTjj,k + AWjjk jjk+1 n+1/2 Xj,k k = kmid+1,...kmax (330) where the subscript jj is used to emphasize the specific column currently to be computed. It can be seen that the set of difference equations now is in the tridiagonal form 50 since the righthandside terms in the equations (329), (3 30) contain only known values from the previous results and the boundary conditions. These values can be computed by PXjk = CPjjkTjjlk + EPjjkTjjk (331) and WX CW...T + EW. T FW.. T j,k = j,k jjl,k j,k jj,k jjjj+l,J (332) The computational algorithm is implemented column by column and the unknown value (Tj,k)n can then be solved by either an iterative or a direct method. In order to do this, equation (329) needs to be assembled into the following matrix form. PD1 PA1 PB2 PD2 PA2 PB3 PD3 PA3 * in PBkdlPDkd1 PAkd1 PBkd PDkd n+1/2 PXl PX2 PX3 Tkd1 PXkd1 j kd PXkd j = 1, 2, 3 .....jmax (333) jj FP,kjj+l,k n T1 T2 T3 51 Similarly, equation (330) can be assembled into the matrix form as seen in (334). WD WA n T n+1/2 n WDkd+1 WAkd+1 n kd+1 n+1/2 kd+1 WBkd+2 WDkd+2 WAkd+2 Tkd+2 WXkd+2 WBkm WDkm WAkm km WXkml kWBmax WDmax j Tkmax j WX kmax j = 1, 2, 3 .....jmax (334) It should be pointed out that the coefficients and/or the right hand side terms in the row marked with symbol "*" in those matrixes needs to be properly modified as well as the limit of either subscript j or subscript k according to the different boundary conditions in the various models considered. For example, along the axis, the symmetry condition requires that aT/ar = 0 and this could be accomplished numerically by equating the temperature To and T2 and by employing a new combining coefficient of PA1* = PA1 + PB1 rather than the original PA1 in the first matrix for evaluating the temperature distribution within the pipe. In addition, if the temperature distributions along certain parts of the outer surface of the wall are given, then the 52 limit of subscript k will be ended with k = kmaxl rather than k = kmax shown in equation (334). The same arguments are also applicable for another subscript j in the axial direction. By using a similar procedures for the second half of each time step, the finite difference algorithm for the kth row then becomes ICP kT + GPj + 4 +FP T+ 1n+1 p Tt1/2 [CPjkkjlkk k+ GPjkkTjkk + FPjkk+lkk n+1 +kk j = 1, 2, ........jmax (335) CWjkk jlkk + GWj,kkj lkk Wjkkj+ = kk / 2j,kk, j = 1, 2, ........jmax (336) where the right hand terms can be computed by pyn+1/2= CP T + HP Tm API T n+1/2 P,kk j,kk j,kkl j,kk j,kk j,kk j,kk+l1 (337) kk (CBjkkTjkk1 + HWjkk jkk Wjkk jkk+l 1 (338) Similarly, to emphasize that the current calculation is in the kth row, we use the symbol kk rather than k in above formula and the matrix form of equation (335) for the kth row can then be written as PG2 PF2 PC3 PG3 PF3 PC4 PG4 PF4 PC * *I n+1/2 PG. PF. 1 Pjml jm1 PCjm PG. jm Tjm1, T. jI y2 PY3 PY4 Pjml PX j jm . n+1/2 k k = 2, 3 .......kmid (339) and the matrix form of equation (336) becomes WG2 WF2 WC3 WG3 WF3 WC4 WG4 WF 4 4 4 4 * In+1/2f WC. WG. WF jm1 jm1 Wjm1 WC.jm DG 3m jm ] T. jm1 T. N+1 WY2 WY3 WY4 WX jm1 k Wjm k = kmid+l, kmid+2.......km (340) N+1/2 k 54 As indicated above, the symbol shows that the coefficients in that row as well as the righthandside terms need to be properly modified corresponding to the different boundary conditions. The matrix terms in equations (329), (330), (331) and (332), (335), (336), (337) and (338) can be computed by 1 APj ( P3 + P4)jk 13  BPjk ( P3 P4). CPj,k 1 P1 P2)jk DP k. = (1 + 2*P4)j k EPk = (1 2*P2)j k FPk= ( P1 + P2)jk GPjk = (1 + 2*P2)j k HPjk = (1 2*P4)jk j,kjk (341) and 1 AW = ( WS + P4)' j,k 2(jW + j,k 1 = (2W3 W4) jk 1 CWjk = ( 2 W1 W2)j DWjk = (1 + 2*W4) jk EWj = (1 2*P2)j k FWjk = ( 2 W1 + W2)j k GWj = (1 + 2*W2)j k HWjk = (1 2*W4).k (342) In order to rewrite the equation in a simple form we define the following linear operator : LY [P]j,k = [BP, DP, AP]j k LY [W]j,k = [BW, DW, AW]j,k (343) LX [P]j,k = [CP, = Iicp GP, FP]j k LX [W]j,k =[CW, GW, FWk]jk and column vectors (TY)j,k = (Tjk' Tjk' Tjk+l (344) (TX)j,k = (Tjlk Tk' Tj+l,k ) Equations (329) which are in the radial direction can be simply written as LY (P] j,k LY [W j,k n+1/2 (TY} j,k n+1/2 (TY} j,k j = 1, 2 .......... k = 1, 2, ......... = PX j,k = WX j,k (345) (346) 56 Equations (333) and (334), which are in the axial direction, then assume the following simple forms: n+1/2 LX [P] j,k n+1/2 LX [W] j,k n+1 (TX) j,k n+1 (TX) j,k n+1/2 = PY j,k n+1/2 = WY j,k (347) (348) j = 1, 2, .......... k = 1, 2, ....... where the right side terms can be estimated by PX j,k = (CP, EP, FP) j,k (TX) j,k j,k j,k = (CW, EW, FW) = (BP, HP, AP) (TX) j,k j,k (349) j,k {TY) j,k WY = (BW, HW, AW) (TY) j,k j,k j,k 57 It is seen, as the result of the ADI "splitting" procedure which has been employed in the algorithm associated with different boundary conditions for the various models, that only a tridiagonal system of linear algebraic equations needs to be solved (i.e., during step 1, the coupled tridiagonal matrix (333), (334) are solved for each jth column of the grid points, while during step 2, the coupled tridiagonal matrix (339) and (340) are then solved for each kth row of grid points). Once the periodic steady solution has been obtained, we can calculate both the tidal displacement and the heat flux at different locations within the pipe. The numerical technique used here to integrate equations (238), (239) and (241) for evaluating the tidal displacement and the heat flux at each specified location can be obtained either by the Trapezoidal rule (with end correction) or Simpson's rule. Both numerical approaches are essentially fourth order methods. Convergence Criteria As is known, once the calculation work has started, the timematching process will be in the loop forever unless a criterion can be derived that indicates when the goal of the current computing work has been reached and further solutionmatching processes do not produce significant increases in accuracy. Such a criterion depends on the purpose of the calculation. If one's goal is to study the 58 transient process, one can indicate a time limit to pause or quit the current computing work. However, if one seeks a final periodic state solution, one has to develop a criterion to test if the solution can be considered acceptable. It seems that the temperature is a good measure of the accuracy of the overall solution process so that the most efficient way is to apply the convergence test on the temperature rather than on the velocity. For the present purpose, two convergence criteria were alternatively used. The choice of the criterion depends on what was the main goal in the particular study case. If the main effort is to observe the distribution of the quasisteady temperature at any phase angle wt within a period, the testing is done by comparing the temperature residual, i.e., by inspecting the averaged temperature difference of each nodal point at thee same wt between adjacent periods. This can be written as SRe (T11 T22 )2 k )/2 Resl =  < el (JMAX)(KMAX) j = 1, 2 .........JMAX k = 1, 2, .......KMAX (350) where Tllj,k is the value of the temperature of node (j,k) at ot in the current period, T22j,k is the value of the ), iste au o h 59 temperature of node (j,k) at the same wt in the previous period, and EI is the convergence parameter. If the goal of the investigation is to examine the effective heat transfer, the convergence criterion is established by computing the residual Z ( 01 2 ) 1/ Res2 = E2 (351) JSEC j = 1,2,......JSEC where JSEC is the number of crosssections where the axial heat flux was examined, E2 is the convergence parameter, and 41, 02 are the crosssection averaged heat flux in the current period and previous period, respectively. The summation is carried out over all sections (JSEC) where the heat flux is computed. The value of the convergence parameter can be determined by a balance of the acceptable solution accuracy and the cost of CPU time. In the current study both el and e2 were taken between 0.0010.01 as the residual of the temperature or the heat fluxes for acceptable value. Grid Generation As is well known, the solution accuracy and efficiency in a large degree depends on the grid system used in the numerical calculation. A good grid is characterized by small skewness, high smoothness and capability of high 60 resolution in the large gradient regions in the physical plane. It has been shown that rapid changes of grid size and highly skewed grids can result in undesirable errors [25]. The success of a numerical simulation of a complex thermal fluid dynamics problem does depend strongly upon the grid system used in the computation. In the present study, the grid mesh was generated with the emphasis on the high resolution capability of its simple geometric boundaries. The technique used here is one of the algebraic schemes which cluster the grid lines near the region desired [34), namely, ASk = AS (1 + e) k1 (352) where AS0 is the minimum specified grid spacing next to the wall or to some inner interfaces within the tube. The parameter e is determined by a NewtonRaphson iteration process so that the sum of the above increments matches the known arc length between k = 1 and k = kmax. The grid networks used in Model 1 and Model 2 are uniformly distributed along the axial direction in both the fluid and the wall; however, in Model 3, most grid points are clustered near the central portion of the pipe so as to gain higher resolution in the region of interest, while along the radial direction a nonuniform grid was employed in all three models. The minimum specified grid spacing ASo used 61 depends on the boundary layer thickness 6, namely, the kinematic viscosity of working fluid and the oscillating frequency. In the present study, it was found that a good choice of this value is ASo = 0.056 for laminar flow cases with a total of about 1520 nodal points distributed along the radius. CHAPTER IV NUMERICAL RESULTS AND DISCUSSION The problem of timedependent enhanced heat conduction subjected to sinusoidal oscillations can now be solved numerically for the boundary conditions appropriate to a long capillary tube according to the various models described in previous chapters. The computational tasks fall into two categories: the first part is a numerical study concerned with the characteristics of the oscillatory pipe flow, which includes an investigation of the velocity profiles, the Lagrangian and the tidal displacement trajectories. The second part of this study includes a thorough investigation of the final periodic state temperature buildup process in oscillating pipe flow, the periodic temperature distribution in the pipe and wall, and a study of the relationship between the enhanced axial heat flux and the tidal displacement. It also includes an investigation of the tuning effect and a comparison of the enhanced axial heat transfer with the corresponding pure molecular axial heat conduction as well as the investigation of the influences of the pipe radius and the wall thickness on the enhanced axial heat transfer. 63 Properties of oscillating laminar pipe flow have been analytically discussed by S. Uchida (36] for several different Wormersley numbers, and these results offer a useful reference for comparing with the present numerical studies. However, in the area of thermal fields associated with enhanced thermal pumping there is little information available for comparison, except for some recent results of Kaviany (13]. The present computational work was carried out on the Vax11/750 computer in the Department of Aerospace Engineering, Mechanics, and Engineering Science, and on the Vax11/780 in the Center for Instructional and Research Computing Activities, University of Florida, Gainesville, Florida. Note that it is very timeconsuming to buildup a final periodic state, for example, if the grid size is 101x22 and the time steps are 2000 per period, it takes 7.5 minutes (CPU) with the VAX11/780 machine or almost 1 hour (CPU) with the VAX11/750 machine to run only one period. It usually takes 2030 periods to reach the final periodic state solution. The process of selecting model size (i.e., total nodal points) is a synthetic balance among the storage requirement, the solution accuracy, and the cost of CPU time for finding an acceptable solution. Such a selection allows each case to be solved with a minimum of expense in computing operation thereby making it possible to do the 64 large number of runs needed to obtain enough data points to plot curves of the desired parameters. The results discussed in this chapter are presented for the purpose of illustrating the effects to be encountered in working with the enhanced thermal pump. These results should be effective in gaining insight into some interesting features of this enhanced heat transfer process. Since the analytic approximations [15, 16, 18, 22] are effective in describing the flow and the heat transfer aspects and to give considerable insight into the problem, they are often used to compare with the numerical results and should give good agreement when applicable. Part 1. Oscillatory Pine Flow Features In order to better understand the mechanism of the enhanced thermal pump, it is necessary to examine the mechanical features of the oscillating pipe flow. As mentioned in previous chapters, the present interest in the enhanced thermal pumping is confined to an investigation of the central part of the slender pipe, thus the flow field can be well approached by a 1D timedependent laminar model which neglects the ends effects. Velocity Profiles Fig. 41 shows the numerically computed time dependent velocity profiles at different phase angles of the exciting pressure when Wormersley numbers are equal to 1, 10, 100, and 1000, respectively. It can be clearly seen 65 0 10 a 100 51 a =1000 Wt = ge s s Is z 20se 1e zIe I Z 0l 3W 33 336e X Fig 41 1D Velocity Profiles in Oscillating Flow for Wormersley Number a = 1, 10, 100 and 1000 66 that the velocity profile at a = 1 presents a quasi parabolic shape at any instant within a period and is in phase with the stimulating pressure gradient. However, at higher frequency, for example, a = 100 and 1000, the velocity profiles can be clearly separated into two regions: in the vicinity of the wall the flow shows a typical thin boundary layer, while in the region far away from the wall the fluid moves as if it were frictionless slug flow. In fact, within this core region the velocity distribution is independent of the distance from the wall [29]. It can be also seen that the phases of the velocity profiles at higher frequencies cases are shifted about w/2 with respect to the 7 1 7 . Fig 42 Velocity Profile (a=5) [29] (R1r) /6 0.2 0.1 0.5 V 0.5 V 2400 2100 1800 (R1r)/6 V 0.5 (R1r) /6 2700 300 3300 360 0.5 0 Fig 43 Magnified View of Velocity Profile Near Wall (Wormersley Number a=10, H20, 6=0.014cm) (Rlr) / T 00 0.1 I 68 stimulating pressure gradient. At the intermediate frequency case (a = 10 of Fig. 41), the slug flow boundary is not so evident, but one can still note a boundary layer near the wall. The same pattern of the velocity profiles associated with moderate Wormersley numbers can also be found in the reference [36] where the velocity profile at Wormersley number a less than 10 was presented (Fig. 42). In order to better see the variation of the timedependent velocity profile within the boundary layer, a set of closer view of the velocity profile with respect to Wormersley number at a = 10 during phase intervals ot = 300 is plotted in Fig. 43. It should be emphasized that the solutions shown here are under the assumption that the secondary velocity in the radial direction is negligible compared to the axial velocity component, namely, the nonlinear inertial terms are not considered in the governing equation (21). This is a reasonable approximation for moderate oscillatory frequency w except near the ends of the pipe. However, at high Wormersley numbers, care should also be taken to avoid violating another assumption, namely, that of incompressibility and the concomitant condition that the oscillating phase does not change between the tube ends. High Wormersley number can be obtained either by increasing the oscillating frequency w or by increasing the pipe diameter for a given fluid. An oscillating flow can be 69 considered incompressible if Axw/2 < 0.05C, where C is the speed of sound in the fluid. For mercury, C = 1360 m/sec, one requires a = 27.1, when Ax = 20 cm and R1 = 0.1 cm. To avoid a appreciable phase difference between the tube ends, one requires that L/C << 2n/w. Both the restrictions are met in the examples to be considered below. The Lagrangian Displacements An alternative interesting representation to the oscillatory velocity field are the Lagrangian displacements of the fluid elements at different radii within the pipe. They have been plotted in Figs. 44 and 45 at time intervals of wt = 300 for Wormersley number a = 0.1, 1.0 and 10. It is noted that since both pipe diameter and the working fluid were fixed in this test, Figs. 44 and 45 represent the relationship between the Lagrangian displacement and the oscillating frequency. The trajectories plotted in Figs. 44 and 45 have been normalized by A/w2, where A = 1/p dp/dxl is the amplitude of the sinusoidal pressure gradient as defined in equation (2 1) and w is the angular velocity. Suffice it here to point out that for the lower frequency case (for example a 0.1 and 1.0) the Lagrangian displacement trajectory shows a foreseeable parabolic pattern at any moment. Nevertheless, the essential distinction between low frequency oscillatory pipe flow and steady HagenPoiseuille flow is that in the former, the Lagrangian displacement trajectories as well as a 0. 1 27(fy a 1.0 = r/R1 .003 if = r/R1 Lagrangian Displacement for a = 0.1 and 1.0 Fig 44 AX  = r/R1 cat=2100 X udt Lagrangian Displacement at a = 10 Fig 45 72 the velocity profiles are periodic so that the fluid particles do not translate axially upon time averaging, while in the later case they will. For intermediate Wormersley number (a = 10), the trajectories of the Lagrangian displacement departs considerably from the standard parabolic shape. This phenomenon can be even more clearly seen in the a = 100 case. Evidently, the higher the oscillating frequency, the thinner the boundary layer (6 = J2T/w). Tidal Displacements Figs. 46 and 47 demonstrate time variation of the crosssection averaged dimensionless Lagrangian displacement at Wormersley numbers in the range from a = 0.1 to 50. The tidal displacement can be obtained by summing the absolute maxima and the absolute minima of these curves. The corresponding nondimensional tidal displacements with respect to Wormersley number from a = 0.1 to a = 12 are listed in table 41 below Table 41 Dimensionless Tidal Displacement at Different Wormersley Numbers a 0.1 1.0 2.0 3.0 4.0 5.0 AX 0.00246 0.14295 0.67303 1.20213 1.41160 1.49041 a 6.0 7.0 8.0 10.0 12.0 AX 1.56294 1.61879 1.66046 1.71900 1.76698 a 0. .05 0.0 .05 .10 Dimensionless Crosssection Averaged Displacement Versus Time a = 0.1 1.0 Xave 2. 1. 0. 1. Dimensionless Crosssection Averaged Displacement Versus Time a 2 50 .10 Fig 46 Xave .15 Fig 47 74 It is noticed that for very small Wormersley number cases (a < 0.5), the crosssection averaged displacement varies like a sinusoidal function with respect to time (the ordinate), and as Wormersley number increases, the cross section averaged displacements are no longer symmetric about the ordinate but rather favor positive values of x. If the Wormersley number is further increased, eventually, the crosssection averaged displacement almost entirely lies on the right half plane. A similar feature can also be seen in Figs. 44 and 45. In fact, this just again shows the existence of the phase shift in the oscillating pipe flow. When the Wormersley number is small the crosssection averaged displacements as well as the Lagrangian displacements show a sinusoidal variation with respect to time and this has w/2 phase lag with respect to the stimulating pressure gradient which is assumed to be a cosine function of wt with zero initial phase angle. However, for higher Wormersley number, the phase lags increase to almost r with respect to the phase of the exciting pressure gradient. Note that the Lagrangian displacements in Figs. 44 and 45 were computed by lining up all the fluid particles on the plane x = 0 at wt = 0, however, by using nonzero velocity at ut = 00 as shown in Fig. 41. This implies that if the phase of the exciting pressure gradient is taken as the base of the measurement, it is generally not possible to assure the same phase to the 75 velocity and Lagrangian displacement as well as the tidal displacement. For example, at phase of the stimulating pressure gradient wt = 00, the corresponding velocity phase may be w/2 and the Lagrangian displacement phase may be w. In fact, this is just as true for very high Wormersley numbers; however, the phase difference with respect to the exciting pressure gradient phase is less than the value shown above for low Wormersley number. One can well see that if the phase lags of the Lagrangian displacement or the tidal displacement is given, one can certainly redraw the diagram shown in Figs. 46, 47 by lining up the phase with itself, and then an exact symmetric pattern of the tidal displacement curve similar to that in the very small Wormersley number case can be obtained. Unfortunately, the phase lags are a function of the Wormersley number and they are not known in advance. In order to verify the numerical method used in this study and to check the ETP code developed, a comparison of the computed dimensionless tidal displacement to the one using the analytic equations (17) and (18) given by Kurzweg [18] for Wormersley number varying from 0.1 to 100 has been calculated and plotted in Fig. 48. The solid line shows the analytic solution obtained by use of equations (1 7) and (18), while the dashed line shows the results with the ETP code developed in this study. The agreement is quite good, particularly when the Wormersley is less than 76 3.0. However, at high Wormersley number the numerical solution shows a very slight deviation from the analytic solution. This deviation is believed due to an inaccurate numerical integration over the crosssection using relatively large time steps (our time steps per period in the calculation were between 10002500). A comparison with using 104 time steps per period for Wormersley number a = 10 was studied and shows some improvement. However, using such small time steps in the present investigation is beyond the capacity of the current VAX computer facility used. The numerical error becomes particularly serious as the oscillating frequency becomes large where the extremely thin boundary layer requires more grid nodal points to resolve the flow variables in the vicinity of the wall. Fig. 48 shows that as the Wormersley number gets large, the dimensionless tidal displacement tends to the limit of 2.0, which agrees with the limit of 1.0 in the analytical solution given by Kurzweg [18] for the reason that the normalization parameter used in [18) is twice as large as that in the present numerical simulation. Fig. 49 shows the required stimulating axial pressure gradient used in the present study for a pipe radius R1 = 0.1 cm and water (v = 0.01 cm2/sec) taken as the working medium versus the dimensional tidal displacement Ax in cm for various Wormersley numbers (namely, oscillating frequency). It is evident from these results that for fixed 77  Numerical Solution  Analytic Solution 1.5 4) 0 r410 0 . 0 H 10 100 101 1D2 Wormersley Number a Fig 48 Relationship Between Dimensionless Tidal Displacement AX and Wormersley Number a 0 1 1 I 1I tI I I i O O O 0  I I 1 0 ito 0) 0 44 01 0) 10 V .0) '4 (uoi) xv 4uaaovTTdsTG IgpTI x M 4P4a 0Q) 0C r $4 >:1 0) 0 *c tP QH CO *4 to S0 . hJ. .1 ( )u4 Irr 79 tidal displacement, the required axial pressure gradient in the large Wormersley number case is orders of magnitude higher than that in the small Wormersley number case. This may eventually put some constraint on the use of very high frequency in the enhanced heat transfer technique. Fortunately, to meet the tuning condition, the required Wormersley number in this case is of order 1. As already discussed in the introduction, in order to gain the benefit of axial heat transfer, the use of large tidal displacement is always preferred. However, such an increase must be limited by the requirement of no convective net mass transfer occurring between two reservoirs and may also be constrained by the ability of the device to withstand the increase in the exciting pressure gradient. Phase Lags We are now in the position to study the phase lags of the Lagrangian displacement in the oscillatory pipe flow. Fig. 410 shows the phase variations (in degrees) along radius for Wormersley number varying from a = 0.1 to 4. Some numerical results are also shown in Tables 42 and 43. All of the data shown in Fig. 410 and the tables have the phase angle measured relative to the exciting pressure gradient. Two features can be seen in Fig. 410; first, in the core portion, the phase lags are almost equal to wr/2 when the Wormersley number is small, while the lags are almost w when the Wormersley number is large, and second, 80 the phase lags vary along radius, especially in the boundary layer. It is such phase lags that allow the existing temperature gradient in the very thin boundary layer of the oscillating pipe flow to act as region of temporary heat storage. It absorbs heat when the temperature of the core Table 42 Nodal point K 1 2 3 4 5 6 7 8 9 10 Phase Lags Along Radius (Working Medium: H20, Ax = 10 cm) a = 0.1 q=r/R1 0.0000 0.1667 0.3333 0.5000 0.5833 0.6667 0.7500 0.8333 0.9167 1.0000 phase 90.63 90.63 90.63 90.63 90.63 90.63 90.63 90.63 90.63 0.00 qr=r/R1 0.0000 0.0366 0.0957 0.1913 0.2593 0.3458 0.4559 0.5958 0.7737 1.0000 phase 102.06 102.06 102.06 101.88 101.70 101.52 101.16 100.62 99.72 0.00 a = 2. q=r/R1 0.0000 0.1095 0.2164 0.2836 0.3620 0.4535 0.5603 0.6849 0.8303 1.0000 phase 134.28 134.10 133.38 132.66 131.76 130.50 128.70 126.00 122.58 0.00 81 Table 42  continued Nodal a = 3 a = 4 a = 7 point K q=r/R1 phase P7=r/R1 phase v,=r/R1 phase 1 0.0000 170.83 0.0000 179.82 0.0000 181.62 2 0.0917 170.47 0.1242 179.28 0.0970 181.62 3 0.1818 169.75 0.2370 178.02 0.1892 181.62 4 0.2703 168.31 0.3395 176.04 0.2769 181.44 5 0.3573 166.51 0.4324 173.70 0.3604 180.72 6 0.4427 164.00 0.5168 171.00 0.4397 179.28 7 0.5267 161.12 0.5935 168.12 0.5152 177.12 8 0.6091 157.88 0.6630 165.24 0.5869 174.06 9 0.6901 154.64 0.7262 162.18 0.6552 170.28 10 0.7696 150.69 0.7835 159.12 0.7201 165.78 11 0.8478 146.73 0.8356 156.06 0.7818 160.74 12 0.9246 142.42 0.8828 153.18 0.8406 155.16 13 1.0000 0.00 0.9257 150.30 0.8964 149.04 14 0.9646 147.60 0.9495 142.56 15 1.0000 0.00 1.0000 0.00 Table 42  continued Nodal Point K 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 a = 10 q7=r/R1 0.0000 0.0884 0.1740 0.2566 0.3364 0.4135 0.4880 0.5599 0.6295 0.6966 0.7616 0.8243 0.8849 0.9434 1.0000 phase 180.54 180.72 181.08 181.80 182.52 182.70 181.80 179.64 175.68 170.46 163.80 156.06 147.24 137.52 0.00 a = 14 a = 20 v,=r/Rl 0.0000 0.0834 0.1647 0.2440 0.3215 0.3970 0.4709 0.5428 0.6130 0.6815 0.7484 0.8136 0.8773 0.9394 1.0000 phase 179.46 179.10 180.18 184.50 190.25 192.41 189.53 183.06 174.07 164.00 153.20 141.70 130.19 118.32 0.00 ,q=r/R1 0.0000 0.1635 0.3032 0.4227 0.5247 0.6120 0.6865 0.7502 0.8046 0.8512 0.8909 0.9249 0.9540 0.9788 1.0000 phase 179.10 180.54 185.93 182.34 186.65 184.14 184.14 181.98 178.38 173.71 167.95 161.12 153.57 146.01 0.00 83 Table 43 Comparison of Phase Lags With Different Working Mediums (AX = 10) a = 1. Phase Lags Nodal Point K 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 q =r/R1 0.0000 0.0161 0.0366 0.0626 0.0957 0.1378 0.1912 0.2593 0.3458 0.4559 0.5957 0.7737 0.0000 a = 5. Water 102.06 102.06 102.06 102.06 102.06 101.88 101.88 101.70 101.52 101.16 100.62 99.72 0.00 Phase Lags Mercury 102.14 102.14 102.14 102.14 102.14 102.14 102.14 102.14 101.78 101.42 101.06 99.98 0.00 qi =r/R1 0.0000 0.0714 0.1429 0.2143 0.2857 0.3571 0.4286 0.5000 0.5714 0.6429 0.7143 0.7857 0.8571 0.9286 1.0000 Water 181.62 181.44 180.90 179.64 178.02 176.04 173.34 170.28 166.50 162.18 157.50 152.10 146.16 139.86 0.00 178.02 175.86 173.37 170.11 166.51 162.20 157.52 152.13 146.37 139.90 0.00 (Z Sa) 00 SQ) 4O4 0 4 > 0 l(l snlpea ssaouosuawjju I I I 0 0 a e as O O a a a O 85 slug flow is higher than that of the boundary layer, and it releases heat as the core temperature is relatively lower. This large temperature gradient enhanced by the existing velocity phase lags allows a large amount heat to be conductively transferred radially within a very short time and subsequently to be transferred axially by a convective coupling. Part 2. The Enhanced Heat Transfer Investigation We have examined some mechanical characteristics in oscillating pipe flow, and compared the computed solution of the velocity field with Uchida's solution. The results for the velocity profiles are in good agreement. Nevertheless, before a detailed examination of the thermal field, it is first necessary to test the current developed ETP code when applying the temperature equations (26),(27). It is seen that the energy equations strongly depend on the velocity distribution and its buildup process, so that one can use analytic periodic velocity state (Eqs. 14 and 15, with no buildup process), and the computed velocity (with buildup process) to verify the correctness of the resulting thermal variables. The enhanced heat flux is a function of both velocity and temperature (Eqs. 234 and 235) and was chosen for a comparison of the analytic and numerical results of the problem. Part A of table 44 shows the results of the computed enhanced axial heat flux as well as the axial conduction heat flux when using the analytic velocity 86 Table 44 The Comparison of Enhanced Heat Flux Using Numerical Velocity with Heat Flux Using Analytical Velocity (Model 3, WaterGlass, Pr = 7.03, a = 3) A. Heat Flux Using Analytic Velocity (cm) (w/cm2oK) 0.9839 0.0230 2.9519 0.1855 4.9194 0.3864 6.8868 0.5810 7.8847 0.6775 9.8390 0.8638 B. Heat Flux 1.0033 0.0234 3.0101 0.1870 5.0164 0.3899 7.0226 0.5790 8.0402 0.6746 10.0328 0.8597 * Model 3, WaterGlass * .......enhanced axial (w/cm 0K) 0.0149 0.0132 0.0098 0.0075 0.0067 0.0055 (w/cm K) 0.0811 0.0724 0.0541 0.0415 0.0369 0.0302 (w/cm4 K) 0.0238 0.0213 0.0159 0.0122 0.0109 0.0089 Using Computed Velocity 0.0149 0.0133 0.0097 0.0076 0.0068 0.0055 heat flux 0.0811 0.0728 0.0533 0.0420 0.0374 0.0302 0.0232 0.0206 0.0155 0.0117 0.0104 0.0085 * Of .....axial heat flux by conduction in fluid * Ow.....axial heat flux by conduction in pipe wall * P = O/Ax2 87 (Eq. 14) for a = 1, Pr = 7.03, pipe radius R1 = 0.1 cm and glass walls of thickness AR = R2 R1 = 0.05 cm (which is almost equal to the boundary layer thickness 6 = 0.047 cm for this case) with Model 3. There were 22 points distributed along the radius (13 points in fluid), and the smallest distance next to the wall is equal to 0.166, while along the axial direction 101 point were used (60 points in the central pipe section). Thirty period runs were considered in order to insure that the final periodic state was achieved. An exception was the Ax = 1 cm case where only 23 period runs were computed. Under the exact same conditions described above, a calculation using the numerical velocity with buildup process was tested. Part B of table 44 gives the axial heat flux results using the corresponding numerical generated velocity profile. A comparison of results shown in Table 44 is quite good. One sees that the disagreement for most terms is less than 0.001. Note for Ax of 10 cm, the ratio O/of is about 15 and gets larger as Ax is increased further. Periodic Temperature Buildup in Thermal Pumping Process The purpose of this facet of the present study was to examine the temperature field buildup pattern and to determine how long before the final periodic temperature state can be reached. The author also intends to verify whether or not there indeed exists a time averaged linear axial temperature distribution within the thermal pump at 