Final Report A STUDY OF FORCED CONVECTION BOILING UNDER REDUCED GRAVITY by Herman Merte, Jr. NASA Grant NAG3-589 Report No. UM-MEAM-91-12 Conducted under: National Aeronautics and Space Administration Lewis Research Center Cleveland, Ohio Grant NAG 3-589 September 1991

n VA~ u,

TABLE OF CONTENTS Abstract......................................... iv List of Figures..................................... v List of Appendices..................................viii 1. INTRODUCTION................................ 1 1.1 General Background............................ 3 1.2 Basic Mechanisms of Forced Convection Boiling............. 5 1.3 Effects of Gravity/Orientation in Forced Convection Boiling........ 10 2. HEAT TRANSFER MODELS AND CORRELATIONS............ 12 3. EXPERIMENTAL APPARATUS....................... 16 4. EXPERIMENTAL PROCEDURES...................... 24 4.1 Spacial Mean Gold Film Heater Surface Temperature........... 26 4.2 Metal Heater Spacial Mean Surface Temperature.............. 26 4.3 Bulk Liquid Temperature.......................... 27 4.4 System Pressure............................. 27 5. EXPERIMENTAL RESULTS..................... 28 5.1 Transient Boiling..................... 28 5.2 Quasi-Steady Boiling............................ 31 6. FUTURE WORK............................. 37 Figures................................... 38 References.................................... 87 Appendices.................................. 99 iii

ABSTRACT This report presents the results of activities conducted over the period 1/2/85 - 12/31/90, in which the study of forced convection boiling under reduced gravity was initiated. The study seeks to improve the understanding of the basic processes that constitute forced convection boiling by removing the buoyancy effects which may mask other phenomena. Specific objectives may also be expressed in terms of the following questions: (1) What effects, if any, will the removal of body forces to the lowest possible levels have on the forced convection boiling heat transfer processes in well-defined and meaningful circumstances? This includes those effects and processes associated with the nucleation or onset of boiling during the transient increase in heater surface temperature, as well as the heat transfer and vapor bubble behaviors with established or steady-state conditions. (2) If such effects are present, what are the boundaries of the relevant parameters such as heat flux, heater surface super-heat, fluid velocity, bulk subcooling, and geometric/orientation relationships within which such effects will be produced? A flow loop was designed and fabricated to permit operation at low velocities and various orientations in the test section. Flat heaters are used with flow parallel to the surface, and include both semi-transparent thin gold films on quartz and copper heaters. The gold films serve simultaneously as heaters and resistance thermometers, and permit visualization from behind the heater surface. Results are presented for both transient and quasi-steady nucleate boiling. The quasi-steady results do not include the full spectrum of orientations possible. The experimental data presented here are the results of the activities of the following Research Assistants: Dr. Jamie S. Ervin, Mr. Longhu Li, and Mr. Kevin M. Kirk. iv

LIST OF FIGURES Figure 1. Influence of buoyancy on temperature distribution in vicinity of heater surface with pool boiling. 2. Forces acting on a growing/collapsing vapor bubble attached to a wall. 3. Example of advantage in use of phase change to transport energy. 4. Early flow loop for the study of forced convection boiling at various orientations and at reduced gravity. 5. Schematic of flow loop system. 6. Test section for heater surfaces. 7. Representation of flow boiling on a flat plate heater. 8. Optical paths for two view photography. 9. Gold film heating surface as simultaneous resistance thermometer. 10. Demonstration of concept for transient heating with a thin film. 11. Thin gold film temperature with onset of heating in pool boiling at a/g = 1. 12. Thin gold film temperature with onset of heating in pool boiling at a/g - 10-5 for 5 seconds. 13. Heat transfer between two communicating semi-infinite solids with uniform plane heat source at the interface. 14. Copper heater surface for quasi-steady operation. 15. Pressure control bellows assembly. 16. Flow loop pressure level control system. 17. Schematic of loop flow control system. 18. Schematic of loop condenser/cooler control system. 19. Schematic of R- 113 degassing unit. 20. Experimental parameters for flow loop operation in earth gravity. 21. Flow gravity boiling heat transfer with various earth gravity orientations. 22. Thin film gold heater temperature in forced convection boiling with horizontal up orientation. v

23. Thin film gold heater temperature in forced convection boiling with vertical orientation. 24. Thin film gold heater temperature in forced convection boiling with horizontal down orientation. 25. Heater surface superheat at incipient boiling for surface Q-16 at 1.7 cm/s. 26. Heater surface superheat at incipient boiling for surface Q-13 at 1.7 cm/s. 27. Heater surface superheat at incipient boiling for surface Q-13 at 4.5 cm/s. 28. Heater surface superheat at incipient boiling for surface Q-16 at 4.5 cm/s. 29. Incipient boiling delay time at 4.5 cm/s for surface Q-13. 30. Incipient boiling delay time at 1.7 cm/s for surface Q-13. 31. Incipient boiling delay time at 4.5 cm/s for surface Q-16. 32. Incipient boiling delay time at 1.7 cm/s for surface Q-16. 33. Representative heat loss calibration - metal heater no. 2 - vertical. 34. Sample forced convection nucleate boiling results. Subcooled, vertical, low velocity, metal heater no. 1. Test no. 62289. 35. Sample forced convection nucleate boiling results. Subcooled, vertical, low velocity. Metal heater no. 2. Test no. 62289. 36. Influence of orientation with forced convection nucleate boiling. Metal heaters. ATsub = 23.2~F, V = 1.7 cm/s. 37. Influence of velocity with forced convection nucleate boiling and vertical up flow (- = 180~) Metal heaters. ATsub = 23.3~F. 38. Influence of velocity with forced convection nucleate boiling and horizontal orientation. Up (- = 90~) and down (- = 270~). Metal heaters. ATsub = 13.1~F. 39. Comparison of metal heater versus gold film quartz heater with otherwise identical conditions: vertical upflow (-0 = 180~), ATsub = 13.2~F; P = 16 psia; V = 4.3 cm/s. 40. Comparison of metal heater versus gold film quartz heater with otherwise identical conditions: Horizontal facing up (& = 90~) and down (-0 = 270~), ATsub = 13.1~F; P = 16 psia; V = 4.3 cm/s. 41. Hysteresis effects with quartz heater at -a = 225~. Test No. 112090. vi

42. Influence of Subcooling for -0 = 0~. Test No. 022591 43. Influence of Subcooling for - = 45~. Test No. 022591. 44. Influence of Subcooling for -a = 90~. Test No. 022591. 45. Influence of Subcooling for -0 = 135~. Test No. 022591. 46. Influence of Subcooling for -& = 180~. Test No. 022591. 47. Influence of Orientation with Low Subcooling. Low Heat Flux. Test No. 112090. 48. Influence of Orientation with Low Subcooling. Test No. 022591. 49. Influence of Orientation on steady flow nucleate boiling. High Subcooling. vii

LIST OF APPENDICES Appendix A. Finite Difference Solution to Mixed Convection in a Rectangular Duct. B. The Velocity Potential of Two Spheres Moving Perpendicular to the Line Joining Their Centers and the Associated Lift Force C. Uncertainty Analysis D. Tabulation of Forced Convection Transient Incipient Boiling Data vii

1. INTRODUCTION This report presents the results of activities conducted over the period 1/2/85 - 12/31/90, in which the study of forced convection boiling under reduced gravity was initiated. The study seeks to improve the understanding of the basic processes that constitute forced convection boiling by removing the buoyancy effects which may mask other phenomena. Specific objectives may also be expressed in terms of the following questions: (1) What effects, if any, will the removal of body forces to the lowest possible levels have on the forced convection boiling heat transfer processes in well-defined and meaningful circumstances? This includes those effects and processes associated with the nucleation or onset of boiling during the transient increase in heater surface temperature, as well as the heat transfer and vapor bubble behaviors with established or steady-state conditions. (2) If such effects are present, what are the boundaries of the relevant parameters such as heat flux, heater surface superheat, fluid velocity, bulk subcooling, and geometric/orientation relationships within which such effects will be produced? The potential effects, implied above, have their roots in observations of nucleate pool boiling under variable gravity perpendicular to the heating surface from high gravity to microgravity to negative gravity [16]. It had been observed that under high gravity conditions the nucleate boiling process is degraded; that is, for a give constant heat flux, the driving potential (heater surface superheat) is increased. For reduced and negative gravity conditions the nucleate boiling process is enhanced; that is, for a given constant heat flux, the driving potential (heater surface superheat) is decreased. These are illustrated schematically in Figure 1, where the observed temperature distributions with pool boiling in a saturated liquid are qualitatively presented. The variation of the buoyancy has an influence not only on the heater surface temperature, but on the boundary layer as well. The research, whose initial results are presented here, involves the determination of the 1

influence of an imposed velocity parallel to the heating surface on the bubble dynamics and on the resulting heater surface temperature and liquid temperature distribution. The enhancement of nucleate pool boiling with reduced gravity is believed to be due to the influence of buoyancy on the size and thickness of the microlayer trapped under a growing vapor bubble. Any residual influence of buoyancy on nucleate boiling in the presence of an imposed bulk liquid velocity, say parallel to the heater surface, depends on the extent to which the microlayer would be affected by the combination of buoyancy and forced convection. This is governed in turn by the various forces acting on the vapor bubble as the dynamic evaporation/condensation processes are taking place. Figure 2 shows the various forces acting on a vapor bubble in nucleate boiling. With forced convection of the bulk liquid parallel to the heater surface two forces are acting in addition to those involved in pool boiling: the drag or liquid shear on and around the bubble as a result of the bulk liquid motion, and the lift generated by the liquid velocity change as it moves around the bubble. A total research program may be subdivided into three sequential phases, each intended to provide the base for the next phase: Phase A: This consists of testing in the laboratory with a flow loop at variable orientations between the flow direction and the gravity vector, using variables of flow rate, heat flux and subcooling with sizes, construction, and orientation of heating surfaces to serve as preliminary models for Phase B. It is expected that such testing will result in data that will provide guidance for the testing program at reduced gravity. Phase B: This would involve tests at reduced gravity using aircraft flying parabolic trajectories, with a portable flow loop and instrumentation resulting from the developments in Phase A. This would serve to determine the parameter boundaries for Phase C more accurately and economically. Phase C: This would involve orbital space flight testing in the shuttle or equivalent vehicle, with long term microgravity and well-defined parameters. The long terms 2

available will permit attaining a uniformity of experimental conditions between the various specific tests. It is foreseen that the test package used, resulting from Phases A and B, would be compact, self contained, and virtually completely automatic. The activities described below include only results from Phase A to date. Although a major emphasis is placed here on experimental measurements to observe the behavior of forced convection boiling under microgravity, appropriate analytical activities are an integral component of the total research program. 1.1 General Background Phase changes with or without forced convection can provide high or low heat transfer rates, depending on the mode, i.e., nucleate or film boiling, evaporation, film or dropwise condensation. The mode is governed by a number of factors such as the degree of superheat or sub-cooling present in the liquid, vapor and container walls, the fluid/solid properties, body forces present, fluid velocities, and system geometry/orientation. Requirements for the proper functioning of equipment and personnel in the space environment of reduced gravity and vacuum introduce unique problems in temperature control, power generation, energy dissipation, the storage, transfer, control and conditioning of fluids, (including cryogenic liquids), and liquid-vapor separation. Boiling in microgravity is fundamentally different from boiling in earth gravity: the buoyancy force which induces liquid and vapor motion in boiling with earth gravity is effectively eliminated in microgravity. Temperature control in certain locations where internal heat generation takes place, either as a result of dissipation or because of the nature of the process, may require that this energy be transported to other locations of the facility, where it can be eliminated by radiation to space. Fig. 3 illustrates two advantages in the use of phase change for the transport of energy in space; not only is the pumping power reduced by a 3

factor of 484 for the conditions assumed, but the mass of the fluid required is reduced by one-half. Certain effects which can be neglected at normal earth gravity, such as surface tension and vapor momentum, can become quite significant at microgravity conditions. Examples of applications in which these effects must be considered are: ullage control in storage containers; mechanisms acting in heat pipes; the effective transfer/flow of saturated or near saturated liquids from one vessel to another. The latter is a particular problem with cryogenic liquids, where the transfer lines must be chilled, resulting in vapor production and two-phase flow. A phenomenon very similar to this in the mechanisms involved, and which gives rise to the fundamental study outlined in this proposal, is the process of forced convection boiling heat transfer under reduced gravity conditions. Applications in space stations are being considered [1], whether for temperature control or for vapor generation itself, as having distinct advantages over passive heat pipes in certain circumstances. Experiments of two-phase heat transfer conducted at earth gravity are generally either in the horizontal or vertical orientation. With the horizontal case separation of phases occurs due to gravity which, together with interfacial shear, can give rise to severe surging or chugging. This behavior may be quite different with significantly reduced gravity conditions. In the vertical orientation, either with upflow or downflow, the body forces accelerating or decelerating the vapor phase relative to the liquid will likewise produce behaviors quite different than in a microgravity environment. Very small temperature differences, whether superheat or subcooling, which may normally be of little importance, can produce significant effects when the processes are diffusion limited, as will be the case under microgravity conditions. Their influences must be understood, anticipated, and given appropriate consideration. Such small temperature differences can arise in large containers subjected to solar heating, for example, even with multi-layer insulation installed, where insulation penetrations are necessary for supports and fill lines. Additional heating can occur with connections to heated engine components. 4

Small temperature differences can also cancel the effectiveness of surface tension control devices such as screens. The effective application of forced convection boiling heat transfer in the microgravity environment of space, then; requires a well-grounded and cogent understanding of the mechanisms involved. Many correlations for a/g = 1 are presently available in the literature and are continuing to appear [e; g., 2-11]. This list is by no means exhaustive, and the mere existence of such a large number means that each has its inadequacies and limitations. Any attempt to. extend these correlations to microgravity conditions, or to modify them using experimental results obtained without adequate consideration of mechanisms of a fundamental nature, as has been proposed [12], will provide results of limited utility. A discussion of the mechanisms in forced convection boiling heat transfer anticipated to be influenced by changes from earth gravity to microgravity will be presented below. A review of early works on pool boiling under reduced gravity is available [13], along with more recent data [14, 15, 16]. 1.2 Basic Mechanisms of Forced Convection Boiling a. Nucleation. The onset of boiling is inherently a transient process in the sense that once having begun, the dynamics of the boiling process will so change the situation that it can only be repeated by beginning anew. A number of studies of nucleation and the inception of boiling under non-forced convection circumstances with reduced gravity have been reported [17 - 20]. The condition at which nucleation takes place essentially depends on the microgeometry of the solid surface, the solid/fluid properties, the surface temperature of the solid, and the temperature distribution in the liquid. The latter two parameters in turn depend upon the imposed heat flux, whether saturated or subcooled conditions prevail, the velocity distribution in the vicinity of the heater surface, and the magnitude of the net buoyancy forces relative to momentum effects associated with the 5

velocity. The ratio of the latter forces constitute the Richardson Number, and is expected to be one of the parameters necessary to quantify the role that microgravity plays in forced convection boiling. The Richardson Number can also be expressed as the ratio of the Grashof Number to the square of the Reynolds Number (Gr/Re2). For Ri>> 1 natural convection dominates, while for Ri<<l natural convection effects can be neglected [113], and for Ri-l the full flow must be considered, and is referred to as mixed convection. In the absence of fluid turbulence this is amenable to computation. An example is included in Appendix A for the geometry used here. It is useful to define a corresponding "two-phase" Richardson Number by replacing the Grashof Number by the Archimedes Number, which is also a ratio of buoyancy to viscous effects, except buoyancy is now in terms of a finite density difference Ap: Ar = x (1) v2 p A "Two-Phase" Richardson Number thus provides a measure of the buoyancy versus flow forces in two-phase flow: Ri(2) = pV2 (2) Results of research on nucleation at standard earth gravity have been reported for both pool boiling [21, 22] and forced convection [23, 24] conditions, while the thickness of the thermal layer at the initiation of nucleate pool boiling has been measured [25]. Once boiling has initiated and reached a steady condition, the nucleation site density becomes an important parameter in the description of the boiling. A reasonable amount of measurements of nucleation site density have been reported for pool boiling [26 - 29], but only one work is know for forced convection boiling [30]. Measurements of the tempera 6

ture distribution in the boundary layer, which influences the nucleation site density, are available for pool boiling only [31, 32], and the interactions taking place between adjacent nucleating sites with pool boiling has been investigated [33]. The nucleation site density is also expected to be dependent on the departure size of the bubbles, and the factors which govern this will be considered below. b. Growth/Collapse. Once a particular nucleating site has become activated, the subsequent rate of growth and possible later departure and/or collapse are dependent on the transient temperature distribution in the vicinity of the bubble interface. The rate of growth affects the bubble frequency and together with the nucleating site density governs the relationship between heating surface superheat and the heat flux for pool boiling [32]. The rate of growth will be influenced by forced convection and reduced gravity only insofar as the temperature distribution is affected. Considerable work, both analytical and experimental, has been reported on the dynamics of vapor bubbles in the liquid bulk and near solid walls with pool boiling [e. g., 34 - 46]. The collapse of cavitation bubbles are mechanistically the same as boiling bubbles [47 - 50], with large collapse rates associated with large subcoolings and large temperature gradients in the immediate vicinities of the bubbles. However, surface tension effects are neglected relative to the dynamic effects. With the slow velocities expected to be utilized with forced convection under the microgravity conditions of space, for energy conservation, it is not anticipated that the large liquid momentum associated with large collapse rates will be present, which result in cavitation damage. Instead, it is intuited that the growth and collapse rates will be relatively small, although still significant in their influence on the heat transfer rates, because the vapor formation arises from the relatively small liquid superheats. In this case the influences of surface tension may very well play a significant role in the heat transfer from the solid surface on which boiling is taking place. Again, a considerable literature exists on this effect with pool boiling [51 - 59], but none with the addition of forced convection. A factor in addition to surface tension which may become significant in the absence of body 7

forces is the momentum effect associated with the density changes of phase change [60, 61]. This can influence the departure size of the vapor bubble as well as its subsequent trajectory, which will be discussed below. One further facet of vapor bubble nucleation and growth as influenced by surface tension should be mentioned here. The superheat that the liquid acquires in the boundary layer adjacent to the heater surface can be considerable, prior to nucleation. It is thus possible for the vapor formed initially to completely envelope the heater surface. With certain configurations such as small wires or cylinders it is possible that subsequent surface tension effects will maintain a stable "pseudo" film boiling process only because of the particular geometry used. It is expected that even if film boiling becomes suppressed to nucleate boiling on a small wire or cylinder, the thermophoretic effects and the resulting heat transfer will be quite different than with flat surfaces. Observations made that pool nucleate boiling is uninfluenced by changes from earth gravity to microgravity [62] are believed due to the large surface tension effects associated with the fine wire used, so that buoyancy is relatively unimportant in either case. The possibility of such effects, together with the fact that a flat surface provides a more well defined orientation for buoyancy and forced convection purposes provides the motivation for using a flat heating surface in the initial studies here. c. Departure Size and Trajectory. The size and trajectory of the vapor bubbles upon departure following growth in the vicinity of the walls will be important factors in establishing the flow pattern taking place in the bulk fluid stream, which can influence the subsequent heat transfer processes taking place as well as the pressure drop. For pool boiling the forces which play a role in the departure process are surface tension, buoyancy, inertia, pressure difference between the inside and outside of the bubble, and the thermophoretic forces resulting from surface tension gradients [63 - 68]. A possible source of error in assessing the departure of vapor bubbles from a solid surface has been pointed out recently [69]. With forced convection, additional forces affecting the departure are 8

shear stresses and lift associated with circulation around the bubble, because of the velocity gradient in the flow field. The various forces acting are illustrated in Fig. 2. An analysis of the lift forces conducted for the case of potential flow is included here as Appendix B. Both analytical and experimental works have been reported [70 - 72], which include the possibility for sliding rather than departing, and with limiting effects taking place at very low velocities. Buoyancy was always present in these experimental works, of course, and a limitation exists in extending such measurements to behavior in microgravity conditions. The trajectory followed by a bubble following departure depends on the dynamics of the growth process, which will be influenced by the degree and distribution of subcooling in the flow stream, along with the fluid velocity gradient. The description of the motion is complicated by virtual mass effects [73 - 76], and by interactions between bubbles and solid walls in the presence of temperature and velocity gradients [77 - 81]. It can be expected that the absence of buoyancy in microgravity will have a significant influence on these interactions. d. Pressure Drop. An extremely rich literature deals with the matter of pressure drop prediction in two-phase flow [e.g., 82 - 91]. In microgravity conditions the pressure drop will be due solely to viscous and to acceleration effects associated with the quality changes with boiling. However, each of these will be dependent upon the size and spacial distributions of the vapor bubbles. As pointed out earlier, these depend upon the departure size and the subsequent growth and trajectories of the vapor bubbles which depend, in turn, on the temperature and velocity gradients. These are expected to be influenced by the removal of gravity forces, for the cases of low velocities to be used. 9

1.3 Effects of Gravity/Orientation in Forced Convection Boiling Little work has been done to investigate the effect of gravity in forced convection boiling, since the liquid momentum is normally assumed to be predominant over buoyancy effects. Some results available are described here. Bubble growth and motion of a single bubble with no gravity has been investigated in an isothermal and superheated fluid with a velocity field present [92], but the test time for these experiments were less than two seconds, and the study ignored the effect of a temperature distribution in the fluid as well as the interaction with other nucleation sites. The relative influence of velocity and buoyancy on the critical heat flux of liquid nitrogen was investigated [93]. The data presented was separated into buoyancy-dependent and buoyancy-independent zones. The inlet velocity required to prevent buoyancy from influencing the critical heat flux was found to be a function of the pressure and subcooling. This result is consistent with expectations, since the magnitude of the buoyancy is directly related to the volume of vapor present which is a function of the system pressure and subcooling. Transient slightly subcooled forced convection boiling experiments were conducted in a drop tower to simulate zero gravity [94]. The liquid velocities used were of the same order-of-magnitude as free convection velocities at earth gravity. It was observed that at zero gravity the vast majority of bubbles remained attached to the surface forming what was called a "bubble boundary layer." This phenomenon was peculiar to zero gravity, since bubbles always separated at earth gravity where the free convection velocities were nearly the same as the liquid velocities in the forced convection boiling. A correlation for the size of bubbles was obtained from a thermal equilibrium analysis and found to be a function of the saturation layer thickness. A few investigators have considered the effect of flow direction on forced convection boiling. The effect of surface orientation on bubble frequencies in nucleate pool boiling of 10

R- 1 was investigated [95], and it was demonstrated that the bubble frequencies increase by increasing the orientation angle. Boiling nitrogen, was studied with upward and downflow [96], and the higher accumulation of void observed in downward flow suggested that different heat transfer coefficients may occur in upflow and downflow. A qualitative comparison of upflow/downflow heat transfer of R- 113 was provided for a Reynold's number range from 1 x 104 to 5 x 104 [97]. These data show that the heat transfer coefficient for upflow is significantly greater than that for downflow in subcooled boiling. A corresponding but smaller difference also exists for saturated boiling. This difference between upflow and downflow is probably due to the fact that in upflow the buoyant and drag forces on a bubble prior to detaching from the surface are additive, while for downflow they are in opposite directions. Thus, it can be expected that bubbles detach from the surface at smaller diameters in upflow than for downflow, resulting in greater microconvection effect and enhanced heat transfer in upflow. In contrast to this work, Ref. [98] reported that for fully developed nucleate boiling with flow velocities of 0.2 and 0.8 m/s the flow direction has a negligible effect on the heat transfer coefficient, even though the Reynold's numbers were lower than those of Ref. [97]. Photographs during boiling at low velocities show a higher frequency of vapor bubble formation in upflow compared to downflow, which was consistent with the observations of Ref. [95]. No obvious effect of flow direction on the heat transfer coefficient for fully developed nucleate boiling was observed in Ref. [99], where the entering liquid Reynold's number was 1.4 x 104. 11

2. HEAT TRANSFER MODELS AND CORRELATIONS Current theory on the mechanisms of heat transfer in nucleate boiling asserts that there is a component of heat transfer due to the rapid flow of heat through a liquid microlayer between the growing vapor bubble and the solid surface as well as a bubbleinduced component due to the enhanced transient conduction which occurs as a result of the pumping action of the vapor bubble. In addition, with forced convection boiling there is also a contribution through single phase convection. There is considerable debate as to the importance of latent heat transport relative to the other mechanisms of heat transfer. Gunther and Kreith [100], using measurements from a photographic study of forced convection boiling, estimated that latent heat transfer accounted for only a small fraction (1 -2%) of the total heat flow. Clark and Rohsenow [101] had similar findings. However, Bankoff and Mikesell [102] argue that if turbulent convective heat transport dominates in the heat flow of condensing bubble surfaces, latent heat may indeed be significant. Zuber [103] postulated that the mechanisms of heat transfer in nucleate boiling differ according to the different two-phase flow regimes. At low heat fluxes the vapor bubbles can be considered as isolated bubbles with no interference from either its predecessor or neighboring bubbles. In this regime, the heat transfer models based upon "bubble agitation" or "bubble pumping" give reasonable results. However, these models give incorrect results at moderate and high heat fluxes where bubble interference does occur. In this region of interference, the dominant heat transfer mechanism is by latent heat transport, according to Zuber. Experimental work by Bilicki [11] also supports the hypothesis that latent heat transport is significant in forced convection boiling. Bilicki argued that if heat transfer in nucleate boiling were enhanced by the bubble-induced turbulence, then the frictional pressure drop and the heat transfer coefficient should change simultaneously in accordance with the analogy between momentum and heat transfer. Yet, it was noted that the pressure 12

drop remained nearly constant during the onset of nucleate boiling, indicating that the bubbles do not act as a stirring device but as thermal sinks for the transport of latent heat. In a heat transfer model given by Del Valle and Kenning [104], an area of influence surrounding each nucleation site is considered, which is repeatedly quenched at the bubble frequency by liquid at the bulk temperature, and the heat transfer directly under the bubbles is modified by microlayer evaporation. In addition, the heater surface area between the areas of influence is assumed to be cooled by single phase convection. It was concluded that the enhanced transient conduction induced by bubble motion was by far the most important mechanism in heat transfer while microlayer evaporation accounted for only 2 -3%, and single phase convection accounted for 5-10% of the total heat transfer. However, there is doubt as to whether the transient conduction model used accurately describes the quenching process, since the model assumes that the bubble-induced convection is sufficiently strong to produce instantaneous replacement of liquid in the quenching process. It has been assumed that the fully developed region of the boiling curve in forced convection boiling coincides with the extrapolation of the pool boiling curve. This assumption appears reasonable, at least for low velocities. However, experiments carried out by Bergles and Roshenow [23] demonstrate that the boiling curve for forced convection boiling is not a simple extrapolation of the pool boiling curve. Their forced convection boiling data does merge into an asymptote at large values of superheat, indicating the invariance of heat flux to velocity and subcooling in fully developed forced convection boiling. However, the slope and intercept of the asymptote differs from their pool boiling curve. Lemmert and Chawla [105] also noted that there was no significant effect of flow velocity or subcooling on the boiling heat transfer coefficient in fully developed forced convection boiling. Owing to the complexity of forced convection boiling, theoretical analyses cannot provide a general equation for boiling heat transfer coefficients for different substances and conditions. Hence, heat transfer calculations require use of empirical correlations, where 13

the correlations will have a relatively large uncertainty and can be used only for restricted cases. In the calculation of heat transfer coefficients, information on the operating conditions, the fluid properties, and the geometry are usually necessary. Chen [106] suggested that the heat transfer coefficient for flow boiling can be expressed as the sum of the heat transfer coefficients of forced convection, hfc, and pool boiling, hb, h = F(hfc) + S(hb) (3) where the forced convection heat transfer is intensified by the factor F (F>1) and the boiling heat transfer is suppressed by the factor S (S<1). In addition to knowledge of the convective heat flux, qfc, and nucleate boiling heat flux, qb, some correlations also require knowledge of the conditions required for the inception of nucleation. Bjorge et al. [107] proposed: q = Jqfc2 + (qb - qbi)2 (4) for subcooled boiling where qbi is the flux on the curve at a predicted superheat for the inception of nucleation, which does not depend on the available cavity sizes. A different superposition scheme proposed in this work obtains better agreement for subcooled boiling: q = qfc 2 + qbi2 1 - { ATsatib/ATsat} 3)2 (5) where the recommended equations for the convective heat flux and boiling heat flux are given in Refs. [108, 109], respectively. The incipient boiling criterion for subcooled 14

conditions is derived by Bjorge [110]. The only empirically determined coefficient needed in the correlation is in the calculation of the nucleate boiling heat flux. Other correlations by Shah [111] and by Gungor and Winterton [112] consist of only the forced convection term, where the boiling effect is included in enhancement factors which are functions of the boiling number, Bo. 15

3. EXPERIMENTAL APPARATUS A small scale forced convection flow loop was designed and fabricated, occupying a total volume of 97 x 75 x 61 cm (38 x 30 x 24 inches), for use with R-1 13 over a range of temperatures from 25~C (72~F) ambient to 600C (140~F) with corresponding saturation pressure variations. An early version of the loop is shown in Fig. 4. As will be described in more detail below, the geometry of the conduit leading to the straightener section and of the preheater and condenser/cooler assemblies were changed. A schematic representation is given in Fig. 5, and provision is made to permit rotation of the assembly through almost 3600 relative to the gravity field vector while under continuous operation. To compensate for changes such as hydrostatic pressure taking place during rotation, the system pressure and temperature at the entrance to the test section are automatically controlled, as is the flow rate. To prevent cavitation at the pump inlet while operating near the saturated liquid state in the test section, heat exchangers are included for subcooling the liquid prior to pumping and flow measurement, followed by heating. The 12VDC centrifugal pump is capable of control over a 10:1 volume flow rate by the pump speed, using the output of a propellertype flowmeter with a microprocessor to control the DC voltage. The outlet of the pump leads to a "T" section, one branch connecting to the pressure control and filling systems while the other branch leads to the heating system. The preheaters in the loop raise the temperature of the R-1 13 to the desired operating temperature. Referring to Fig. 5, subcooled R-1 13 leaving the condenser/cooler is pumped into preheater #1, a counter-flow heat exchanger. The heat exchanger raises the R- 113 temperature to within about 40C of the set point temperature. Preheater #2 then heats the R113 to the desired temperature level at the inlet of the test section, as indicated by a resistance thermometer at that location. Preheater #1 consists of a one pass multiple tube heat exchanger, with 35 stainless steel tubes, 0.953 cm OD x 0.699 cm ID x 33.02 cm length (3/8 in OD x 0.275 in ID x 13 inches). Preheater #2 is identical to each of the three heat exchangers used as the 16

condenser/cooler, with these four (4) heat exchangers being similar to Preheater #1 except that each contains only 13 tubes instead of 35. A 2 kw a.c. heating unit in series with an Emerson (Model SA55CXJAR-4814) 1/3 hp pump provides hot water to the shell side of the heat exchangers with R1 13 flowing within the tubes. A temperature controller maintains the desired hot water temperature by sensing the temperature of the water exiting the heater, and the R-113 outlet temperature from preheater #2 is controlled manually by a flow valve in the series with the pump. It is planned that this will also be automated in the future. The test section is shown in Fig. 6 and consists of a rectangular flow channel 4-1/8" wide, 14" long, with four different possible heights (1/8", 1/4", 1//2", 1"). Only test sections with heights of 1/8", 1/2" and 1" have been fabricated to date, which permit varying the bulk flow velocity by a total factor of over thirty (30) with the existing pump. However, the test section with a height of 1/8" results in difficulties in observation from the side, and its utility may be limited. The maximum Reynold's possible at present with R113 in the 1" test section is about 5000. For this initial basic study, the rectangular flow channel in Fig. 6 provides for the use of flat heater surfaces, eliminating complications of surface tension effects associated with curved surfaces, and also provides a more well-defined flow field in the vicinity of the heater surface than would be possible were tubing or cylinders used. Additionally, the orientation between the surface and the gravity vector is more well-defined regardless of whether the gravity field is a residual one in space or earth gravity. The configuration shown in Fig. 6 permits the simultaneous use of three (3) pairs of identical heaters if the study of upstream or downstream interactions between boiling surfaces is desired, or the use of three (3) different pairs of heaters, either simultaneously or independently, without requiring the draining, disassembly and refilling of the loop system to change surfaces. This has been found to be a time consuming process. The use of heaters in opposing pairs within the channel permits simultaneous operation with 17

opposite orientation relative to the gravity vector and, where desired in the future, the study of the interaction between the boiling boundary layers. Further, the availability of three opposed pairs of openings in the test section permits the insertion of diagnostic devices into the flow stream if desired, such as local fluid velocity probes (Hot wire, LDV) or temperature traverses. The duct is a welded assembly of 1/2" aluminum plates with the various quartz windows and heater substrates held in place by bolted flanges with appropriate gaskets. A 10" long upstream section provides for smoothing and transition of the duct dimensions from the preceding 900 turn in the loop. Measurements of the velocity profiles in both the 1" wide and 1/8" wide test sections were made, using the hydrogen bubble technique in water. At the relatively low velocities used, in Reynold's number ranges from 300 to 1700, a maximum pump efficiency of 10% was measured, which is typical for these small pump sizes. Additionally, it was found that the velocity profile in the test section was quite unsymmetrical, owing to the attempts to make the loop compact in length. After a number of experimental trials involving relocation of the flowmeter, redesign of the flow turning section, which included the installation of flow turning vanes preceded by a large cross section flow duct filled with small chrome plated metal spheres to break up the liquid jet issuing from the tubing used, an acceptably symmetrical velocity distribution was obtained. The heated surface itself is rectangular, 19.1 mm x 38.1 mm (3/4" x 1-1/2"), in size. This size is large relative to the maximum anticipated bubble sizes, again to eliminate any geometrical dependency, and is small enough in an absolute sense so as to keep electrical heater power requirements within a manageable level. A representation of the flat heater surface relative to the flow is shown in Fig, 7. The heated surface is mounted on a circular substrate which can be rotated in its own plane, so that the heated length or aspect ratio in the flow direction can be changed conveniently by a factor or two. This will have the effect of changing the thermal boundary layer thickness for given levels of flow velocity and heat flux. Two different types of test surfaces are presently used, mounted as opposing pairs to 18

permit simultaneous or consecutive operation at a/g = +1 and a/g = -1. The first pair of heaters in the flow direction of Fig. 6 consists of semi-transparent gold coatings (400 A) on quartz to provide both steady and transient nucleate boiling, using the surface coating simultaneously as a resistance heater and as a resistance thermometer, and permit simultaneous high speed photography of the process from the side and through the heating surface as illustrated in Fig. 8. This permits obtaining data on the departure size and trajectory of the bubbles, along with the nucleation site density and frequency of bubble departures. A sketch of this heater is given in Fig. 9, and shows the means by which the current carrying and the potential lead electrical connections are carried through the quartz surface without introducing any impediments to the fluid flow. By passing a D.C. current through the thin film and measuring the voltage drop and current with sufficient accuracy the heat flux input and the mean instantaneous surface temperature can be measured, following appropriate calibration of the electrical resistance versus temperature. Accuracies of +1~F have been attained with reasonable precautions, and its reliability and durability have been thoroughly tested. The heater can be used in several ways, referring to Fig. 10. If the fluid remains motionless, the temperature distribution in both the quartz substrate and the fluid can be computed from the measured power input, using classical techniques, and the measured surface temperature can be compared with the computed surface temperature. This is demonstrated in Figures 11 and 12 for a constant and uniform heat flux input in pool boiling, with the measurements following the 1-Dimensional computations up to the onset of natural convection at a/g = 1 in Fig. 11 and up to the onset of boiling at a/g = 10-5 in Fig. 12. Once the fluid has been set in motion, whether by natural or forced convection or by boiling, or by both, the transient measurements of the thin gold film temperature and the power input as in the later time intervals in Figures 11 and 12 permit the computation of the mean heat flux to or from the substrate, and hence to the fluid. A limitation of the temperature measurement with the heater surface in Fig. 9 is that only the integrated mean 19

surface temperature is measured. With forced convection over the surface and a uniform imposed heat flux it can be anticipated that a temperature variation will arise in the flow direction as the thermal boundary layer develops. Under microgravity conditions and in the absence of boiling and forced convection in the fluid, the division of energy between the substrate and fluid can be computed using the well-known solution for transient conduction heat transfer in two semi-infinite solids having a plane heat source at their interface. For constant properties and constant plane heat generation rate at the interface, the division of heat between the two materials remains constant. The equations describing the transient interfacial temperature and division of heat transfer rate for a uniform constant heat generation rate at the interface are presented in Fig. 13, along with the properties for quartz (Fused polycrystalline), pyrex, BK-7, R-113 and the resulting division of heat transfer rates between the R-113 and each of the substrates. For steady-state conditions the fraction of heat input transferred to the R-1 13 with forced convection and/or boiling can be determined once the steady heat loss through the substrate is known as a function of the interfacial and the surrounding substrate holder temperatures. This is obtained by calibration. Fig. 14 shows a copper heater surface with the same size as the gold film heater of Fig. 9, and is intended to provide a metallic substrate to the fluid more representative of engineering-type surfaces. Because of the large heat capacity this can be used only for quasi-steady operation. It is gold plated to provide the same heater-fluid surface energy relationship as with the gold-coated quartz surface. It also has the same external configuration as the quartz heater, and thus is interchangeable in the test section. A film heater is compressed against the underside of the copper body, which is then encapsulated within a stainless steel housing. A copper foil 0.001 inch thick is then soldered across the entire upper machined surface to eliminate the crevices which otherwise serve as false 20

nucleating sites. Heat losses through the underside are obtained by calibration, as a function of the heater surface temperature. A Heise Model 623 pressure transducer is used for the measurement and control of pressure at the inlet to the test section. This has a maximum pressure of 50 psi and a sensitivity to 0.0025 psi. An uncertainty in pressure of.025 psi corresponds to an uncertainty of approximately 0. 1~F in the saturation temperature for R- 113. Temperature control at the inlet to the test section is maintained to within 0. 1F of the set point, which means that pressure control must be to within.025 psi to remain within 0.10F of the saturation temperature. For pressure calibration, a Ruska Pressure Calibration system was used. The air piston portion of this system is capable of resolving from 0.0001 psi at a level of 2 psi, to 0.005 psi at a level of 600 psi. This was used to calibrate both a precision Heise Bourdon tube gage and the pressure transducer. The Bourdon tube pressure gage is used for inplace periodic check calibrations of the pressure transducer, using the filling connection to the test vessel. A 7" diameter stainless bellows is installed in an appropriate housing for the pressure control within the test section, and is shown in Fig. 15. Both visual and electrical indications of the bellows position are provided. This, together with the two pneumatic solenoid valves and a modulating proportional valve constitute the mechanical components of the pressure control system, which controls the steady state pressure to within +.025 psi. The system is shown in Fig. 16. The flow control system illustrated in Fig. 17 serves to maintain the flow rate constant, as the system flow resistance varies due to changes in void fraction with boiling, or as the voltage output of the power supply changes. The turbine flowmeter, manufactured by the Halliburton Co., is rated for 0.3-3.0 GPM and has been calibrated over the full range, using water. 21

A schematic of the loop condenser/cooler control system, used to automatically maintain the desired degree of subcooling at the pump inlet, is given in Fig. 18. The condenser system also condenses any R-113 vapor which may exit the test section, and consists of two heat exchangers and a secondary pump. Heat exchanger #1 consists of three identical shell and tube heat exchangers in series, with each well baffled shell 5.7 cm OD (2.25 inches) containing 13 tubes 0.95 cm OD x 0.70 cm ID x 33.02 cm in length (3/8 in OD x 0.275 in ID x 13 in length), and transfers heat from the R-113 in the main flow loop to water circulating through heat excnanger #1. Heat exchanger #2, which consists of a plexiglass shell and copper tubing formed into a spiral, transfers heat from heat exchanger #1 to cooling water supplied from the building. A degassing unit and storage tanks for the R-1 13, including a pressure gage on each, have been fabricated from 304 s.s. A typical assembly is shown in Fig. 19. The concept involved is a combination of distillation at room temperature, leaving behind high boiling point components such as oils and solids, and freezing of the R-113 on the fins using Liquid Nitrogen within the inner vessel. By maintaining a sufficiently low pressure, the air components, except for water vapor, remain in the gaseous state and are removed by the vacuum pump. A molecular sieve is installed prior to the freezing vessel to absorb the water vapor. A filling system, mounted on a portable cart, has been fabricated for use with R-1 13, and consists of (1) a vacuum pump, including a trap, for removing the air prior to filling, (2) a flexible connection to the R-113 source tank, mounted on the cart, (3) a Heise compound Bourdon tube-type pressure gage covering the range 15 psig to 25 psig, with minor divisions of 0.05 psi, easily readable to.025 psi, (4) a flexible connection to the test vessel, and (5) associated valves and fittings. Fig. 20 demonstrates the variety of parameters possible with the flow loop operating in earth gravity alone. Only a relatively few test results have been obtained to date, to be presented in the next section. Fig. 21 presents the definition of the orientation angle 22

between the flow direction parallel to the heater surface and the buoyancy vector, used in these results. For an angle 8 = 0~, the flow is downward, opposite to the net buoyancy acting upward. 9 = 90o means that buoyancy is perpendicular and away from the heating surface. 9 = 1800 means buoyancy is in the same direction as the flow, while - = 270~ means buoyancy is acting perpendicular into the heating surface keeping any heated liquid or vapor formed in the vicinity of the heating surface. 23

4. EXPERIMENTAL PROCEDURES Because of the high solubility of air in R-1 13, precautions were taken in the filling process to insure that no contact was permitted between the degassed R- 113 and air. A filling cart held a vacuum pump, a valve system and associated piping, and a pressure gage. The vacuum pump was used for the evacuation of the test loop and the connecting line. The piping connected the test loop to the tank holding the degassed R- 113, and the valve system permitted the evacuation of the loop and piping without opening the valve of the R-1 13 container. These were evacuated with a vacuum pump for eight hours to remove air and water vapor. The pressure gage provided an approximate vacuum indication before the R-113 was allowed to flow from the storage container. In addition, this filling cart supported the elevated storage tank of the degassed R-1 13 which provided the hydrostatic pressure for filling. The storage container was irradiated with heat lamps to maintain a positive pressure relative to the atmosphere. When the pressure gage mounted on the storage container indicated a gage pressure of 34.47 kPa (5 psig), the valve between the storage container and the stainless steel lines was opened and R-1 13 flowed into the loop. The system was flushed with R-1 13 vapor several times before filling with the liquid, to further remove any possible residual gases. To avoid contamination of the degassed R- 113 by atmospheric air, the pressure in the loop was always maintained above the ambient, either by the bellows actuated pressure system or by heating. For each test, a single point calibration of the thin film heater surface was performed before heating the loop to the desired operating temperature. As discussed earlier, the linear temperature-resistance relationship changed with time, with observable changes occurring after about a week. However, the slope of the temperature-resistance relationship remained constant, and hence, only a single calibration point was necessary for the determination of the new temperature-resistance relationship. Prior to the heating of the loop, the surface temperature of the heater was known with certainty, since the entire loop 24

was at the known ambient temperature, determined with calibrated thermocouples. This known surface temperature was determined simultaneously with the heater surface resistance for the calibration point data. The forced convection boiling loop was used to examine the effects of an externally forced velocity on incipient boiling with transient heating and various heater orientations, as well as on steady-state boiling. Consequently, both the gold film heaters with an imposed step increase in heat flux and the copper heaters with steady power input were used. The bulk liquid velocities selected were sufficiently low so as to not mask buoyant effects. Velocities of 1.7 and 4.5 cm/s produced corresponding Re Numbers of 508 and 1343 in the test section. For the non-boiling conditions, the maximum Richardson Number is estimated to be 60, computed from Equation (2) so that buoyancy effects may be expected to dominate. The experiments were initiated after the forced convection boiling loop had reached steady-state conditions of temperature, pressure, and flow rate. All tests were conducted at the nominal liquid pressure of 120.66 kPa (17.5 psia), with fluctuations of less than 0.7 kPa (0.1 psia). The bulk liquid temperature was maintained either at 400C (1040F), which resulted in a nominal subcooling of 12~C (22~F), or at 46~C (115~F), which resulted in a nominal subcooling of 70C (13~F) for the imposed system pressure. Fifteen hours were required to reach steady-state conditions at the low subcooling level; four to five hours were necessary for the high subcooling level. Two orientations of the test section were used for the incipient boiling tests here: vertical and horizontal. The flow loop was rotated such that the heater surfaces were parallel to the earth gravity field (in the vertical orientation), and the buoyancy force was in the direction of fluid flow. The second test surface orientation was obtained by rotating the loop such that earth gravity was perpendicular t the heating surfaces. As a result, one of the surfaces was horizontal up, and the opposing surface was horizontal down, with 25

buoyancy aiding forced convection or hindering it, respectively. For tests conducted with steady-state boiling, a number of other orientations were used in addition to these. 4.1 Spatial Mean Gold Film Surface Temperature The spatial mean thin gold film heater surface temperature, Tw, is determined from the measured test surface resistance, Rw, the single point calibration resistance, RC, with dT its corresponding calibration temperature, TC, and the slope of the calibration curve, R, in the following equation: Tw = T + d (Rw -Rc) (6) For the constant heat flux tests, Rw, Rc, TC, and d were determined with representative values and uncertainties of 2.9985 +~0.0008 Q, 2.6266 + 0.0002 Q2, 20.00 + 0.010C, and 214.83 + 0.23~C/Q, respectively, as given in Appendix C.1. As a result of the uncertainties in these quantities, the heater surface temperatures were determined with an uncertainty of + 1.0C, as described in Appendix C-1. 4.2 Metal Heater Spatial Mean Surface Temperature The metal heater surface temperature was determined from the chromel-constantan thermocouple imbedded in the copper heating block near the heater surface, as shown in Figure 14. The error in taking the temperature at the location of the thermocouple to be equal to the surface temperature was estimated to be less than 0.0020C, by assuming one dimensional conduction in the copper block and negligible contact resistance between the soldered copper foil and the copper block. An error of this magnitude is less than the 26

uncertainty of + 0.050C associated with the thermocouple temperature measurement and, as a result, was neglected. 4.3 Bulk Liquid Temperature The liquid bulk temperature was measured by calibrated chromel-constantan thermocouples, with an uncertainty of + 0.050C (~ 0. 1~F). A discussion of the estimate of this uncertainty-is in Appendix C.2. 4.4 System Pressure The liquid pressure at the inlet to the test section was measured by calibrated transducers. The uncertainty of the pressure measurement was less than the desired uncertainty of + 0.172 kPa (~ 0.025 psi). The details of the estimate of the uncertainty are given in Appendix C.3. 27

5. EXPERIMENTAL RESULTS Experimental results obtained during the study period are presented below. The transient heating results are presented first, followed by the quasi-steady heating results. 5.1 Transient Boiling The transient boiling experiments were performed using a step in heater surface heat flux. A step change in heat flux is the most elementary form of imposed heat flux variation possible, as all other time varying imposed surface heat fluxes can be constructed from a series of steps in heat flux. Since the resistance of the gold film heater surface is computed from the measured current and voltage drop across the entire heater surface, the resulting calculated surface temperature is a spatially averaged temperature. The transient heater surface temperature for a representative test with pool boiling was shown in Fig. 11, along with a tabulation of the test conditions, and applies to the case where the heater surface is facing horizontal up in earth gravity. The onset of natural convection appears as an irregularity in the temperature-time plot. The time from the energization of the heater to the onset of natural convection is designated as tnc and is identified by the departure of the heater surface temperature from the one dimensional semiinfinite media transient conduction solution, and also by the observed onset of fluid motion as characterized by a wave-like disturbance recorded photographically. The next significant event in Fig. 11 following the onset of natural convection is incipient boiling, as indicated. Incipient boiling, or the onset of boiling, is defined as the appearance of the first vapor bubble on the heating surface. In some tests, this incipient boiling takes place at the point when the mean heater surface temperature reaches a maximum, while in other tests, the onset of boiling occurs prior to this point. For tests in which the latter took place, the level of the ensuing maximum heater surface temperature depended on the manner in which the boiling propagated across the heated surface. If 28

incipient boiling occurred as an almost explosive event over the entire heater surface, then the temperature associated with boiling incipience was the maximum surface temperature. If incipient boiling occurred at a location near a corner of the heater, for example, the heater would become cooled locally while the remainder of the heater surface continued to rise in temperature. The heater surface temperature measurement is a spatially averaged quantity, as described earlier, and the measurement would continue to rise until the boiling spread sufficiently to produce a subsequent decrease in the mean value. This boiling propagation will be defined in more precise terms later. Following the onset of boiling, a quasi-steady boiling region is noted in Fig. 11. In this domain of the temperature-time plot, boiling has spread across the entire heating surface, and the quasi-steady boiling temperature level is less than the maximum heater surface temperature, but above the saturation temperature for the liquid at the system pressure, as expected. Figures 22 - 24 present the measurements of the transient heater surface temperature, which is a spatial average, for experiments conducted with forced convection. The independent test variables included three heat flux levels, two subcoolings, and two flow velocities, in addition to the orientations described below. High speed video recordings of the incipient boiling across and through the heating surface were made. With both surfaces positioned opposite each other in the forced convection boiling loop, surface Q-13 was heated in the horizontal down position and surface Q-16 was heated in the horizontal up position. In the vertical flow experiments, where the forced convection boiling loop was rotated 90~, the R-113 flowed upward such that buoyancy assisted the externally forced flow. The data are tabulated in Appendix D. To illustrate the effect of orientation on the heater surface temperature, heater surface temperatures resulting from a nominal q"T of 4 W/cm2 and with similar initial conditions were used in Figs. 22 - 24. In Fig. 22, for the heater surface in the horizontal up orientation, the heater surface temperature decreased to a lower quasi-steady boiling level 29

once boiling had propagated over the entire surface, as also occurred in the case of pool boiling in Fig. 11 with the heater surface in the horizontal up orientation. In Fig. 23 with the surface in the vertical position, the heater surface temperature dropped to a quasi-steady boiling level after incipient boiling, which was considerably lower than that for the horizontal up orientation. This was expected for the heater surface in the vertical orientation, since the buoyancy force and the inertial force acted together to remove the vapor, which otherwise had an insulating effect, from the heater surfaces. As in the pool boiling vessel with the heater surface in the horizontal down position, a single vapor bubble covered the heater surface in the forced convection boiling loop with the heater surface horizontal down, even at the highest flow rate here, 4.5 cm/s. As a result, the heater surface temperature was unsteady and continually increasing, as shown in Fig. 24. It may also be noted that the measured surface temperature followed the one dimensional conduction prediction up to incipient boiling. Approximate adherence to the one dimensional conduction prediction is a consequence of 83 percent of the input energy going into the quartz substrate. The effect of orientation in earth gravity on the heater surface temperature at a q"T of 4 W/cm2 as shown in Figures 22 through 24 is representative of the effect of the heater surface orientation on the measured heater surface temperature at the other heat flux levels. Figures 25 through 28 present the heater surface superheats at the onset of boiling as a function of q"T for the various orientations, subcoolings, and fluid velocities. Figures 25 and 26 are for the low velocity, 1.7 cm/s. The heater surface superheats lie in a band between 10 and 400C for surface Q-16 in Figure 25, positioned horizontal up, and between 15 and 55~C for surface Q-13 in Figure 26, positioned horizontal down. As in the case of pool boiling, the time delay between the onset of boiling and the spread of boiling is most prevalent near a q"T of 8 W/cm2 for surface Q-13. The largest heater surface superheats at the onset of boiling are for the low subcooling tests for a given heater surface orientation. The largest surface superheat at a velocity of 1.7 cm/s occurred for Q-13 in the vertical 30

position in Figure 26. In Figures 27 and 28 for the high velocity of 4.5 cm/s, the greatest heater surface superheats are for the low subcoolings in the horizontal down position. In Figure 28 for surface Q-16 at the fluid velocity of 4.5 cm/s, the heater surface superheats show no time delay between incipient boiling and boiling spread, and the superheats lie in a band between 25 and 400C. From examination of Figures 25 through 28 for the vertical orientation, it may be concluded that the bulk velocity effects are small relative to buoyancy induced convection. Figures 29 and 30 present the delay time for the onset of boiling as a function of q"T with two subcooling levels with the horizontal down and vertical orientations for surface Q13, for flow velocities of 4.5 cm/s and 1.7 cm/s, respectively. For the same subcooling level, the horizontal down surfaces tend to have smaller boiling inception delay times than the vertical heater surfaces, attributed to the assistance of buoyancy in the vertical orientation. The range in delay time before incipient boiling is nearly the same at both flow velocities. A time interval between incipient boiling and boiling spread, as was the case in pool boiling, occurred at the q"T of 8 W/cm2. As also was the case observed with pool boiling, the delay time to boiling inception decreased with increasing q"T, as shown in Figures 29 and 30. Figures 31 and 32 parallel Figures 29 and 30 for surface Q-16, but with the horizontal up in place of the horizontal down surface orientation. For the same subcooling, the horizontal up surfaces have larger boiling inception delay times than the vertical heater surfaces, and the range in delay time before the onset of boiling is essentially the same at both flow velocities. 5.2 Quasi-Steady Boiling The early measurements of the forced convection boiling heat transfer behavior in the flow loop described here were made in two groups, designated A and B here, depending on the orientations of the heater surface relative to earth gravity, and are so presented 31

below. In all cases the liquid velocities are quite low, so that buoyancy indeed plays a role, and the liquid velocity vector is parallel to the heater surface as illustrated in Fig. 7. Thus, the angle between the flow velocity and gravity is changed as the loop is rotated, and is defined in Fig. 21. Group A consists of testing conducted with flow in the vertical up direction, so that the flow velocity enhances buoyancy with 9 = 180~ in Fig. 21. Opposing sets of heaters, including both the metal heaters of Fig. 14 and the quartz-gold film heaters of Fig. 9, were installed in the test section of Fig. 6. Since these heaters were fabricated so as to be identical, operation in the horizontal orientation would simultaneously provide results for - = 900 and e = 2700 in Fig. 21, which are also included in Group A. However, as will be demonstrated, the behavior of the two surfaces presumed to be identical were not the same, so that results of changes between the vertical and the two horizontal orientations could only be compared in a qualitative manner. The differences in the changes were nevertheless quite dramatic. Subsequent to the testing which constituted Group A, the loop was modified so that it could be rotated through a total of almost 3600. In this way the same heater surface can now be subjected to all possible orientations between the flow direction and earth gravity. Results of these tests to date are presented as Group B below, and include only the use of the quartz-gold film heaters. In using the metal heater of Fig. 14 the determination of the boiling heat flux from the measured power input requires that the heat loss from the lateral sides and the rear be estimated. Referring to the representative quantities plotted in Fig. 33, from the measured power input qT-in with non-boiling convection on the fluid side, the computed heat transfer to the fluid, qFc, together with the computed heat transfer to the air qA using wellestablished correlations, the heat loss qL is determined by the difference. Since the primary path of the loss cannot be identified precisely, an attempt was made to relate this loss to the average of the measured front and rear surface temperatures, Tw and Twb, respectively, 32

and to what is deemed to be the main sink for the heat losses, the test section housing whose temperature is close to the fluid temperature Tf. These are used to define the abscissa in Fig. 33. The ordinate on the left is expressed in terms of the power or heat transfer for the system of Fig. 14, while that on the right is expressed in terms of the heat flux at the boiling heat transfer surface area. It should be noted that the heat loss to the air out of the back side is negligible, and that the heat loss qL is virtually a linear function of the temperature difference defined above. Although the losses here appear to be a major part of the measured power input, the temperature differences are quite large, whereas with nucleate boiling taking place these differences will be considerably smaller, and the heat losses will be represented by the values on the left portion of the curves. Such calibration curves must be generated for each surface and for each orientation. 5.2. 1 Group A Results: Metal and quartz heaters; 0 = 900, 1800, 270~. Initial nucleate boiling results are shown in Figures 34 and 35 for the conditions given at the top of each figure. The numbers accompanying each data point represent the sequences in which the data were taken. The two sets of curves reproduce each other surprisingly well, and such should not normally be anticipated. Noteworthy are the large heater surface superheats prior to the onset of nucleate boiling, the lack of hysteresis effects in the nucleate boiling curve, and the distinctive cessation of boiling at heater surface superheats considerably lower than the onset. The non-boiling convection data at the end coincide with that at the beginning. The relatively large increase in heat flux once boiling begins occurs with a small increase in power, and is a result of the considerable decrease in heat loss accompanying the decrease in heater surface temperature. Subsequent to these tests the power supplies were modified to provide heat flux levels up to 20 W/cm2. Fig. 36 shows the influence of orientation changes between vertical and horizontal-up and horizontal-down at the low velocity used. The non-boiling measurements for the three 33

orientations are indistinguishable from one another. In the vertical-upflow configuration, the two surfaces show a displacement of heater surface temperature by about 2 1/2~F at any given heat flux, to be contrasted with the data in Figs. 34 and 35, which were virtually identical. Rotation to horizontal-up (0 = 900) produces an enhancement in performance, attributed to a decrease in the superheated boundary layer thickness because buoyancy is acting perpendicular to the flow direction. On the other hand, rotation to horizontal-down (9 = 270~) results in dryout of the surface, since buoyancy holds the vapor generated against the heater surface, and the velocity is too low to produce shear stresses sufficient to move the vapor away. The associated heat transfer rate is lower than the single phase liquid convection because of the lower thermal conductivity of the vapor, even with the reasonably high degree of inlet liquid subcooling. The influence of two low velocities (1.7 and 4.3 c/ms) are demonstrated in Fig. 37 for the vertical orientation with liquid upflow, and in Fig. 38 for the horizontal orientation, facing both up and down. The change in heat flux with velocity in the non-boiling domain is clearly seen with the vertical up flow in Fig. 37, and with the horizontal orientation facing upward in Fig. 38. This is evidence that mixed convection, the combination of single phase free and forced convection is taking place. On the other hand, when the horizontal orientation is changed to facing down, also in Fig. 38, boiling and dryout again occurs, with no observable difference in behavior for the two velocities used. These two low velocities likewise have no influence on the nucleate boiling process in the horizontal orientation facing upward (9 = 90~), which is indicative that buoyancy effects are dominating the process here. This is to be contrasted with nucleate boiling in the vertical upflow orientation in Fig. 37. For the low velocity of V = 1.7 c/ms the two heater surfaces differ by about almost 3~F in superheat, owing to somewhat different nucleating characteristics. When the flow velocity is increased to 4.3 cm/s, these now take on an identical behavior, and with a larger slope. 34

Comparisons between the two types of heater surfaces used are presented in Fig. 39 for the vertical orientation with upflow (9 = 180~) and in Fig. 40 for the horizontal orientation facing upward (9 = 900) and downward (8 = 270~). For the vertical orientation in Fig. 39, differences between the opposing pairs of heaters are noted, with constant heater surface superheat difference of about 70F for the quartz heaters and 20F for the metal heaters. For this same velocity, the two metal heaters had identical surface superheats in the vertical orientation in Fig. 37, for which a liquid subcooling of 23.30F was present, to be contrasted with 13.20F in Fig. 39. This demonstrates that liquid subcooling can play a significant role. In spite of the larger differences between the quartz and metal heaters, the slope of heat flux-heater surface superheats are quite similar, with the metal surface superheats displaced to lower values, as is well known owing to its "better" nucleating characteristic. The negative slope of the quartz heaters at the lower levels of heat flux, with the data obtained with increasing heat flux, are manifestations of hysteresis effects, which are magnified considerably with smooth surfaces such as is present with a polished quartz surface. Stable nucleating sites are less likely to be present, and the spreading of the boiling process over the heater surface must be against both buoyancy and the flow direction. Hysteresis effects will be presented in somewhat more detail in the following section. As must be the case with the vertical orientation, the behavior of all four surfaces is identical in Fig. 39 in the non-boiling domain. Displacement of the heater surface superheat persists between the quartz and metal in the horizontal orientation facing upward (9 = 900) in Fig. 40. A smaller hysteresis is evident with the quartz heater in this orientation, compared to the vertical one of Fig. 39, since the increase or spreading in the number of nucleation sites no longer needs to act against buoyancy, but only against the flow velocity. Only one data point was obtained with boiling with the quartz surface in the horizontal facing down orientation (8 = 270~), since the gold film used at this time could not withstand the relatively high surface 35

temperatures encountered. This point lies quite close to the data obtained with the metal heater, as it should since surface characteristics have no influence in the vapor bound or dry-out circumstance. 5.2.2. Group B Results: Quartz heater, 8 = variable. The results presented below were obtained over a wider range of orientations than those in the previous section, with only the gold film quartz heater used. Since the spectrum of orientation angles covered is as yet incomplete, detailed interpretations of the results would be premature and will not be given here. Fig. 41 presents the complete data obtained with both increasing and decreasing heat flux, and it is noted that the behavior for these two cases is identical in the steady state with a sufficiently high level of heat flux, 6 W/cm2 in this case. Figs. 42 - 46 demonstrate the influence of bulk liquid subcoolings of 40F and 20~F for orientation angles varying from 0~ to 1800, from flow directly counter to buoyancy to flow parallel to buoyancy. Data such as Figs. 42 - 46 are combined in Figs. 47 - 49 to illustrate the influences of orientation with the low velocity of about 4.3 cm/s used. Figs. 47 and 48 apply for low levels of subcooling, while Fig. 49 applies for the high subcooling. The increasing numbers on the plots correspond to increasing magnitudes of the orientation angle. All of the results in Figs. 42 - 49 were obtained with increasing heat flux. Fig. 47 includes only lower levels of heat flux, with orientation angles in the "inverted" positions, while Fig. 48 has higher levels of heat flux, with orientation angles in the upward-facing positions. Fig. 49 covers a wider range of orientations and levels of heat flux with the high subcooling of 200F. 36

6. FUTURE WORK a. Examination of the steady state nucleate boiling results with low velocities, in section 5.2 above, makes it clear that in conducting such experiments in earth gravity little can be discovered about the potential behavior for corresponding conditions in long term microgravity. Experimental work in both short and long-term microgravity will be required to understand the basic elements that constitute forced convection nucleate boiling at the smaller flow rates with low Reynolds Numbers, and should include heat flux levels up to the critical values. b. A number of extensionsof the initial work presented here can be carried out in earth gravity: i. Operation at near saturated liquid at the entrance to the test section, in which it is anticipated that the influences of orientation will be more dramatic. ii. Extend the levels of heat flux up to the critical or "bur-out" level as functions of orientation, subcooling, and low-level fluid velocities. iii. Study the influence of small variations from the horizontal down orientation (9 = 2700) with various small velocities to establish the conditions necessary to remove the vapor from the heater surface. 37

Tsat a/g X I NW ftm lb 1 T I THERMAL BOUNDARY LAYERS "\ ~,> 1 \- -—,I /-10 i //////////,/////// //// //////////////// HEATING a/g= 1 a/g < 1 SURFACE q" = CONSTANT Fig. 1. Influence of Buoyancy on Temperature Distribution in Vicinity of Heater Surface with Pool Boiling. 38

Lift ) Liquid-Vapor Interface Shear V (Bulk Liquid) Evaporation/Condensation (Molecular Momentum) Licuid V (interface) VaDor Liquid Pressure (Bo yancy) * Liquid-Vapor Surface Tension (Marangoni Convection) Wake Drag a Pressure * Liquid Visccsity Heater Surface ~~ Liquid-Solid-Vapor Surface Tension Fig. 2. Forces Acting on a Growing/Collapsing Vapor Bubble Attached to a Wall. 39

Thermal Management in Space (Change in pumping power requirements using phase change) m gout (Radiation to space) Ap - Pin Ip l liquid heat capacity: latent heat: q 7mHC C\T =.106kg/s rmLH - = 6.81 10-3kg/s hfg 7hHC or ~. = 15.56 rhLH Power: AP m f V2 L in laminar flow 1 fV P m h. Ap _ m2. L PHC ( mnc2 LHC 15.562 22 = 484 PLE mrLH LLH / Also: Mass of fluid required is reduced by 1/2!!! Fig. 3. Example of Advantage in use of Phase Change to Transport Energy. 40

Pressurizer Assembly Test Section uonaenser / cooier Assembly Pump Assembly Fig. 4. Early Flow Loop for the Study of Forced Convection Roiling of Various Orientations and at Reduced Gravity. 41

Heater Surftfces i^ u Flow CoVrve Valve Fig. 5. Schematic of Flow Loop System.

Viewing, Windows N 4Swe I I 4,5 Fluid Inlet Fig. 6. Test Section for Heater Surfaces.

To PO Tsai g V(y) T(y) Y, X0 Fig. 7. Representation of Flow Boiling on a Flat Plate Heater. Fig. 7. Representation of Flow Boiling on a Flat Plate Heater. 44

M1 L1 i,t/I\ M6 ap — _- __ ---- -__- - - - Test section in loop 'T, r% A I IS-1 M20 II M5 \ oO oo TS-2 F --- L2 I L3 I \ M4,' \ \ M3 NOAH SEDCAEAO Note: M1,M2 & M3 --- Cold Mirrors NOVA HI SPEED CAMERA OR HIGH FRAMING RATE VIDEO Fig. 8. Optical Paths for Two View Photography. 45

Gold Flm Surface ~nmmm I 2 5 4 7A A 4 PARTS I Teflo Suppor Ring 2 Qunz 3 Qumnt Leas 1000 A Gold Lyer 4 Semi-Tranparent Gold Boiling Sorface IS X 0.75" X- 400 A 5 Potamu Lids 6 "0' Rig Sea 7 Korvr Lads I Iwww — --- - ----- -* L 5 3 t_ Fig. 9. Gold Film Heating Surface as Simultaneous Resistance Thermometer. 46

Tsat TEST FLUID.(Tsub) X T ----t 400 A TH (t) I Il *' r (~> - -11 -QUAti I -.0UBSTRATt= - - - - - ~ 1 - I - * I I~~~~~~~~~~~~~~~ - TRANSPARENT GOLD FILM (HEATER + RESISTANCE THERMOMETER) - TRANSPARENT GOLD FILM (HEATER + RESISTANCE THERMOME TER) AIR Fig. 10. Demonstration of Concept for Transient Heating with a Thin Film. 47

160 '1 V - - ID Transient Surface Q-5 Conduction Prediction a/g = +1 q"T = 5.7 W/cm2 140 - ATsub = 2.470 C P= 117.1 kPa Ts= 51.940 C Ti = 49.470 C 6.894 Run # 79 C PAUT0429.605 / Observed o 120 Incipient <Sriia~~ / o Boiling- v I, /, t,.*, E0 100- Eu \ Qi-y Quasi-Steady 80- \ Boiling Observed Onset of Natural Convection tnc = 1.2 sec 60 -402 0 2 4 6 8 Time (Sec) Fig. 11. Thin Gold Film Temperature with Onset of Heating in Pool Boiling at a/g = 1. 48

160 140 - ci O e4 O U 0 o 5 C.) - o ci [0 am cn 120 -100 - ID Transient Conduction Prediction Surface Q-5 a/g 10-s q"T = 5.8 W/cm2 ATsub = 2.89~ C P = 110.11 kPa Ts= 50.0~ C Ti= 47.11~ C Run # M-17-16 PBMT1115.605 /. ~ ---- Impact /\ " ""...!" Observed /. 'Incipient Boiling t/ = 3.083 sec 80 60 40 i II. 0 2 I 4 Time (Sec) 6 8 Fig. 11. Thin Gold Film Temperature with Onset of Heating in Pool Boiling at a/g - 10-5 for 5 seconds. 49

Ti ) - To = 2(asa9)1/2 qT frl/2(k al/./2+ksa, /2) Ti/2' ( ka T(X, ) = To t' =2 CO ~It r tl/2 qI qz s k%( as ) 1/2 k a s 2Z: a = k PC qa q" T =_ 1 s Z Fig. 13. Heat Transfer Between Two Communicating Semi-infinite Solids with Uniform Plane Heat Source at the Interface. 50

A ITXM 1. TEKUMXAL STRIP SCREW 2. TERMINAL STRIP 3T.C.NUT-SEALiO 4. EALIO BODY M. FM HEATE L. HEAT HOLDOI WN SCREW 7. HATER HOLDOOWN BIACZOT L G.OM' C OP!P' POFLOLDUD IN PEACE 9. VACUUM VtENrTT 1l. STYCAST L AE THEUMR CUPLE 2 COP IIHEATE BOMY 1I. TULON SUIIPPOTBRA 14. TLON LOCATTO PIN 1s. Sa HEATE' ENCLDOSU A CTITT'fA1 VYTF W Fig. 14. Copper Heater Surface for Quasi-steady Operation. 51

Fig. 15. Pressure Control Bellows Assembly. 52

PRESSURE TRANSDUCER POWER AND SIGNAL \ SOLENOID VALYES Prpeortlena YVlve Mdlre YVlve PRESSURE TRANSDUCER HEISE GAUGE MODEL 623 (or Seera-Hodel 270) 0- SPSIA 0-5 YDC Fig. 16. Flow Loop Pressure Level Control System. 53

Flowrate Flowmeter set oPint ( I Power Pump Out rr Ampontroller Drive Pump-~ E+ r Motor 4.- -, I Fig. 17. Schematic of Loop Flow Control System.

mIo Main Loop P R1 13 1 Flowmeter r 3 Pass Counterflow Heat Exchanger */ 1 Counterflo Heat Exchanger Water From Building Suppl y ressurizer 4 Command ends On RTD Feedbnck From Main Pump I nlet Water nw *2 To Drain Fig. 18. Schematic of Loop Condenser/Cooler Control System. 55

Freezng Fins Watr Vapor Absorber Sight Windows Source Vessel Non-Condensable Gas Removal Cold Trap Vacuum Pump R-113 Reciver Vessel tig. iy. Schematic ot K-113 Degassing Unit. 56

Forced-Convection Nucleate Boiling Nucleation & Steady-State Incipient Boiling Boiling Transient Quasi-Steady Pool Forced Convection Independent Variables Fixed Fluid Fixed Pressure R-113 (near Ambient) r Semi-transparent Gold on Quartz Heater Surfaces < Gold on Metal (Copper) Variable AT subcooling Variable Velocity (velocity distribution) Variable Heat Flux Variable Orientation Dependent Variables Conditions for incipient boiling Vapor bubble motion Mean heat transfer behavior Fig. 20. Experimental Parameters for Flow Loop Operation in Earth Gravity. 57

Definition of Orientation Angle '//////7//////7///////////////////////// Inlet Flow lI g m 9=90 q /// ^//////////// q" Orientations used in experiments =180 = 135 8 =90 Fig. 21. Flow Gravity Boiling Heat Transfer with Various Earth Gravity Orientations. 58

n1 A mg p U O Cu h., C U 130 120 110 100 -90 -80 -70 4 1D Conduction Prediction Surface Q-16 ~ ' a/g = + q"r = 4.7 W/cm' ATsub = 12.6 OC P = 120.45 kPa T,= 52.8 OC Ti = 40.2 C Velocity = 4.5 cm FAUT0607.412 -incipient Boiling Quasi-Steady Boiling 60 50 40 I a I I a I a L/S I — ~~~~~~~~ I --- 0 *I I T i I 2 4 6 8 Time (Sec) I 10 12 14 Fig. 22. Thin Film Gold Heater Temperature in Forced Convection Boiling with Horizontal Up Orientation. 59

140 al I c Ce a 3 b. m (u I. 0 lu2 cn 130 - 120 110 -100 -90 - 80 - 70 60 - 4 50 -40 0 a 4 I ID Conduction Prediction [ / Incipient / Boiling / * "* * 4..* r4 4 — Quasi%/,e* / rBoilin /. Surface Q-13 I.*- eag v q"T = 4.2 W/cm2 / ' ATsub =12.0 C P = 119.56 kPa *i T =52.6 C | ~ Ti=40.6~C Velocity = 4.5 cm/s FAVT0516.400 i~~~~~~~~~~~~~~ - 1 Steady g I __ _ i, I 0 2 4 I ' - ' I 6 8 tO 12 lr 14 Time (Sec) Fig. 23. Thin Film Gold Heater Temperature in Forced Convection Boiling with Vertical Orientation 60

1 A III 130 120 -110 - U Q 100 o a too u 90 0. E U U c8 S 80 cn V) / ~ ee ~ -— incipieicnt Boiling ll..! / ID Conduction Prediction 14 14 I 70 60 50 Surface Q-13 a/g = - q"r = 5.4 W/cm2 aTsub = 12.6 ~C P = 120.45 kPa Ts= 52.8 ~C Ti = 40.2 ~C Velocity = 4.5 cm/s FADT0607.423 40 I F' a' a a r I 0 I * 1 - - - I t * I 2 4 6 Time. 8 10 12 (Sec) Fig. 24. Thin Film Gold Heater Temperature in Forced Convection Boiling with Horizontal Down Orientation. 61

80 70 J 60 -uO.0 = 40 U I 30 u c2 20 - Subcooling (Deg C) a 7 a/g= v o 12 + 7 a/g = + x 12 Vertical ' — X a x x $ X lw x 10 - 0 - + x v I I ' I I ~ ' I ' i I 0 1 2 3 4 5 6 Total Heat Flux. q"T 7 (W/cm2) 8 9 10 Fig. 25. Heater Surface Superheat at Incipient Boiling for Surface Q-16 at 1.7 cm/s. 62

80 70-! c O C4 u a w I m o C. U cz 0 2i 60 -50 - 40 30 -20 m IN I a 0 4 0 4 4 Subcooling (Dcg C) a 7 a/g= v o 12 * 7 a/g = -I * 12 Initial Boiling i mmmm Boiling Propagation Vertical 0 0 0 I I I I o I 1 I I II II II bL I 10 - I - I I I I...... I V I I I ' ' Il 0 1 2 3 Total 4 Heat i a ~ I 5 Flux, 6 7 q"T (W~cm2) 8 9 10 Fig. 26. Heater Surface Superheat at Incipient Boiling for Surface Q-13 at 1.7 cm 63

80 I -, 70 - Boiling Propagation I Ml' Initial Boiling Subcooling (Deg C) a 7 a/g = v 0 12,1 * 7 a/g = -1 * 12 - 60 -0,0 |, 40 -0 30 -0 20 $.M mC i a a.. 4 0 0 a a / - Vertical a 10 0 $ I 1 I I I 4( 0 I _,. I I I. * 0 1 2 3 Total 4 Heat I ' i I * 5 6 7 8 9 10 Flux, q"T (W/cm2) Fig. 27. Heater Surface Superheat at Incipient Boiling for Surface Q-13 at 4.5 cm/s. 64

80 70 G 60- a/g= +1 3 x Ij Vertical 40 '= 40- \ - +. 30 I20 10 -O ' r I " ' I " ' " ' w ~ T I * I.. c2 30 0 1 2 3 4 5 6 7 8 Total Heat Flux. qT, (W/cm2) Fig. 28. Heater Surface at Incipient Boiling for Surface Q-16 at 4.5 cm/s. Fig. 28. Heater Surface at Incipient Boiling for Surface Q-16 at 4.5 cm/s. 9 10 65

100 I 9 6' Boiling Propagation Initial Boiling 10 1 rM u 2 I E,;.1 msm 1 Ot 0 0 0 0 0 a Vertical 9I Subcooling (Deg C 7 'a/g = v o 12 0 7i a/g = -1 0 12 Test Surface Q-13.) I N I I I I.01 Z....... 0 2 4 6 Total Heat Flux, q"T, (W/cm2) 8 1C Fig. 29. Incipient Boiling Delay Time at 4.5 cm/s for Surface Q-13. 66

100 10 c0 I —, * [-, El 1 1 01 0 2 4 6 8 Total Heat Flux. q"T, (W/cm2) Fig. 30. Incipient Boiling Delay Time at 1.7 cm/s for Surface Q-13. 10 67

100 E 1 q Subcooling (Deg C) l oo.s 10 ' 0 12 Test Surface Q-16 01 0 2 4 6 8 It Total Hat Flux q", (W/cmVertical Fig. 31. Incipient Boiling Delay Time at 4.5 cm/s for Surface 16. 68 Test Surface Q-16 0 2 4 6 8 i. Total Heat Flux, q"T, (W/cm2) Fig. 31. Incipient Boiling Delay Time at 4.5 cm/s for Surface Q-16. 68

0 Vertical 10 t 0 a, 1- Subcooling (Deg C) E 1' '| n 7 7 a/g=v v 0 12.1' E 7 a/g = +1 0 12 Test Surface Q-16.01 O 0 2 4 6 8 10 0 Total Heat Flux, q"T, (W/cm2) Fig. 32. Incipient Boiling Delay Time at 1.7 cm/s for Surface Q-16. 69

40 30 -20' cl 10 cr Cl u cr -10 0 20 40 60 [(T + T )/2 - Tf ]F 80 Fig. 33. Representative Heat LossCalibration - Metal Heater No. 2 - Vertical. 70

P=17.91psia. Tsat=128.6*F., Tsubl18-F, v=1.72cm/s. -*- Boiling nd after 16:r 5 15 4 14 The number indicate the sequence of the est dat 13 3 1 ~~ 17 18 2 19 20 22 21 11 10 24 5 -1 I ~ '~ I 3 -20 -10 0 10 20 30 40 50 ATwsup (K) Fig. 34. Sample Forced Convection Nucleate Boiling Results. Subcooled, Vertical, Low Velocity, Metal Heater No. 1. Test No. 62289. 71

P=17.91psia. Tsat128.6-F, &Tsubl18.F, v=i.72cm/s 6. -r i cr _- Before boiling 16 5 - *- Boiling end aft 15 The n Im r indicate the sequnc of the tat dat 4~ 3 - 14 12(13) 17 18 19 20 22 11 2423 24l 10 9 0 - 8 6 7 ~ 1 2 3 4 -1.1 I u -20 I I - I ' I I -10 0 10 20 30 40 50 ATwsup (K) Fig. 35. Sample Forced Convection Nucleate Boiling Results. Subcooled, Vertical, Low Velocity. Metal Heater No. 2. Test No. 62289. 72

15 i -20 0 20 40 Twup ( 60 Fig. 36. Influence of Orientation with Forced Convection Nucleate Boiling. Metal Heaters. ATsub = 23.2~F, V = 1.7 cm/s. 73

15 i i * Hvl Sfl nb 1.7c o HvlSflblg * HvlS nb * Hvl Si2 blg - L -Surface 2 5 i Nonboiling va4.3cm/s Surface 1&2 -> Nonboiling v= 4.3 cm/s 4.3 cm/v=1.7cm/s -20 0 40 60 Twsup (K) Fig. 37. Influence of Velocity with Forced Convection Nucleate Boiling and Vertical Up Flow (= 180~). Metal Heaters. ATsub = 23.3~F. 74

1 1 - aLA.LIJ &1 kI 10 * HvlHUBIg ~ * Hvl HD Nb '< a/g=+l, Htl up, surface 1 v=1.7 & 4.3 cm/s 9 = 90~ o Hvl HD Big Lvl HUNb A Lvl HU Big * LvlD Nb o LvlHDBIg j4 Nonboiling v=4.3cm/s, = 90~ a/g=-l, Htl Dwn,~ c e~surface 2 ~A Nonboiling v=1.7&4.3cm/s v=1.7cm/s = 270~ 2- *9- = 900 0 -20 0 20 40 60 80 100 TwsupK) Fig. 38. Influence of velocity with forced convection nucleate boiling and horizontal orientation. Up (&= 900) and down (O= 2700). Metal heaters. ATub = 13.1~F. 75

15 10 5 -20 0 20 40 60 Twsup (K) Fig. 39. Comparison of metal heater versus gold film quartz heater with otherwise identical conditions: vertical up flow V (t= 180~); ATsub = 13.2~F; P = 16 psia; V = 4.3 cm/s. 76

10 a/g=+l, Metal heater Nonboiling, Quartz a/g=+l Metal a/g=+~ 5, ~~ ~~i i ~boiling, a metal Boiling a/g=-l Quartz -20 0 20 40 60 Fig. 40. Comparison of metal heater versus gold film quatrz heater with otherwise identical conditions: Horizontal facing up (4O= 90~) and down (4= 270~); ATsub = 13.1~F; P = 16 psia; V = 4.3 cm/s. S ) 80 vsup(K 77

8 - 6 - N -3:3 -I 0 -1 *rn 11 Fluid: R-113 Surface Q23 Flowrate = 4.3 C/ys Subcooling = 5 F Inlet temp = 117 F - = 225 degrees a Increasing flux * Decreasing flux 4 - w * a U 0 E a a a * 0 U 0 10 20 30 40 50 60 ATsup (F) Fig. 41. Hysteresis Effects with quartz heater at 4 = 225~. Test No. 112090. 78

Heat Flux vs. Superheat 12 10 0 0 0 E 0 x UL. 0 0 6 - Angle = 0 degrees Inlet Temperature = 120.0 Flow Rate = 4.3 cm/s Gold Heater Surface (Q23) 0 ATsub = 4 F * ATsub= 20 F 0 0 (1 3: 4 - * E 2 0 0 *; w — * a 0 EX [X 3 0 I 10 20 30 40 50 60 AT superheat (F) Fig. 42. Influence of Subcooling for — = 0~. Test No. 022591. 79

Heat Flux vs. Superheat 14 12 10 CJ E U 31, x ) Q) 8 - -t a * a3 * X0 Q~ I Angle = 45 degrees Inlet Temperature = 120.0 Flow Rate = 4.3 cm/s Gold Heater Surface (Q23) 0 ATsub = 4 F * ATsub = 20F 6 - 4 2 -0 - 9 0 0 U U 0 0 a 0 a U 0 1 10 2 20 30 4 40 50 60 AT superheat (F) Fig. 43. Influence of Subcooling for-4 = 45~. Test No. 022591. 80

Heat Flux vs. Superheat 14 12 10 Angle = 90 degrees Inlet Temperature = 120.0 Flow Rate = 4.3 cm/s Gold Heater Surface (Q23) * ATsub = 4 F * ATsub = 20 F E E 0 x 3 lU.. 8 6 Ic =D 4 2 0 10 20 30 40 50 60 70 AT superheat (F) Fig. 44. Influence of Subcooling for4 = 90~. Test No. 022591 81

Heat Flux vs. Superheat 15 -10 CJ < E o x 0 =. -r 0. 0 0 Angle = 135 degrees Inlet Temperature = 120.0 F Flow Rate = 4.3 cm/s Gold Heater Surface (Q23) 0 ATsub = 4 F * ATsub = 20 F 5 * 13 * 0 a * r 'IE * 0 10 20 30 40 40 50 0 60 AT superheat (F) Fig. 45. Influence of Subcooling for 9 = 135~. Test No. 022591. 82

Heat Flux vs. Superheat 15 -10 - I 0 Angle = 180 degrees Inlet Temperature = 120.0 F Flow Rate = 4.3 cm/s Gold Heater Surface (Q23) 0 ATsub=4F * ATsub= 20 F < cq E x 3 t. * 0 0 0 * a -) I 5 - * 0 * 0 * a ~ la 0a * % * 0 * ^0 =- - a 0 10 20 30 40 5 5o 60 70 AT superheat (F) Fig. 46. Influence of Subcooling for4 = 180~. Test No. 022591. 83

'h 8 N E x 0 4 - w I / I I / 0 / 0* o = 285 0=0 Fluid: R113 Surface Q23 Flowrate = 4. / c'/ Subcooling = 6 F Inlet temp = 117 F X 0=315 a0 I * mG = 225 2 0 - * o1 * * o a EX v U * * m a 0 2 20 40 60 80 ATsuperheat (F) Fig. 47. Influence of Orientation with Low Subcooling. Low Heat Flux. Test No. 112090. 84

Heat Flux vs. Superheat at Varying Orientations 14 12 10 N E x X 8 I * U 11 01 0 Inlet Temperature = 120.0 Subcooling = 4.0 F Flow Rate = 4.3 cm/s Gold Heater Surface (Q23) Angle = 0 Angle = 45 Angle = 90 Angle = 135 Angle = 180. 6 0 3) I U 0 / 0. m 0 * 4 m U 0 2 0 I E O U* 0 II [ *0 J" o U. 1 ~ ~ U * * I. - I 30 I I 40 20 0I 50 6 60 70 AT superheat (F) Fig. 48. Influence of Orientation with Low Subcooling. Test No. 022591. 85

14 12 10 < - 5. 3 x LL CZ 0) Inlet temperature = 120.0 F Inlet velocity = 4.3 cm/s Subcooling = 20.0 F Gold Heater Surface (Q23) ( a Angle=0 () * Angle=45 () Angle 90 * ~ Angle 135 ( * Angle=180 a Angle= 225 8 6 4 2 0 0 10 20 30 40 50 60 AT superheat (F) Fig. 49. Influence of orientation on steady flow nucleate boiling. High subcooling. 86

REREENCES 1. Ollendorf, D., "Thermal Control for the 1990's," Technology for Space Astrophysics Conference: The Next 30 years Sponsored by AIAA/SPIE/OSA, pp. 107-111, Danbury, Conn., Oct. 4-6, 1982. 2. Chen, John, "Correlation for Boiling Heat Transfer to Saturated Fluids in Convective Flow," Ind. Eng. Chem. Process Design Develop, Vol. 5, No. 3,.pp. 322-329, July 1966. 3. Bogdanov, F. E., "Effect of the Flow Velocity of a Vapor-Liquid Mixture of Coolant, and the Vapor Content, on Surface Heat-Transfer Coefficient in Boiling of Water Inside Tubes," Soviet Atomic Energy 29, 6, 1229-1232 (Dec. 1970) (Trans. Atomnaya Energiya 29, 6, 454-456, Dec. 1970). 4. Shekriladze, I. G., "Analogies Between Heat and Momentum Transfer in a Forced Liquid Flow With Surface Boiling in a Channel," J. Eng. Phys. 229, 4, 1249-1253 (Nov. 1976) Translation of Inzhenerno-Fizicheskii Zhurnal 29, 4, 631-636 (Oct. 1975). 5. Zaletnev, A.F., "Model of Heat Transfer Associated with the Surface Boiling of Liquid in Tubes," J. Eng. Phys. 31, 3, 1005-1009 (June 1977) Translation of Inzhenerno-Fizicheskii Zhurnal 31, 3, 396-401 (Sept. 1976). 6. Campolunghi, F., et al, "Subcooled Boiling Correlation for the Thermal Design of Commercial Steam Generator," Thermotecnica 31, 4, 185-194 (Apr. 1977). 7. Bennett, D. L. and Chen, J. C., "Forced Convective Boiling in Vertical Tubes for Saturated Pure Components and Binery Mixtures," AIChE J., 26, 454-461 (1980). 8. Stephan, K., and Auracher, H., "Correlations for Nucleate Boiling Heat Transfer in Forced Convection," Int. J. of Heat and Mass Transfer, Vol. 24, No. 1 pp. 99 -107, Jan. 1981. 87

9. Bjorge, R. W., Hall, G. R., and Rohsenow, W. M., "Correlation of Forced Convection Boiling Heat Transfer Data," Int. J. of Heat and Mass Transfer, Vol. 25, No. 6, pp. 753-758, June, 1982. 10. Dhar, P. L. and Jain, V. K., "On Evaluation of Correlations for Prediction of Heat Transfer Coefficient in Nucleate Flow Boiling," Int. J. of Heat and Mass Transfer, Vol. 25, No. 8, pp. 1250-1251, Aug. 1982. 11. Bilicki, Z., "Latent Heat Transport in Forced Boiling Flow," Int. J. of Heat Mass Transfer, Vol. 26, No. 4, pp. 559-566, April 1983. 12. Bradshaw, R. D. and King, C. D., "Conceptual Design for Spacelab Two-Phase Flow Experiments," NASA CR-135327, Dec. 1977. 13. Siegel, R., "Effects of Reduced Gravity on Heat Transfer," Hartnett, J. P. & Irvine, T. F., Jr., edited by Advances in Heat Transfer, Vol. 4, New York, Academic Press, pp. 143-228 (1967). 14. Merte, H., Littles, J. W., & Oker, E., "Boiling Heat transfer to LN2 and LH2: Influence of Surface Orientation and Reduced Body Forces," Proc. XIIIth International Congress of Refrigeration, Commission 1, Paper 1.62, Washington, D.C., Aug. 27-Sept. 3, 1971. 15. Merte, H., Hastings, L., "Zero Gravity Boiling Heat Transfer," Proc. Cryogenic Workshop Conference, NASA, Huntsville, Alabama, pp. 277-298, March 29-30, 1972. 16. Merte, H., "Nucleate Pool Boiling in Variable Gravity," in "Low-Gravity Fluid Dynamics and Transport Phenomena," Vol. 130. Progress in Astronautics and Aeronautics, Edited by Jean N. Koster and Robert L. Sani, AIAA, 1990. 17. Coeling, K. J., Merte, H., "Incipient and Nucleate Boiling of Liquid Hydrogen," Trans, ASME, J. Eng. for Industry, 91B, No. 2, pp. 513, May 1969. 88

18. Merte, H., Jr., and Littles, J. W., "Zero Gravity Incipient Boiling Heat Transfer," Proc. Space Transportation System Propulsion Technology Conference, Vol. 4 (Cryogens), pp. 1312-1348, NASA, Huntsville, Alabama, April 6 & 7, 1971. 19. Oker, E., Merte, H.. "A Study of Transient Effects Leading Up to Nucleate Boiling," Proc. 6th International Heat Transfer Conference, Vol. I, paper PB-5, pp. 139, Toronto, Aug. 7-11, 1978. 20. Nghiem, L., Merte, H., Winter, E. R. F., Beer, H., "Prediction of Transient Inception of Boiling in Terms of a Heterogeneous Nucleation Theory," J. Heat Trans., ASME, Vol. 103, No. 1, pp. 69-73, Feb. 1981. 21. Anderson, D. L. J., Judd, R. L., Merte, H., Jr., "Site Activation Phenomena in Saturated Nucleate Boiling", Proc. Symposium, "The Role of Nucleation in Boiling and Cavitation," ASME Paper 70-HT-14, ASME Joint Fluid Mechanics and Heat Transfer Conference, May 25-27, 1970, Detroit, Michigan. 22. Shoukri, M. & Judd, R. L., "Nucleation Site Activation in Saturated Boiling," J. Heat Transfer, Trans. ASME, 97C, #1, pp. 93-98, Feb. 1975. 23. Bergles, A. E., Rohsenow, W. M., "The Determination of Forced Convection Surface - Boiling Heat Transfer," J. Heat Transfer, Trans, ASME, Vol. 86, pp. 365, 1964. 24. Sato, T., and Matsumura, H., "On the Conditions of Incipient Subcooled-Boiling with Forced Convection," Bull. JSME 7, 26, 392-398, May 1964. 25. Grant, I. D. R. & Patten, T. D., "Thickness of the Thermal Layer at the Initiation of Nucleate Pool Boiling," Proc. Inst. Mech. Engrs., Vol. 180, Pt. 3C, pp. 124 -134, 1965-66. 26. Gaertner, R. F. & Westwater, J. W., "Populations of Active Sites in Nucleate Boiling Heat Transfer," CEP Symp. Ser., Vol. 46, No. 30, pp. 39, 1960. 27. Gaertner, R. F., "Distribution of Active Sites in Nucleate Boiling of Liquids," CEP Symp. Ser., Vol. 59, No. 41, pp. 52, 1963. 89

28. Singh, A., Mikic, B. B., & Rohsenow, W. M., "Active Sites in Boiling," J. Heat Transfer, Trans. ASME, 29, Ser. C, #3, pp. 401-406, Aug. 1976. 29. Kocamustafaogullari, G., & Ishii, M., "Interfacial Area and Nucleation Site Density in Boiling Systems," Int. J. Heat Mass Transfer, Vol. 26, No. 9, pp. 1377-1387, Sept. 1983. 30. Treshchev, G. G., "The Number of Vapor-Formation Centers in Surface Boiling," in Convective Heat Transfer in Two-Phase and One Phase Flows, edited by V. M. Borishonskii and I. T. Pallev, Transl. from Russian by Israel Program for Scientific Translations (1969). 31. Marchs, B. D. and Dropkin, D., "Measured Temperature Profiles within the Superheated Boundary Layer Above a Horizontal Surface in Saturated Nucleate Pool Boiling of Water, " Trans. ASME 87 C (J. Heat Transfer) 3, 333-341, August 1965. 32. Judd, R. L. & Merte, H., Jr., "Influence of Acceleration on Subcooled Nucleate Pool Boiling," Proc. 4th International Heat Transfer Conference, Vol. VI, paper B. 8.7, Paris, Aug 31-Sept. 5, 1970. 33. Judd, R. L., Lavadas, C. H., "The Nature of Nucleation Side Interaction," Journal of Heat Transfer, Transactions of the ASME, 102, 3, pp. 461-464 (Aug. 1980). 34. Dougherty, D. E., Rubin H. H., "The Growth and Collapse of Vapor Bubbles on a Boiling Surface," In: Proceedings of the 1963 Heat Transfer and Fluid Mechanics Institute, Stanford Univ. Press, pp. 222-235, 1963. 35. Skinner, L. A., "Vapor Bubble Growth on a Semi-Infinite Solid," Zeitschrift fur Angewandte Mathematik und Physik, 19, 6, 833-843 (1968). 36. Madhavan, S., Mesler, R., Lawrence, Kan., USA, "A Study of Vapor Bubble Growth on Surfaces," Paper B 2.6, Proceedings of Fourth International Heat Transfer Conference, Paris, 1970. 90

37. Guy, T. B., Ledwidge, T. J., "Numerical Approach to Non-Spherical Vapour Bubble Dynamics," Int. J. Heat Mass Transfer, 16, pp. 2393-2406 (1973). 38. Akiyama, M., "Bubble Collapse in Subcooled Boiling," Bulletin of the JSME 16, 93, 570-575 (Mar. 1973). 39. Shima, A. and Nakajima, K., "The Collapse of Non-Hemispherical Bubble Attached to a Solid Wall," J. of Fluid Mechanics, pp. 369-391, Vol. 80, Part 2 (April 1977). 40. Chesters, A. K., "Modes of Bubble'Growth in the Slow Formation Regime of Nucleate Pool Boiling," Int. J. Multi-Phase Flow 4, 3, 279-302 (Aug. 1978). 41. Winters, W. S., and Merte, H., "Experiments and Nonequilibrium Analysis of Pipe Blowdown," Nuclear Science and Engineering: 69, 411-429 (1979). 42. Shima, A. and Sato, Y., "The Collapse of a Bubble Attached to a Solid Wall," Ingenieur-Archiv, 48, 2, 85-95, (Apr. 1979). 43. Burow, P., "The Computation of Bubble Growth on Heating Surfaces at Pool Boiling Conditions as a Numerical Solution of the Conservation Equations," Warme und Stoffubertragung, Vol. 13, No. 4, pp. 253-264, Apr. 1980. 44. Tokarev, V. M., Bogovin, A. A., "Experimental Study of Dynamic Effects in Bubble Boiling of Liquids," J. of Heat Transfer, 34, 7, 1420-1423 (June 1980). 45. Koffman, L. D., Plesset, M. S., "Experimental Observations of the Microlayer in Vapor Bubble Growth on a Heated Solid," J. Heat Transfer, Vol. 105, pp. 625, Aug. 1983. 46. Witze, C. P., Schrock, V. E. and Chambre, P. L., "Flow About a Growing Sphere in Contact with a Plane Surface," Int. J. Heat Mass Transfer, Vol. 11, No. 11, pp. 1637-52, Nov. 1968. 47. Levkovskii, Yu. L., "Collapse of Spherical Gas-Filled Bubble Near Boundaries," Soviet Physics Acoustics, 20, 1, 36-38 (July/Aug. 1974) (Translation of Akusticheskii Zhurnal 20, 62-66 (Jan./Feb. 1974). 91

48. Lauterborn, W. & Bolle, H., "Experimental Investigations of Cavitation-Bubble Collapse in the Neighborhood of a Solid Boundary," J. Fluid Mech. 72, Part 2, 25, pp. 391-399, Nov. 1975. 49. Lenoir, M., "Numerical Calculation of the Implosion of a Cavitation Bubble Near a Solid Boundary or a Free Surface," J. de Mecanique, 15, 5, 725-751 (1976). 50. Plesset, M. S. & Chapman, R. B., "Collapse of an Initially Spherical Vapour Cavity in the Neighborhood of a Solid Boundary," J. Fluid Mech., Vol. 47, Pt. 2, pp. 283-290, 31 May 1971. 51. Kenning, D. B. R. & Cooper, M. G., "Flow Patterns Near Nuclei and the Initiation of Boiling During Forced Convection Heat Transfer," Proc. Instn. Mech. Engrs., Vol. 180, Pt. 3C, pp. 112-123, 1965-66. 52. Statnikov, Yu. G., "Microstreaming About a Gas Bubble in a Liquid," Soviet Physics-Acoustics, 13, 3, 398-399 (Jan/Mar. 1968). (Translation of Akusticheskii Zhurnal 13, 3, 464-466, (July/Sept. 1967). 53. Bieniasz, B., "Wetting Tensions Acting on a Vapor Bubble in Boiling Process," (in Polish), Mechanika Teoretyczna i Stosowana 9, 4, 529-540 (1971). 54. Kao, Y. S. & Kenning, D. B. R., "Thermocapillary Flow Near a Hemispherical Bubble on a Heated Wall," J. Fluid Mech. 53, 4, 715-735 + Figs. (June 1972). 55. Gaddis, E. S., "The Effects of Liquid Motion Induced by Phase Change and Thermocapillarity on the Thermal Equilibrium of a Vapor Bubble," (in English), Int. J. Heat Mass Transfer 15,, 11, 2241-2250 (Nov. 1972). 56. McGrew, J. L., Rehm, T. L. & Griskey, R. G., "The Effect of Temperature Induced Surface Tension Gradients on Bubble Mechanics," Applied Scientific Research 29, 3, 195-210 (June 1974). 57. Babskiy, V. G. & Sklovskaya, I. L., "Hydrodynamics in Weak Force Fields - Onset of Steady Thermocapillary Convection in a Spherical Fluid Film Under 92

Conditions of Weightlessness," NASA Technical Translation, NASA TTF-16668, 14 pp. (Dec. 1975). 58. Baranenko, V. I., Chichkan, L. A., "Thermocapillary Convection in the Boiling of Various Fluids," Heat Transfer-Soviet Research, 12, 2, pp. 40-44, (Mar/Apr 1980). 59. Cowley, S. J., and Davis, S. H., "Viscous Thermocapillary Convection at High Marangoni Number," J. of Fluid Mech., 135, pp. 175, Oct. 1983. 60. Kaznovskii, S. P. et al, "Reactive Forces Exerted on an Expanding Vapor Bubble by Evaporating Molecules," High Temperature, 14, 5, 894-900 (Mar. 1977). 61. Blangetti, F., and Naushahi, M. K., "Influence of Mass Transfer on the Momentum Transfer in Condensation and Evaporation Phenomena," I.J.H.M.T. 23, #12, pp. 1694, Dec. 1980. 62. Weinzierl, A., Straub, J., "Nucleate Pool Boiling in Microgravity Environment," Proc. of 7th Int. Heat Transfer Conf. Sept. 6-10, 1982, Munich. 63. Fedotkin, I. M., Konstantinov, S. M. & Tereshchenki, A.A., "Separation of a Vapor Bubble and Calculation of the Critical Bubble Radius," J. Eng. Phys. 24, 5, 589-592 (Feb. 1975). 64. Kirichenko, Yu, A., "Evaluation of the Conditions of Vapor Bubble Separation During Nucleate Boiling," J. Engineering Physics, 25, 1, 811-817 (Mar. 1975). 65. Malcotsis, G., "The Mechanism of Bubble Detachment from a Wall at Zero and Negative Gravity," J. Fluid Mech. 77, 2, 313-320 (Sept. 1976). 66. Kutateladze, S. S., and Gogonin, I.., "Growth Rate and Detachment Diameter of a Vapor Bubble in Free Convection Boiling of a Saturated Liquid," High Temperature, 17, 4, 667-671 (Jan. 1980). 67. Mitrovic, J., "Das Abreissen von Dampfblasen an festen Heizflachen," Int. J. Heat Mass Transfer, Vol. 26, No. 7, pp. 955-964, July 1983. 93

68. Prinsnyakov, V. F., "Breakoff of Vapor Bubbles from a Heating Surface," Journal of Engineering Physics 19, 5, 1435-1441, June 1973). 69. Cooper, M. G., "The 'Mirage' in Boiling," Int. J. Heat Mass Transfer, Vol. 26, NO. 7, pp. 1088-1089, July 1983. 70. Koumoutsos, N., Moissis, R. & Spyridonon, A., "A Study of Bubble Departure in Forced-Convection Boiling," J. Heat Transfer, Trans. ASME, Ser. C90, 2, 223 -230 (May 1968). 71. Anderson, G. H., Minns, D. E., London, England, "Nucleate Boiling in a Flowing Liquid," Paper B4.1, Proceedings of Fourth International Heat Transfer Conference, Paris, 1970. 72. Winterton, R. H. S., "Flow Boiling: Prediction of Bubble Departure," International Journal of Heat and Mass Transfer, Vol. 27, No. 8, pp. 1422, August, 1984. 73. Jansen, F., "Influence of Condensation on the Upward Motion of Steam Bubbles," Compt. Rend. Acad. Sci., Paris 258, 25, 6061-6064, June 1964. 74. Moelem, D. & Sideman, S., "The Effect of Motion on Bubble Collapse," Int. J. Heat Mass Transfer, 16, #12, pp. 2321-2329, Dec. 1973. 75. Srebnyuk, S. M., Gorban, V. A., "Concerning the Induced Mass of a System of Bubbles," Fluids Mechanics-Soviet Research 8, 3, 93-102 (May/June 1979). 76. Drew, D., Cheng, L., et al, "The Analysis of Virtual Mass Effects in Two-Phase Flow," Int. J. of Multiphase Flow, 5, 4, 233-242, Aug. 1979. 77. Bratukhin, Yu. K., Evdokimova, O. A., "Motion of Gas Bubbles in a Nonuniformly Heated Liquid," Fluid Dynamics 14, 5, 679-681 (Mar. 1980). (Translated by Consultants Bureau). 78. Kuznetsov, V. M., Lugovtsov, B. A. and Sher, E. I., "On the Motion of Gas Bubbles in a Fluid Under the Effect of a Temperature Gradient" (in 94

Russian),PMTF: Zhurnal Prikladnoi Mekhaniki i teknicheskoi Fiziki, No. 1, 124 -126 (1966). 79. Van Wijngaarden, L., "Hydrodynamic Interaction Between Gas Bubbles in Liquid," J. Fluid Mech. 77, Part 1, 9, pp. 27-44, Sept. 1976. 80. Bhaga, D., Weber, M. E., "In Line Interaction of a Pair of Bubbles in a Viscous Liquid," Chemical Engineering Science, 35, 12, pp. 2467-2474 (1980). 81. Lee, S. H. and Le Al, L. G., "Motion of a Sphere in the Presence of a Plane Interface. Part 2, An Exact Solution in Bipolar Co-ordinates," J. Fluid Mech., Vol. 98 (1) 193-224 (May 1980). 82. Hatch, M. R. and Jacobs, R. B., "Prediction of Pressure-Drop in Two-Phase Single-Component Fluid Flow," AIChE J. 8, 1, 18-25, March 1962. 83. Kotake, S., "Heat Transfer and Skin Friction of a Phase-Changing Interface of Gas-Liquid Laminar Flow, Int. J. Heat Mass Transfer, 16, #12, pp. 2165-2176, Dec. 1973. 84. Davis, M. R., "The Determination of Wall Friction for Vertical and Horizontal Two-Phase Bubbly Flows," J. Fluids Engineering, Trans. ASME, 61, #2, pp. 173-179, June 1974. 85. Wisman, R., "Analytical Pressure Drop Correlation for Adiabatic Vertical TwoPhase Flow," Applied Scientific Research 30, 5, 367-380 (March 1975). 86. Kopalinsky, E. M. & Bryant, R. A. A., "Friction Coefficients for Bubbly TwoPhase Flow in Horizontal Pipes," AIChE J 22, 82, (1976). 87. Moussalli, G. & Chawla, J. M., "Void Fraction and Pressure Drop in Bubble Flow," (German) Forchung im Ingenieurwesen 42, 5, 149-153 (1976). 88. Friedel, L., "Pressure Drop in Flow of Gas/Vapor-Liquid Mixtures in Pipes," Chemie Ingenieur Technik, 50, 3, 167-180 (Mar. 1978). 89 Ishii, M. and Zuber, N., "Drag Coefficient and Relative Velocity in Bubbly Droplet or Particulate Flows," AIChE J. 25, 5, 843-855 (Sept. 1979). 95

90. de Socio, L. M., Gaffuri, G., "Analytical Results on Two-Phase Flows in Ducts," Warme und Stoffubertragung, Vol. 14, No. 3, pp. 183-188, Oct. 1980. 91. Van der Welle, R., "A Numerical Procedure for Vertical One-Dimensional TwoPhase Flow," Nuclear Engineering and Design, Vol. 65, No. 1, pp. 113-130, May 1981. 92. Cooper, M. G., Mori, K., and Stone, C. R., "Behavior of Vapour Bubbles Growing at a Wall with Forced Flow," Int. J. Heat Mass Transfer, Vol. 26, No. 10, pp. 1489-1507, 1983. 93. Papell, S. S., "Buoyancy Effects on Forced-Convective Boiling," Paper 67-HT63, ASME, August 1967. 94. Cochran, Thomas H., "Forced-Convection Boiling Near Inception in Zero Gravity," NASA Report No. TN D-5612, 1970. 95. Chen, L. T., and Chuang, W. C., "The Effect of Orientation on R 1 Bubble Frequencies in Nucleate Pool Boiling," Letters Heat Mass Transfer, 6, pp. 429 -438, 1979. 96. Simoneau, R. J. and Simon, F. F., "A Visual Study of Velocity and Buoyancy Effects on Boiling Nitrogen," NASA TN-D-3354, Sept. 1966. 97. Thorsen, R. S., Dobran, F., and Alcorta, J. A., "A Comparative Study of Vertical Upflow and Downflow in a Uniformly Heated Boiling Liquid," Proc. 5th Int. Heat Transfer Conf.., Tokyo, Paper B.4.3, pp. 155-159, 1974. 98. Bartolini R., Guglielmini, G., and Nannei, E., "Experimental Study on Nucleate Boiling of Water in Vertical Upflow and Downflow," Int. J. Multiphase Flow, Vol. 9(2), pp. 161-165, 1983. 99. Hahne, E., Shen, N., and Spindler, K., "Fully Developed Nucleate Boiling in Upflow and Downflow," Int. J. Heat Mass Transfer, Vol. 32, No. 10, pp. 1799 -1808, 1989. 96

100. Gunther, F. C., and Kreith, F., "A Photographic Study of Bubble Formation in Heat Transfer to Subcooled Water," Progress Report No. 4-120, J. P. L. California Institute of Technology, Pasadena, California, 1950. 101. Clark, J. A., and Rohsenow, W. M., "A Study of the Mechanism of Boiling Heat Transfer," Trans. ASME, 73, p. 609, 1951. 102. Bankoff, S. G., and Mikesell, R. D., "Bubble Growth Rates in Highly Subcooled Nucleate Boiling," Preprint 5 presented at AIChE-ASME Heat Transfer Conference, 1958, Chem. Engng. Prog. Symp. Series No. 29, pp. 95-109, 1959. 103. Zuber, Novak, "Recent Trends in Boiling Heat Transfer Research, Part I: Nucleate Pool Boiling," Applied Mechanics Reviews, Vol. 17, No. 9, pp. 663-672, (Sept. 1964). 104. Del Valle M., Victor H., and Kenning, D. B. R., "Subcooled Flow Boiling at High Heat Flux," Int. J. Heat Mass Transfer, Vol. 28, No. 10, pp. 1907-1920, 1985. 105. Lemmert, M. and Chawla, J. M., "Influence of Flow Velocity on Surface Boiling Heat Transfer Coefficient," Heat Transfer in Boiling, edited by E. Hahne and U. Grigull, Academic Press, pp. 237-247, 1977. 106. Chen, J. C., "A Correlation for Boiling Heat Transfer to Saturated Fluids in Convective Flow," Ind. Engng. Chem Process Des. Dev., 5(3), pp. 322-329, 1966. 107. Bjorge, W., Hall, R., and Rohsenow, M., "Correlation of Forced Convection Boiling Heat Transfer Data," Int. J. Heat Mass Transfer, Vol. 25, No. 6, pp. 753 -757, 1982. 108. Colburn, A. P., Trans. Am. Inst. Chem. Engng., Vol. 29, p. 174, 1933. 109. Mikic, B. B., and Rohsenow, W, M., "A New Correlation of Pool Boiling Data Including the Effect of Heating Surface Characteristics," Trans ASME, Series C., J. Heat Transfer, Vol. 91, p. 245, 1969. 97

110. Bjorge, R. W., "Recent Developments in Boiling Heat Transfer," SM Thesis in Mechanical Engineering, MIT, 1979. 111. Shah, M. M., "Chart Correlation for Saturated Boiling Heat Transfer; Equations and Further Study," ASHRAE Trans., Vol. 88(1), pp. 185-196, 1982. 112. Gungor, K. E., and Winterton, R. H. S., "Simplified General Correlation for Saturated Flow Boiling and Comparison of Correlation with Data," Chem. Engng. Res. Des. Vol. 65, pp. 148-156, 1987. 113. Arpaci, V. S., and Larsen, P. S., "Convection Heat Transfer," Prentice-Hall, Inc., Englewood Cliffs, N. J., 1984. 114. Kline, S. J., and McClintock, F. A., "The Description of Uncertainties in Single Sample Experiments," Mechanical Engineering, Vol. 75, 1953. 98

APPENDICES Appendix A. Finite Difference Solution to Mixed Convection in a Rectangular Duct. B. The Velocity Potential of Two Spheres Moving Perpendicular to the Line Joining Their Centers and the Associated Lift Force C. Uncertainty Analysis D. Tabulation of Forced Convection Transient Incipient Boiling Data 99

Appendix A Finite Difference Solution to Mixed Convection in a Rectangular Duct by Kevin M. Kirk December 19, 1990 A two-dimensional transient finite difference solution is presented for the problem of laminar mixed convection in a rectangular duct. The incoming flow to the duct, which can be specified as being either uniform flow or fully developed poiseuille flow, passes over a heated patch subjected to a constant heat flux. The problem is solved for varying orientations but with particular attention given to vertical upflow, vertical downflow, and horizontal flow.

Table of Contents I. Introduction II. Problem Formulation A. Dimensional Governing Equations B. Boundary and Initial Conditions C. Non-Dimensionalization of Equations III. Method of Solution IV. Numerical Results and Discussion V. Program Listing VI. References 2

I. Introduction The presence of gravity can give rise to significant buoyancy forces in low speed flows where behavior predicted solely from the consideration of forced convection can lead to erroneous results. Furthermore, the direction of gravity relative to the heating surface can have significant influence on the importance of gravity in natural convection. Nyugen et al. (1979) studied mixed convection for the case of horizontal flow between parallel plates with isothermal walls and nonuniform velocity and temperature profiles. Their results demonstrated that heat transfer above a heated place was enhanced while that below a heated plate decreased with respect to purely forced flows. Other analytical investigations studying the effect of buoyancy on fully developed laminar flow between parallel planes (Gill and Del Casal, 1962) have showed that a separated recirculating flow can occur at the upper wall. Kennedy and Zebib (1982) presented numerical and experimental work on the effects of free convection on forced laminar flow with a local heat source on the lower wall. Their results demonstrated that a recirculation region adjacent to the top wall and above the heat source was generated due to the effects of natural convection. A numerical study by Davis (1976) examining the effect of orientation on mixed convection over a heated patch of constant temperature predicted recirculation for the case of vertical downflow. The present numerical study was done with the assumption of constant heat flux on a small heated patch. The heat transfer characteristics and flow patterns were studied for four different orientations where particular attention is given to the occurrence of flow recirculation since it is well-known that an inflection in velocity can give rise to instabilities. A study of mixed convection also has relevance in the study of incipient forced convection boiling for the case of low speed flows. For example, when a sufficiently high heat flux is applied to a surface, boiling will commence when the surface has locally attained the required superheat for boiling inception. Knowledge of the temperature distribution in the vicinity of the heater surface is critical to understanding the initial nucleation and subsequent spreading phenomenon. Hence, if buoyancy is assumed significant, a transient mixed convection solution is necessary and such a solution is given numerically in this work for laminar flow. 3

II. Problem Formulation A. Dimensional Governing Equations The geometry of the problem is illustrated in Figure 1 with the inlet flow being directed in the +x direction and horizontal flow with the surface facing upwards occuring when 0 = 900. < b '->b Inlet Flow L Y.......T....x. 9 $ T T T T ///////////////////////// q"G =90~ Figure 1 - Two dimensional problem geometry The standard problem of two-dimensional mixed convection in a rectangular field containing a Newtonian fluid for which the Boussinesq approximation may be applied can be described by the following dimensional equations: Continuity au av=o ax ay (1) Momentum in the x-direction au au au _ a U a2 a+u" + va- x1 g(T- To) cosO+ I L x2 + y 21 at ax ay poax ~ ' po ax2 ay2J (2) 4

Momentum in the y-direction av Dv av ~3v 1 a2vP + UT +V y +v- Po= L + gP(T - To) sine + 2 + —2 at ax ay PoaY Poa}x2 ay2 (3) Energy aT aT aT a2T a2T1 t a+x a y ax2 ay2J (4) Alternatively, the problem can be formulated in terms of the stream function and vorticity. This approach ensures automatic satisfaction of the continuity equation, eliminates the pressure gradient term and reduces the number of differential equations to be solved by one. The stream function, A, and vorticity, co, are defined as U "av V ay 3x av au _ Ja2 +E axay Lax2 ay2 (5) Differentiating equation (2) with respect to y and equation (3) with respect to x and, then, subtracting (3) from (2) yields the following equation after making appropriate substitutions for vorticity and the stream function ao aaoa ( vO o avaO aT 9\ La2 a2eo + a + ac - -ava = g3(-sine + -i-cos4+ 4 - + 1 at ayax axay ax ay polx2 ay2j (6) After substituting the stream function into the energy equation for the velocity, the energy equation becomes aT aT nT avaT a2T a2T1 at ayax axay lax2 ay2j (7) Equations (5), (6), and (7) form a set of three coupled non-linear differential equations to be solved numerically with dependent variables T, o, and V. 5

B. Boundary and Initial Conditions From the governing equations, it is observed that one condition for each of the variables, T, o, and V, are needed on each boundary as well as an initial condition for each variable. Boundary conditions on V at y=O and y=L are obtained from the no slip and no penetration conditions. That is, on a solid boundary, IV = constant, and "/l = 0, where '1 is the direction normal to the solid surface. The boundary condition for V at the inlet (i.e. x=O) is obtained from specifying either poiseuille or uniform flow as the inlet flow condition. At the exit, it is assumed that the curvature of the streamlines is negligible. Hence, at the exit plane — (Xexit,y,t) = 0 Dx2 Boundary conditions for vorticity, o, must be expressed in terms of the stream function. For uniform flow at the inlet, x=0, vorticity is allowed to diffuse upstream. Thus, applying equation (5) at the inlet we get a2w ax2 For the inlet condition of poiseuille flow, the boundary condition at x=0 is obtained simply by applying equation (5) at the boundary. Along the wall, y=0, the flow is assumed to be essential parallel, hence, aV/x =0 and a/ax =0. From these assumptions, the following second-order Woods formula can be obtained through application of equations (5) and (6) at an insulated wall for steady flow 3(Ni,1- Vi,2) Oi,= 2'-i,2 + (Ay)2 Finally, at the exit plane since the curvature of the streamlines is considered negligible, equation (5) reduces to a2N ay2 6

The boundary conditions for temperature are a constant inlet temperature, adiabatic wall surfaces, a constant heat flux applied at the heater surface, and negligible curvature of the isotherms at the exit. The initial and boundary conditions are summarized as follows: at t=O T(x,y,O) = T, for uniform flow s(x,y,O) = U,y o(x,y,0) = 0 for poiseuille flow V(x,y,O) = 6U4j2y 3-L2] o(x,y,0) - -y = - - at x=O T(O,y,t) = T. for poiseuille flow [y2 y31 V(O,y,t) = 6U r 3L2 c(0O,y,t) = -I- -1 --- \ax2 ay2/ for uniform flow 4(0,y,t) = U.y as2 O(0,y,t) = - 'x2 at x=xexit ay X2(exit,y,t) = 0 ax2 a2T -(xexit,y,t) = 0 ax2 (O(Xexit.y,t) = 2 ay2 7

at y=O x(x,O,t) = 0 (x,O,t) = 0 ay to(x,0,t) =- 3( i +2) 2 i,2 (Ay)2 on heater surface -IqT (x,Ot) ay kf at y=L x,L;Lt) = 0 ay W(xL,t) = U.L co(x,L,t) = - -0 + M- ) 2 i,M-1 (Ay)2 C. Non-Dimensionalization of Equations It is usual to put the equations into dimensionless form to reduce the number of independent parameters and, hence, the computational effort. Non-dimensionalization also offers opportunities for enhanced physical understanding through the representation of the relative forces into dimensionless parameters. In order to make the equations (5) thru (7) dimensionless, the following dimensionless variables are defined where the characteristic velocity is the mean flow velocity, UO, and the characteristic length is the heater length, 1. =x=A y =y/ A o = to(WU4) T - To, = V/(U.1) 0 = t (q"l/k) 12 Substituting the dimensionless variables into the governing equations, the following dimensionless equations are obtained le '+ aw = Ri _os9 +a 9 ine + 1 co () -+ a - ---— x - = Ri cos 0+ -sinO + I a2i ( ~Pe 3tBBy \9B /~ Re \2l yl2j (8) 8

Lao + la ao av ao a Pe ar ay ax Ay Pe \a2 ay2L (9) co = -~ +! a\X2 ay2/ (10) where Ri =Gr = gl (q"Vk) Re2 U2 U.. Re Pe = 09 xa The dimensionless initial conditions become at t=O (x,y,O) = O for uniform flow V(x,y,0) = y o(CXy0) = 0 for poiseuille flow =[21 Y312 M2L 3L2 -and the dimensionless boundary conditions are illustrated in Figure and the dimensionless boundary conditions are illustrated in Figure 2 9

For poiseuille flow: ( o,yt) - a- LaX2 as2J For uniform flow: (t a2= ~ ao2 Vo(,yt) = (',L/lt) = L/l (,LlI) = - i +3(/i l' i, M-1) L = iM,1/-1 + 23(^n (ar) Alb,\~,\-,~\\~\~\~\\\\ T L i -> w(ri.7,,t) = - - — >,.,t) ='o - ae k-1 -1 o(~,O,lt) = - i,2. 3((Vi,W - 'i,2) =~ - ~r2= +(A$)... V<, (A o V(rrO,t) = O a) Vorticity and stream function boundary conditions ae ~_,,\,,,,,,,,,,,,,,,a,,,v, (o,y,t) = o L l - a2e 4-> 2exi,,t,) = o -> a.2 k-1 -' t,\\\\\~\\\ ae o) = IltJlOot) = 0 a.7 ae,-<O,) = -1 ay b) Temperature boundary conditions Figure 2- Dimensionless Boundary Conditions 10

III. Method of Solution The second order derivatives were approximated using central differences. However, approximating the convection terms in the governing equations with central differences suffers from the well known Reynold's cell number restriction. This instability may be overcome with the use of upwind differencing, instead of central differencing, since diagonal dominance is ensured. However, it has been recognized (Davis, 1974) that the first order upwind difference approximation can generate significant truncation errors which may give rise to a "false diffusion," especially in the case of high Reynolds number flow where physical diffusion is weak. Nonetheless to assure convergence over a wide range of Reynolds number an upwind differencing scheme was used where the local stream direction was determined from the the stream function values around the point in question. The equations were solved using an alternating direction implicit method. The method employs two difference equations which are used in turn over successive time steps, each of duration At/2. The first equation is implicit only in the x-direction, while the second equation is implicit only in the y-direction. The advantage of, this method is that this technique requires only the solution of a tridiagonal matrix which affords a straightforward solution. A stability analysis was performed experimentally through modification of the mesh size and time step until convergence was achieved for a given set of dimensionless parameters. 11

IV. Results and Discussion Representative contour plots are given of the isotherms and the stream function in Figures 3 thru 9 for different orientations and values of Ri, and Re. In all the figures it is assumed that Pr=l, hence, Re = Pe. Low values of Re were used in the representative contour plots in order to accentuate the effects of buoyancy. It is seen from the contour plots of vertical upflow in Figure 6 that the buoyancy acts to accelerate the flow near the heated surface thereby enhancing heat transfer. The acceleration of the flow is illustrated by the increased density of the streamlines near the heating surface, and is accompanied by an increase in the mean Nusselt number across the heater as shown in Figure 10. In the contour plot of vertical downflow, shown in Figure 4, the flow is decelerated due to the buoyant force acting up on the heated fluid near the heater surface, and the mean Nusselt number across the surface decreases. If the Richardson number is sufficiently large in vertical downflow, the numerical solution predicts flow reversal and circulation where the point at which the onset of circulation occurs is shown in Figure 10. Inlet Poiseuille Flow Vertical Orientation Pr = 1 5 Re=10 Vertical Down 4 * Re=30 Vertical Down z U A. ----" " 0 Re=100 Vertical Down i 3~i Kr*"-\ | Re=10 VerticalUp * Re=30 Vertical Up 0 3 Re =100 Vertical Up 2 _ Onset of Circulation 0 10 20 30 40 Ri Figure 10 - Mean Nu vs Ri for Vertical Flow 12

Without circulation, the steady-state maximum temperature in vertical downflow occurs at the downstream edge of the heater surface, however, when circulation does occur, the steady-state maximum temperature is located at the upstream edge of the heater surface. It is interesting to note that for a range of Re from 0 to 100, it is shown that the value of Ri at which the onset of circulation occurs seems to be solely a function of Re for a given Pr, and is defined below as the critical Richardson number, Ricr. In other words, as Re increases, less buoyancy is required for the onset of recirculation. This trend has also been predicted by Davis (1976). This relationship is shown graphically in Figure 11, and is given by the equation Ricr = 48.98Re"0.368 or in terms of the Grashof number Grcr = 48.98Re1'632 100:. 10 U 1 ------ --------------- ------ ---- -------------- ------- -------------------- --- ---------- --- - --------------- ----------- -------- ------- ------------ -------- ------ -------------- ---- ------- - ----------- --------- -- ---------- ---------- - ---- --- ------- ----- - -------------- ------------- --------- -------- ------- ----- ----- -- Vertical Downflow Inlet Poiseuille Flow Pr= 1 * Critical Ri 10 Re 100 Figure 11 - Critical Ri vs Re for Vertical Flow In horizontal flow with the surface facing up, the streamlines initially bend toward the heated surface as shown in Figure 8 indicating that the flow 13

upstream of the surface is drawn toward the surface. However, at near the mid-point of the heater surface the flow begins to move away from the surface due to the influence of buoyancy. Conversely, in horizontal flow with the surface facing down, illustrated in Figure 10, the streamlines initially bend away from the surface and, again, at near the mid-point of the heater the flow is drawn back toward the heater side due to buoyancy effects. As expected, the highest temperature on the heater surface in horizontal flow was always at the downstream edge of the surface as seen in the contour plots of the isotherms, Figures 7 and 9. Furthermore, a higher maximum temperature was attained with the surface facing down than with the surface facing up in horizontal flow with identical values of Ri, Re, and Pe. The higher temperature is due to the degradation of the heat transfer with the surface facing down as illustrated by the decrease in the Nusselt number with increasing Richardson number in Figure 12. On the other hand, buoyancy enhanced heat transfer when the surface was facing up in horizontal flow which is also demonstrated in Figure 12. Recirculating flow is also predicted for the case of horizontal flow with the surface facing downward for large values of Ri, however, whether or not this represents reality remains a question left to experimentation. 6 Inlet Poiseuille Flow Horizontal Orientation Pr 1 5 z I 0 Re=30 Horiz. Up ~ |C~t 1/ | * Re=50 Horiz. Up: 1 ' l "Y/,"a|i Re=1l00Horiz.Up ~ Re=30 Horiz. Down * Re=50 Horiz. Down 3 a gAd -r Re=100 Horiz. Down * I; - - -- - -- O nset of C irculation 2 0 25 50 75 100 125 150 175 200 Ri Figure 12 - Mean Nu vs Ri for Horizontal Flow 14

The effect of buoyancy on heat transfer at orientations other than horizontal and vertical orientations is given in Figure 13 for a constant Reynold's number. ~~~~~~~~~~5 - ~4,5 3 Re = 100 4 8 Pe =100.^"^ ^ ^ ^ ^ ^ ^ 21 Nu-0 2 * Nu-45 3 \ 3 3 NU-90 ~. \4 o Nu-135 Z~3 2 2~~~ \ |1~5 * Nu-180 ~~z 2 - ~~~~~~~\ 6 a Nu-225 \ 7 A Nu-270 8 a Nu-80and238 1 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 RI Figure 13 - Mean Nu vs. Ri at Varying Orientations As 0 increases from 00 (vertical downflow), there is a considerable decrease in the tendency for the buoyant force to degrade the heat transfer. This trend continues until 0 equals 41,i where Ri has little influence on Nu and can be approximated by a horizontal line as seen on Fig. 13. By interpolation, the angle of {1 has been determined as approximately 80 degrees. As 0 increases past (1 the buoyant forces no longer degrade the heat transfer from the plate but instead help to enhance it by increasing Nu. The greatest value for Nu occurs between 0 = 135 and 0 = 180 since components of buoyant forces and forced convection are in the same directions for upward flow. As 0 rotates from 180 to 270 degrees there once again begins a degradation of heat transfer as the buoyant forces become more directed toward the heating surface. This ultimately results in warm fluid remaining closer to the heating surface and a decrease in heat transfer. 15

Ri =40 Re = 10 Pe = 10 r -- -- 7 - r - r -I- -I - - r -I- -- 7 - r - r -- - 7 - r - r -I- 7- 7 - r - L - I I I I I I I I I I I I I i I I I I I I I I I I I L _I I_ I I 1 _ I I_ I I _ 1 _ I I I I I _ I I _ I I I I I I_ _ I i 1 I I I I I I I I I I I! I I I I I I I I I -!- t- -- - - - - - — 1 — -I- — 1 — -. -- - - -- -4I- - -- 4 - 4 — — I- — 1- - - _I I I I I I I I I I I I I I I I I I I I I / I I I I! I I I I I I I I I / 1 I I I I I I I I I I I I I I, 1- - - t- I -I I -I - I- I- -II I I I - - I I I I I I I I I I 1 / 1 I I I I I I I I I I I I I I I I I I I \ I I I I I ~ a x ~... —r" i i i i i '"+.,. i t I I I I I I I rI- - -- - + -- - t- - 1 t r e- - i- -I- - - — I- t- - - r- - -I - - - II /_//I.1C Figure 3 - Dimensionless Isotherms Vertical Downflow

Ri = 40 Re = 10 Pe = 10 I '- -1 i - I- - - I I - -7 I - m — ' -T - I Ii I- - - I 1l I I I I I I I I I I I I I I I I I I I I I I I I 1. I I I I I I I I I I I I I I I I I I I I I I I I I- - -T -r - I I I I I I r- I I I 1I1 iI I-I I I I I I I I 1 I l I I 01. i I I III I I I T I I I I I ItL Figure 4 - Dimensionless Stream Function Vertical Downflow

Ri = 40, Re = 10 Pe = 10. r --- —r-7 T r- — n —T7 — r- — i —7-T7r- -I-i7-rT r —I —r-i L —1 —1- T-r- -— T ---1I. -1- -— T — ] --- — '1- — T-] — - r ----1 — q ---T - -— 1 I I I t - t - r-I I I I I- I I I I I I I I I I I I -— l —l-4-4 —L_ ---- - -4-_- -___-_J --- -- 4 —t —l —4-4 ---__- J___4__4 I -! - --- - - - -- - - - --- - - - - - -- - - - --- - -- -- -- - - I I I I I I I I I I I I I I I I I I I I I I I I.OE L I_ I 1_ _ __ _1 L I_ I IJ I I I I I Figure 5 - Dimensionless Isotherms Vertical Upflow

Ri = 40 Re = 10 Pe = 10 I - I I I I I I I I - I - I I I I I I - I I I I I- I t I _ I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I _ I _ I _ I I I I _ I I __I I _ 1I _ I I L I I I I I I I ' - I- I I I I I I I I I I I I I I I I I I I I I I T I I I I I I I I. I I I I I I I I I II I! I I I — LI —I- - -- - L' - — I- - I — i~~~~~~~~~~- -. _ _ _ - I I- I I - I I I I I I T.-I- - i _ I I I I I I I I I L -1-.-1- -+-F- I _ I _I _I _ __I I _ I _ I I I 1 I_ I I L_ I I I I F I I I Stream Function -— I ---I-1 -' ---.- -1 —I -I-r —I —I —I-'-T-1 -I I I I I I I I I I I I I I I I I I I I I I ~ I I I I _ I _ I _ I I I I! I I I I I I I I I I I I I I I I I I I I I I I I I I Figure 6 - Dimensionless Stream Function Vertical Upf low

Ri = 100 Re =30 Pe = 30 r-I — - - - r --— I- -r-r --- -7-r-r-r -- - — T — r —I —I- - -TI I I I I I I I I I I I I I I I I I I I I I I I --- - - 4- --- I- -4 - 4 — -. -. - 4 — -- I4 -4- ----—. 4- 4- - I I I I I t II I I I I I I I I I I I I I I I I I I I I I I I I I II -I ---- -I- ---- I- t- - -- - -1- --- 1- - I- --- - --- - t — t- -1 - - -- I i I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I -I -I I I I I I — I I I I I I I I I I I I I I I I I I I I I I i I I I I - -- I — + —4 ---------- — + —I --- ----- ----- I —H- - ---— I-H —+ — I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I -I-~ ~ ~ ~ ~ ~ ~~- - - F-.14.16.18 - -.12.10.08- - I I I I I I I I I I I I I I I I I I I I I I Figure 7 _- Dimensionless Isotherms Horizontal Flow with Surface Facing Up I I I I I I I I I I I I I I I I II I I I I I I I I I I I I L t.14.16.18 I.12.10.08 Figure 7 -Dimensionless Isotherms Horizontal Flow with Surface Facing Up

Ri = 100 Re =30 Pe = 30 - 7 -r - - - -FT - 7 - - m - TT - — 7 - T1 I I I I I I I I I I I I I I I I I l L-. LI ---_ — l-L _ —_.. -L -— 1 ---l —I-._-_-_ ---I ---- _ ---. L I I i i I I I I I i I i, I I I I I I I I I I I I I I I I I I I I L-" —l.J.-. |._|- - - -. - -L - -L- -.- - -- -—. - -- - L- — L -.- 4 ir- - -I-7 - - iI- -- 7 - r — -I- -— \ —\ —' —I — ---- -1.- -7 ^-I - -I ~ t ~ T "t-I -I I- - r - I I II I- I I r - I —I —i- I I I I I I __ t i ___ i i l- 1 -1 r ' 1 i l 1 rl l l _: —i____' ----!-\ —___l:___i iTl __l —!4 Figure 8 - Dimensionless Stream Function Horizontal Flow with Surface Facing Up

I I I I I - -- -I - - I I I I I - T - - I - -- - -1 -- -- - T - - F F IL__l_ I I I I I- - -I- - 4 - - I I I -I- - - - - IP I _ "X1 oQ I D: 0U o H H C H. r' 0 U):3 H. o r 0 ~ v3 D 0l U) 0 U ) H~ )t Cr) (D 1- 3 C) 1-+. 0z) O L- - -I- - -I I --- I- -- I I r - - - - n - L - -I- - -t 1 — -I- - 4 - I I I -- - I - -4 -- - I- - -I - 1 I 1 L 1. 1 _ — _ 4 -- - I I 1 — I —II I I I I I I I I -1- I- +- -I I I _ I -l I I I I I I I I I I I I _ +- - _ I I I I I I I I 1- -L I I - T- -r I I ---- I 11 I _ _ I _ _I I I I -I- - 4 - - - I I I I t I I I I -I ~ I I_ _ _ _ I I I I -I - - - - I I I -I- -t -- - I I I -I- - 7 -- I I I -I - I I _ _ I _ _ I I -I- - 4 - - - -I- -t - - t I I I -I - - - - I I I I _ _ _ _ L I I I -I- - - - 1 --I- - -- I I I -I- - 4 - - F - -II - -I- -I - -II I - -II I I I P - - H. I I I I C) l - 4 I I D I _ J I I - -1 I 0 - -I I I I I I -I- - 4 I I - -I- -II - -I- -1 - - -1 I -I I - 1 I -4 - I I -I- - - I I - - I - -I- -1 -I _ 1 _ - L I I — I — -- I I I I I III l I I I -I- - 4 - -I- - I I I -I - - - - - - -II I. I-. I — I I I -1 -4 I - -I I o o o

Ri = 100 Re =30 Pe = 30 m_-m - -1 m LJ^J.-1-7- rn7 mr —7 T — m -— 7-rn F- F1 I I 1_ I i 1 I I I 1 II I I I I I I I I L I I I I _ I I I! I I I I i I I I I_ I I I I I I I I I I I I I -I I I I I I I I I I I I - -I —41- ( -- - 4- IT I -f - I- -I- -1-( - - -— I-4-4 —1 - I - I I I I I I I I I I I I I I I i! I I I I I I I I i I I I I I I I I I I I I I I I I I I I I I I rL. -A - - L- - L -I ----. - A' - -' - L -- -- - 4 — - - - L 1. - I I I I I I I__ _ I i _ _ __ _ _ _ i i i I I I I LI I I I I I I I I I 1 I _ 1 I _ L I I I I I I I I I I I I I I I I I I I I I I I I - I I I I _ I _II I. I I I I I I I _ I I I I 1 I I_ ~ 1 1 I I I I 1 I I I I I I I I I I III I I I I I I I I I I I I I I I I I I I I I I I II I I I 1 I I I I I 1 I I I I 1 I 1 I I I I I 1 t I I _X _ _ I I I I I I I I I I I I I I I I I I I I I I I I I I I_ i I I I I _ _ - I, 1 I I I _ I I I _ I I Figure 10 - Dimensionless Stream Function Horizontal Flow with Surface Facing Down

V. Program Listing C PROGRAM TRANSADI C C BY KEVIN M KIRK THE GREAT C C This program solves for the transient temperature and velocities C fields in the presence of combined natural and forced C convection in a rectangular duct. The problem is solved C for the case of either slug or poiseuille flow over a heated patch C subjected to a constant heat flux. C C The program has the capability to solve C the problem for varying orientations. That is, gravity can be C perpendicular or parallel to the heating surface. C C DECLARATION OF VARIABLES REAL*8 PSI(150,150),W(150,150),TH(150,150),A,B,C,DX REAL*8 Q,UINF,G,TC,DIFFU,MU,RI,RE,PE,L,PR,UP,UM,VP,VM REAL*8 E1 (10000),PSIOLD(150,150),DIFF1,DIFF2,TINF,NU,NUT REAL*8 F,DT,A 1(10000),B 1(10000),C1 (10000),D 1 (10000),NUX REAL*8 TOLD(150,150),ALPHA,PI INTEGER N,Nl,N2,M,IJ,SELECT,ITER,ORIENT,COUNT,NSTEP INTEGER NUM,K LOGICAL CONVERG C C PSI=DIMENSIONLESS STREAM FUNCTION C W = DIMENSIONLESS VORTICITY C TH = DIMENSIONLESS TEMPERATURE C DX = MESH SIZE C TINF = TEMPERATURE OF FLUID AT INLET C Q = IMPOSED HEAT FLUX ON SURFACE C UNIF = INLET FLUID VELOCITY C G = GRAVITY C TC = THERMAL CONDUCTIVITY C DIFFU = THERMAL DIFFUSIVITY C MU = KINEMATIC VISCOSITY C RI = RICHARDSON NUMBER C RE = REYNOLD'S NUMBER C PE = PECLET NUMBER C N = NUMBER OF NODES IN X DIRECTION 24

C M = NUMBER OF NODES IN Y DIRECTION C A = LENGTH OF CHANNEL C B = HEIGHTH OF CHANNEL C C = ENTRY LENGTH C ALPHA=ANGLE OF INCLINATION WITH VERTICAL C C Define default properties for freon (R-113) C DIFFU=4.5E-6 MU=4.63E-7 TC=0.066 PI=3.14159 C C Initialize counters and time. Open output files. C NCOUNT=0 TIME=0.0 OPEN(UNIT=4,FILE='theta.dat',STATUS='NEW') OPEN(UNIT=5,FILE='stream.dat',STATUS='NEW') C C Input dimensions of duct and time step C WRITE(*,200) READ(*,*)A,B,C WRITE(*,214) READ(*,*)L,DX WRITE(*,21 1) READ(*,*)DT,STOP,PSTEP NPRINT=NINT(PSTEP/DT) NSTEP=NINT(0.2/DX) N=NINT(A/DX)+1 M=NINT(B/DX)+1 WRITE(4,300) WRITE(5,295) N1=NINT(C/DX)+1 N2=NINT((C+L)/DX)+1 NUM=(N-2)*(M-2) I2=(N1+N2)/2 C C CONVERT TO METERS C A=A/100.0 B=B/100.0 25

C=C/100.0 L=L/100.0 DX=DX/100.0 WRITE(*,202) WRITE(*,203) WRITE(*,204) READ(*,*)SELECT WRITE(*,208) READ(*,*)ALPHA WRITE(*,209) READ(*,*)NTYPE WRITE(*,212) WRITE(*,213) READ(*,*)MODE IF (SELECT.EQ.2) THEN WRITE(*,210) READ(*,*)TINF,BETA,Q WRITE(*,220) READ(*,*)UINF,G C WRITE(*,230) C READ(*,*)TC,DIFFU C WRITE(*,240) C READ(*,*)MU C C MAKE UNITS CONSISTENT IN SI UNITS C Q=Q*(100.0**2) UINF=UINF/100.0 G=G*9.8 RI=(G*Q*(L**2))*BETA TEMP=((UINF**2)*TC) RI=RI/TEMP RE=UINF*L/MU PR=MU/DIFFU PE=RE*PR WRITE(*,215)RI,RE,PE ELSE WRITE(*,250) READ(*,*)RI,RE,PE ENDIF WRITE(*,256)N,N1,N2,M,NSTEP C C SET UP FILE FOR PLOTTING PROGRAM 26

C WRITE(4,257)((N-1 )/NSTEP)+ 1,((M-1 )/NSTEP)+1 WRITE(5,257)((N- 1)/NSTEP)+ 1,((M- 1 )/NSTEP)+1 WRITE(4,259)(((I-1 )*DX* 100),I=1,N,NSTEP) WRITE(5,259)(((I- 1 )*DX* 100),I=1,N,NSTEP) WRITE(4,259)(((I- 1 )*DX* 100),I= 1,M,NSTEP) WRITE(5,259)(((I-1 )*DX* 100),I=1,M,NSTEP) PAUSE DX=DX/L C C SET INITIAL CONDITIONS FOR ALL NODES C ASSUME POISEUILLE OR SLUG FLOW, AND UNIFORM INITIAL C TEMPERATURE C DO 10,I=1,N DO 20,J=1,M Y=(J-1)*DX IF (NTYPE.EQ.1) THEN PSI(I,J)=Y W(I,J)=0.0 ELSE S=L/B PSI(I,J)=6.0*((S*(Y**)2)/20)-((S**2)*(Y**3)/3.0)) W(I,J)=-6.0*(S-(2.0*Y*S**2)) ENDIF TH(I,J)=0.0 20 CONTINUE 10 CONTINUE C C INCREMENT TIME AND COUNTER FOR PRINTING OUTPUT C TO SCREEN C 25 TIME=TIME+DT NCOUNT=NCOUNT+1 C C UPDATE BOUNDARY CONDITIONS C C AT INLET ASSUME UNIFORM OR POISEUILLE FLOW, AND UNIFORM TEMPERATURE C 1=1 IF (NTYPE.EQ.1) THEN DO 30J=1,M 27

TEMP=2*PSI(I,J)-5 *PSI(I+ 1,J)+4*PSI(I+2,J)-PSI(I+3,J) W(I,J)=(- 1.0/(DX**2))*TEMP 30 CONTINUE ELSE DO 31,J=2,M-1 TEMP 1=2*PSI(I,J)-5 *PSI(I+1,J)+4*PSI(I+2,J)-PSI(I+3,J) TEMP2=PSI(I,J+1 )+PSI(I,J-1 )-2*PSI(I,J) W(IJ)=(-1.0/(DX**2))*(TEMP1+TEMP2) 31 CONTINUE TEMP=2*PSI(I, 1)-5*PSI(I,2)+4*PSI(I,3)-PSI(I,4) W(I,1)=(-1.0/(DX**2))*TEMP TEMP=2*PSI(I,M)-5*PSI(I,M-1)+4*PSI(I,M-2)-PSI(I,M-3) W(I,M)=(-1.0/(DX** 2))*TEMP ENDIF C C AT Y=0, assume adiabatic surfaces everywhere C except on heater surface. Use Woods approximation for C vorticity on the surface. C J=1 DO 40,I=2,N1-1 TH(I,J)=( 1.0/3.0)*(4*TH(I,J+1 )-TH(I,J+2)) W(I,J)=-0.5*W(I,J+ 1 )-(3*(PSI(I,J+1 )-PSI(I,J))/(DX**2)) 40 CONTINUE DO 50,I=N1,N2 TH(I,J)=(1.0/3.0)*(4*TH(I,J+1 )-TH(I,J+2)+2.0*DX) W(I,J)=-0.5 *W(I,J+ 1 )-(3 *(PSI(I,J+ 1 )-PSI(I,J))/(DX * *2)) 50 CONTINUE DO 60,I=N2+1,N TH(I,J)=( 1.0/3.0) *(4*TH(I,J+1 )-TH(I,J+2)) W(I,J)=-0.5*W(I,J+1 )-(3 *(PSI(I,J+1 )-PSI(I,J))/(DX**2)) 60 CONTINUE C C At exit, assume constant temperature gradient, constant C curvature of streamlines C I=N DO 69,J=2,M-1 TH(I,J)=0.5 *(5 *TH(I-1,J)-4*TH(I-2,J)+TH(I-3,J)) W(I,J)=(-1.0/(DX**2))*(PSI(I,J-1 )-2*PSI(I,J)+PSI(I,J+ 1)) 69 CONTINUE C C AT Y=B 28

C J=M DO 70,I=2,N TH(I,J)=(1.0/3.0)*(4*TH(I,J-1)-TH(I,J-2)) W(IJ)=-0.5 *W(I,J-1 )-(3*(PSI(I,J-1 )-PSI(I,J))/(DX**2)) 70 CONTINUE C C PRINT OUT TEMPERATURE AND STREAM FUNCTION TO THE SCREEN C AT PRESELECTED TIME INTERVALS C IF ((NCOUNT.EQ.NPRINT).OR.(TIME.LE.DT)) THEN WRITE(*,261)TIME DO 32,J=M,1,-NSTEP WRITE(*,270)(PSI(I,J),I=1,N,NSTEP) 32 CONTINUE WRITE(*,201) DO 34,J=M,1,-NSTEP WRITE(*,280)(TH(I,J),I=1,N,NSTEP) 34 CONTINUE WRITE(*,290) READ(*,*)KEV IF (KEV.EQ.1) THEN GOTO 125 ENDIF NCOUNT=0 ENDIF C C C USE ALTERNATING IMPLICIT DIRECTION SCHEME TO SOLVE FOR TEMPERATURES C AT INTERIOR NODES. START WITH THE X-DIRECTION BEING IMPLICIT C AND USE UPWINDING SCHEME BASED ON LOCAL VELOCITY DIRECTION TO C AID CONVERGENCE C DO 80,J=2,M-1 DO 85,I=2,N-1 K=(J-2)*(N-2)+I-1 UP=0.5*(PSI(I,J+1 )-PSI(J-1 )+PSI(I+1,J+1 )-PSI(I+1,J-1)) UM=0.5 *(PSI(I,J+1 )-PSI(I,J-1 )+PSI(I-1,J+1 )-PSI(I-1,J- 1)) VP=-0.5*(PSI(I+1,J)-PSI(I-1,J)+PSI(I+1,J+ 1 )-PSI(I- 1,J+1 )) VM=-0.5*(PSI(I+1,J)-PSI(I-1,J)+PSI(I+1,J-1 )-PSI(I- 1,J- 1)) Al(K)=-PE*((UM+ABS(UM))/(2.0*DX**2)+(1.0/(PE*DX**2))) TEMP=(UP+ABS(UP)-UM+ABS(UM))/(2.0*DX**2) 29

B 1 (K)=PE*((2.0/(PE*DT))+TEMP+(2.0/(PE*DX**2))) Cl(K)=PE*((UP-ABS(UP))/(2.0*DX**2)-( 1.0/(PE*DX**2))) TEMP1 =(VP-ABS(VP))*TH(IJ+1 )+(VP+ABS(VP)-VM+AB S(VM))*TH(I,J) TEMP1 =TEMP1 -(VM+ABS(VM))*TH(I,J-1) TEMP =(TEMP1)/(2.0*DX**2) TEMP3=TH(I,J+ 1 )+TH(I,- 1 )-2.0*TH(I,J) TEMP3=( 1.0/(PE*DX**2))*TEMP3 D1(K)=PE*((2.0*TH(I,J)/(PE*DT))-TEMP 1 +TEMP3) IF (I.EQ.2) THEN A1(K)=0.0 ENDIF IF (I.EQ.(N-1)) THEN D (K)=D 1(K)-C1(K)*TH(I+1,J) C1(K)=0.0 ENDIF 85 CONTINUE 80 CONTINUE C C SOLVE FOR TEMPERATURES AT HALF TIME STEP C CALL TRID(NUM,A 1,B 1,C1,D 1,E1) DO 81,J=2,M-1 DO 86,I=2,N-1 K=(J-2)*(N-2)+I-1 TH(I,J)=E1(K) 86 CONTINUE 81 CONTINUE C C NOW MAKE THE Y-DIRECTION IMPLICIT AND SOLVE FOR TEMPERATURE C DO 82,I=2,N-1 DO 87,J=2,M-1 K=(I-2)*(M-2)+J-1 UP=0.5*(PSI(I,J+1 )-PSI(I,J-1 )+PSI(I+1,J+1)-PSI(I+1,J-1)) UM=0.5*(PSI(I,J+ 1 )-PSI(I,J-1)+PSI(I-1,J+ 1 )-PSI(I- 1,J-1)) VP=-0.5*(PSI(I+1,J)-PSI(I-1,J)+PSI(I+1,J+1)-PSI(I-1,J+1)) VM=-0.5*(PSI(I+1,J)-PSI(I-1,J)+PSI(I+ 1,J-1)-PSI(I- 1,J-1)) A (K)=-PE*((VM+ABS(VM))/(2.0*DX**2)+(1.0/(PE*DX**2))) TEMP=(VP+ABS(VP)-VM+ABS(VM))/(2.0*DX**2) B 1(K)=PE*((2.0/(PE*DT))+TEMP+(2.0/(PE*DX**2))) Cl(K)=PE*((VP-ABS(VP))/(2.0*DX**2)-(1.0/(PE*DX**2))) TEMP 1=(UP-ABS(UP))*TH(I+ 1,1 )+(UP+ABS(UP)-UM+ABS(UM))*TH(I,J) TEMP =TEMP 1 -(UM+ABS(UM))*TH(I-1,J) 30

TEMPI =(TEMP1 )/(2.0*DX**2) TEMP3=TH(I+ 1,J)+TH(I- 1,J)-2.0*TH(I,J) TEMP3=( 1.0/(PE*DX**2))*TEMP3 D1(K)=PE*((2.0*TH(I,J)/(PE*DT))-TEMP 1 +TEMP3) IF (J.EQ.2) THEN IF ((I.GE.N1).AND.(I.LE.N2)) THEN B1(K)=B 1(K)+(4.0/3.0)*A1(K) C1(K)=C1(K)-(1.0/3.0)*A1(K) D1(K)=D1(K)-A1(K)*(2.0*DX/3.0) A1(K)=0.0 ELSE B 1 (K)=B 1(K)+(4.0/3.0)*A 1(K) C1(K)=C1(K)-(1.0/3.0)*A1(K) A1(K)=0.0 ENDIF ENDIF IF (J.EQ.(M-1)) THEN B 1 (K)=B 1(K)+(4.0/3.0)*C1 (K) A1(K)=A1(K)-(1.0/3.0)*C1(K) C1(K)=0.0 ENDIF 87 CONTINUE 82 CONTINUE CALL TRID(NUM,A1,B1,C1,D1,E1) DO 88,I=2,N-1 DO 89,J=2,M-1 K=(I-2)*(M-2)+J- 1 TH(I,J)=E1(K) 89 CONTINUE 88 CONTINUE C C USE ALTERNATING IMPLICIT DIRECTION SCHEME TO SOLVE FOR VORTICITY C AND START WITH THE X-DIRECTION BEING IMPLICIT C DO 90J=2,M-1 DO 100,I=2,N-1 K=(J-2)*(N-2)+I- 1 UP=0.5 *(PSI(I,J+1 )-PSI(I,J- 1 )+PSI(I+1,J+1 )-PSI(I+ 1,J- 1)) UM=0.5*(PSI(I,J+1)-PSI(I,J-)+PSI(I-1,J+1)-PSI(I- 1,J- 1)) VP=-0.5*.(PSI(I+ 1,J)-PSI(I-1,J)+PSI(I+1,J+1 )-PSI(I-1,J+ 1)) VM=-0.5*(PSI(I+1,J)-PSI(I-1,J)+PSI(I+1,J-1 )-PSI(I-1,J- 1)) Al(K)=-PE*((UM+ABS(UM))/(2.0*DX**2)+(1.0/(RE*DX**2))) 31

TEMP=(UP+ABS(UP)-UM+ABS(UM))/(2.0*DX* *2) B 1(K)=PE*((2.0/(PE*DT))+TEMP+(2.0/(RE*DX**2))) C (K)=PE*((UP-ABS(UP))/(2.0*DX**2)-(1.0/(RE*DX**2))) TEMP1 =(VP-ABS(VP))*W(I,J+1 )+(VP+ABS(VP)-VM+ABS(VM))*W(I,J) TEMP1=TEMP1 -(VM+ABS(VM))*W(IJ- 1) TEMPI =(TEMP1 )/(2.0*DX**2) TEMP2=(TH(I+1,J)-TH(I-1,J))*SIN(ALPHA*PI/180.0) TEMP2=TEMP2-(TH(I,J+1 )-TH(I,J-1 ))*COS(ALPHA*PI/180.0) TEMP2=-(RI/(2*DX))*TEMP2 TEMP3=W(I,J+ 1 )+W(I,J- 1 )-2.0*W(I,J) TEMP3=(1.0/(RE*DX**2))*TEMP3 D1 (K)=PE*((2.0*W(I,J)/(PE*DT))-TEMP1 +TEMP2+TEMP3) IF (I.EQ.2) THEN D1(K)=D 1 (K)-A 1 (K)*W(I- 1,J) A1(K)=0.0 ENDIF IF (I.EQ.(N-1)) THEN D1 (K)=D 1(K)-C 1(K)*W(I+ 1,J) C1(K)=0.0 ENDIF 100 CONTINUE 90 CONTINUE CALL TRID(NUM,A 1,B 1,C 1,D 1,E1) DO 72,J=2,M-1 DO 73,I=2,N-1 K=(J-2)*(N-2)+I-1 W(I,J)=E1(K) 73 CONTINUE 72 CONTINUE C C MAKE THE Y-DIRECTION IMPLICIT NOW AND SOLVE FOR VORTICITY C DO 91,I=2,N-1 DO 101,J=2,M-1 K=(I-2)*(M-2)+J-1 UP=0.5*(PSI(I,J+1 )-PSI(I,J-1 )+PSI(I+1,J+1 )-PSI(I+1,J-1)) UM=0.5* (PSI(I,J+1 )-PSI(I,J-1 )+PSI(I- 1,J+ 1 )-PSI(I-1,J-1 )) VP=-0.5 *(PSI(I+1,J)-PSI(I- 1,J)+PSI(I+1,J+ 1 )-PSI(I- 1,J+1)) VM=-0.5*(PSI(I+1,J)-PSI(I-1,J)+PSI(I+1,J-1 )-PSI(I-1,J- 1)) Al(K)=-PE*((VM+ABS(VM))/(2.0*DX**2)+(1.0/(RE*DX**2))) TEMP=(VP+ABS(VP)-VM+ABS(VM))/(2.0*DX**2) B 1(K)=PE*((2.0/(PE*DT))+TEMP+(2.0/(RE*DX**2))) Cl(K)=PE*((VP-ABS(VP))/(2.0*DX**2)-(1.0/(RE*DX**2))) 32

TEMP1=(UP-ABS(UP))*W(I+1,1)+(UP+ABS(UP)-UM+ABS(UM))*W(I,J) TEMP =TEMP1 -(UM+ABS(UM))*W(I-1,J) TEMP1=(TEMP1 )/(2.0*DX**2) TEMP2=(TH(I+1,J)-TH(I-1,J))*SIN(ALPHA) TEMP2=TEMP2-(TH(I,J+1)-TH(I,J-1))*COS(ALPHA) TEMP2=-(RI/(2*DX))*TEMP2 TEMP3=W(I+1,J)+W(I- 1 J)-2.0*W(I,J) TEMP3=(1.0/(RE*DX**2))*TEMP3 Dl(K)=PE*((2.0*W(I,J)/(PE*DT))-TEMP1 +TEMP2+TEMP3) IF (J.EQ.2) THEN D1(K)=D 1(K)-A1(K)*W(I,J-1) A1(K)=0.0 ENDIF IF (J.EQ.(M-1)) THEN D1(K)=D 1(K)-C 1(K)*W(I,J+1 ) C1(K)=0.0 ENDIF 101 CONTINUE 91 CONTINUE CALL TRID(NUM,A 1,B1,C1,D 1,E1) DO 92,I=2,N-1 DO 93,J=2,M-1 K=(I-2)*(M-2)+J-1 W(I,J)=E1(K) 93 CONTINUE 92 CONTINUE C C SOLVE FOR STREAM FUNCTION THROUGH ITERATION C 95 DIFF2=0.0001 DIF=0.001 DO 102,I=2,N-1 DO 103,J=2,M-1 PSIOLD(I,J)=PSI(I,J) TEMP=PSI(I+ 1,J)+PSI(I-1,J)+PSI(I,J+ 1 )+PSI(I,J- 1) PSI(I,J)=0.25*((W(I,J)*DX**2)+TEMP) DIFF1=ABS(PSI(I,J)-PSIOLD(I,J)) IF(DIFF1.GT.DIFF2) THEN DIFF2=DIFF1 DIF=DIFF2/PSI(I,J) ENDIF 103 CONTINUE 102 CONTINUE 33

I=N DO 104,J=2,M-1 PSIOLD(I,J)=PSI(I,J) PSI(I,J)=0.5 *(5 *PSI(I- 1,J)-4*PSI(I-2,J)+PSI(I-3,J)) DIFF1=ABS(PSI(I,J)-PSIOLD(IJ)) IF (DIFF1.GT.DIFF2) THEN DIFF2=DIFF1 DIF=DIFF2/PSI(I,J) ENDIF 104 CONTINUE C C ITERATE UNTIL STREAM FUNCTION CHANGE IS LESS THAN C HALF A PERCENT C IF (DIF.GT.0.005) THEN GOTO 95 ENDIF C C CHECK IF SYSTEM HAS REACHED STEADY STATE C IF (MODE.EQ.2) THEN DO 304,J=1,M DO 301,I=2,N IF (TH(I,J).GT.0.000001) THEN DIFF1 =AB S(TH(I,J)-TOLD(I,J)) DIFF1 =DIFF1/TH(I,J) IF (DIFF1.GT.(10.0*DT)) THEN DO 302,J1=1,M DO 303,I1=2,N TOLD(I 1,J1)=TH(I1,J1) 303 CONTINUE 302 CONTINUE GOTO 25 ENDIF ENDIF 301 CONTINUE 304 CONTINUE GOTO 125 ENDIF C GOTO 25 C C OUTPUT RESULTS 34

C 125 WRITE(*,260)RI,RE,PE,TIME DO 110,J=M,1,-NSTEP WRITE(*,270)(PSI(I,J),I=1,N,NSTEP) 110 CONTINUE DO 111,J=1,M,NSTEP WRITE(5,270)(PSI(I,J),I=1,N,NSTEP) 111 CONTINUE WRITE(*,201) DO 120J=M,1,-NSTEP WRITE(*,280)(TH(I,J),I= 1,N,NSTEP) 120 CONTINUE J=1 NUT=0.0 DO 113,K=N1+1,N2 NUX=((K-N 1)*DX)/TH(K,J) NUT=NUT+NUX 113 CONTINUE NU=NUT/(N2-N1) PRINT*,'AVERAGE NU = ',NU DO 121,J=1,M,NSTEP WRITE(4,280)(TH(I,J),I= 1,N,NSTEP) 121 CONTINUE PAUSE close(4) close(5) C C 200 FORMAT(' ENTER CHANNEL LENGTH, HEIGTH, ENTRY LENGTH IN CM') 201 FORMAT('') 202 FORMAT(' CHOOSE EITHER 1 OR 2:') 203 FORMAT(' 1. ENTER DIMENSIONLESS PARAMETERS (RI,RE,PE)') 204 FORMAT(' 2. ENTER DIMENSIONAL QUANTITIES') 205 FORMAT(' ENTER N,M,N1,N2,NSTEP') 208 FORMAT(' ENTER ANGLE OF INCLINATION WITH VERTICAL') 209 FORMAT(' ENTER FLOW TYPE: 1) SLUG FLOW 2) POUISELLE FLOW') 210 FORMAT(' ENTER INLET TEMP(K),BETA (1/K), HEAT FLUX (W/cm2)') 211 FORMAT(' ENTER DT,FINAL TIME,PRINT TIME (1 S =.0001)') 212 FORMAT(' ENTER MODE: 1) TRANSIENT) 213 FORMAT(' 2) STEADY-STATE') 214 FORMAT(' ENTER HEATER LENGTH, DX IN CM') 215 FORMAT(' RI= ',F9.2,' RE= ',F9.1,' PE= ',F9.1) 220 FORMAT(' ENTER INLET VELOCITY (cm/s), GRAVITY (g)') 35

230 FORMAT(' ENTER THERMAL CONDUCT (W/(m K)),DIFFUSITY (mA2/s)') 240 FORMAT(' ENTER KINEMATIC VISCOSITY (m^2/s)') 250 FORMAT(' ENTER RI,RE,PE') 254 FORMAT(' SOLUTION DOES NOT CONVERG!!!') 256 FORMAT(' N= ',13,' N1= ',13,' N2= ',13,' M= ',I3,'STEP= ',13) 257 FORMAT(I3,I3) 259 FORMAT(100F6.2) 260 FORMAT(' RI=',F9.4,' RE=',F7.1,' PE=',F7.3,' TIME=',F7.3) 261 FORMAT(' TIME =',F9.4) 270 FORMAT(100F7,3) 280 FORMAT(100F6.3) 290 FORMAT(' CONTINUE? (1=NO)') 295 FORMAT(' ISOX DATA FILE:STREAM LINES') 300 FORMAT(' ISOX DATA FILE:TEMP LINES') END C C TRIDIAGONAL MATRIX SOLVER C SUBROUTINE TRID(N,D 1,D2,D3,R,F) REAL*8 D1(10000),D2(10000),D3(10000),R(10000) REAL*8 F(10000) INTEGER N,I DO 15,I=2,N D2(I)=D2(I)-(D3(I-1 )*D 1 (I)/D2(I-1)) R(I)=R(I)-(D 1 (I)*R(I- 1 )/D2(I- 1)) 15 CONTINUE C C USE BACKWARD SUBSTITUTION TO SOLVE FOR F(X) C F(N)=R(N)/D2(N) DO 25,I=N-1,1,-1 F(I)=(R(I)-D3(I)*F(I+ 1 ))/D2(I) 25 CONTINUE RETURN 36

VI. References Arpaci, Vedat S., and Larsen, Poul S., Convection Heat Transfer, PrenticeHall, Inc., 1984. Carnahan, Brice, Luther, H.A., and Wilkes, James 0., 1969, Applied Numerical Methods, John Wiley & Sons, New York. Davis, Graham de Vahl, 1986, Finite Difference Methods for Natural and Mixed Convection in Enclosures, Proceeding of the Eighth International Heat Transfer Conference, Vol. I, pp. 101-109. Davis, Graham de Vahl, and Mallinson, G.D., 1976, An Evaluation of Upwind and Central Difference Approximations by a Study of Recirculating Flow, Computers and Fluids, Vol. 4, pp. 29-43. Gill, W.N., and Del Casal, E., 1962, A Theoretical Investigation of Natural Convection Effects in Forced Horizontal Flows, AIChe Journal, 8, pp. 513-518. Hatton, A.P, and Woolley, N.H., 1971, Laminar Combined Natural and Forced Convection in a Rectangular Field, Cl19/71, in Heat and Mass Transfer by Combined Forced and Natural Convection, pp. 47-53. Kennedy, K.J., and Zebib, A., 1982, Combined Forced and Fee Convection between Parallel Plates, Paper 82-IHTC-152, Proc. 7th Int. Heat Transfer Conf. Heat Transfer Conf., Hemisphere, New York. Nyugen, T.V., Maclaine-cross, I.L., and Davis, Graham de Vahl, 1979, Combined Forced and Free Convection Between Parallel Plates, Numerical Methods in Thermal Problems, pp. 269-278. Pine Ridge Press, U.K. 37

Appendix B To National Aeronautics and Space Administration Lewis Research Center Study for Forced Convection Boiling Under Reduced Gravity Grant NAG 3-589 THE VELOCITY POTENTIAL OF TWO SPHERES MOVING PERPENDICULARLY TO THE LINE JOINING THEIR CENTERS AND THE ASSOCIATED LIFT FORCE Longhu Li Herman Merte University of Michigan Ann Arbor, Michigan August, 1991

1 THE VELOCITY POTENTIAL OF TWO SPHERES MOVING PERPENDICULARLY TO THE LINE JOINING THEIR CENTERS AND THE ASSOCIATED LIFT FORCE I. Introduction Vapor bubbles generated on a flat heating surface in forced convection boiling, aside from being subjected to the influences of forces involved in pool boiling conditions, are also subjected to the pressure gradient in the bulk flow field and to the lift produced due to the pressure distribution caused by the velocity redistribution. Drew (1987) used a numerical method to compute the lift force for small spheres moving in a shearing fluid near a wall, but the result is not convenient to use analytically and does not include the condition when the sphere touches the wall. An alternative to the numerical solution is a potential flow solution involving the use of a potential flow model to approximate the real situation when a bubble is on or near a wall with bulk liquid flow parallel to the wall in forced convection boiling. The simplest potential flow model which best approximates the real situation is that of a spherical bubble attaches to or near a wall, with an ideal fluid sweeping over this infinite plane. Since this potential flow model is equivalent to two spheres in contact or separated, and moving perpendicularly to the line joining their centers in an infinite ideal fluid, it is natural that the solution should be pursued along this line. A number of classical solutions to this problem exist. The solution of two spheres moving in a perfect fluid was first attempted by Stokes (1843), then by Herr Bjerknes (1863), and Hicks (1879). However, it was not until 1887 1 that Basset(1887) and Herman(1887) solved the problem for two spheres moving perpendicularly to the line joining their centers. Even though they were the first ones to use the first order associated Legendre polynomials to solve the problem(Basset used the

2 integral expression and Herman used the derivative expression of the polynomial), they did not solve the problem completely; instead, they used only the first degree of the polynomial in the final expression for the velocity potential, neglecting all higher degree terms, most likely due to the lack of computational means. This gave the impression that the method was good only for approximate solutions, and much effort was subsequently devoted to seeking more accurate solutions. The problem was essentially solved by Basset (1887). Miloh (1977) used the proposition proposed and proved by Basset(1887), which later appeared in Hobson(1931). The solution of Miloh for the subject problem is identical to that obtained by using Basset's method directly. The only difference is in the solution technique: in the formulation of Miloh the potential coefficients are solved from an infinite set of equations, while in the formulation of Basset the potential coefficients are solved by adding a set of infinite summation series. It will be demonstrated in the following sections that the two formulations are equivalent. Basset was interested in the kinetic energy associated with the motion of the spheres, so no results were given relative to the motion interaction forces between the two spheres. Miloh did give this force, but with insufficient accuracy. The solution was based on solving an infinite set of equations, and the residual normal velocity produced by the velocity potential on the spherical surface can be a check on the accuracy of the solution. This check was not conducted by Miloh; it was simply assumed that the solution satisfied the prespecified boundary conditions, which was not the case. Miloh had the advantage of knowing the form of the potential function, and therefore was able to set out to solve the coefficients in that function directly; in Basset's time the potential function was completely unknown, and instead of using the geometric imaging method of single dipoles, as was done by prior workers, he was the first person to successfully use the function imaging method to lead to the form of the potential function in the spherical coordinate system.

3 II. The construction of the velocity potential The following is an extension of the work of Basset. The configuration and coordinate systems are shown in Fig. 1. Two spheres of radii rbl and rb2 are moving in an ideal stationary liquid with respective uniform velocities U1 and U2 in the negative X direction. Although the generality of U1 * U2 will be retained for some time, the particular case of interest is where U1 = U2, rbl = rb2, and the X direction is perpendicular to the line joining the centers of the spheres as shown. For clarity, the notation rbl and rb2 will be retained, although it is understood that rbl = rb2. _ X 4A - Ic ry^ x Fig. 1. Two spheres moving perpendicularly to the line joining their centers Spherical coordinates r, 0, ), and r',0', * are fixed at the centers of the spheres rbl and rb2 at 01 and 02, respectively. The Cartesian coordinate system X,Y, Z has the origin at the mid-point between the two spheres, with the X axis perpendicular to this line and the angle 0 measured in the X-Y plane. Lengendre polynomials with degree n and order m, whose origins are at 01 and 02, respectively, are defined by the following relations:

4 or Pn =M(sinO')m f ( + 2 -1 cosa)nm(sina)2m da pm = M(sin0')m (c + t'2 -1 cosa)nm(sina)2m da P = M(sin)m J (sina)2mda ( + &-1 Tcos(a)n+m+l In P'm = M(sinO')m (sina)2mda ( + ' 2-1 cosa)n+m+.................... (1)....................(1') where M(n+m)! (n-m)!1 3*...(2m-1)7 ' = cos O' Basset was the first one who proposed and proved the following propositions, which relate the associated Legendre polynomials PMn and P'm:,m n rm" (n+m)! m (n+m+l)! r m rn+l (n-m)!cn+ml 2m! (2m+1)! c Pm + rtn+i(n-m)!cn+m+l 2m? Am (2m+1)! c Pm+ + (-1)k (n+m+k)! r P^ +)....................(2) (2m+k)! c+ m+k when r<c, and (.1)n-mp n r n rm (nc+m)!,m (n+m+l)! r',m rn+l (n-m)! cn+m+l 2m! Pm+ (2m+l)! c pm+l + (n+m+k)! r'k,m +(2m+k)! (-) Pm+k +..*)....................(3) when r'<c Since velocity potentials are linear, the total velocity potential 0 at any points outside the spheres is the sum of the potentials Olcosqp produced by the sphere of radius rb centered at 01 and 02cosq produced by the sphere of radius rb2 centered at 02. For the sake

5 of simplicity, here 01 and 02 are used to denote the amplitude of the two velocity potentials, i.e. b= (D1 + 42) cos(................. (4) Once 01 is known, <2 can be determined directly from geometric symmetry. (1 must be constructed in such a manner that it satisfies the boundary conditions on the surfaces of the two spheres: ~ b1 1 ar I bl U l.................... (5) a<i I 0....................(6) ar r —rb2 Beginning with the well-known potential for a single sphere moving in an infinite quiescent pool of incompressible ideal liquid: 3 rbl 1 1 1 2r P....................(7) which satisfies the boundary condition given by Eq. (5). Equation (3) is used to transform Eq. (7) into the velocity potential in the vicinity of the sphere centered at 02: 3 0o rbl r l ri =c u X i(8) 3 i = ci' P i'....................(8) 2c3 i=l c1 Eq. (8), however, does not satisfy the boundary condition Eq. (6). To satisfy Eq. (6), Eq. (9) is added to Eq. (8): 3 2i+l, rb. i rb2 Pi U l1 i +i i-....................(9) 1 2c3 i 1 i+ 1,i 1 01 in the vicinity of the sphere 02 is now given by: 3 0 2i+1,i rbl i i rb2. (ri = U1 r,'-+ i-...................(10) 2,3 + c

6 Although Eq. (10) now will satisfy the boundary condition Eq. (6), it will no longer satisfy boundary condition Eq. (5). To demonstrate this, using Eq. (2), Eq. (9) may be transformed into: (9)= -— rb 23 i 2c im 2i+l i rb2 1 (-'l(i+j)! rJ i i+l ci (i-1)!ci+2 j= (j+1)! cJ-1 P '2 2i+ )1 rlo- rb3lU 2c3 i=l.2 2i+1 i rb2 (i+l)!c2i+l (-1- (i+j)! rJ (j+1)! cJ- PJ 3 - j 00 00 rb Ul Z 2c3 j=l i=l 2i+1 (-1)J' rJ l i2(i+j)! rb2 (j+l)! c- Pi (i+l)! c2i+l....................(11) upon adding Eq. (11) to Eq. (7), the velocity potential t1 in the vicinity of sphere 01 is now 3 rbl O1 = UI2r2 3, o P1 rbl U1 i v3 2i+1 (-1)J1 ri 1 i2(i+j)! rb2 (j+l)! cj-l J (i+l)! c2i+l....................(12) Obviously, the sum of the two terms on the right hand side of Eq. (12) cannot satisfy Eq. (5). To satisfy the boundary condition given by Eq. (5), Eq. 3 o o 2j+l 1 2i+1 rblU (-lJ-lj rb Pi i2(i+j)! b2 2c3 j i (j+l)(+)! - r+l (i+)! c2i+l (13) is added to Eq. (12),....................(13) The velocity potential in the vicinity of the sphere at 01 is now 3 3 o* rbl 1 rblU1 ~1 =-r PI + 2~3 j=l 0i i=l 2j+l 1 2i+l (-li'1 (j rbl P i2(i+j)! rb2 (j+1)! j+i rj+l Cj-1 (i+l)! c2i+l...................(14) By summing Eqs. (7), (9) and (13), the potential function produced by the motion of sphere 01, valid at all locations and still satisfying Eq. (5) but not Eq. (6), is:

7 3 rbl 1 U1 Pl + ^3^ 3 00 rbl; 13 2c3i~ u ^y 2i+l,1 i rb2 p i+1 i-i ri+l C r 3, 00 00 rblU1 v 2c j=l i=l 2j+l 1 2i+1 (.)-lj rbl Pj i2(i+j) rb2 (j+l)(j+)! c-1 r1 (i+1)! c2...1.. (15) (i+1)(j+1)! eJ-1 rJ+1 (i+l)! c2i+f 'c-......................... The above procedure can now be expanded, alternately satisfying Eqs. (5) and (6), resulting in Eq. (16), the general potential function for sphere 01 valid at all locations. The convergence of this function is shown by the convergence of the normal velocity, i.e., its derivative, at the sphere surface in section IV below. 3 bl 1 '1 = U1 Pl + 3 rbl 12c3 00 2i+1 P1 { 2 i r b2 Pi + i+ ci- ri+l i=lI j=0l j=l i=l 2j+l 1 2i+l (,bl Pj l 2(i+j)! r b2 (j+1)(j+1)! cJ-1 r+l1 (i+l)! c2i+1 21+l,1 2j+l 2i+1 k rb2 P j2(j+k)! rbl i2(i+j)! rb2 il (k+1)(k+1)! ck-l rk+l ((j+1)!)2 c2j+ (i+l)! c2i+1 + 00 k=l 00 j=l - 00 00 21+1 1 2k+l 2j+l +V. 1 1- 1 r bl Pl k2(k+l)!rb2 j2(j+k)! rbl + (- 1)1- I(+1)(1+1)! 1 I+ ((k+.1).1 1=1 k=l j (1+1)(1+1)! cl- r ((k+)!)c2k+l (j+1)!)2c2+ 2i+1 i2(i+j)! b2 +...(16) (i+ 1)! c2 +1....................(16) (i+l)! c2i+l By symmetry 02 can be written directly: 3 rb2 n,2 = v2 pI 2r'2 00 00 j=l i=l 00 00 4 -+ i:i k=l jl i 00 00 1=1 k=l 1 I - 3 rb2 2 2c3 ~00 ~ 2i+l 1 (I ( 1)i-. Ii rbp + i+1 i-i i+l i=l c r 2j+l,1 2i+l j b2 p j i2(i+j)! rbl (j+l)(j+l)! cJ-1 rJ+l (i+l)! c2i+l + 2k+l 1 2j+l 2i+l (-lk-i k rbl Pk j2(j+k)! b2 i2(i+j)! bl =1 (k+l)(k+l)! ck-l r+l ((j+1)!)2 c2j+l (i+l)! c2i+l 00 i=l 0 21+1,1 2k+l 2j+l 1 rb2 P k2(k+ l)! rbl j2(j+k)! b2 i=l (1+1)(1+1)! Cl-l r'l+l ((k+1)!)2c2k+l ((j+1)!)2 c2j+

8 2i+l i2(i+j)! rbl (i+ 1)! c2i+l....................(17) 1 Substituting Eqs. (16) and (17) into (4), and collecting terms that have the same pn results in: 4D= U1 'cos(P....................(18) where ' is the amplitude for the potential of unit velocity in the direction of U1, given by 3 3 rbl 1 U2 rb2,1 = 2r2 2r 2r2 Pi +UI 2r,2 1 ~~ 1,1 Yf(n 1 A +2 Pn B,n+2 n + ( n-1 AnIrbl n1 + Bnrb2 rtn+l n=l....................(19) An's and Bn's in (19) are given by the following infinite summation series 3 n-I 3 oo n-l 2i+1 U2 rb2 n rb l rb n rbl i2(i+n)! rb2 An= U1 23c n+1 n -' ' U1 2c3 n+1 cn-T cn-l(n+l)!(i+1)!c2i+ 3 00 U2 rb2E U1 2c3 j= 3 0o 2c3 k=l 00 j=l n-l 2j+l 2i+1 n rbl j2(j+n)! rb2 i2(i+j)! rb n+1 cn-l(n+l)!(j+l)! c2j+l (j+)!(i+l)! c2i+l + n-l 2k+l 2j+l n bl k2(k+n)! rb2 j2(j+k)! rbl n+l cn-1 (n+1)!(k+1)! c2k+l (k+l)!(j+l)! c2j+l +.................... (20) 3 0~ n-I 2i+l rb2 n rb2 i2(i+n)! rbl 1 2c3 n+l cn-l-(n+l)!(i+1)! 2i+l 2i+1 i2(i+j)! rb2 (j+ 1)!(i+ 1)! c2i+l 3 n-l rbl n rb2 U. Bn 2c3 n+1 cn-l U 3 oo n-l 2j+l 2i+l '2c3 ji il n+l cn- (n+l)!j+l)!c2j+l (j+1)!(i+1)!c2i+l + 3 oo oo oo n-l 2k+l 2j+l U2 rb 2 2 n rb2 k2(k+n)! rbl j2(+k)! b2 U12c3 k= j= i= n+l cn-1 (n+l)!(k+l)! c2k+1 (k+l)!(j+1)! c2j+ 2i+1 i2(i+)! rb (21) (j+l)!(i+l)! c2i+1 +.................... or n-l n rbl n n+l n-l 3 U2 rb2 U1 2c3 3 00 rbl 2-c3 i=l 2i+1 i2(i+n)! rb2 (n+l)!(i+ 1)! c2i+

9 3 ooj r 2j+1 2i+l U2 r2brb 2 +U2 yB x j2(j+n)! rb2 i2(i+j)! rbl U1 2c3 j= i=l (n+l)!(j+l)! c2j+l (j+l)!(i+l)!e2i+l + 3 ~ o0 00 2k+l 2j+1 +r1, V k2(k+n)! rb2 jj+k)! rbl +2-c3^ = 1 (n+l)!(k+l)! c2k+1 (k+l)!(j+l)! C2j+ k=l j=l i=l 2i+l i2(i+j)! rb2 * (j+l)(i+ )!c2i+ +......................20') n-i 3 3 - 2i+1 rnb2 rb 2 r3 b i2(i+n)! bl rbl + U+ n n+l cn-l 2c3 U+ 2c3 (n+l)!(i+l)!c2i+l + r3 ~ o 2j+1 2i+1 rbl j2(+n)! rbl i2(i+j)! rb2 + ---+ +2c3 j= i= (n+l)!(j+l)!c2j+1 (j+l)!(i+l)!c2i+l + 3 oo 0 2k+l 2j+l U2 rb2 Z k2(k+n)! rbl j(j+k)! b2 Ui 2c3 k=l j=l (n+l)!(k+l)! c2k+l (k+l)!(j+l)! c2j+l 2i+1 i2(i+j)! r2i +....................(21') (j+ 1)!(i+ 1)! e2i+i In (20) and (21), U1 and U2 are the parallel velocity components of sphere 1 and sphere 2, respectively. They may be identically the same or different in either or both direction and magnitude, and the positive direction of U1 is in the negative X direction. In the case that either the direction or magnitude or both of one velocity is different from that of the other, the above result is correct only at the moment when the two spheres are moving in such a way that these two velocity components are perpendicular to the line joining their centers. A^'s and Bn's in (20) -and (21) can be evaluated one at a time. However from (20) and (21) it can be seen that the result of a lower summation order term in An could be used in the evaluation of a next higher summation order term in Bn and vice versa. It is obvious that it would be more advantageous to evaluate these simultaneously in groups (with the upper limit of N prespecified) in order to reduce computational effort. These have been computed below for two identical spheres, with U1= U2, and for the cases of various center distances.

10 For the case of two identical spheres touching and moving with the same velocity, U1=U2, rbl== rb and c=2r Eq. (19) reduces to: l b2 b b 1b I b,1 "n+2 P n' zt' =P p + Anrb {(-1)n'l rn l +I n l }.................... (22) =rb 2r'2 I 1rb (()l ).(22) where 0i 00 1 1 n 1 n 1 n i2(i+n)! 1 An 24 ~ ' {2n-'1 n+ 2n-l n+il (n+ )!(i+ 1)! 22i+1 + 1 n Y j2(j+n)! 1 i2(i+j)! 1 2`n-ln+j1 I2. (n+l)!(j+1)! 22j+ (j+ 1)!(i+1)! 22i+1 00 00 00 1 n X Y k2(k+n)! 1 j2(j+k)! 1 i2(i+)! 1 (23) 2+1 +1 (n+)!(k+)! 22k+ (k+l)!(j+)! 22j+22i (j+l)!(i+l)!22i+l + "....................(23) or 00 A n l E i2(i+n)! 1 + j2(j+n)! 1 i2(i+j)! 1 ~1 i 1 (n+l)!(j+l)! 22j+1 (j+l)!(i+l)! 22i+1 00 00 00 + Y X k2(k+n)! 1 j2(j+k)! 1 i2(i+j)! 1 k^1 1 j-Ii I (n+l)!(k+l)! 22k+1 (k+l)!(j+1)! 22j+l (j+l)!(i+l)! 22i+1 +................... (23') The An's in Eq. (22) were evaluated by using Eq. (23). Table 1 gives the first 16 terms for the velocity potential with different numbers of terms. The An's were computed for N up to 2000. The relative error for each term in Eq. (23) was set as 10-17, and the relative error for An's was set as 10'11. Eq. (23) was also used to evaluate the An's for N up to 30, with relative error set as 104. Since each term in Eq. (23) uses the results of the previous term no relative error was set for each term in Eq. (23) in this case; instead each term is computed for the index from 1 to 30. The results are identical to those given by

11 Miloh (1977). This reveals why the values given by Miloh were smaller than what is obtained as the total number of terms is increased. Even though N was chosen to yield a maximum relative error of ~10i4 between successive approximations in the computations by Miloh, the total relative error is larger tr. a ~10-4, as can be seen in the values of the 16th term for n=30 and n=1000 in Table 1. For the case of two identical spheres moving in the same direction with the same velocity and not touching, the coefficients are computed by using Eq. (20') for center distances CIRb = 2.000001,2.000002, 2.000003, 2.000004,2.000005,2.000007,2.00001,2.000015, 2.00002, 2.00003, 2.00005, 2.00007, 2.0001, 2.0002, 2.0003;2.0005,2.0007, 2.001, 2.0015, 2.002, 2.003, 2.005,2.007,2.01, 2.015, 2.02, 2.02244,2.03,2.05, 2.07, 2.10, 2.20, 2.3,2.4,2.5, 2.6, 2.683, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.415, 3.5, 3.683, 4.0, 4.45, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0,20.0 and 22.0. When the center distance between the two spheres is large(C/Rb 22.02), the computation of the coefficients using Eq. (20') converges very rapidly. The first 16 terms of the coefficient in the velocity potential for a few selected center distances are given in Table 2. The results for the remaining cases listed are available, but not given here. Table 1. The first 16 coefficients A vs N(c/r = 2.000) n b n\N 30 100 200 500 1000 2000 1.04032383.04034059.04034122.04034132.04034133.04034133 2.02998921.03001425.03001519.03001535.03001536.03001536 3.01952214.01955302.01955418.01955437.01955438.01955438 4.01258031.01261584.01261716.01261738.01261740.01261740 5.00830055.00834003.00834151.00834175.00834177.00834177 6.00567546.00571841.00572002.00572028.00572030.00572030 7.00403218.00407826.00407998.00408027.00408029.00408029 8.00296983.00301876.00302059.00302089.00302092.00302092 9.00225711.00230868.00231061.00231093.00231096.00231096 10.00176079.00181483.00181685.00181718.00181720.00181721 11.00140299.00145934.00146145.00146180.00146182.00146183

12 12.00113708 13.00093425 14.00077610 15.00065050 16.00054916.00119562.00099487.00083872.00071504.00061556.00119781.00099714.00084106.00071745.00061804.00119817.00099751.00084145.00071785.00061845.00119819.00099754.00084148.00071788.00061848.00119820.00099754.00084148.00071789.00061848 Table 2. n\ — =2.000000 2.000001 rb The first 16 coefficients A vs c/rb n b 1.04034133 2.03001536 3.01955438 4.01261740 5.00834177 6.00572030 7.00408029 8.00302092 9.00231096 10.00181721 11.00146183 12.00119820 13.00099754 14.00084148 15.00071789 16.00061848.04341190.03001519.01955422.01261725.00834162.00572016.00408016.00302080.00231084.00181709.00146171.00119809.00099744.00084138.00071779.00061838 2.000010.04339905.03001369.01955274.01261584.00834030.00571892.00407898.00301967.00230976.00181606.00146071.00119712.00099650.00084047.00071690.00061753 2.000100 2.001000 2.010000 2.100000 2.200000.04032728.02999903.01953837.01260231.00832765.00570706.00406780.00300908.00229968.00180643.00145150.00118827.00098798.00083225.00070896.00060984.04020703.02986126.01940547.01247927.00821460.00560291.00397134.00291927.00221566.00172752.00137714.00111801.00092144.00076911.00064893.00055267.03915365.03195988.02664009.02869903.02153822.01679388.01833209.01245318.00902192.01153133.00696264.00464970.00738535.00388762.00237215.00487594.00219639.00121337.00333046.00126340.00062634.00235087.00074188.00032747.00170883.00044490.00017376.00127349.00027221.00009364.00096886.00016959.00005125.00074966.00010734.00002846.00058816.00006886.00001602.00046681.00004468.00000913.00037413.00002926.00000525.00030236.00001931.00000305 III. The transformation of Equations (20) and (21) into an infinite set of equations As was pointed out in the introduction, the formulation of the potential coefficients by an infinite series for each coefficient as given by Eqs. (20) and (21) is the same as th( formulation of Miloh, in which the coefficients are given as an infinite set of equations This will be demonstrated now. First rearrange Eq. (19) as:

13 00 1,1 = I,(n_)n-_ I n+2 Pn n+2 n ' ( )n-' A'n rb 1 r +l- + B'n n+b2 r n+ n=l r....................(24) where A'n = 2 (n-1)+An U2 B'n = 2u 8(n-1) + Bn....................(25) 12 and 8(n-1) is the Kronecker delta function which is one(l) for n=1, and zero(O) otherwise. Next rearrange Eq. (25) and substitute XA = and E =B U2 - C C U-[ n-l 3 o 1 A,e3 n A A'n= 2 8(n-1) +(n+l)(n+T)! {A2X (n+l)! + 2 m=l m2(m+n)! 2m+l (m+)! + 00 +,2AB l m=l 3 r 1 cc m=l i= m2(m+n)! 2m+l i2(i+m)! 2i+l (m+ 1)! (m+ 1)!(i+ 1)! i=1 o n m2(m+n)! 2m+l j2(j+m)! 2j+l i2(i+)! i+.. 1 (im+1)! A (j+l)!(i+ )! 1 i=l'' 1 A+nl1 (n+1)! 2 )(n+ )(n+l )! (2 X(n+l)!+ m=l m-1 (m+n)! m+2 mAB (m-i)! B (m+l)(m+l)! 3 *['t (m+l)!+ 00 1=1 i2(i+m)! 2i+l (i+1)! A + 3 2 *Ji2(j+m)! j=l i=l 1.2j+l i2(i+j)! A (j+l)!(i+ 1)! i+1 +...]....................(26) Since 23 (n+1)! = 8(m-1) (m 1)! + substitute Eq. (27) into Eq. (26): n-l At 1 2 n-1 _lA ( em+2 (m+n)! 2 =(n+l)nn+l)! {2 ( m-)XB (-)!................ (27) X (m+n)! m+2 (m-1)! K m=l

14 rm-1 ( B (m+ l)(m+ l)! 3 00 A 2 3 IT(m+l)!+| XBL i=l i2(m+i)! 2i+l (i+)! AA + (i+1)! 3 00 +2 J=1 i=l j2(m+j)! (j+l)!.2j+1 i2j+i)! 2i+l i A (j+l)!(i+l)! ** n-l nx'n-1 1 flA 1 8(n- )+ )! -2 (n+l)(n'+l)! (m+n)! (m-l)! m-1 3 mAn A m+2 (E(m-1) + (m1)(m)! [ (m+1)!+ )i 2 (m+ )(M+M2m+l)! + 00 i=l i2(m+i)! (i+l)! 3 00 2i+l A) AA +T L J=l i0 i=l j2(m+j)! (j+1)! A2j+1 i2(+i)! 2i+1 A (j+l)!(i+l)! *....................(28) rm-l mk' g m2(l mB B'm 2 8(m-l)+ (m+l)(m+l)! 3 { (m+l)! + 2 00 e 2 X i 1 i2(m+i)! 2i+l (i+1)! X"A + 3 0o XAV +T 2 j=1 j2(m+j)! (j+l)!.2j+1 i2 +i)! 2i+l + A (j+l)!(i+l)! i....................(29) Substitute Eq. (29) into Eq. (28), yields: n-l 1 n' A A'n 8(n-1) + )! 2 (n+ 1)(n+1)! 00 B (m+n)! m+2 m=lm (m-)! B....................(30).(30) In a similar manner it could be demonstrated that n-l Bn=|8n-l)+ B'n = 2 (n-1) +(n+l)(n+)! 00 A' (m+n)! rm+2 __ (nm (m-l)! "A m=l....................(31).(31) These are equivalent to the results given by Miloh. IV. The evaluation of the residual normal velocity produced by Eq. (19) at the Lt the surface of the spheres The adequacy of the velocity potential constructed above, or others constructed byicted by previous researchers, can only be checked by the evaluation of the residual normal normal

15 velocity, which should be zero on the spherical surface. Since the procedure for evaluating the normal velocity on the two spheres are the same, it would be sufficient to evaluate the residual normal velocity on either one. The residual normal velocity will be evaluated here on the sphere with center at 01. For any given 0 the normal velocity varies with the angle qp, and is a maximum when p=0~ or p=180~, that is, in the plane parallel to the direction of movement of the spheres containing both centers. It is sufficient to compute the residual normal velocity for qp=0~, with 0 varying from 00 to 1800 since the ones for qp=1800 are the same in magnitude as those for qp=0~ and opposite in sign. vI I I'....................(32) nr r=r bl ar Ir=rbl or 3,1 vI 1 U2 a rb2p I 1 v r r1 =r b '+ ) + (n+l) (-1)n An + +U r 2r'2 bl nn P n+2,1 a rb 2 p n ar n+1 rn+.................. (33) +B (33) nar n+l rtn+l r=rb The recursive scheme given by Hobson(1931) ( p290, Eq. 164) (n-m) pn () = (2n-1) l pn- () -(n-l+m)p...................(34) together with the initial values P = sin0 P0 =0...(35) P =0....................( is used to evaluate pn in Eq. (33). Using the definition for p'n given by Eq. (1) taking m=l and dividing it by rn+l, we have:,1 s Pn n(n+l) n rn+l r2n+1 sinO'r ('r'+ '2 - rcosa)n sin2ada................(36) Substituting the relations in Eq. (37) into Eq. (36), Eq. (38) is obtained.

16 rsinO = r'sinO' r' L'= c + riA r'2 = r2 + c2 +2 crcosO....................(37),1 1 1n n(n+l)rsinO P n = ---n+2 ---- IZ (pr+ c + ]j~2 -1 r cosa)"n l sin2ada c(r2+c2 +2 c r ) 2.................... (38) The left side of Eq. (38) is a function of r and 0, and can be readily differentiated with respect to r. Performing this operation and evaluating at r=rb, yields: 2 n+2 12- - "-1 - a rb2 P ) rb2 sin( r'bl bl ar'n+1 Inr- )~ 2i In-l} Or n+l rn+l bl +2 (2 c r bin+ 2 c 2 n 1+ c2 c rbl rIr rbl rbl bl...................(39) where n n(ntl) c In= n(n+) (+L+ +2-1 cosa)n-' sin2a da x (1rbl.................. (40) It can be proved that In can be computed by the following recursive relation I=n1 c2 c I= - (2n-1) (c — + g) In- n (1 --- 2 )I.2...................(41) rbl r bI n-i ~~~~~~~~bi1 with the initial values Io=0 11=1...................(42) 2 1 c C n+ - Since the evaluation of (1 + - +2- A ) 2 in Eq. (39) is difficult or impossible when n 2 rbl rbl is large, Eq.(39) is rearranged as:

17 c2 c n +2 1) n c - 1 a fb2 P n'b sinr. blbl ar n+1 rn+l b - n+2 n+1 rbl Ku= In -.+......... (44) birb 2 3 (1 + 2 + 2- ) 2 rb2 rb l KI can be computed by the following recursive procedure 2 rb -=. Kn= - ( (2n-1) ( + ) Kn - n Kn-2 } (n+l) (1 + 2 + 2-A ) rbl bl with the initial values of K being given by _, 1 c2 c (1+ 2 +2C K)2 r bl KO=O 0................... (46) K0 = (46) The residual normal velocity at the sphere surface produced by the second term and the series in Eq. (19) are plotted for two identical touching spheres for the cases of 30, 100, 200, 500, 1000 and 2000 terms in the velocity potential in Figures 2-7 respectively. The normal residual velocity computed by using the first 30, 100, 200, 500 and 1000 terms of the velocity potential having 2000 terms are plotted in Figs 2b through 6b. From Figs 2 to 7 and 2b to 6b it can be seen clearly that as the number of terms in the velocity potential increases the amplitude of the residual normal velocity at the sphere surface diminishes and the region having higher amplitude narrows. In Fig. 7 the largest amplitude is about 2% and the region of high amplitude is about 2 degrees. Comparing Figures 2 through 6 with 2b through 6b, one can see that the normal residual velocities in the latter group are considerably larger than those in the former

18 group. Therefore, if a velocity potential with fewer terms is desired, the coefficients in that potential should be computed specifically for that potential and it should not be constructed by simply truncating a potential that has more terms. The residual normal velocity at the sphere surface produced by the second term and the series in Eq. (19) are plotted for two identical spheres very near to each other, using 2000 terms in the velocity potential, for various separation distances, in Figures 8 - 13. When the center distance between the two spheres is greater than 2.00002, the constructed velocity potentials produce a practically zero normal residual velocity, and therefore are not plotted. It should also be noted from Figs 2 through 13 that the normal velocity as computed from Eq. (19) converges and therefore Eq. (19) converges. V. The accuracy check of the recursive schemes of Equations (34) and (41) It can be seen that in each of the two recursive procedures Equations (34) and (41) and in Equation (39) two terms are being subtracted one from the other. If in any step the two terms have the same magnitudes and with the first few digits being the same that many significant digits would have been lost in that step. Since nearly 2000 such steps were performed, it is therefore prudent to ask how accurate the results given by the recursive schemes are. To verify the value of the associated Lengendre polynomials given by the recursive scheme Eq. (34), the following direct scheme was tried: n(n+) sinJ ( + 1 osa)n lsin2a da n=n~l) ( s~ i (n-1)(n-G2) _ 32 (n -i2 1)p. I I4) (n)..(n-5( 1)2 1 1 =n(n+l) ({. n-2 2! ' 2 4 4!.ln (.12 1) 24 6+' (n-1)...(n-2k) ln_2_ 1 3 2k-1 1n+(2k)! p -1(p2 1)K2 4 " 2k 2(k+l) +...sin................. (47)

19 When the order of n is large, some terms in Eq. (47) will have a large magnitude and cause a big roundoff error. For example, when n=1000, the largest terms in Eq. (47) have the magnitude of 1089: with double precision only 16 significant digits are available, and the roundoff error is of the order of 1073. Eq. (47) therefore can not be used for large orders of n. To facilitate numerical computation, the complex integrals of the first order Associated Lengendre polynomials as defined by Eq. (3) with m=l are transformed into real integrals as given by the following Pnn(n+l) sinO (A + 2-1 cosan lsin2a da x nn-I n(n+) sin (2sin2a + cos2a ) 2 (cosco + i sino )sin2a da C ItOr'~I 2 n-I =n(n+l) sin0 J (2sin2a + cos2a) 2 cos sin2a da..................(48) where sinO (o = (n-l) Arctan(c-s cosa) = (n-l) Arctan( tanO cosa)...................(49) for cosO >0, and sine o)=(n-l)(e +Arctan(-Co cosa)) =(n-1)(x+Arctan(tanO cosa))................... (49') for cos0<0. Special care must be exercised for 0=90. Examination of Eq. (47) reveals that 0=900: 0 when n is even { (-1 n-2 n-4 3 1 (50) n-i ii n 4 2 whennisodd n-1 n-3. 2 Simpson's composite rule was first used to evaluate Eq. (48), and the result was not satisfactory for higher orders of p1, Romberg's algorithm was then used, and was much

20 superior to Simpson's composite rule. The results of the numerical integration show that the recursive scheme Eq. (34) gives at least 12 accurate digits. To verify the accuracy of the result of Eq. (39) with the recursive scheme of Eq. (41), Eq. (39) is first rearranged in the following manner so that numerical overflow can be avoided during the computation: n+2,1 n+2 a rb2 p n { (( 2 arn+l In+1 r=r 5 - 1n1 J r r r.bi. - r2 rbl rbl (1+ - + 2 ----)2 D rbl rbl c rbl nbl......(51) where n Jn - - ( + r +-1 cosa)nl1 sin2a da rt(l+c 2 c )n- rbl 2 rbl rbl................... (52) Since Jn is a complex integral, Jn is transformed to a real integral to facilitate the numerical computation: (A + C )2 + sin2 cos2a nJn f (cs I-1 --— } 2 (coso' + i sinto') sin2a da (1 + C )2 + sin2 cos2a n-l 2 |_____ _ 2 2bl 2 K {1 - )} 2 coso' sin2a da ' ( 1 + ( C )2 + 2CL, )2 rbl rbl w'= (n-l) Arctan( sin cosa ) cos0 + c rbl...................(53)...................(54) where Again, Romberg's Algorithm was used to compute the integral Jn and is used in Eq. (39). The result shows that the recursive procedure gives at least 12 accurate digits.

21 VL The interacting force on a sphere The interacting force on the spheres could be determined by integrating the local pressure on the sphere surface around either one of the two spheres, since the force on the two would be equal in magnitude and opposite in direction. However, since the basic problem is unsteady and the expression for pressure would involve a time dependent term, it would be convenient to add a uniform flow in the X direction in the velocity potential as given by Eq. (4) and thereby translate the unsteady problem of two spheres moving in the negative X direction into a steady problem of two fixed spheres in a uniform flow field. = U1 r sinO cos(p +...................(55) Since the lift force exerted on a sphere due to the interaction between the sphere and a plane is the same as the interaction force between two spheres when the sphere center-toplane distance equals half the center-to-center distance between the two spheres, henceforth the term "lift force" will be used indiscriminately for both cases. The velocity components in the r, 0 and qp directions can be determined by differentiating Eq.(55). Since the velocity potential given by Eq. (55) is not exact the r direction velocity on the sphere may not be zero. However, if it is assumed to be zero then the expression for r direction velocity could be used to reduce the expression for the velocities in the 0 and qp directions. The local pressure on the sphere surface is found from the Bernoulli Equation, and then integrated numerically around the entire surface to determine the net force on the sphere. Alternatively, the lift force could also be computed by using the equation given by Miloh, 00 F = -2r pr U2 n(n+2)2An An+l = -21 rU2f...................(56) n=l where f is the lifting force coefficient and is defined as

22 f= Z n(n+2)2 A A+l...................(57) n=1 The lifting coefficients for two identical touching spheres with different number of terms in the velocity potential as computed by using Eq. (57) are listed in Table 3 and plotted in Fig. 14. Table 3. Lifting coefficients f for touching spheres with various number of terms n in the velocity potential N f N f 10.4373060 20.4921019 30.5143508 60.5404147 100.5527179 200.5633335 250.5657013 500.5708787 1000.5738523 2000.5755495 0.6 0.5 0.4 0.3 02 0.1 0.0 i.'..........,......M....... 1 10 100 1000 N 10000 Fig. 14 The lift coefficient (f) for two touching spheres or a sphere touching a plane against the numbers of terms in the velocity potential (N) Fig. 14 shows that as N increases f increases, but the increase levels off as n becomes large and f has a finite limit as N goes infinity. This limit can be estimated by using the data given in Table 3. It can be seen from Table 3 that each time n doubles the increase in f is less than 2/3 the increase it made previously, i.e. f(4N)- f2N) (f2N) N)) f4N) - f(2N) < 3 (f(2N) -f(N))..................(58)

23 where f(N) is the lifting coefficient f for a velocity potential having N terms. It is reasonable to assume that Eq. (58) will hold for N>2000, and the upper limit of f is then given as: 2 4 f < f(2000)+ (f(2000)-f(1000) ) +- (f(2000)-f(1000))+......................(59) Upon solving the geometric series in Eq. (59), we have f<3f(2000) - 2f(1000)...................(60) or f<5789421...................(60') 37 which is approximately 3. It may be seen noted in Table 3 that the value of f for N=2000 is more than eleven per cent (11%) larger than that for N=30, while its value for N=2000 is less than.3% larger than that for N=1000. For two identical spheres with varying center distances, f computed by using Eq. (57) are plotted in Fig. 15, and some selected values and the corresponding center distance c are listed in Table 4 below. Table 4. Some selected lifting coefficients f for various center distance c c f c f 2.000000.5755125 2.000001.5742419 2.000002.5732937 2.000005.5713256 2.000007.5703674 2.000010.5691914 2000030.5642572 2.000070.5585942 2.000100.5555539 2.000500.5347092 2.001000.5206195 2.005000.4685720 2.010000.4345427 2.050000.3187479 2.100000.2526425 2.500000.0903779 3.000000.0396850 4.000000.0119794 10.00000.0003003 22.00000.0000128

24 Table 4 shows that when the center distance between two identical spheres is 10 times their radius the interaction force between the two is almost negligible. Expressing this in another way, when the sphere center is 5 times its radius from a plane the sphere almost does not "see" the presence of the plane, and the interaction between them is negligible. VII. The Lift on a Hemisphere with the Base on an infinite plane parallel to the fluid velocity The velocity potential for this case is the same as a fixed sphere in an infinite ideal fluid flow field: Urb sinO U = U r sinO cosp +cosp................... (61) 2r2 The tangential velocities in the 0 and qp direction can be obtained by differentiating the velocity potential as given by Eq. (61) with respect to 0 and Ap: ~1 X, U 3 vo I r =UcoSOcos o = UcosO c osp................... (62) b- r as I-b 2.. 1 _ U 3 v = I Usin - 2 sinp = -2 Usin...................(6 P b rsinO" acp r=-Usin2.(63) The fluid pressure on the surface of the hemisphere in the absence of earth gravity or when gravity is parallel to the plane to which the hemisphere base is attached is determined from the Bernoulli equation as: PV2 pi )2) P I P v N )2) PIrrb= Po P 2 = P ob' ((vo r=rb =Po -8Pl (cos20 cos2p + sin2)................... (64) Integrating the pressure distribution across the entire surface of the hemisphere and noting that the differential force is in the opposite direction to the surface normal, the Z-direction component of the force on the hemisphere is:

25 2 2n Fz= -j P I r cos A= rb Po+ 8 Pl U2 (cos2 cos2 +sin2p)r cos9 sinO d(p dO 0 1 =rb Po+ g PlU2 x rb |(2 + 1) d =-rp+ r pIu2.................. (65) o The first term on the right in Eq. (65) represents the force when the pressure is uniform in the fluid and there is no flow, and the second term represents the lift force induced by the motion of the fluid around the hemisphere. The lift force coefficient in this 27 case is therefore 2. In the case where the direction of gravity is opposite to that of the base plane normal, the pressure on the surface of the hemisphere is: Pitr = Po - pi U2 (cos2O cos2q + sin2p) - pi gZ r P1U (cos2Ocos2q =PO - p U2 (cos20 cos2V + sin2 ) - pl g rb sine................... (66) where PO is the pressure at the base plane, and the force component in the direction of the plane normal is: Fe 2r 2A ' 3 27 2 rz= - b + rb2 rbplU2.................... (67) b 0 3 r1Plg+ 72r'p1 The first two terms on the right in (64) are the forces exerted on the hemisphere by the hydrostatic pressure and the third term on the right represents the dynamic force on the hemisphere as induced by. the fluid flow. When the gravity is in the same direction as the base plane normal, the pressure on the hemisphere surface is: P lr= Po ' 8 Pl U2 (cos20 cos2q) + sin2 ) + p1 gZ =PO -8 p U2 (cos2O cos2p + sin2p ) + p g rb sinO......... (68) and the force on the hemisphere in the direction of the plane normal is: Fz= XrP _- 'rb l 2...................(69) -rb 0.3 ii32 1..

26 Similar to Eq. (67), the first two terms on the right in Eq. (69) are the hydrostatic force, and the third term is the dynamic force on the hemisphere. It can be seen from Eqs (65), (67) and (69) that the dynamic force or the lift on the hemisphere induced by the fluid motion is independent of the base plane orientation and the 27 lift coefficient for a hemisphere in a uniform flow field is 2. VIL RFERENCES Basset, A. B., "On the Motion of Two Spheres in a Liquid, and Allied Problems", Proceedings, Mathematical Society, Vol. 18, 1887, pp.369-377. Drew, Donald A., " The lift force on a small sphere in the presence of a wall", Chemical Engineering Science, Vol. 43, No. 4, 1988,pp.769-773. Endo, Dyuro, "The Force on Two Spheres Placed in Uniform Flow", Proc. PhysMath. Soc. Japan,Vol. 20, 1938,pp.667-703. Herman, R. A., " On the motion of two spheres in fluid and allied problems", Quartly Journal of Pure and Applied Mathematics, Vol. xxii, 1887, pp. 204-262. Hicks W. M., "On the motion of two spheres in a fluid", Philosophical Transactions, Royal Society of London, Vol. 171, 1880, pp. 455-492. Hobson, E. W., "Theory of spherical and ellipsoidal harmonics", Cambridge University Press, 1931. Kawaguti, Mitutosi, "The Flow of a Perfect Fluid around Two Moving Bodies", J. of The Physical Society of Japan, Vol. 19, No. 8, August, 1964, pp.1409-1415. Love, John D., "Dielectric sphere-sphere and sphere-plane problems in electrostatics", Quartly Journal of Mech. Appl. Math., Vol. 28, Pt. 4, 1975, pp.449-471.

27 Miloh,T., "Hydrodynamics of deformable contiguous spherical shapes in an incompressible inviscid fluid", J. of Engineering Mathematics, Vol. 11, No. 4, October, 1977, pp. 349-372. Stokes, "On some cases of Fluid Motion", Cambridge Philosophical Transactions, Vol.8, 1843. Voinov, O. V., "On the motion of two spheres in a perfect fluid", J. of Applied Mathematics and Mechanics, Vol. 33, No. 4, 1969, pp. 638-646. Witze, C. P., Schrock,V. E. and Chambre, P. L., "Flow about a growing sphere in contact with a plane surface", Int. J. Heat Mass Transfer, Vol. 11, 1968, pp. 1637-1652.

25. 20.~_..0 --- —- L-____..__ L.-L_____.._ 15.-0-_-~____ L-L ~ L - - O.O_10.0~- -__L- -L________ LLL___L,,,,__ I I/\I 5.LOL - _-__-____L-________L —_ - L --- —15.0_ __L_ L.....L_ _ _. L.._...._ __ -10.0L - -- -- L,,,, _ _ __ _ L_______ -20.0- -------— L_-_______L-__ -____L_- - - - 4 I I, -25.0 2 o5 140.0-or_.0 —/......L.......170.0......./___18_ Fig. 2 RESIDUE VELOCITY(IN PER CENT OF STREAM VELOCITY) AGAINST DEGREE FOR THE PARAMETER(NUMBER OF TERMS IN VELOCITY POTENTIAL) N-30

25. Fig.2b RESIDUAL VELOCITY(in percent of stream velocity) AGAINST DEGRE USING THE FIRST 30 TERMS OF A POTENTIAL HAVING 2000 TERMS

I 25.0. — L ~ — - - - -- - --- - -, --- - - -.C.......... L -L_........_ 20 15.C )..., _ _ I _ __ __ L - - - - - - - -- L- _- - - - - - - L - - - - - - - - - 10. OL___ --- —L --- - -— L._-_- L —. - II 7 5. _ --- —--------- -L ----- i. I I I I i I i i i.I i I i i I 0. 0. I~~ ~~rL \ \ I ~I -5.o L _ — L- -- -L --- — L- I 10. L - - - ------ ---- L L — ---- -15. - - - - -L L ---- ----— L -_ ------- -1 -20.0 O - - -- - - L -_ _- - - L - - - - — L -. I II I -j.C I I r~~~~~~~~~ I I7 -. 18 1 4. U I V. U I 0U. U I V. V Fig. 3 RESIDUE VELOCITY(IN PER CENT OF STREAM VELOCITY) AGAINST DEGREE FOR THE PARAMETER(NUMBER OF TERMS IN VELOCITY POTENTIAL) N=100

_ __ 25. 20 15 10 I I i Fig.3b RESIDUAL VELOCITY(in percent of stream velocity) AGAINST DEGRE' USING THE FIRST 100 TERMS OF A POTENTIAL HAVING 2000 TERMS.~~~~~~~~~~~~~~~~~~~~~~~~~~ _,,.

25., 20.o-..-_ __ __ L -.-__ __ L -_ _-_ _ L _I I I 15.0 ______ L-_____ L _ _____ L_ _ — _ — 10.0_ _ _ _ _ L _-_ _ _L - - -- L -. --- 5.O - - - - - - - - L -_ ___ L I I O.C -5. 0 — -- - — L _ _ _ -_ --- — L - -- __ _ __ L 10.0 -- ___ _ L.._ -_ - - - L - - - - - - -— L 15.0 - -- - - L - - ~ ~ ___ L _ I -20.0L --—. -- L - - - - - - - - - L - - Il n1 r n. 4 n 1 n 1 I V. V* I V. U I OU. U I /V. U I E Fig. 4 RESIDUE VELOCITY(IN PER CENT OF STREAM VELOCITY) AGAINST DEGRE' FOR THE PARAMETER(NUMBER OF TERMS IN VELOCITY POTENTIAL) N-200

i 25.0 I I 20.CL- -- -- L -- _ - - - L - - L - - - - - - - - -: I. 15._0 ____ L -_-_-_._.L -.- - -— L_ --- 10.0 - - - - - - - -- L - - - - - - - L. --- L 5.( 0. OL 0 - I, I I r L I I \ I \-ii;nftFtiC\(jrti i I I i I -5. - - - - - - - - L,,- -. --- L -- - -10.0 - - - - — L - - -- - L- - - - ---- L I I1 -15.0O_ — L - _.-U.-..L. — L-_,,, - ------ i I I -20.OL - - - - - - L. —.L - - --- - - -L - - - - - - - - - C I I -25. I _ _ ~~~~~~~~~rI 1 I__ I 1 4. U i. U 1 oU. U I /V.V I Fig.4b RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGRE USING THE FIRST 200 TERMS OF A POTENTIAL HAVING 2000 TERM. _ _ -

25. A.1 Vr _ I 20.O- -_ _-_ _ L _-_-_ _ —L. -- - L — _ -- 15.0 10.C 5.( I i OL 0 -0 -I. -------- L ----- -- - - - - - — 'L- -- - - -- I - L — - - - - - - - - - - - - - - L- _ -- - I I I ii i I L - I - - - I I I -- - - - L -- - - I.I I 11 - --- -.41 i I i I I Ii I I I I11 -- - - - - - - L - - - - - - - L - - -- - - - - L -10.0 — - - - - - - - -L - -- - - - - - -L --- —--— L --- —---- I I I I I I I I I -15.0 -------- L L - - - - - L ----L __ _ -- I I I I I I I I I -20.OL - - - - - - - -L - L -_ L L - ---- I 0 I I I~~~~~~ -25. 18 140. 1 DU. U IOU. U. I /V. V I Fig. 5 RESIDUE VELOCITY(IN PER CENT OF STREAM VELOCITY) AGAINST DEGREE FOR THE PARAMETER(NUMBER OF TERMS IN VELOCITY POTENTIAL) N=500,,.,,,,...,.,,.,,,,~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

25. I I iI i I I I I i I 20.0-~._ _...._ L. ~... L~ ~.... L_ L.. _ - 15.0- ~... L-_........L~ ---.L-_.O.oL _ _ _ _ _ _ - - - L _ - - - - - - - - - 0. ^v\ vN v I I ' * I I I 10.0(. -. -. - - - -—. L — --. - -— L - -----— L- ---- - - -10....... L- _ L-20.0- ~ - -.~ _ __ L -~ - L- - L — - - - - - 25.01 140.0 150.0 1b0.0 10.. Fi9.5b RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGRE USING THE FIRST 500 TERMS OF A POTENTIAL HAVING 2000 TERM

V 20.0 _ _ _ _. — - L - -- -- - ---- L -. _ ---,.- L I I 15.0 ---- - L -— _...L.- --- - L ------ - 10.C I I, ~-~ --- —. ---- ~L~I I -.- -__....__ L -.-.-_ - - - - L 5.0L - I I I i I I 'I - - - - - - - - L - I I I I L I I I a^,IAAI I 0. 0 ir.\ I I I I I I - - --- - — 1 -ft-ow I WONM Myp -5.- - - - - - -L ~ - L- - - - - L - --- ----- -1 0.O O --- —--— L - ---- -- - L ------ - -L — --- --- I I I I I I I -i5.0- -.-.. - L -. -. -. - - L -- - -... - - L -20. C -,f e. > - - - - - - - -- - - L -- ~-. ---L_ _ -2 3. 1, % ~ I,. l r I qA-r, 1 7n r 1 4V. U I 3Vu. U 1 oU. V I /U. Fig. 6 RESIDUE VELOCITY(IN PER CENT OF STREAM VELOCITY) AGAINST DEGRE FOR THE PARAMETER(NUMBER OF TERMS IN VELOCITY POTENTIAL) N=1000

25.0 20.0 _ L_........, L - L_ — 15.~C L -- L ~..... L 10.0 L- -L-LI I I 0.0 —_-lO.OC _, _ ______ L _____ L _ - - - - - - L- - -5.C 0, ___L-L - L _ ___ L -20.O0 L....._...L._.L L- -. --- L- - 1540.O I 0.0 Fig.6b RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGRE USING THE FIRST 1000 TERMS OF A POTENTIAL HAVING 2000 TERM

25. A VI 20.0 — L_ _ __ _ _ L_ - __ - - 15.0 L. --..________ L _ - - _ _ _ _ -- L 10o. ---_ ___ L _ _______L- _ _ -— L -- 5.0 ______ _ L _________ L -. _-. --- L 0.O n,. i. ll I.... I 11., II,1 -5.0 - _ _ _ _ - - - - - L I 10._C ___ __ L -____ __L — ___ __ L _- -- - 15.0 - __ __ _L __,,,_, _ L - - - L — --. --- - - I I I I I I L........._ _L....._ __ _ _ -20.0 -- _ __L ___ I _l40. - I n n II -& - 7 I V. v I JV. V I OV. V I /V. V Fig. 7 RESIDUAL VELOCITY(in percent of stream velocity) AGAINST DEGRE' FOR THE PARAMETER(NUMBER OF TERMS IN VELOCITY POTENTI'AL) N=2000

I i i I I i L 25.0 I - 20. - LU_..._ _... L.U U _ L_I I I 15.0L - - _ _ L......._.L.....-L _ _ __ 1I I I 10.C 0_____L_ _..L — _- _.... L-_ _ ___ ___L________ — 5.C I L I - - - - - - L - _- _- - - - L_ L 0.( 1.111i '11 I I I -5. - - _ _ _ L -_ _ _ L I I I -10. —I I I I I I I I I I I I I 1 1 -20.OL --- -.. L. --- - - L-_ _ _- _ -L_ - - j -2 A. I n n- 1 en n --- —- 1 n IV. V I J V. IOU. V I /V. Fig. 8 RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGRE FOR THE PARAMETER OF C/Rb=2.000001 and N=2000

i i i L 25.0 r I 20.0- _. ___. __ L _ __ __ L --- -L - - --- - I I. 15.0- ___.___ L _..._____ L.. ---- L - ----- - I I I I I I 10.0 _ ~ _.. -.I I I 10.0 - - - - - - - - _ L - - - - - _ _._ L - - _ - - - - _ L _ _-_ - - - - - _ 5.( O. 0 0 - - - - - - - - - ~- - - - - - - - -I~ L -- - - - - - - - - L L....All I -5.0- - --- — L-L-L -- -- I I I I I I I I I I I -1.0C -- - - L -- - -... -- -L- -L - - - - -- - - - 1 I. - I - - - L - -20.O --- —--- L — _-____ L -_ -- — L- - C I n I n n A 1 7 n I TV. V I; V. U I OU. U I /U. 1E Fig. 9 RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGREr FOR THE PARAMETER OF C/Rb-2.000002 ond N=2300

_ 25 n * I 20.C F. --- _L-o L - _ - -- L 15.0 --- - -- L -... L --.. L 10.0~L~L~- I I I I I I I I lO.C _ _ _ __ L.... - _. _._._.L - _ __L___ 5.( D I L I I I - - - - - - - - - L 0.< I I 0 J, -5. L ~_ ~_~_ _ L~ _ _ L _ - _ _ - _- - -10.C - - —.-.- L_ L -.. L.-L- __ — -- -15.~ _ _ _ ___L~~ _ _ L __ -L --- 0 0~ L- - - - - - L -- - -- - L - - - - - - - - 140. U 1 0. 0 160.0 Fig.10 RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGRE FOR THE PARAMETER OF C/Rb=2.000003 and N=2200 - -- - -

I 25. _ I v _ 20.0L- - - - - -- L - - - - - - - - - L - c --- -- 15.0 { _ 'I I -_ -_- - - - L - -_- - - - - - - L -_- - { _ _ 10.0 - - - - - - - - - L.. — - L - - I I - -- - - L - - -- - I 0 - - - - - - L - I I.ol 1~~~~~~~~~~~~~ I -w —4 -5. I _ _,. _ I _ _ " I _- - L _ _ _ _ _ - - L _ - - - - - - L -- I.10.0_ __ --- —--— L_ _ _____ L_ -- ---- -— _ -L -15. C0 - - - - - - - L -20.0C - - - - L I I I I I I I I I I i I I -1 -25 0I I I o400. I U. U 1 OU. U I /V. Fig.11 RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGRF FOR THE PARAMETER OF C/Rb=2.000005 ond N=2500

25. A u_ vr 20.0 -. -. -.L - -. --.- L -------— L - -- I. 15.0L- L --- L ------- 10.C 5.( C I I L D - - - - - - - - -~ L - - - - - - - L I I - - - - - - - - - L - L L _ _ __ _- _- - - L I I I I I I -5.( _ _ - -. - - L - - - -- - L - - - - L - -- ---- I I.........-L. —. --- —-L.. -10. -- --. -- L 15.0L --- — — L- ------ L — L --- —-L --- -20. 0 O I — -— _ - - - - L - - -- - - -- L - - i -25. I LI A.-~ r n:.-_ -- 1 r A 1 "7n n A 1.r 1 4. U I. I OU. U IIV * Fig.12 RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGRE' FOR THE PARAMETER OF C/Rb=2.000007 and N=2300

25 n. V 20.C L I -- L --- - - - --- L - -_- - - - - L ---- 15.0L -- --- ---— LL — ----— L --- 10.( o-. - - - - - - - - - L - - - - I _ _ _ _ L - - _ _ _ _ _ _ - L - I I I. I I I. -- -- -.L-_ -- - - - - - - LI I i I 5 I I I i C - - - - - - - - - - L - - - I I I I I - -5. C IILI I ~- - - - - - - - - - I II I 10.0 - - - - -- - L - -__ _ _ L -.L_. --- 15.C- - - - — L ----. _-L — ---- -- - L -- I -20.0 -25.0 A - _ _ _ _ _ _ _ _ - L - - - _ - - - - - L -- - - - - - - - L I I I I n A I A ri rl I 61 11 n, Lo r~~~~~~~~~l r,% ~ -70 fl 1 U. 1:. V I OU. I V.V Fig.13 RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGR FOR THE PARAMETER OF C/Rb-2.00001 and N-2200

Ii f 0.6 I1l 2 I. II0.2 I -g- - - - - - - - - - - — I I I I I I \ I I I 0 I I I I0 \ 1. I I I O.1-_ - -.... —... - --.. --- —--- I I l\ II I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I\ II I I I I 100 10.0 10.0 10.0 10.0 10.0 RbFiO. 15 LIFT COEFFICIENT (f) AGAINST CENTER DISTANCE (C/2Rb —.0, I I I I I I I I I I I I I I I I I I I I I I I I I I I 0. u _. 10 -- I

Appendix C. UNCERTAINTY ANALYSIS The uncertainties associated with the calculation of the gold film temperature, thermocouple temperature, system pressure, and total heat flux are estimated following the procedure of Kline et al. [114]. According to this procedure, the uncertainty ur in computing r is given by: Ur = ar a' Z1 ~Za zl I + A UZr J^2^ +... + Uzn )U azn znJ J (C.1) where r is a function of zl, z2,..., Zn and uzl, Uz2, u., Uz, are the uncertainties in each of these values, respectively. C. 1 Gold Film Temperature The CR7X was used for calibration of the gold film heater surfaces and for the imposed constant heat flux tests. The uncertainty in the heater surface resistance, which ultimately results in the heater surface temperature, was less for calibration of the heater surfaces than for actual power tests. From Equation (C.1) the uncertainty in the gold film resistance is expressed as: URW Rw (- W 2 (+ R A 2h J f '21 VW +j ^Vs + Rsh (C.2) Values of the quantities in equation from a representative calibration are: uw = 0.6gV UVsh = 0.6gV Vw = 120mV Vsh = 50mV URSh = 0.001l R= sh 3.760gO These values substituted into Equation (C.2) result in an uncertainty of 0.03%. For power runs, the above values are different and result in a greater uncertainty. Representative values for a power run are: 1

uyV = 0.05mV uVsh = 0.6gV V = 120mV Vh = 50mV URsh = 0.0000020 Rsh = 0.016314Q with a resulting uncertainty in the heater surface resistance of 0.11%. Calibration data with the largest scatter would be expected to provide a conservative estimate of the uncertainty in the slope. From the calibration data of surface Q-2, the corresponding temperature and resistance data with the largest deviation from the fit line are 54.85~C at 2.8801 Q2 and 83.59~C at 2.0171 KQ. An estimate for the uncertainty in the slope of the linear relationship between the surface temperature and the resistance can be found by using Equation (C.1) together with the definition of this slope and the above values for temperatures and resistances, respectively. The uncertainty in the slope for these values is + 0.23~C/i. The uncertainty in the surface temperature measurement can now be estimated with a known estimate of the uncertainty in the heater surface resistance and of the slope in Equation (6). Carrying out the indicated operations of Equation (C. 1) on the equation for the surface heater temperature, Equation (6), the following expression for uncertainty in the temperature measurement results: ( 2 2 2M1/2 UTw = (UTc) + (um) + (u(Rw-RC.3) Representative values for the quantities in equation are: T 0050C = 0.0035 um = 0.23~/c.Q UR 0.003Q UC =.0f 2

These values result in an estimated uncertainty in the surface temperature measurement of about + 1 ~C using the CR7X data acquisition system. C.2 Thermocouple Temperature The thermocouple spool from which all thermocouples used here were fabricated was calibrated against a platinum resistance thermometer which had an uncertainty in temperature of + 0.001~C, and a calibration equation for the spool was determined from a fourth order polynomial least square fit of the measured thermocouple voltage and the determined platinum resistance thermometer temperature. In terms of the uncertainty in the CR7X voltage measurement, which was + 0.6 plV, the uncertainty in the temperature determination is estimated to be better than + 0.05~C ( 0. 1~F). C.3 System Pressure There was little discrepancy, on the order of + 0.005 psi between the Ruska pressure determination, with any uncertainty in pressure of + 0.0003 psi and the calculated pressure resulting from the calibration equations. The uncertainty in pressure resulting from the use of these equations was then taken as + 0.005 psi. From Equation C.1: up=(. VPI (C.4) where P is the calculated pressure from the calibration equation and V is the pressure transducer voltage signal. For a typical run the following values were obtained: 3

UV = mV = 9.8608 - 2.7352 x 10-2 V Dv - P = 16psi V = 1.6V The uncertainty in pressure at this level is estimated from equation to be 0.06%. C.4 Total Heat Flux The expression for q"T was given earlier in Equation (7). Using Equation C.1 the expression for the uncertainty in q"T is: UqT (( + VT URsh~ + 1/2:VT + Vsh J RshJ A (C.5) In the experiments using the CR7X the following were typical values: UVT = 0.3pxV VT= 5V UVh = 0.3 gV Vsh = 0.08 V URSh = 0.000002 Q Rsh = 0.016314 Q UA = 0.08 cm2 A = 7.25 cm2 With Equation C.5 and these values, the estimated uncertainty in the calculated q"T is 1%. 4

Appendix D. Tabulation of Forced Convection Transient Incipient Boiling Data Run No. Fl F2 F3 F4 F5 F6 F7 F8 F9 F10 F I F12 F13 F14 F15 F16 F17 F18 F19 F20 F21 F22 F23 F24 F25 F26 F27 F28 Filename p (kPa) T. (~C) 40 40 39.9 ATsub (~C) 12.5 12.5 12.6 qt (W/cm2) 2.50+0.03 4.30+0.03 6.75+0.06 FAVT0515.212 FAVT0515.412 FAVT0515.712 a/g 1* tspr (sec) (sec) V V V V V Tw* (OC) Tspr (~C) Flow Rate (1/min) Surface No. FAVT0516.212 FAVT0516.412 FAVT0516.712 FAVTO527.213 FAVTO527.413 FAVT0527.807 FAVT0527.707 FAVT0527.707 FAVT0527.407 FAVT0527.207 FAVT0527.207 FAVT0530.807 FAVT0530.407 FAVT0530.407 FAVTO530.207 FAVTO530.212 FAVTO530.207 FAVTO530.207 FAVTO530.807 FAVTO530.812 FAVTO530.412 FAVTO530.212 FAV'r0530.2 12 119.49 119.56 120.93 119.83 119.83 119.76 119.83 119.83 119.76 119.9 119.76 119.83 119.83 119.83 119.83 119.83 119.83 119.83 119.83 119.83 119.83 119.83 119.83 40.6 40.6 40.7 45.6 45.6 45.6 45.5 45.6 45.6 45.6 45.6 45 45 45 45 44.8 44.9 44.9 44.9 40.3 40.3 40.3 40.3 12 12 11.8 7 7.1 7 7.1 7 6.98 7 7 7 7.6 7.6 7.6 7.8 7.7 7.7 7.7 12.3 12.3 12.3 12.3 2.38:0.02 4.21:0.03 6.60+0.03 2.49+0.03 4.51+0.02 7.95+0.05 7.51+0.06 7.50+0.04 4.37+0.03 2.60+0.04 2.60+0.03 7.83+0.06 4.60+0.05 4.61+0.02 2.66+0.03 2.50+0.20 2.58+0.02 4.60+0.02 7.63+0.04 7.51+0.05 4.41+0.03 2.58+0.01 2.58+0.03 V 4.164 V 6.968 V 1.352 V 31.146 V 4.825 V 1.525 V 0.182 V 0.303 V 3.182 V 9.722 V 7.412 V 1.674 V 6.152 V 3.346 V 14.322 V 19.408 V 6.824 V 2.48 V 0.369 V 1.606 V 3.971 V 8.152 V 17.272 4.176 6.974 1.36 31.156 4.83 1.53 1.996 1.657 3.196 9.744 7.424 1.677 6.154 3.354 14.336 19.42 6.832 2.49 0.378 1.611 3.981 8.16 17.824 62 84 80 82 85 91 53 75 91 88 83 90 89 80 78 80 85 88 87 94 87 80 77 62 84 80 82 85 91 102 102 91 88 83 90 89 80 78 80 85 88 87 94 87 80 77 6846.15 6846.15 2641.62 2641.622 2631.729 2602.05 2641.622 2641.62 2641.62 2641.62 2631.73 6737.32 6737.32 6737.32 6737.32 6737.32 6737.32 6737.32 6737.32 2641.62 2641.62 2641.62 2641.62 Q13 Q13 Q13 Q16 Q16 Q13 Q13 Q13 Q13 Q13 Q16 Q16 Q16 Q16 Q16 Q13 Q13 Q13 Q13 Q13 Q13 Q16

Appendix D. Continued Run No. F29 F30 F31 F32 F33 F34 F35 F36 F37 F38 F39 F40 F41 F42 F43 F44 F45 F46 F47 F48 F49 F50 F51 F52 F53 F54 F55 F56 F57 F58 Filcnamc FAVT0530.412 FAVTO530.812 FAUT0607.812 FAUT0607.412 FAUT0607.212 FADT0607.212 FADT0607.412 FADT0607.412 FADT0607.812 FADT0608.813 FADT0608.413 FADT0608.213 FAUT0608.213 FAUTO608.413 FAUT0608.813 FAUT0608.213 FAUT0608.807 FAUT0608.813 FAUT0608.413 FAUT0608.213 FADT0608.213 FADT0608.213 FADT0608.813 FAUT0610.808 FAUT0610.408 FA UT0610.208 FADTO610.408 FADTO610.408 FADT0610.808 FADT0610.208 P (kPa) 119.83 119.83 120.38 120.45 120.45 120.45 120.45 120.45 120.45 120.58 120.58 120.58 120.58 120.58 120.58 120.58 120.63 121.03 121.03 121.12 121.18 121.2 121.2 120.72 120.76 120.76 120.95 120.95 120.95 120.76 T.0 (~C) 40.3 40.3 40.2 40.2 40.2 40.2 40.2 40.2 40.2 39.4 39.4 39.4 39.4 39.4 39.4 39.4 45.3 45.6 45.6 45.7 45.5 45.5 45.4 44.7 44.7 44.8 44.8 44.8 44.7 44.8 ATsub (~C) 12.3 12.3 12.6 12.6 12.6 12.6 12.6 12.6 12.6 13.4 13.4 13.4 13.4 13.4 13.4 13.5 7.6 7.3 7.3 7.3 7.3 7.5 7.6 8.2 8.2 8.2 8.1 8.1 8.2 8.1 q" (W/cm2) 4.48+0.03 7.84+0.03 8.11+0.03 4.67+0.04 2.86+0.02 3.01+0.02 5.45+0.05 5.4+0.1 8.7+0.1 7.79+0.06 4.64+0.04 2.67+0.03 2.57+0.04 4.57+0.06 8.48+0.10 8.49+0.04 7.90+0.06 7.86+0.05 4.47+0.04 2.60+0.01 2.70+0.02 4.37+0.05 7.93+0.07 7.95+0.07 4.60+0.04 2.67+0.03 4.68+0.04 4.70+0.04 7.89+0.05 2.58+.04 a/g 1* (sec) V 4.842 V 1.282 1 1.18 1 8.515 1 32.024 -1 8.81 - 1 2.7525 -1 2.108 - 1 0.039 -1 0.16 -1 2.131 - 1 6.298 1 46.76 1 13.968 1 1.777 1 1.714 1 1.912 1 1.769 1 11.906 1 33.568 -1 3.15 -1 3.464 - 1 0.205 1 1.831 1 6.448 1 27.552 -1 8.039 - 1 6.276 - 1 0.447 - 1 55.22 Ispr (sec) 4.85 1.286 1.207 8.524 32.034 8.81 2.7525 2.108 0.538 0.676 2.131 6.298 46.78 13.974 1.784 1.721 1.912 1.776 11.912 33.584 3.15 3.464 0.205 1.836 6.458 27.588 8.041 6.279 0.44 7 55.22 Tw* (~C) 83 86 83 87 78 90 96 91 42 52 79 75 73 86 83 78 87 88 89 76 75 94 66 92 85 88 111 107 74 109 Tspr (~C) 83 87 83 87 78 90 96 91 67 76 79 75 73 86 65 61 75 74 62 76 75 94 66 92 90 88 111 107 74 109 low Rate (1/min) 2691.09 2671.3 6638.39 6658.18 6658.18 6658.18 6658.18 6658.18 6658.18 2671.3 2661.41 2661.41 2661.41 2671.3 2671.3 2671.3 2701 2701 2701 2701 2701 2701 2700.98 6658.2 6658.2 6658.2 6658.2 6658.2 6658.2 6658.2 Surface No. Q16 Q16 Q16 Q16 Q16 Q13 Q13 Q13 Q13 Q13 Q13 Q13 Q16 Q16 Q16 Q16 Q16 Q16 Q16 Q16 Q13 Q13 Q13 Q16 Q16 Q16 Q13 Q13 Q13 Q13 CD C0 -ro-o I, 01,