Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Oil flow in the oilwell tube annulus of vertical bearing assemblies Piao, Yinghu 1997

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1997-196399.pdf [ 8.74MB ]
Metadata
JSON: 831-1.0080938.json
JSON-LD: 831-1.0080938-ld.json
RDF/XML (Pretty): 831-1.0080938-rdf.xml
RDF/JSON: 831-1.0080938-rdf.json
Turtle: 831-1.0080938-turtle.txt
N-Triples: 831-1.0080938-rdf-ntriples.txt
Original Record: 831-1.0080938-source.json
Full Text
831-1.0080938-fulltext.txt
Citation
831-1.0080938.ris

Full Text

OIL FLOW IN THE OILWELL TUBE ANNULUS OF VERTICAL BEARING ASSEMBLIES By Yinghu Piao B.Sc. & M.Sc. (Thermo-Energy Engineering) Tianjin University, China M.A.Sc. (Mechanical Engineering) University of British Columbia A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF D O C T O R OF PHILOSOPHY in T H E FACULTY O F G R A D U A T E STUDIES DEPARTMENT OF MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA 1997 © Yinghu Piao, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mechanical Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1Z1 Date: Abstract Oil leakage from the oilwell tube annulus of self-contained vertical bearing assemblies has been of practical industrial concern for many years. Earlier observations of bearing behavior indicated that this leakage was associated with the complicated physical pro-cesses that describe the flow of oil in the bearing's reservoir. In this latest study, we seek to obtain a better understanding of this oil flow, particularly in the annular clearance space surrounding the oilwell tube. This will be done by means of numerical simulation and experimental flow visualization. A test rig was designed and built for the purpose of simulating the bearing's oilwell tube annulus. A major feature of this rig was to provide visual access to the annu-lar clearance space, and also to the region beneath the rotating runner where strong secondary flow effects are known to exist. A light sheet visualization technique, using micro air bubbles as the tracer, was the primary method for tracing the secondary flow pathlines. In particular, the effect of runner speed on these pathlines was investigated. Appearances of the air bubbles in the light sheet (radial-axial plane) were timed and analyzed in order to quantitatively examine the velocity of the oil flow, and to compare the experimental results with numerical data. The oilwell tube has an axisymmetric geometry. Accordingly, a three dimensional CFD code for laminar, axisymmetric flow with a free surface (A3D-CFD code) was de-veloped for the purpose of analyzing the flows in the oilwell tube. This code was verified using results from the experimental study, and also using data obtained from classical analytical studies of free surface shapes, velocity distributions and flow patterns. The experimental and analytical data were in good agreement with the computational results. ii Using the A3D-CFD code, oil flow patterns in the oilwell tube annulus of the test rig were computed. This study showed that the primary (circumferential) swirling flow was well stratified inside the oilwell tube annulus, while a strong secondary flow beneath the runner was in the form of a vortex. The strength of this secondary vortex (in terms of the local secondary flow vorticity u>ec at the core of the secondary vortex) was also studied. The non-dimensional vorticity (U>«C/VR) was correlated with the Reynolds number (Ren). A numerical technique was developed for the purpose of tracing micro air bubbles in the oil flow field computed from the A3D-CFD code. By considering buoyancy effects due to the gravitational and centrifugal forces, the pathline of an air bubble was determined. Various parameters that affect the path of the air bubble, such as bubble size, oil viscosity and runner speed, were included in this study. The computed pathlines of the air bubble were in reasonably good agreement with the experimental results. iii Table of Contents Abstract ii List of Tables vii List of Figures viii Nomenclature xii Acknowledgement xvi 1 Introduction 1 2 Experimental Apparatus 13 2.1 Simulated oilwell tube test rig 14 2.2 Experimental apparatus 14 2.3 Experimental procedure 19 2.4 Optical distortion 20 2.5 Air bubble size 22 2.6 Dynamics of an air bubble in purely rotating flow 24 2.7 Air bubble size and oil static pressure 25 3 Experimental Results and Discussions 27 3.1 Image processing and restoration 27 3.2 Time sequence 29 3.3 Pathlines at different locations 32 iv 3.4 Pathlines at different runner speeds 35 4 Computational Fluid Dynamics Code 37 4.1 2D CFD code: SOLA-VOF 37 4.2 The A3D-CFD code . 38 4.2.1 Governing equations 39 4.2.2 Finite difference discretization 40 4.2.3 Free surface cells and boundary conditions 43 4.2.4 Wall boundary conditions 45 4.2.5 Obstacle boundaries 47 4.2.6 Algorithm 48 4.2.7 Shape of the free surface 52 4.2.8 Secondary flow stream function 54 4.2.9 Secondary flow vorticity 55 4.3 Verification of the A3D-CFD code 55 4.3.1 Flow near a disk rotating in a fluid at rest 55 4.3.2 Taylor vortices 58 4.3.3 Shape of free surface in a rotating annulus 61 5 Computational Modeling of Flow 67 5.1 Modeling of flow in the simulated bearing assembly 67 5.2 Modeling of flow in the simulated oilwell tube annulus 70 6 Computational Results and Analysis 73 6.1 Test oil properties 73 6.2 Flow in the simulated vertical bearing assembly 74 6.3 Computed flow field in the oilwell tube annulus 79 v 6.4 Pathline of an air bubble 81 6.4.1 Effect of bubble size 92 6.4.2 Effect of viscosity 94 6.4.3 Effect of runner speed 101 6.5 Secondary vorticity at the vortex core 106 7 Concluding Remarks and Future Work 108 7.1 Concluding remarks 108 7.2 Future work 110 Bibliography 111 Appendices 115 A Image restoration of optical distortion in the experiment 115 B Dynamic analysis of a micro air bubble in a purely rotating flow 119 C Air bubble size affected by the oil static pressure 125 D Subroutines in Figure 4.23 128 E Free surface cell and its side ratios 131 F Dimensional analysis of shape of free surface 133 G Analysis on isobaric lines in a rotating annulus 135 H Computation of trace of an air bubble in the computed flow field 138 vi List of Tables 4.1 Summary of wall boundary conditions used in the CFD code 48 6.2 Measured kinematic viscosity coefficients of the test oil 73 6.3 Time t of reappearances of air bubbles in the (r-z) plane 94 6.4 Time t for reappearance of the air bubble in the (r-z) plane 98 6.5 Time t for reappearance of the air bubble in the (r-z) plane 105 6.6 Calculated secondary vorticity togc (rad/s) at core of the secondary vortex. 106 vii List of Figures 1.1 A schematic of vertical bearing assembly 2 1.2 A schematic of a simulated guide bearing assembly 3 1.3 A schematic of concentric rotating cylinders 5 2.4 Simulated oilwell tube test facility 15 2.5 Section of the simulated oilwell tube test rig 16 2.6 A schematic of the experimental apparatus 17 2.7 Optical distortion (scale units: mm) 21 2.8 Air bubble sizes in flow visualization 23 2.9 Air bubble sizes in still oil 23 3.10 An example of image processing: (a) before processing; (b) after processing. 28 3.11 An example of image restoring: (a) before restoration; (b) after restoration. 29 3.12 The first appearance of air bubbles (U>R = 3.14 rad/s) 30 3.13 The second appearance of the air bubbles (WJJ = 3.14 rad/s) 31 3.14 The third appearance of the air bubbles (LOR = 3.14 rad/s) 31 3.15 Secondary pathline of the air bubbles (t = 59 seconds, WR = 3.14 rad/s). 32 3.16 Secondary pathlines for different radial locations of the air bubble probe (u)R = 4.19 rad/s) 33 3.17 An image of the air bubbles as the tip of the probe was close to the oilwell tube wall (wB = 3.14 rad/s) 34 3.18 Secondary streamlines at different runner speeds: (a) UR = 3.14 rad/s; (b) U>R = 4.19 rad/s; (c) U>R — 5.23 rad/s 35 viii 4.19 An illustration of the staggered grid 40 4.20 A finite difference cell and its neighboring cells 42 4.21 A free surface cell 44 4.22 Meshed computational domain and wall boundaries 47 4.23 Computational flow chart of the modified A3D-CFD code 50 4.24 Free surface representation 53 4.25 Location of ip in a finite difference cell 54 4.26 Velocity distributions close to the center of the disk 57 4.27 Computational and analytical velocity distribution of a rotating disk. . . 57 4.28 Computed Taylor vortices: (a) stream function; (b) secondary velocity vectors; (c) tangential velocity contours. 59 4.29 Schematic of test facility for free surface shape measurement 62 4.30 Computed and measured free surface shapes 63 4.31 Computed and measured height differences of free surface 63 4.32 Non-dimensional shapes of free surface 65 5.33 Modeling of the simulated bearing assembly 69 5.34 Modeling of the simulated oilwell tube annulus 71 6.35 Kinematic viscosity of the test oil versus temperature : 74 6.36 Velocity distributions in the simulated bearing assembly (WR = 15.7 rad/s, v = 1.82 x 10 - 4 m2/s): (a) secondary velocity vectors; (b) primary velocity contours 75 6.37 Static pressure distribution beneath the runner 77 6.38 Computed flow field (ReB=711): (a) secondary flow vectors; (b) tangential velocity contours 80 6.39 Computed secondary flow streamlines (Re#=711) 81 ix 6.40 Close-up views of computed flow in flow visualization region (Re#=711): (a) secondary flow vectors; (b) tangential velocity contours 82 6.41 Computed location and time history of an oil particle (Refl=711) 84 6.42 Computed pathline of an oil particle (Refl=711) 84 6.43 Computed location of the air bubble, and its time history (ReJR=711; da = 0.1 mm) 86 6.44 Computed paths of the air bubble and an oil particle (Refl=711) 87 6.45 Appearances of air bubbles in a flow visualization test (Ren=711) 88 6.46 Computed pathline of the air bubble after the "distortion" (Re.R=711; da = 0.1 mm). . 90 6.47 Computed pathlines ("distorted") of the air bubble at different starting points (Refl=949; da = 0.1 mm) 91 6.48 Computed pathline ("distorted") of the air bubble (da = 0.1 mm) at dif-ferent runner speeds: (a) ReH = 711; (b) Re f i = 949; (c) ReR = 1185. . . 92 6.49 Location of air bubble and the effect of bubble size (Re.R = 711) 93 6.50 Pathlines affected by air bubble sizes (Re^ = 711) 94 6.51 Computed secondary flow fields at a constant runner speed O>R = 3.14 rad/s: (a) i/=1.235 xl0- 4 m 2 /s; (b) i/=1.09 xl0~ 4m 2/s; (c) i/=0.75 xlO" 4 m2/s; (d) u=0A0 xlO" 4 m2/s 96 6.52 Computed velocity profiles at r=150.2 mm (u;^  = 3.14 rad/s; da = 0.1 mm). 97 6.53 Computed velocity profiles at r=124.0 mm (a;^  = 3.14 rad/s; da = 0.1 mm). 99 6.54 Locations of the air bubble as affected by oil viscosity (WR — 3.14 rad/s; da = 0.1 mm). . . 100 6.55 Pathlines of the air bubble as affected by viscosity (a;^  = 3.14 rad/s; da = 0.1 mm). . 100 6.56 Computed velocity profiles at r=150.2 mm (z/=7.5 x 10~5 m 2 /s ; <ia=0.1 mm) 102 6.57 Computed velocity profiles at r=124.0 mm (i/=7.5 x 1 0 - 5 m 2 /s ; da=0.1 mm) 103 6.58 Locations of the air bubble as affected by runner speed (u = 7.5 x 1 0 - 5 m 2 /s; da = 0.1 mm) 104 6.59 Pathlines of the air bubble as affected by runner speed (y = 7.5 x 1 0 - 5 m 2 /s ; da = 0.1 mm) 105 6.60 Secondary vorticity at the core of the vortex vs. Re# 107 A.61 Illustration of the true scales and its distorted image 115 A.62 Optical compression ratio Cr vs. rd 117 A.63 Graphic expression of the optical correction function 118 E.64 A free surface cell 132 E . 65 Three different cases for a free surface cell (NF(ij) = 3; TANTH(i j) > 0). 132 F . 66 Illustration of the rotating annulus 133 G. 67 The annular radius ratio coefficient K vs. the annular radius ratio <f>. . . 137 H . 68 Interpolation of velocity at an arbitrary point P 141 H.69 New location of the air bubble after the time interval At 142 xi Nomenclature A constant a acceleration B constant C radial clearance Cr optical compression D drag d diameter F Fraction of volume of fluid (Eq. 4.9), also force / side ratio in a free surface cell Fr Froude number g acceleration of gravity; also acceleration of other body force H height h height i , j cell numbers in r-, z-directions cell numbers neighboring i,j io,jo largest cell numbers K annular radius ratio coefficient (Eq. G.70) / length m mass n number counter p pressure R radius xii r radial (horizontal) distance Re Reynolds number r, z, 6 radial, axial, circumferential in cylindrical coordinates Ta Taylor number t time U terminal velocity component in r direction u,v,w velocity components in r-, z-, 8 directions V terminal velocity; also terminal velocity in z direction v velocity; also velocity component in z direction w velocity component in 8 direction z axial (vertical) distance a over-relaxation parameter 8 small number 9 circumferential angle fi absolute viscosity coefficient v kinematic viscosity coefficient p density a surface tension coefficient r shear stress ip stream function <j> annular radius ratio (R2/Ri) u angular velocity; also vorticity J7, V functional relationships xiii subscript: 1 inner cylinder of an annulus 2 outer cylinder of an annulus a air bubble b buoyancy c core of vortex; also centrifugal d optical distortion; also drag force E east e east f free surface i inner cylinder of an annulus max maximum min minimum N north n north; also normal to the wall o outer cylinder of an annulus R the runner; also radius r radial direction; also relative s south s secondary; also south tc transient characteristic (Eq. B.49) W west w wall boundary; also west e tangential xiv superscript: n time counter at previous computational cycle n + 1 time counter at present computational cycle * non-dimensional parameter xv Acknowledgement I owe an enormous debt of gratitude to my supervisors, Prof. K. R. Brockwell and Prof. S. M. Calisal, for their help, encouragement and patience through all stages of this study. Their guidance and expert supervision have been invaluable in finishing this study successfully. My gratitude also goes to Mr. J. Ferguson of General Electric of Canada for his generous assistance in designing and manufacturing the test facility for this project. I wish to express my sincere appreciation to Dr. I. S. Gartshore and Dr. R. L. Evans for their insightful suggestions and discussions. My special thank goes to Dr. M. J. Hinchey of MUN for his kind help of making the SOLA-VOF code available. My grateful appreciation also goes to Mr. D. Kleinbub, Mr. C. Gerber, Mr. T. Vanderhoek and Mr. C. Trickett of the National Research Council of Canada, Mr. S. Oshika, Mr. T. Besic, Mr. A. Schreinders and Mr. A. Steeves of the Mechanical Engineering Department for their technical assistance. I would like to specially thank my school friends, especially Dr. S. Duan and Dr. P. He, for their valuable academic discussions and help on the experimental and numerical aspects of this study. I would also like to take this opportunity to thank my families for their persistent faith in me and my friends for their various kinds of support. Financial support for this research by the Natural Sciences and Engineering Research Council of Canada, General Electric of Canada, National Research Council of Canada is gratefully acknowledged. xvi Chapter 1 Introduction In large vertical rotating machinery, such as a hydro-generator, the massive rotating shaft is supported axially by a pivoted shoe thrust bearing and located radially by a pivoted shoe journal bearing. Both the thrust and guide bearings are usually contained in a bath of oil, which also contains an integral cooler to maintain safe bearing operating temperatures. In this self-contained bearing assembly, proper oil flow is essential to maintain proper lubrication of, and to dissipate heat from, the bearing surfaces. A schematic of a vertical bearing assembly is shown in Figure 1.1. The oilwell tube forms the inner wall of the oil bath and would normally prevent leakage of oil from the bath into the machine. Historically, previous designs of oilwell tube have been successful in preventing such leakage. However, in recent years, serious leakage problems have arisen as a result of increases in both the size and operating speed of the machines. Oil leakage from the oilwell tube is a practical problem that has been of some con-siderable concern to hydro-generator manufacturers for a number of years. In 1975, Berggren and Engelhart [3] of General Electric, Schenectady, NY, USA, first reported oil leakage from a large vertical motor bearing reservoir. In 1978, General Electric Canada also noticed oil leakage during shop testing of a large vertical guide bearing. This leak-age seems to be associated with a sudden loss of flow stability of the oil in the oilwell tube annulus, resulting in a serious loss of oil from the bearing and contamination of the machine underneath the bearing. Usually, oil aeration is associated with this process [5], 1 Chapter 1. Introduction 2 Guide Bearing Cooling Coils o 'o o o o o o o Figure 1.1: A schematic of vertical bearing assembly. and may ultimately result in damage to the bearing surfaces as a result of the change in oil properties. The cost associated with terminating machine operation for maintenance and replacement of spare parts is very high. Growing concern about this oil leakage problem initiated a joint research between the National Research Council of Canada and General Electric of Canada in 1988. GE Canada built a test facility to visually observe the oil flow in the oilwell tube of a bearing assembly. The oilwell tube was made from transparent plexiglass, and enabled visual access to the oil flow in the oilwell tube annulus from the bottom of the rig. A para-metric study on the flow instability of the oilwell tube of vertical bearing assemblies was completed in 1992 [5]. A schematic of the test rig used in this project is shown in Chapter 1. Introduction 3 Figure 1.2. Supporting frame Pseudo-guide A bearing Figure 1.2: A schematic of a simulated guide bearing assembly. From this experimental study, it is now possible to describe the physical processes associated with the oil flow in the oilwell tube annulus. At a constant rotational speed, oil in the oilwell tube annulus primarily swirls around the axis of the shaft as a result of the viscous shear stresses exerted by the contacting surfaces of the rotating runner. At the same time, oil flows from the oilwell tube annulus through a "pseudo-guide bearing" clearance into the oil reservoir and then returns to the oilwell tube annulus through a series of radial holes located in the base of the bearing assembly. The oil in the reservoir is cooled by a heat exchanger (cooling coil). This recirculation of oil is essential for removing heat from the bearing surfaces and also to maintain effective lubrication of the bearing. Chapter 1. Introduction 4 At low rotational speeds, the oil is clear and transparent and the free surface takes on a distinctive profile as illustrated in Figure 1.2. As the speed increases, oil in the reservoir becomes aerated, which passes into the oilwell tube annulus through the radial recirculation holes. When the speed is further increased beyond a certain limit, the oil-air interface in the oilwell tube annulus suddenly degenerates into a violent sloshing mixture of air and oil, which may lead to oil leakage over the top of the tube. Also, it is noted that oil moves into the annulus from the recirculation holes and returns to the oil reservoir through the guide bearing clearance as a result of flow recirculation. Due to limited visual and lighting accesses to the flow in the oilwell tube annulus, only primary (swirling) flow could be observed through the transparent oilwell tube. However, this provided some preliminary understanding of the problem. It is clear that the flow in the oilwell tube annulus is a very complicated process, and is closely associated with the runner speed. Also, while this study could not identify the primary cause of the flow instability, some qualitative relationships in relation to the onset of flow instability could be identified [5j. • Not a clear relationship between the tube annulus clearance and the onset of flow instability (critical runner speed). • Oil viscosity has little influence on the critical runner speed. • Still oil level has a clear influence. The higher the still oil level, the higher the critical speed. • Isolating the upper part of the oilwell tube annulus from the lower part (at the bottom of the runner) results in a significant increase in the critical runner speed. • From non-dimensional criteria, it does not seem to be possible to predict the critical runner speed. Chapter 1. Introduction 5 These earlier studies suggest that the problem is not simply related to just the swirling flow (with free surface) in the oilwell tube annulus, even though the flow instability seems to first appear at the free surface. The onset of flow instability is greatly affected by the boundary condition at the bottom of the oilwell tube annulus. In the present study, we seek to obtain a better understanding of this oil flow, par-ticularly in the annular clearance space surrounding the oilwell tube. Basically, the geometry of the oilwell tube annulus is formed by two concentric cylin-ders as in the case of the classic fluid dynamics problem where flow between two concentric rotating cylinders is considered. However, this problem is complicated by the presence of boundary conditions at the bottom of the annulus and an oil-air interface, namely the free surface. The classical concentric rotating cylinder problem usually assumes in-finitely long cylinders with the annulus filled with a homogeneous fluid. Figure 1.3 shows a schematic of concentric rotating cylinders. / / n ; ' :::::::::::::: / ' •*•':*• 'M ' J: / / z — " " <r ^ * / /* / v i Figure 1.3: A schematic of concentric rotating cylinders. A classic solution for this flow satisfies the Navier-Stokes equations in an infinite Chapter 1. Introduction 6 annulus and the no-slip condition at the cylinder walls. The steady solution to the flow is simply given by: ve = Ar + — (1.1) r where A - ~ ^ B - - r i r l ( u o - U i ) ' o t o x ?*i and r0 are the inner and outer radii of the annulus, respectively and u>,- and u)Q are the angular speeds of the inner and outer cylinders, respectively. This simple physical system has attracted the interest of scientists for well over a century. The attraction to such a simple system is basically twofold. One, this simple system can create very complicated physical phenomena at the point where the fluid flow loses its stability. Two, because of the simple geometry of the cylinders, there is a good chance of analytically or numerically modeling the complicated unstable flow processes [2]-In an experiment on a stationary outer cylinder and a rotating inner cylinder, Couette [9] observed that the torque exerted by the fluid on the outer cylinder was proportional to the velocity of the inner cylinder. However, when the velocity exceeded a certain critical value, the torque increased at a greater rate than the velocity. This was attributed to a transition of the fluid flow pattern from a steady laminar motion into a somewhat more complicated flow. This was one of the earlier works that detected flow instability in such concentric rotating cylinders. Rayleigh [33] developed a theory for inviscid rotational instability, assuming axisym-metric geometry. Stability of flow is ensured if ^-{rv6)2 = 2r*wuz > 0, (1.2) dr Chapter 1. Introduction 7 where ui = v$/r is the angular velocity and is the local vorticity. Rayleigh's criterion asserts that the axisymmetric rotating flow is stable when the rotation u> and the vorticity uz have the same sign. Applying the steady-state solution from Equation 1.1 to Rayleigh's criterion, the local vorticity u>z = 2A and -!L(rv6)2 = 4r3uA. dr It should be noted that Rayleigh's criterion is solely dependent upon the sign of the constant A. Using Rayleigh's criterion, the following two cases, which are of practical interest, are examined. Case (a) : Inner cylinder is rotating, while outer cylinder is stationary; [ Case (b) : Outer cylinder is rotating, while inner cylinder is stationary. Case (a) yields a negative value of A which leads to unstable flow according to Rayleigh's criterion, no matter how small the rotational speed is. Case (b) yields a positive value of A which leads to stable flow according to Rayleigh's criterion, no matter how large the rotational speed is. This inviscid theory only indicates the tendency of flow towards instability. However, viscosity does have some effect on the flow in the annular space. A viscous instability analysis was performed by Taylor [36] for a small-gap problem. Taylor found out that as the cylinder speed exceeded a certain critical value, characterized as the onset of flow instability, the secondary flow took the form of cellular toroidal vortices spaced regularly along the axis of the cylinder. These are the so called Taylor vortices. The critical Taylor number for the onset of such an instability is given as, Chapter 1. Introduction 8 Without the small-gap assumption, the energy stability analysis [20] predicts that circular Couette flow will be stable if v{rl~rl) (log r- 0/ri) 2' This criterion from the energy stability theory depends on the magnitude of the angular velocity \ui0 — u>i\. It does not distinguish case (a) and case (b), as mentioned above. Nevertheless, it does predict the Taylor vortex type instability quite well. In other words, the energy stability theory holds when u>j > u?a. The unsymmetric characteristic of the circular Couette flow to the interchanges of u>i and a>0 was observed by Mallock [24] and Couette [9]. Later, Coles [8] systematically dis-cussed the unsymmetric behavior and distinguished the transition as being two separate motions. The first kind of transition is by spectral evolution and is a characteristic of the motion when the inner cylinder has a larger angular velocity than the outer one, i.e. case (a). The second kind of transition, called catastrophic transition, is characteristic of the motion when the outer cylinder has a larger angular velocity than the inner one, i.e. case (b). The first kind of transition is what has been described earlier as the onset of Taylor vortices. It should be noted that in this first kind of transition, as the Taylor number increases, the vortices begin to develop circumferential waves but remain laminar. Tur-bulence ensues at Ta = 160,000 as given by Coles [8], which is much higher than the critical Taylor number that predicts the onset of Taylor vortices. In other words, for case (a), the critical angular speed (u>i\Turbulence) f ° r turbulent transition is about ten times the critical angular speed Baylor) f ° r the onset of Taylor vortices. According to Rayleigh's criterion, case (b) is stable. But it is clear from experimen-tation that all laminar flows are conditionally stable and become unstable at a finite limit. With the inner cylinder stationary and the outer one rotating, flow in the annulus Chapter 1. Introduction 9 becomes unstable as the speed of the outer cylinder exceeds a certain limit even though the critical speed might be much higher than that for the first kind of transition. This leads to the second kind of transition, defined by Coles [8]. As the outer cylinder speed is slowly increased, the flow remains laminar and steady up to a relatively high speed. As the speed is further increased and eventually exceeds a certain limit, patterns of al-ternating laminar and turbulent flow will occur with almost explosive suddenness. This intermittent turbulent flow turns into fully turbulence as the speed is increased further. In the first kind of transition, the critical speed for the turbulent transition is much lower than that for the second kind of transition. Wendt [39] investigated laminar to turbulent transition in circular Couette flow. By interpreting Wendt's results for the case r0/ri = 1.176, the critical Reynolds number for the laminar to turbulent transition is given by 2 Re; = ^ i - « 2,000 (when co0 = 0) and 2 Re0 - « 100,000 (when = 0) It should be noted that based on the earlier description, the critical Reynolds number for the onset of Taylor vortices in such a case will be of the order of 200. Most of the work on circular Couette flow neglects end effects by assuming cylinders of infinite length, even though, in reality, end effects cannot be eliminated for a finite cylinder. Nevertheless, even such a simple geometry gave rise to complex flows. There are relatively few studies on the end effects which complicates the problem further. One such study was conducted by Cole [7], who looked at the influence of the end effect on Taylor vortex formation and transition. It was found that the end effect did not have significant influence on the onset of Taylor vortices, as compared to the onset of the wavy vortices at a higher speed. Another study by Burkhalter and Koschmieder [6] Chapter 1. Introduction 10 on a similar problem also examined end effects on the vortex structure. Wild et al [41] investigated the end effect on the axial distribution of shear stress, thus the torque on the cylinder, and found that the higher the aspect ratio (rotor length to radial gap width), the lower the averaged azimuthal shear stress. Another relevant problem is the endwall effect on rotating flow in such a cylindrical container that the top endwall, bottom endwall as well as the, cylinder sidewall are all independent components and may rotate independently [4, 11, 15]. As the aspect ratio of the height to the radius of the cylindrical container increases, the problem changes from one which is rotating disk dominant to one which is rotating cylinder dominant. This particular geometry somewhat resembles the geometry underneath the bottom of the runner of the present study. Referring to the simulated vertical bearing assembly shown in Figure 1.2, it can be seen that the runner, which rotates in the oil bath, basically separates the oil into two regions, the oilwell tube annulus and oil reservoir. Comparing to the classical circu-lar Couette flow, the oil in the oilwell tube annulus may experience the second kind of transition while the oil in the oil reservoir (in the annulus between the runner and the perforated baffle if the later is present) may experience the first kind of transition. How-ever, this present research will mainly focus on the flow in the oilwell tube annulus where the oil leakage problem is of particular interest. The oilwell tube annulus is basically formed by the outer surface of the stationary oilwell tube and the inner surface of the rotating runner. Although there are some geometric similarities to the classic circular Couette flow problem, the oilwell tube annulus also has some significant differences. For example, the oilwell tube annulus has finite axial length, and a very complex geometry in the region close to the bottom of the runner. Also, there is an air-oil interface, namely the free surface, in the annulus. These characteristics of the oilwell tube annulus result in complex flow phenomena. Ferguson [12] examined Chapter 1. Introduction 11 the "end effects" resulting from the complex geometry at the bottom of the annulus and concluded that these conditions would have a large influence on the behavior of the oil in the annulus. This was confirmed by a separate experimental study, by Brockwell and Kleinbub [5]. In a vertical bearing assembly, the air-oil interface may drastically deform and entrain air into the oil at higher operating speeds. This type of oil aeration in a bearing assembly is usually characterized by the occurrence of a creamy yellow oil-air mixture. Because of the recirculation of oil in the bearing assembly, this milk-like mixture is transported into other regions of the bearing assembly, until eventually all oil in the bearing assembly reaches an equilibrium state and is equally aerated. In this state, this aerated oil may be treated as a homogeneous mixture. The aerated oil may contain as much as 15% by volume of dispersed air [13], which implies that the density of the homogeneous mixture may decrease by as much as 15%, compared to that the un-aerated oil. One practical difficulty remains in obtaining a satisfactory definition for the viscosity of the two-phase mixture, although it is known that, in general, the dynamic viscosity of the mixture is less than that of the oil. Both Isbin et al. [19] and Duker et al. [10] have attempted to model the viscosity of such two-phase mixtures. In a parametric study on a simulated vertical bearing assembly, various oilwell tube annulus clearances and runner geometries were tested. Kaiser [21] presented the ex-perimental results in a non-dimensional form and attempted to link the flow instability to a critical criterion, i.e. Froude number or Taylor number. The suitability of this type of correlation was questioned by Brockwell and Kleinbub [5]. This is in accordance with a separate experimental study on a thrust bearing test rig reported by Monk [26]. These reports concluded that there was no correlation between the critical speed and the Reynolds number, the Taylor number or the Froude number. Furthermore, no over-all general conclusions could be drawn since the appearance and the speed at which Chapter 1. Introduction 12 the instability occurred seemed to vary according to changes made to the oilwell tube geometry. In this type of vertical bearing assembly, as the annular clearance of the oilwell tube is changed (usually by increasing or decreasing the bore of the runner), other geometric dimensions are usually kept unchanged. Therefore geometric similarity is already vio-lated. Furthermore, the complex geometry at the bottom end of the oilwell tube annulus may not be well represented by one single significant length (i.e. the radial clearance of the annulus). In this sense, to find a unified criterion based on the annular clearance that describes the onset of instability does not seem to be possible. From a parametric point of view, previous studies on the flow instability in vertical bearing assemblies have given some insight as to the complexity of the problem. This complexity can definitely be attributed to the complex flows that are set up as a result of the geometry (both in the tube and below it). Therefore, it is very important to know the flow patterns as a whole, and velocity distributions in particular, that are present in the bearing assemblies. The oil flow in the oilwell tube annulus is the main focus of the present study. To achieve this objective, a new test rig with a transparent runner was built to enable observation of both the primary and secondary flow in the oilwell tube annulus. Furthermore, a CFD code for 3D axisymmetric flow has been developed to simulate the oil flow in the bearing assemblies and, in particular, in the oilwell tube annulus. Chapter 2 Experimental Apparatus The main purpose of the experiments is to better understand the flow patterns in the vertical oilwell tube annulus of a vertical self-contained bearing assembly, particularly where there is the potential for serious leakage problems to take place at high operating speeds. The experimental results are later compared with computed results from the axisymmetric 3D CFD code developed for this study. To maintain continuity of the contents, the earlier test facility, the simulated vertical guide bearing assembly as shown in Figure 1.2, is also described. The oilwell tube of this.earlier test rig was made of a transparent plexiglass to provide visual access to the flow in the oilwell tube annulus. Basically, flow visualization was limited to the primary swirling flow in the annulus. A series of parametric studies was conducted, and some physical processes relating to the oil were observed. During operation, oil flows from the oilwell tube annulus through a pseudo-guide bearing clearance into the oil reservoir and then returns to the oilwell tube annulus through a series of radial holes. The oil is cooled by a heat exchanger (cooling coil) in the reservoir. At low operating speeds, the oil is clear and transparent and the free surface takes on a distinctive profile. But as the operating speed increases, oil in the reservoir becomes aerated and passes into the oilwell tube through the recirculation holes. When the operating speed is further increased beyond a certain limit, the oil-air interface in the oilwell tube suddenly degenerates into a violent sloshing mixture of air and the oil, which may lead to leakage of oil over the top of the tube. 13 Chapter 2. Experimental Apparatus 14 From the earlier study, it was deduced that the transition from a stable to an unstable state of flow is a very complicated process. Accordingly, it was felt that this current study was needed to concentrate on, and understand, the flow patterns in the oilwell tube annulus where the flow transition initially takes place. This decision has led to the design and construction of a second test rig which incorporates a simulated oilwell tube. 2.1 Simulated oilwell tube test rig The purpose of the simulated oilwell tube test rig is to focus attention on the three dimensional characteristics of the oil flow in the oilwell tube annulus. This test facility, which is shown in Figure 2.4, was vertically mounted on the laboratory wall. The test section consists of a stationary oilwell tube unit and a rotating runner unit with the stationary oilwell tube unit suspended from the top frame. The rotating runner unit was mounted on the top end of the shaft and driven by a variable speed motor through a gear-belt speed reduction unit. Details of the test section are shown in Figure 2.5. The runner cylinder, which forms the outside wall of the oilwell tube annulus, is made of transparent plexiglass and enables the flow in the oilwell tube to be visually observed. A separate cylinder, containing oil underneath the runner, is also made of transparent plexiglass, where secondary flow is relatively strong. The stationary radial vanes are designed to reduce the bulk swirling motion of the fluid caused by the rotating disk upon which the top end of the shaft was mounted. 2.2 Experimental apparatus In addition to the primary swirling flow, the secondary vortex flow may play an important role in terms of the flow transition from a stable to an unstable state and is Chapter 2. Experimental Apparatus 15 Figure 2.4: Simulated oilwell tube test facility. Chapter 2. Experimental Apparatus 16 Runner Pseudo-guide bearing Vanes Figure 2.5: Section of the simulated oilwell tube test rig. also of interest. A schematic of the experimental apparatus is shown in Figure 2.6. The video camera was set about 37 centimeters from the center of the image plane to the camera lens. A close-up lens was used to magnify the image. Light sheet visualization was a commonly used flow visualization technique, which was first attempted in high-speed gas flow by McGregor [25] in 1961. Using light sheet generated by laser, Wei et al [38] successfully captured the Gortler vortex in Taylor-Couette fluid flow. This optical technique is usually used together with a tracer method by injection of foreign particles such as aluminum powders [1], fluorescent dye and hy-drogen bubbles [35]. In the present experiment, a vertical light sheet was projected into Chapter 2. Experimental Apparatus 17 Pressure Gauges Valve 1 Figure 2.6: A schematic of the experimental apparatus. Chapter 2. Experimental Apparatus 18 the oilwell tube annulus region, while micro air bubbles were injected into the oil flow as tracers. As a micro air bubble passed through the plane illuminated by this sheet of light, it scattered the light, which could be visually observed. A Panasonic Super VHS video camera was used to record the image of the bright air bubbles. The recorded images on a video tape were digitized with a software called CAPTURE and the digitized im-ages were further processed with a software called IMG WORKS on a Silicone Graphics computer afterwards. Due to swirling feature of the flow, the same air bubble would reappear in the illuminated plane periodically. Also, by continuously adding tracers at the same location, a stream of bubbles would be released. Reappearance of a bubble in the illuminated plane would be at a different location in the presence of secondary flow. After a sufficient lapse of time, many reappearances of the bubbles were completed. These locations of the reappearances would he on a secondary pathline of the bubbles in the r-z plane. By observing development of these bubbles in the illuminated plane, information on the secondary and primary flows can be obtained. The sheet of light was produced by a slide projector. A slide was made such that light was only allowed through a vertical thin slit, resulting in a sheet of light that was about 1.5 mm thick. The air bubble injection probe was made from a glass capillary tube of 2.0 mm OD and 1.16 mm ID. The capillary tube was heated at its center with a electric coil and was pulled at the two ends until it broke into two pieces, resulting in two cone shaped sections with very thin tips. Each of these two pieces could be used to make an air bubble probe. Since the outer diameter (OD) of the cone shaped section gradually reduced towards the tip end (so did the inner diameter), by properly cutting off the tip end, the OD of the tip could be chosen. After some experimental trials, the OD of the tip, measured at about 20 - 30 pm using a microscope, was used in the experiments. The ID was not known, but would be even smaller. This probe was located at the lower part of the oilwell tube annulus and was arranged so that it could be slid in and out radially. Chapter 2. Experimental Apparatus 19 Therefore, air bubbles could be released at different radial locations. The runner was driven by a motor-reducer assembly. The motor was an Elektrim a.c. motor, rated at 5.0 horsepower at an output speed of 3420 rpm. The reducer was a Sew-Eurodrive R60P reducer with a fixed speed reduction ratio of 5.16. The belt pulley transmission system provided additional speed reduction of 2.0. Therefore, the total speed reduction ratio from the motor to the rotating runner was 10.32. The motor speed was governed by a Parajust a.c. motor speed controller, at a speed range of approximately 10~100% of its rated speed, providing a runner speed that ranged between 33 and 330 rpm. 2.3 Experimental procedure Initially, the air bubble probe was moved to its correct location. Subsequently, valve-1 was opened to check the pressure in the compressed air cylinder and throttling valve-2 was slightly opened to ensure that air was discharging out at a very low rate before connecting the air supply to the air bubble probe. Valve-3 was opened before submerging the probe in the oil to avoid the back flow of the oil into the probe. Then the oil was poured into the test chamber from the center cylinder to a preset level. After immersion, the air bubble probe was checked to ensure that the air was flowing smoothly into the oil. At this stage, it was important to keep the air flowing until the testing commenced, since the probe could be easily clogged by back flow of oil into the orifice, when the back end of the probe was open to the ambient pressure without back pressure or valve-3 was closed for a long period of time, say, a few hours or more. When this happened, the probe was very difficult to clean, and, in most cases, a new probe was required. Before commencing the test, the light source and video camera were both positioned in their correct locations, and valve-3 was closed. The oil chamber was checked for air Chapter 2. Experimental Apparatus 20 bubbles that accumulated on the underside surface of the runner after a previous test. It was noted that these small bubbles could join and form larger bubbles resided on the underside surface. These air bubbles were removed with a brush. At this time, the motor could be started and adjusted to a correct speed. The speed of the runner was monitored with a tachometer directly pointing at the shaft. There was a few minutes wait to allow both the runner and oil flow reach a steady state. Valve-3 was opened to inject air into the oil flow. Video-recording was conducted in a dark environment, after switching on the slide projector and illuminating the region that was of interest. The stream of air bubbles passing through the illuminated plane was then recorded on video cassette tape. This experimental procedure was repeated for several runner speeds and different locations of the air bubble probe. 2.4 Optical distortion The dimensions of the fluid domain of the recorded image are 43 mm and 12.7 mm in the radial and axial directions, respectively. A significant problem associated with the present flow visualization study is the optical distortion of the recorded image due to both the curvature of the transparent rotating cylinder and the positioning of the video camera. According to the experimental set-up, the image will be distorted in the radial (horizontal) direction, but not in axial (vertical) direction. An example of the distorted image is shown in Figure 2.7. In the figure, a meshed board was placed in the oil at the same location that the images of flow visualization were taken. The meshed board was divided by cells of 5 mm by 5 mm squares originally. After optical distortion, the cells were non-linearly compressed in the radial direction such that the radial dimension of the fluid domain (43 mm) reads about 26.5 mm (see Figure 2.7). Scales in bright (white) lines were overlaid Chapter 2. Experimental Apparatus 21 Figure 2.7: Optical distortion (scale units: mm). on the recorded image. The radial axis was scaled 1:1 to the axial scale and the origin of the coordinates was at the lower left-hand corner of the fluid domain shown in the image. Accordingly, the radial coordinate rj of the distorted image needs to be corrected to the actual radial coordinate r. The distortion may be described by a rather complicated mathematical procedure. But here a simple physical method is used. Details are shown in Appendix A. From Eq. A.44, the radial coordinate 7V in the distorted image is converted to the actual radial coordinate r, r = 0.000332 r3d + 0.004114 r\ + 1.295 rd, where both r and rd are in millimeters. It should be noted that in the above formulation, both r and r<£ originate from the lower left-hand corner of the fluid domain in the image, not from the axis of the rotating body. Chapter 2. Experimental Apparatus 22 2.5 Air bubble size Air bubble size is a very important parameter in terms of its dynamics in the rotating oil flow as analyzed in Appendix B. In the flow visualization, the air bubbles were quite small and their approximate sizes were evaluated from a photographic method. A close-up lens was attached to the video camera to magnify the image, and a scale with a 1 mm increment was placed in the oilwell tube annulus. The experimental procedure used was as follows. The motor was started and the runner speed was adjusted within the experimental range. After a sufficient amount of air bubbles was present in the oil, the rig was stopped and there was a pause to allow for the oil to stop moving. The air bubbles would gradually move upward as a result of the buoyancy force due to the gravity. Bubble sizes were then measured using the scale. A wire of 0.37 mm diameter was also placed vertically in the oil so that the horizontal width of the bubble could be determined. A photograph showing the relative sizes of air bubbles as the runner rotates (U>R = 3.14 ~ 5.19 rad/s) is shown in Figure 2.8. The oil temperature was 27 °C, corresponding to a viscosity of v = 0.75 x 1 0 - 4 m2/s. Figure 2.8 indicates that air bubble sizes were generally around, or less than, 0.1 mm in diameter. However, it should be noted that air bubble sizes in still oil were much larger (about 0.3 mm in diameter), as shown in Figure 2.9. Again, the diameter of the vertical wire shown in Figure 2.9 was 0.37 mm. The reason for this difference in size is that when the oil is flowing across the tip of the air bubble probe, the motion tends to break the air stream from the probe into smaller bubbles. In the experiment, the air bubble size and the air flow rate through the air bubble probe can not be deliberately chosen since it depends on many factors such as shape of the air bubble probe, bore at the probe tip, back pressure of air of the probe as well as oil flow properties. However, the back pressure and the bore of the probe can be deliberately Chapter 2. Experimental Apparatus 23 Figure 2.9: Air bubble sizes in still oil. Chapter 2. Experimental Apparatus 24 chosen. As discussed in the following section, a small air bubble makes a better tracer, which implies that small bore of the probe is desirable. However, in the experiment, if the air bubbles are made too small, which usually corresponds to a small air flow rate, the air bubbles, illuminated by the light sheet, may not be bright enough to be properly observed, especially in the less transparent oil (compared with water). Therefore, proper probe size and the back pressure of the air were determined by experimental trials. 2.6 Dynamics of an air bubble in purely rotating flow As described in Appendix B, the air bubble gains both radial and axial velocity as a result of gravitational and "centrifugal" buoyancy forces as it follow the main swirling flow. Obviously, in order to have the air bubble accurately trace the oil flow, the ra-dial and axial velocity components must be small compared to the tangential velocity component. Based on the solutions (Eq. B.47, Eq. B.52 and Eq. B.54) given in Appendix B, the velocity components of the air bubble may be written as As time t —>• oo, the three velocity components approach to their terminal velocities, ua ->• Ua = Ylv r ' Va Va = 12i/' WA —»• Wa = W. Chapter 2. Experimental Apparatus 25 The transient characteristic of the air bubble in the oil flow can be evaluated by deter-mining the time t taken for the air bubble to reach 99 percent of the terminal velocities from rest, which is given by Equation B.49, _4.6d» The time ttc is proportional to the area of the bubble and is inversely proportional to the kinematic viscosity of the oil. Therefore, this suggests that a small air bubble is desirable from the point of view of fast response of the air bubble to the dynamic change of the oil flow. For an air bubble of diameter da = 1.0 x 10 - 4 m in oil with kinematic viscosity v = 7.5 x 10 - 5 m2/s, the resulting time is ttc = 2.6 x 10~5 second. Hence, this shows that the air bubble responds very quickly to dynamic changes in the oil flow. Eventually, the air bubble will follow the primary swirling flow, although it will also gain radial and axial velocity components. The magnitude of these secondary velocity components is proportional to the term (d\jv). Therefore, a small air bubble in a viscous fluid such as oil will make a good flow tracer. Also, the bubbles should be placed in a region where the tangential velocity is sufficiently high compared to Ua and Va. Therefore, the bubbles should not be placed too close to a stationary wall or object. 2.7 Air bubble size and oil static pressure Air bubble size is an important factor for its dynamic behavior in the rotating oil flow. However, air bubble size in the oil varies with the pressure acting on its surface by the oil. Static pressure of the oil is stratified and changes from place to place. As an air bubble travels in the oil flow, it also experiences a pressure change applied by the oil, hence, the size change of the air bubble. In Appendix C, possible size change of an air bubble in the oil flow is discussed, which Chapter 2. Experimental Apparatus 26 concludes that the possible relative change of an air bubble is ^ x 100% = -0.042% da at the runner speed UIR = 5.19 rad/s. The minus sign implies that oil pressure increase leads to a decrease in the air bubble size. At lower runner speeds, the pressure change of oil will be less, which leads to smaller value of (Ada/da). From the calculation, it can be concluded that a change in air bubble size due to the oil static pressure in the experimental conditions can be neglected. Chapter 3 Experimental Results and Discussions The purpose of the experimental study was to visually observe secondary flow in the oilwell tube annulus, particularly in the region underneath the bottom of the runner, where the secondary flow was believed to be strongest. 3.1 Image processing and restoration The image recorded on the Super VHS video tape was selectively digitized with the image software CAPTURE provided by Silicon Graphics computer. The development of the micro air bubbles in the flow visualization tests can be clearly replayed on a television monitor from the recorded tape. However, the digitized image could be developed from only one single frame of the images. Each single frame of image had quite poor quality in terms of visual presentation, even though the quality of the images during continuous play looked much better. Therefore, the digitized images were further processed with the software IMAGWORKS also provided by Silicon Graphics computer. Images before and after the image processing are shown in Figure 3.10. In the image, the left bright vertical line was the intersection of the light sheet and the outer surface of the stationary oilwell tube cylinder, while the right bright vertical line was the intersection of the light sheet and the inner surface of the transparent rotating cylinder. The sheet of light was projected from the right side of the image and the axis of the concentric cylinders was somewhere to the left of the image from the view point of the camera. The area contained by the two bright lines is the fluid domain. The relative 27 Chapter 3. Experimental Results and Discussions 28 (a) (b) Figure 3.10: A n example of image processing: (a) before processing; (b) after processing. location of the image in the test section can be referred to in Figure 2.6. In the fluid domain, bright spots in the image are considered to be the air bubbles illuminated by the light sheet. It can be seen that these bright spots form a closed path that is considered to be the secondary flow pathline of the air bubbles. As mentioned earlier, the pathline shown is optically distorted and the distortion can be restored using the correction function, Eq. A.44. A n example of restoring the secondary pathline using this correction function is shown in Figure 3.11. Figure 3.11(a) was derived from Figure 3.10. To extract the pathline from the photo image point by point is a lengthy business. Therefore, in later comparisons between the experimental and computational pathlines, instead of restoring the experimental pathlines, the computed pathlines are deliberately distorted by simulating the optical distortion inherent in the experiments. Chapter 3. Experimental Results and Discussions 29 (a) (b) Figure 3.11: An example of image restoring: (a) before restoration; (b) after restoration. 3.2 Time sequence Ideally, as a tracing particle is injected into the swirling oil flow, it travels with the oil flow and reappears in the r-z plane illuminated by the light sheet. Due to the secondary flow, the tracing particle may reappear at a different location in the r-z plane, but still on the same secondary streamline of the previous appearances. If the magnitude of the secondary flow is small in comparison to the primary flow, the neighboring spots of the reappearances should be close to each other. After a sufficient time period, with continuously adding air bubbles into the oil stream, the illuminated points will indicate the secondary flow streamline. Also, from the time interval between each neighboring reappearance, the average velocity of the secondary flow can be calculated. However, the experimental scenario was quite different. Let's examine the time history of the air bubble stream as it appears in the illuminated r-z plane. First, take the instant that the air bubble stream first appears in the illuminated plane as a reference, and set the starting time t = 0. The first appearance of the air bubble stream in the illuminated plane was quite concentrated at one bright spot. A picture taken just after the air bubble stream first appears in the illuminated plane is shown in Figure 3.12. Chapter 3. Experimental Results and Discussions 30 Figure 3.12: The first appearance of air bubbles ( U R = 3.14 rad/s). The air bubble stream travels with the swirling oil flow, and circulates back to the illuminated plane. The second appearance occurred at t = 3.90 seconds and was less concentrated at the "spot". Figure 3.13 shows a picture just after the second appearance of the bubbles. The third appearance occurred at t — 7.65 seconds. After the third appearance, the bright "spots" were hardly distinguishable. Figure 3.14 shows a picture that was taken just after the third appearance of the bubbles. In most cases, because neighboring spots were quite far apart, it was impossible to determine the average velocity of the secondary flow. However, after a sufficient period of time, the bright spots would form an enclosed loop, which indicated the secondary pathline of the air bubbles. A secondary pathline traced by the air bubbles is shown in Figure 3.15. Chapter 3. Experimental Results and Discussions 31 Figure 3.14: The third appearance of the air bubbles (U>R = 3.14 rad/s). Chapter 3. Experimental Results and Discussions 32 Figure 3.15: Secondary pathline of the air bubbles (t — 59 seconds, U>R — 3.14 rad/s). 3.3 Pathlines at different locations In the experiment, the position of the air bubble probe could be changed, thus permit-ting the air bubbles to be released at different radial locations. At one particular location of the probe, the secondary pathline of the air bubbles could be traced. Accordingly, a series of locations of the air bubble probe were tested, and images for these different radial locations of the probe are shown in Figure 3.16. The radial location of the air bubble probe in stationary oil could be measured directly. However, because the thin tip of the probe was quite flexible and would bend in the flow of the oil, the actual radial location of the probe was not known. Nevertheless, what was of concern was the location of the air bubbles appearances in the illuminated plane, rather than the location of the probe. In the experiments, the radial location of the probe was adjusted before each test, when the oil was still. Starting with, a distance of 10 mm from the oilwell tube wall, the tip location of the air bubble probe was increased in 5 mm Chapter 3. Experimental Results and Discussions 33 Figure 3.16: Secondary pathlines for different radial locations of the air bubble probe (ufR = 4.19 rad/s). Chapter 3. Experimental Results and Discussions 34 increment radially between one test and another. In Figure 3.16, figure (a) corresponds to a tip location of 15 mm from the oilwell tube wall (in still oil). Following figures ((b) to (f)) are for 5 mm incremental steps of tip location. It was observed that if the probe was too close to the oilwell tube wall, the air bubbles appeared in a rather random fashion in the image. A n image recorded when the tip of the probe was close to the oilwell tube wall (10 mm from the tube wall in still oil) is shown in Figure 3.17. Figure 3.17: A n image of the air bubbles as the tip of the probe was close to the oilwell tube wall (w« = 3.14 rad/s). When the tip of the probe was less than 10 mm from the oilwell tube wall, the air bubbles could hardly be observed in the illuminated plane. This might be attributed to the fact that the tangential velocity of the oil was relatively small in this region, where, as discussed earlier, air bubbles would not make good tracers of flow. Chapter 3. Experimental Results and Discussions 35 3.4 Pathlines at different runner speeds In the experiment, different runner speeds were tested at each location of the air bubble probe. The lowest runner speed was UJR = 3.14 rad/s and was increased in an increment of 1.05 rad/s. Figure 3.18 shows images at different runner speeds at one location of the air bubble probe. (c) Figure 3.18: Secondary streamlines at different runner speeds: (a) OJR = 3.14 rad/s; (b) U R = 4.19 rad/s; (c) OJR = 5.23 rad/s. It can be seen that, at lower runner speeds, the air bubbles traced a clear secondary Chapter 3. Experimental Results and Discussions 36 flow pathline. However, as the runner speed increased, the pathline became more blurred. This may be attributed to fact that the air bubbles move faster at higher runner speeds. Thus, the time it takes for a bubble to pass through the thin light sheet is shorter, resulting in a smaller exposure time of the bubble in the illuminated plane. Another possibility might be an increase in the randomness of the flow at higher speeds. Chapter 4 Computational Fluid Dynamics Code 4.1 2D C F D code: S O L A - V O F Initially, the two dimensional CFD code - SOLA-VOF, developed by the Los Alamos Scientific Laboratory [27], was employed for this study. The SOLA-VOF code solves for two-dimensional transient laminar fluid flow with free boundaries and uses the concept of fractional volume of fluid (VOF) [17] to treat arbitrary free boundaries. SOLA-VOF uses an Eulerian finite-difference method and employs rectangular cells with the capability of generating variable cell sizes. In its original form, SOLA-VOF solves the r- and z-momentum equations together with the continuity equation. For incompressible Newtonian flow, the governing equations in cylindrical coordinates may be written as: The r-momentum (radial) equation: (4.4) The z-momentum (axial) equation: (4.5) The continuity equation: (4.6) Convective time derivative in 2D cylindrical coordinates: •V) = u— + v d_ dz (4.7) 37 Chapter 4. Computational Fluid Dynamics Code 38 Laplacian operator in 2D cylindrical coordinates: v = i£(4U£ ,4 , ) r dr \ dr I dz2 For free fluid boundaries (free surfaces), special boundary conditions need to be imposed. However, the free surface, which is developing with time for a transient flow, can not be automatically tracked from solving the above governing equations only. Special numerical treatment is required to track the free surface. Marker and Cell method was first proposed by Harlow and Welch [17]. In this method, a set of marker particles is used to indicate fluid configuration, thus the free surface. A more flexible and efficient method for treating complicated free boundary configurations, Volume of Fluid (VOF) method, was reported by Hirt and Nichols [18]. In the VOF method, the fractional volume of fluid (F) is defined such that when the cell is fully occupied by fluid, F = 1; and when the cell is empty, F = 0. When 0 < F < 1, the cell contains a free surface and the value of F represents the fractional volume of the cell occupied by the fluid. The free boundary is traced by updating the F function. The F function is governed by the equation, dF dF dF n 1r- + u— + v— = 0. 4.9 Ot Or oz In cylindrical coordinates, the 2D SOLA-VOF code assumes that there is no flow (zero velocity) in the circumferential direction. 4.2 The A 3 D - C F D code In the present problem, primary flow is in the circumferential direction as a result of the rotating shaft, and is the source of the secondary flow in the r-z plane. Therefore, the 2D CFD code, in its original form, is not appropriate for the three dimensional flow problem associated with the oilwell tube test rig. For the present study, a 3D Chapter 4. Computational Fluid Dynamics Code 39 axisymmetric flow CFD (A3D-CFD) code has been developed [30, 31, 32] and is described as follows. 4.2.1 Governing equations The one prominent feature of most rotating machinery, and the oilwell tube annulus of a vertical bearing assembly, is the axisymmetric geometry. Usually, the assumption of axisymmetric flow leads to significant simplifications of the numerical solution procedure. In the following development of the A3D-CFD code described below, an axisymmetric flow is assumed, and the development is based on the above mentioned 2D code. For 3D axisymmetric flow, the governing equations, in a cylindrical coordinate system with the coordinates r, z and 0, take the forms (radial) r-momentum equation: * + ( V . V ) „ - V = - I | + * + „ ( V . - £ ) , (4.10) (axial) z-momentum equation: ^ + (V.V)t, = --p^z+9z + U (V2*) ' ( 4 > 1 1 ) (circumferential) -^momentum equation: ^ + (V-V)w + ^-uw = u (v2w - ^ , (4.12) and the continuity equation: I | . M + | . w = 0 . ( 4 . 1 3 , For an axisymmetric three dimensional flow, the convective time derivative (V«V) and the Laplacian operator V 2 keep the same form as the 2D code. To solve these governing equations for a particular flow, the wall boundary conditions must be defined. If there is a free surface involved, special conditions also need to be imposed at the interface. Chapter 4. Computational Fluid Dynamics Code 40 4.2.2 Finite difference discretization The Eulerian finite-difference method is used to discretize the governing equations, based on rectangular cells in r-z plane. The 2D rectangular cell corresponds to a rectan-gular ring in the 3D cylindrical coordinate system due to the axisymmetric assumption. Staggered grid is used for both velocity components u and v. Therefore, for a rectangular cell, velocity components u and v are calculated for points that he on the sides of the cell; while other scalars (i.e. p, F etc.) are calculated at the center of the cell. The velocity component, w, is also calculated at the center of the cell. If not specially mentioned, "i" and "j" in the formulation refer to the discretized cell numbers in the r and z directions, respectively. A typical staggered grid is shown in Figure 4.19. * Al-(i) • • i v(ij) ) - -u(i-- • — T -1J) (i,j) u( » - • — 4 v(i,j-1) * - • AZ 0) I f i i i i Figure 4.19: An illustration of the staggered grid. The discretized form of the r- and z-momentum equations are not presented here, Chapter 4. Computational Fluid Dynamics Code 41 but can be found in reference [27]. The -^momentum equation can be rewritten as dw dluw) dlvw) 2 1 / 1 d{r2Tr6) drze dt dr dz where the shear stresses are given by dr dz (4.14) TT8 = P 8_ dr (=) TZ8 = P dw dz (4.15) (4.16) In this formulation, the viscous dissipation term in Eq. 4.12, is rewritten in terms of the viscous shear stresses. Such a formulation facilitates the imposition of boundary conditions when the free surface is present. A finite difference cell and its neighboring cells are shown in Figure 4.20. Cap-ital letters "E" (east), "W" (west), "N" (north) and "S" (south) are the centers of the neighboring cells and the small letters "e"," w"," n" and "s" represent the four sides of the cell For a homogeneous full fluid cell, the finite-difference approximation of the -^momentum equation is u>+(hj) ~ w{hJ) , u(i,j)we - u(i-l,j)ww At Ar(i) | v(ij)wn-v(i,j-l)wa w(ij)[u(ij) + u(i-lj)] Az(j) r(i) 1 ([r27>o]e ~ V2rre}w [rz6\n - [rzg]/ p \ [r[i)YAr{i) * Az(j) t (4.17) where At = t+ — t is the time step and w+ is the velocity at time t+. In the formulation, velocity at each side of the cell is unknown since velocity w is calculated at the center of each cell in the CFD code. By linear interpolation, velocity at a side of cell is obtained from the velocity at the cell and the velocity at the neighboring Chapter 4. Computational Fluid Dynamics Code 42 zN zG) z.J Ar(i-1) r, i W Ar(i) T i i n w (U) -r • Ar(i+1) -H r(i) T AzG+1) AzG) AzG-1) 1 Figure 4.20: A finite difference cell and its neighboring cells. cell that has the same common side. Velocities at the four sides of the cell are given by Ar(i + 1) Ar(i) J \Ar(i) Ar{i + 1) f ww w(^J)\ ( 1 -l WR, W, Ar(t - 1) + Ar(t) ) [ Ar(i) + Ar(i - 1) J ' W N + ^ 7 # U ^ + ^ ^ T V 1 , (4-20) kAz(j + l) Az(j)J \Az(j) Az(j + 1), ws , w(z,j)\ / 1 1) -1 v Az(j - 1) + Az(j) ) \Az(j) + Az(j - 1) J ' ( 4 > 2 1 ) From Eq. 4.15 and Eq. 4.16, the shear stress terms are formulated as \rK 1 - 2 ^ r e (^B _ . . L "Je Ar(i) + Ar{i + l)\rE r(i) )' { ' ' \ 2 1 2V rw ( _ (A OI) L r 7 V*L Ar(i) + Ar(t - 1) ^ r(i) r w ) ' ^ ' ; Chapter 4. Computational Fluid Dynamics Code 43 M » = A,(,-) + A . f j + l ) ( " * - "<'•>'»• ( 4 - 2 4 ) ^ ' • ' A ^ l + M J - l ) ^ — ( 4 - 2 5 ) However, for a cell containing a free surface, additional procedures have to be taken to reflect the free surface boundary conditions such as zero shear stress condition on the free surface. 4.2.3 Free surface cells and boundary conditions The shape of the surface in a free surface cell is represented by the following three parameters from the CFD code: F(i,j) represents the fraction of fluid in a cell, NF(i,j) indicates which side of a cell the fluid is on, TANTH(i j) shows the slope of the free surface with respect to the side of the cell. If the shape of the surface in a free surface cell is approximated to a straight line, the orientation and location of the free surface line in the cell can be defined by the above three parameters. Figure 4.21 illustrates the free surface in a free surface cell. Shear stresses in the ^-direction are given by Equations 4.15 and 4.16. The classical free surface condition states that there are constant pressure and zero shear stresses at the free surface. The pressure term appears only in the r- and z-momentum equations. The constant pressure condition at the free surface is imposed in the 2D solver [27]. The zero shear stress condition implies that the shear stresses are only transferred within the fluid and not at the free surface. Therefore, shear forces applied to each side of the cell are dependent on how much fluid is in contact with each side. Here, four new parameters are defined so as to describe the ratio of fluid contacting each side. The four side ratios Chapter 4. Computational Fluid Dynamics Code 44 w Ie AZ Ar Figure 4.21: A free surface cell. are denned as Az Az Ar Ar where I represents the length of the contacting fluid on the corresponding side of the cell. For example, fe = 0.85, fw = 0.43, /„ = 0 and /„ = 1 in Figure 4.21. The method used for calculating each of the four parameters is presented in Appendix E. With these four side ratios, the -^momentum equation may be reformulated for those cells containing a free surface. It should be noted that in a free surface cell the resident mass of fluid is decreased as a result of the fact that the cell is only partially occupied by the fluid. If the fluid density is p and the air density is negligible, then an equivalent density for the free surface cell is pF(i,j), where F(i,j) is the fractional volume of fluid of the cell as described earlier. Chapter 4. Computational Fluid Dynamics Code 45 Therefore, for a free surface cell, the RHS of Equation 4.17 can be modified to reflect the fact that the cell is partially empty. RHS of Eq. 4.17 = —L- (fe [ryk~fw}r'Tr6L + M%zAi!fi) . (4.27) PF(i,j) \ [r(z)]2Ar(z) Az(j) ) K ) For a fluid cell F(i,j) — 1, the four side ratios are all equal to unity and the above equation becomes identical to that which is used for the fully occupied cell. Therefore, the formulation for the free surface cell is also valid for a full fluid cell. 4.2.4 Wall boundary conditions For a specific problem, the wall boundary conditions need to be specified to solve the -^momentum equation. Shear stresses in the ^-direction at the wall boundary due to viscous fluid flow depend on the orientation of the wall. This is because of the non-interchangeable feature between the radial and axial coordinates in the cylindrical coordi-nate system. Wall shear stresses for a rectangular domain in the r-z plane are illustrated as follows. In 3D r-z-9 coordinates, the rectangular domain represents a rectangular ring or a cylinder because of the axisymmetric assumption. For a rigid no-slip wall, the wall shear stress is equal and opposite to the fluid shear stress at the wall. When the wall is parallel to the axis (the axial coordinate z) at r = rw, the wall shear stress in general terms is d / w I wall = -f1 •w-(~) orn \ r J (4.28) r„=0 where rn is measured from the wall and normal to the wall. When the wall is facing the same direction as r (i.e. Wall (W) in Figure 4.22), then rn = r — rw and d (w rw = -p dr \r J (4.29) r=rw Chapter 4. Computational Fluid Dynamics Code 46 If the wall is facing the opposite direction to r (i.e. Wall (E) in Figure 4.22), then rn — —(r — rw) and the minus sign in the RHS of Equation 4.29 should be replaced with a positive sign. When the wall is perpendicular to the axis (the axial coordinate z) at z = zw, the wall shear stress in general terms is T~w = - Tzolwall = -r1 where z n is measured from the wall and normal to the wall. When z n is in the same direction as z (i.e. Wall (S) in Figure 4.22), then zn = z — zw and (4.31) If z n is in the opposite direction to z (i.e. Wall (N) in Figure 4.22), then zn = — (z — zw) and the minus sign in the RHS of Equation 4.31 should be replaced with a "+" sign. For a rigid free-slip wall, rw = 0. In the finite difference formulation, the meshed fluid domain is enclosed by fictitious cells that enable the wall boundary conditions to be specified. The computational domain is shown in Figure 4.22. Orientations of the four walls are represented by East, West, North and South. For a rigid no-slip wall, rw ^ 0, but w at the wall is equal to 0. In the CFD program, shear stress at the wall applied from the fluid in the fictitious cell is specified to be equal, but in the opposite direction, to that applied from the fluid in its neighboring fluid cell on the opposite side of the wall, since the size of the fictitious cell is assigned to be the same as its neighboring cell in the fluid domain. Eventually, the wall conditions are fulfilled by specifying the velocity at the fictitious cell according to the velocity at the neighboring cell. For a rigid free-slip wall, TW = 0. In the CFD program, the angular velocity (w/r) in dw z _ = 0 (4.30) Tw = -p dw dz Chapter 4. Computational Fluid Dynamics Code 47 Jo Wall (W) 2 1 1 2 Fictitious Cells Wall (N) y Wall (E) Wall (S) lO Figure 4.22: Meshed computational domain and wall boundaries. the fictitious cell is assumed to be the same as that in the neighboring fluid cell on the opposite side of the wall. Boundary conditions applied at the wall boundaries, used in the CFD program, are summarized in Table 4.1. For an explanation of the symbols used in the table, refer to Figure 4.22. Other wall boundary conditions such as continuative boundary (zero pressure gradient across the imaginary wall boundary) and pressure boundary (specifying the pressure in the fictitious cells) can be imposed in a similar manner. 4.2.5 Obstacle boundaries Besides the wall boundaries, solid bodies inside the wall boundaries can be set as obstacle boundaries. In the CFD code, those cells that defines an internal obstacle are blocked out. Velocity components, u and v, on the four sides of an obstacle cell are normally set to zero. However, in the program, an obstacle cell is allowed to have a Chapter 4. Computational Fluid Dynamics Code 48 Free-slip No-slip Wall (W) w(l,j) _ w{2J) r(l) - r(2) w[l,j) _ w(2,j) r(l) " r(2) Wall (E) w(i0J) _ w(i0-l,j) r{i0) r(i0-l) Miod) _ w(* 0-i,i) r(io) r(i0-l) Wall (N) w(hJo) = w(i,j0-l) w(i,jo) = -w(i,j0-l) Wall (S) w(i, 1) = w(i, 2) w(i, 1) = — w(i, 2) Table 4.1: Summary of wall boundary conditions used in the CFD code. finite tangential velocity. Therefore, in the A3D-CFD code, not only stationary, but also rotating solid bodies, can be modeled. 4.2.6 Algorithm c Compared to the 2D governing equations, the 3D r-momentum equation (Eq. 4.10) has an additional term (—w2/r) on the LHS of the equation. This additional term accounts for the centrifugal acceleration of the swirling flow, acting as a body force. The z-momentum equation is unchanged and the ^ -momentum equation is added to the algorithm. In the ^ -momentum equation, velocity components u and v are present in the convective terms. In the r-momentum equation, term (—w2/r) requires a solution for w. Thus, the three velocity components u, v and w are coupled and a simultaneous solution is needed to solve the full 3D velocity field. Chapter 4. Computational Fluid Dynamics Code 49 From the comparison of the governing equations, one can see that the 2D code SOLA-V O F can be used to form part of the 3D axisymmetric flow solver. In general, the 2D code can solve for u and v, provided that the tangential velocity, w, is specified either from the initial setup or from a previous iteration. Computed values of u and v can then be put back into the ^-momentum equation and, by solving the ^-momentum equation, a new value of w may be obtained. This new value of w is then put back into the 2D solver and the iteration procedure is repeated until a convergence in w is obtained. The flow chart showing the computational routine for solving a 3D axisymmetric flow is presented in Figure 4.23. The solid blocks represent the original 2D S O L A - V O F code, and the dash-dot blocks represent the new portions of the program. When the program starts, it reads the input data provided by the user. Then subrou-tine M E S H S E T generates computational mesh from the input data. Subroutine S E T U P initializes constants and initial velocity components u and v at time t — 0. In modeling of the present problem, the initial velocity components are all set as zero. Subroutine B C is then called, which sets the wall boundary conditions and free surface conditions. Subroutine V F C O N F computes the value F (fractional volume of fluid) for each cell. From the initial values of u and v, subroutine W T H E T A sets initial wall boundary con-ditions for w and computes w. Then, subroutine P E T A C A L determines the slope of a free surface and the side of a free surface cell on which the fluid is primarily sitting on. Flow field at time t = 0 is obtained now. Subroutine O U T P U T can be called to export a data file of the computed results. From the data file, the flow field can be viewed using software T E C P L O T . Subroutine S T R E A M computes the secondary flow stream function from the computed secondary velocity components u and v. To advance time step, subroutine D E L T A D J computes the allowable time step At and also adjusts the time step according to number of iterations. The program then checks i f the preset time limit is reached. If not, the time is advanced by one step. Chapter 4. Computational Fluid Dynamics Code 50 START Read Input Data MESHSET SETUP t = cycle = 0 BC Rotating B.C. Wb.c=^(t) TILDE H I BC I PRESSIT Have Pressures Converged ? I YES g r =w 2 / r j NO VFCONV WTHETA ! BC I YES Have w(i,j) Converged ? 1 PETACAL NO OUTPUT i .1 DELTADJ |——j STREAM ^ Time to Stop ?) YES t = t + At cycle = cycle +1 — * ^ S T O P ^ Figure 4.23: Computational flow chart of the modified A3D-CFD code. Chapter 4. Computational Fluid Dynamics Code 51 It should be noted that the rotating boundary condition is imposed by assigning values of w for the rotating boundary cells. The values of w can be set as a function of time t. In modeling of the present problem, the values of w are assigned as such that the angular speed of the rotating boundary is set as zero at t = 0 and is gradually increased with time to a constant speed set by the user. As the time advances, a steady flow field can be obtained for a constant speed of the rotating boundary. At the new time level, variables u, v, w and p are to be determined from the initial conditions or from the previous time level. Using the w values from the previous time level, subroutine TILDE computes an explicit solution for u and v for the new time level. Then, subroutine PRESSIT iterates u, v and p such that mass is conserved in each cell of the mesh. The program checks if pressures are convergent or mass is conserved. Iteration is performed until the convergence is achieved. However, these values of u and v are computed for the values of w from the previous time level. Using these values of u and v, subroutine WTHETA computes w again. This value of w, at the first iteration, is likely to be quite different from the value w for the previous time level, which is considered that the solution for w is not convergent. Then, the newly computed values of w are put into the r-momentum equation and the r- and z-momentum equations are solved again. From the new values of u and v, subroutine WTHETA computes w. The iterative procedure continues until a solution for w is convergent. Then, the velocity field (u, v, w) at the new time level is obtained. Details of the iterative solution procedure for the 2D code are given by Nichols et al. [27]. The original 2D code was written in subroutine format. Most of these subroutines are kept intact or modified only slightly in the A3D-CFD code. Each subroutine shown in Figure 4.23 is listed in Appendix D, with a brief description of its major function. It should be noted that only the major subroutines of the program are listed in Appendix D. In evaluating the new value of w for the next iteration, a relaxation factor a is Chapter 4. Computational Fluid Dynamics Code 52 introduced. The new value w for the next iteration is given by w = w + a(w — w) (4.32) where w is the value of w obtained from the current iteration, while w is the initial value of the current iteration. Hence, the contents of the parentheses represent the change in iy as a result of the current iteration. In general, when the change in w produced by the current iteration is small enough, then w is considered to be convergent. Underrelaxation (a < 1) is a technique often used to avoid divergence in the iterative solution of strongly nonlinear equation [28]. In particular, the convergence criterion for the solution to w in the code is defined such that if E & . X & W , ; ) ! ( 4 ' 3 3 ) the solution is considered to be convergent, io and jo are the numbers of the discretized cells in the r and z directions, respectively. The LHS of the equation is the ratio of the average difference of the velocities between the current and previous iterations to the average velocity obtained from the current iteration. When the ratio is less than a user defined small number (<J), the solution is considered to be convergent. Apart from these main features of the A3D-CFD code, some other minor features were also added to the code for the purpose of this study. These features are discussed in the following sections. 4.2.7 Shape of the free surface This part of work was carried out solely for the purpose of plotting the shape of the free surface using a software package called TECPLOT. As described earlier, the orientation of a free surface in a free surface cell is represented by the following three parameters: F(i,j), NF(ij) and TANTH(i j) . As long as these three parameters are known Chapter 4. Computational Fluid Dynamics Code 53 for a given cell the free surface, represented by a line in that cell, can be located. The general shape of the free surface is represented by all the lines of the free surface cells. The idea of plotting the free surface using TECPLOT is that the free surface is plotted by smoothly connecting the mid-points of those lines of the free surface cells. The procedure is detailed as follows. First, the coordinates (R, Z) of the mid-point of the line is calculated for each cell containing the free surface. Second, a constant value of 0.5 (in the subroutine output) is assigned to each mid-point (R, Z). Third, another constant value of 0.0 (in the subroutine output) is assigned to cells not containing a free surface. Fourth, using TECPLOT, a contour line of value 0.5 connects all the mid-points, which represents the shape of the free surface. The method of representing the free surface is illustrated in Figure 4.24, in which the mid-points are represented as dots. Presentation of shape of the free surface is also included in subroutine output. Free Surface Z -( U ) o .» R " r Figure 4.24: Free surface representation. Chapter 4. Computational Fluid Dynamics Code 54 4.2.8 Secondary flow stream function With the A3D-CFD code providing solutions to the three momentum equations, the secondary flow stream function ip(r, z), which is also of interest for this project, may also be determined. In cylindrical coordinates, the stream function is defined as: 16V 10^ u = v = - — . r oz r or (4.34) In a finite difference cell, the stream function is calculated at the top-right corner of the rectangular cell as shown in Figure 4.25. k i r < ( i . j ) o -\|/(i-1J-1) 4 _u(ij) . W - 1 ) 1 I ' v(i,j-1) Figure 4.25: Location of if) in a finite difference cell. Using the finite difference method, Equation 4.34 may be approximated as: [r(t) + 0.5Ar(t)] Az(j')' r(i) Ar(i) (4.35) or Ar(i)' Az(j) + v(i,j — l) r(i) Ar(t). (4.36) Chapter 4. Computational Fluid Dynamics Code 55 Boundary conditions for the stream function ^ may be denned such that ip = 0 at all boundary walls and at the free surface. In a finite difference cell, the free surface boundary condition is applied by setting tp(i,j) = 0 if cell is empty (no fluid). This part of the work is included in subroutine stream. 4.2.9 Secondary flow vorticity As with the stream function, the secondary vorticity function is also of interest. Secondary vorticity (in the 9 direction) is defined as: du dv ,M „„. - = &-*• ( 4- 3 7 ) In a finite difference scheme, the vorticity function u>«(z,j) is also calculated at the top-right corner of a rectangular cell The finite difference approximation of the vorticity function is given by .v + ! ) - « ( « , j) _ v{i + l,j)-v(i,j) . . e { h J ) 0.5[Az(j) + Az(j + 1)] 0.5[Ar(») + A r ( t + 1)]' For non-fluid cells the boundary condition for <JJ$ is set to zero. The computation of the secondary vorticity function can be found in subroutine vorticity. 4.3 Verification of the A 3 D - C F D code To determine the rehabihty of the A3D-CFD code, computational results from the code are compared with some existing analytical solutions and experimental results. Code verification includes the following cases. 4.3.1 Flow near a disk rotating in a fluid at rest In this section, three velocity components of a homogeneous fluid are examined. The sample case is the flow near a disk of infinite diameter, rotating about an axis Chapter 4. Computational Fluid Dynamics Code 56 perpendicular to its plane with a uniform angular velocity (u>) in a fluid otherwise at rest. The analytical solution to the problem is given in [34]. Velocity distribution is represented in terms of the following assumed functions: F { ( ) = J L ; H(0 = -L.; 0(0 = -ru> y/v(jj ruj where Due to the axisymmetric nature of the flow, only half of the disk is modeled, with the axis assigned as the left, wall boundary of the domain. Solid free-slip condition is applied at this wall boundary. In the computation, the disk has to be of a finite diameter, and in order for the computational solution to be comparable to the analytical one, the ratio of the radius of the disk to the thickness of the layer of fluid carried with the disk must be large. Consequently, the dimensions of the computational domain are 0.25 m x 0.025 m in the radial and axial directions, respectively. The disk is located at the bottom of the computational domain and the top and right-hand walls are continuative boundaries. The continuative boundaries allow the fluid flow in and out the boundaries freely. Computed velocity distributions close to the center of the disk are shown in Figure 4.26. As stated earlier, the analytical solution assumes a disk of infinite diameter. In the computational study, it is appropriate, therefore, to examine the flow close to the axis where the flow is less affected by the finite treatment of the right-hand boundary. The velocity distribution at r = 0.016m in terms of F((),H(() and G(() are compared with the analytical solution as shown in Figure 4.27. Agreement between the computational and analytical results is reasonable, except that the magnitude of the velocity components from the computation are somewhat smaller than those from the analytical solution. Discrepancies between these two results can probably be attributed to the fact that the disk diameter had to be treated as finite Chapter 4. Computational Fluid Dynamics Code 57 0.030 0.020 Z(m) 0.010 h 0.000 h , free-slip boundary continuative boundary • i i i i i i i i i i i i i i w (m/s) 0.000 0.010 0.020 0.030 r(m) 0.040 Figure 4.26: Velocity distributions close to the center of the disk. Figure 4.27: Computational and analytical velocity distribution of a rotating disk. Chapter 4. Computational Fluid Dynamics Code 58 in the computation, which would result in reduced centrifugal pumping effect close to the right-hand wall. 4.3.2 Taylor vortices The Taylor vortex is a well-known phenomena of flow instability between two concen-tric rotating cylinders [36]. To explore the capability of the A3D-CFD code to predict such a vortex flow, calculations have been performed for a case when the outer cylinder is rotating at 1.5 times the speed of the inner cylinder, but in an opposite direction (1^2/^1 = —1.5). The computational domain is chosen such that the radii of the inner and outer cylinders are 7*1 = 38.0 mm and r2 = 40.35 mm, respectively, and the vertical (axial) height is 30 mm. The ratio of the axial height to the radial clearance (r2 — ri) is 12.8. The wall boundaries are specified such that the left and right-hand walls are rigid no-slip and the top and bottom walls are rigid free-slip. It should be noted that the axial length of the annulus is relatively short in comparison to the case considered by Taylor [36]. However, the free-slip condition for the top and bottom walls ensures that zero shear stress is imposed on the fluid from these walls (earlier experimenters constructed long cylinders to minimize the shear stress effect from the ends of cylinder). In the computation, one effect from the short length of the annulus is that the solid free slip wall condition does not permit axial flow across the walls. Computed Taylor vortices for a case when uii/v = 5.5 x 107 (rad/m2) are shown in Figure 4.28. In figure (a), the dashed lines represent the stream functions with negative values, where the direction of the streamline is counter-clockwise, while the solid lines represent the opposite to the dashed lines. In figure (c), the dashed lines represent the tangential velocities with negative values, which are in the opposite direction to the velocities of the solid lines. It should be noted that only about a quarter of the axial (z) Chapter 4. Computational Fluid Dynamics Code 59 computational domain is presented in this figure. z I J (a) (b) (c) Figure 4.28: Computed Taylor vortices: (a) stream function; (b) secondary velocity vectors; (c) tangential velocity contours. Two columns of vortex arrays with an inner column of strong vortices and outer column of weaker vortices (as Taylor predicted) are shown in Figure 4.28(a). Compared to Taylor's analysis, the shape of stream functions are somewhat "distorted" and the cores of vortices in the inner column are not equally spaced axially, as shown. The tangential velocity is significantly distorted after the onset of the instability as shown in Figure 4.28(c). The stable state solution to the tangential velocity given by Eq.1.1 is not valid after the onset of the flow instability. It should be noted that in the Taylor's analysis [36], the convective terms involving the disturbed velocity components were neglected for the sake of simplicity, and may lead to a solution which does not include some of convective interactions from the three Chapter 4. Computational Fluid Dynamics Code 60 velocity components. However, the computational results have clearly shown that there are convective interactions among the velocity components. Such interactions were also indicated in reference [41]. As mentioned earlier, the computation shows the core of the vortices as being shifted (in comparison to the equally spaced vortices of the Taylor's results). However, the computed vortices are axially spaced regularly in terms of a counter-rotating vortex pair and the axial distance between two neighboring vortex pairs is a constant. In order to compare to the Taylor's results, the axial distance between two neighboring vortex pairs is compared to twice the axial core distance of neighboring vortices from the Taylor's results. In this case, the computed axial distance of a neighboring vortex pair is about 3.0 mm, which is close to twice the spacing (2 x 1.35 mm) obtained by Taylor. This discrepancy may be attributed to the limited axial dimension of the computational domain. It should be noted that in the computation, the Taylor vortices appear as pairs along the axial direction. Regardless of the axial dimension of the domain, the total number of the computed vortices will be an even number along the axial direction. In the computation (not shown), swirling motion is dominant with minimal random secondary flow at low rotational speeds. As the speed increases, and eventually exceeds a certain value (LOI/U = 4.0 x 107 rad/m 2 ), vortices start forming at both the top and bottom ends of the annulus and develop toward the middle section of the annulus. These vortices are very weak and tangential velocity maintains the profile of the stable state solution (independent of z). However, as the speed further increases and the value of OJX/V reaches 4.5 x 107 (rad/m 2), the vortices become stronger and the tangential velocity profile is noticeably deformed from its stable state solution. Prom the computation, the instability occurs at about — = 4.0 ~ 4.5 x 107 rad/m 2 . v Chapter 4. Computational Fluid Dynamics Code 61 From interpretation of Taylor's experimental observation, the instability first appears at — = 4.3 x 107 rad/m2. v Therefore, the computation using the A3D-CFD code predicted the critical speed of the instability reasonably well in this case, compared with Taylor's observation. 4.3.3 Shape of free surface in a rotat ing annulus In this section, free surface shapes are examined. As a sample case, consider a vertical annulus, partially filled with fluid, which is formed by two concentric cylinders with the inner one stationary (u>i = 0) and the outer one rotating at a speed of u>2. Experiments were conducted on such an annulus, with the objective of observing the shape of the free surface with different rotational speeds of the outer cylinder. A schematic of the experimental test facility is shown in Figure 4.29. Density and viscosity of the working fluid are p = 876 kg/m3 and p = 6.7 x 10 - 2 Pa-s. Computed and measured free surface shapes, for speeds of u>2 = 8.37 rad/s and u>2 — 16.75 rad/s, axe shown in Figure 4.30. Hollow circles indicate actual measured free surface locations. With regard to the shape of the free surface, correlation between the experimental and computational results were excellent. Based on the analysis given in Appendix F, a dimensionless free surface shape may be determined after introducing the following non-dimensional parameters (refer to Appendix F): AH AH* = Ah* = AR' Hr) AH' 9n- g , . r-Rj. Chapter 4. Computational Fluid Dynamics Code 62 Support Frame \ ^ — Rotating cylinder Figure 4.29: Schematic of test facility for free surface shape measurement. Non-dimensionalized height differences of the free surface are shown in Figure 4.31. The non-dimensional axis gR represents the ratio of the centrifugal acceleration at the inner surface of the rotating cylinder (r = R2) to the gravitational acceleration. Computed free surface height differences AH* are in a good agreement with the experimental results for different non-dimensional radial acceleration gR. Furthermore, the height difference AH* is, approximately proportional to, but slightly curved against the non-dimensional acceleration gR. More precisely, (d{AH*Y { d{gR) t slightly increases with gR. Chapter 4. Computational Fluid Dynamics Code 63 0 . 1 5 r z(m) 0 . 1 0 0 . 0 5 h 0 . 0 0 r-0 . 0 0 0 . 0 2 0 . 0 4 0 . 0 6 0 . 0 8 r(m) (a) co, = 8.37 rad/s (b)oo2= 16.75 rad/s Figure 4.30: Computed and measured free surface shapes. 2 . 5 rr < x < -y y' • y y y' y' Experiment Computat ion Isobar -Solid body V " y y y y y y -y . y y' y' y y' - s y' s' /' y' s' y ^ y A R / R 2 = 0.508 v 2 / g R 2 = 5.61 xlO" 6 s y' y ^ 0 . 0 0 . 5 1.0 1.5 2 . 0 2 . 5 3 . 0 3.5 C02R 2/g Figure 4.31: Computed and measured height differences of free surface. Chapter 4. Computational Fluid Dynamics Code 64 It might be interesting to examine the shape of the isobaric line when the annulus is infinitely long and filled with fluid. The annulus radius ratio <f> = R2/R1 = 2.034. From Eq. G.70 in Appendix G, K = 0.4364. Therefore, the height difference of an isobar across the annular clearance is given by ^ | = 0 .4364^. (4.39) This is represented by the dashed line in Figure 4.31. When the gravitational force is dominant (gR is very small), AH is also very small. In such a case, the free surface level change can be approximated by the above equation. However, as gR increases, the height difference AH of the free surface increases at a greater rate than that of the isobaric line. Therefore, as gR increases, the free surface condition tends to increase the height difference AH* in comparison to this isobaric line. At a large value of gR, the centrifugal acceleration is dominant. As gR increases, the oil approaches to the solid body rotation with the rotating cylinder (assume that the cylinder is infinitely long as to permit the large deformation of the free surface). For the limit of the solid body rotation, the angular velocity of the oil a> = o>2, and >l rR2 wl AH =— [ 2rdr = ^(Rl-Rl). g JRi 2gK 2 1J Therefore, the non-dimensional height difference AH* for solid body rotation may be written as: * § - I ( l + 1 ) ^ = 0 .7458^ . (4.40) A i l 2 V <t>) 9 9 V ' This equation is shown in Fig. 4.31 as dash-dot line. Non-dimensional free surface shapes, from both the experiment and computational sources, are shown in Figure 4.32. The shape of free surface h(r) is non-dimensionalized by the height difference of free surface from the experiment, which also makes the com-putational profile comparable to the experimental profile. Chapter 4. Computational Fluid Dynamics Code 65 x < 0.0 0.2 0.4 0.6 0.8 (r-R,)/AR Figure 4.32: Non-dimensional shapes of free surface. It is shown that the free surface shape is well predicted by the computational model. From the experimental results, it can been seen that the non-dimensional shape Ah* is affected by gR. The computational results predict the same trend of the effect as the experimental ones, even though there are some discrepancies in magnitude. In general, the shape of the free surface in the experiment is well predicted by the A3D-CFD code. Now, let's examine the shape of an isobaric line. Eq. G.67 gives the shape of the isobaric line, in which z — z\ = h(r). Eq. G.67 divided by Eq. 4.39 yields h AH 0.4364fl2A# A 2 2 ( r 2 - J R 2 ) + 2ABln(^- ) This equation gives the isobaric line and is shown in Figure 4.32 as indicated. At a small value of gR, the free surface profile approaches to the isobaric line. At a large value of gR, the free surface profile approaches to the solid body rotation Chapter 4. Computational Fluid Dynamics Code 66 (at a speed of u}2) as described earlier. The free surface shape for solid body rotation, can be expressed as h u,2(r2-R2) u>2R2 [(r/fiQ 2 - 1] ( 4 4 1 ) (r/R,)2 - 1 <f>(<f> - 1) AR 2ARg 2g Equation 4.41 divided by Equation 4.40 yields h _ r2-R\ _ (r/R^2 - 1 AH ~ R\ - R\ ~ <t>2-\ ' This equation gives the free surface shape for the solid body rotation and is shown in Figure 4.32 as indicated. Chapter 5 Computational Modeling of Flow One of the prominent features of a vertical self-contained bearing assembly is its axisymmetric geometry. Accordingly, the A3D-CFD code has been developed for the purpose of studying, numerically, the oil flow in a vertical bearing assembly, with par-ticular reference to the oilwell tube annulus. The numerical aspect of the present study consists of the following two stages. • modeling of the oil flow in the simulated bearing assembly, to obtain a general idea of the flow patterns in the self-contained bearing assembly; • modeling of the oil flow in the oilwell tube annulus, to obtain details of the oil flow in the annulus and its boundary region. 5.1 Modeling of flow in the simulated bearing assembly The simulated vertical bearing assembly test facility, shown in Fig. 1.2, was built in 1988 for the purpose of making a visual observation of the oil flow in the oilwell tube annulus. Generally, the test rig enabled the bulk swirling motion of the oil in the oilwell tube to be observed. Also, an extensive experimental parametric study was conducted and general physical processes of fluid in the oilwell tube annulus were observed [5]. Because of limited visual access, flow visualization was limited to qualitative obser-vation of the bulk swirling flow of oil in the oilwell tube annulus. However, details of 67 Chapter 5. Computational Modeling of Flow 68 oil flow in the bearing assembly are also of interest from both qualitative and quantita-tive points of view. The purpose of the numerical modeling is therefore useful from two standpoints: (1) to obtain a better understanding of the general oil flow patterns in the bearing assembly, and (2) to compare these computational results with the experimental observations. The test rig is designed so that oil returns from the oil reservoir to the oilwell tube annulus via sixteen equally spaced 12.5 mm diameter holes. In the CFD modeling, these sixteen holes are modeled as a channel such that the wetted area of the slot at the oilwell tube side is the same as the total wetted area of the sixteen holes. From this assumption, the height of the slot is determined. It should be noted that, in the model, a match of the wetted area may not produce a head loss in the channel which is the same as the head loss of the sixteen radial holes. Furthermore, flow patterns in the model, especially close to the channel, would be different from the experiment. However, this simplification is still expected to serve the main purpose of this section of the study, which is to qualitatively understand the flow patterns in the bearing assembly. A perforated baffle fitted in the reservoir of the experimental set-up is modeled as a slotted cylinder, which is concentric to the axis of the bearing assembly. The cooling coil in the experimental set-up is not included in the model. Since the oil velocity is very small in the region around the cooling coil, the effect of this simplification on the oil flow in other regions should be minimal. A no-slip boundary condition is assumed for all wall boundaries. The solid bodies in the domain are modeled by setting the cells covering the bodies as obstacle cells. In the code, velocity components u and v on faces of obstacle cells are set to zero. The value assigned to the velocity component w of the obstacle cells is dependent on the nature of the solid body. That is, for stationary obstacles, w is set to zero and for rotating obstacle cells, i.e. the runner, w is set to a value which is equal to the tangential velocity of the rotating body. This rotating boundary condition can be given as a function of time t, Chapter 5. Computational Modeling of Flow 69 in accordance with real start-up situation. The gradual increase of rotating boundary velocity is also advantageous in terms of obtaining a convergent solution compared to a one-step-increase of the velocity. The numerical model of the bearing assembly and the associated computational mesh are presented in Figure 5.33. At the walls, a no-slip boundary condition is assumed. 0.12 h N Runner 0.08 Oilwell tube 0.04 0.00 0.08 Perforated baffle Oil level Oil reservoir Pseudo-guide bearing Oil recirculation channel 0.14 0.10 0.12 r(m) (a) (b) Figure 5.33: Modeling of the simulated bearing assembly. As the cell size —> 0, the finite difference formulation of a partial differential equation Chapter 5. Computational Modeling of Flow 70 (i.e. Navier-Stokes equation) approaches to the exact solution of the PDE. Therefore, small cell size is desirable to minimize the discretization error. However, the total number of cells of the mesh is usually limited by the memory and speed of a computer. As a compromise, a mesh refinement technique is often used in a viscous region such as a boundary layer to obtain an adequate description of the flow [29]. The SOLA-VOF code has an automatic mesh generator [27]. The idea is to build up a mesh by stacking together a series of submeshes. For a variable submesh, the mesh generator expands the cells starting with the minimum cell size (specified by the user) and increases the cell sizes quadratically based on the domain and the number of the cells (specified by the user) in the submesh. The current model of the simulated guide bearing assembly employs a variable mesh with a 50 x 94 (radial x axial) grid which is shown in Figure 5.33 (b). In order to model the small flow geometries of the bearing assembly, i.e. the radial clearance space of the pseudo-guide bearing and the oil recirculation channel, finer grids are employed in these regions. Also the secondary flow underneath the runner region is expected to be strong. Therefore the grids in this region are refined. The perforated baffle is coarsely modeled as uniformly distributed rings. 5.2 Modeling of flow in the simulated oilwell tube annulus Numerical modeling of the oil flow in the oilwell tube annulus is the principal area of focus of this study. The modeling is based on the simulated oilwell tube annulus test rig shown in Fig. 2.5. It should be noted that the vanes present in the bottom part of the rig are designed to minimize the effect of the rotating disk (directly connected to the shaft) on the flow in the test head. In the experiments, it was observed the oil flow in this region was effectively confined with minimal influence on the flow in the test head. To simplify numerical modeling, only the test head (upper part of the rig) is modeled. A schematic Chapter 5. Computational Modeling of Flow 71 of the numerically modeled oilwell tube annulus and its associated computational mesh are shown in Figure 5.34. The fluid domain is composed of a 50 x 74 variable mesh as <g 0.12 N Runner Oil level 0.08 Oilwell tube 0.04 0.00 h I , Pseudo-guide -V bearing 0.12 0.14 0.16 r(m) (a) (b) Figure 5.34: Modeling of the simulated oilwell tube annulus. shown in Figure 5.34(b). The mesh grids are refined in the region underneath the runner and the pseudo-guide bearing clearance space. A no-slip boundary condition is applied to the four wall boundaries. It should be noted that to simulate the experimental set-up Chapter 5. Computational Modeling of Flow 72 the fax right-hand wall below the bottom of the runner is modeled as a rotating boundary (refer to the computed results in Fig. 6.38(b)). Grid dependency test showed that grid-independent results were obtained for the present computational mesh. The test was conducted at earlier stage and could not be reported due to the loss of files caused by a computer damage. Chapter 6 Computational Results and Analysis 6.1 Test oil properties The fluid used in the experiment is a commonly used lubricant oil (ISO VG46). Density of this oil is p = 876 kg/m3 and its kinematic viscosity is measured with a Saybolt viscometer. To measure the viscosity of the oil, the test oil is poured into a cylindrical container, heated to a pre-assigned temperature and kept constant at that temperature in an oil bath. To measure the viscosity in Saybolt Universal Seconds (SUS), the time taken for 60 nulhlitre of the oil to flow through a capillary tube (located at the base of the cylindrical tube) into a separate receptacle is recorded. The time in seconds is then converted to a kinematic viscosity in centistokes using a conversion graph. The viscosity of the oil measured at various temperatures is shown in Table 6.2. A plot Temperature (°C) 19 25 31 40 50 SUS (s) 448.3 360.4 265.0 178.3 119.6 Viscosity (cSt*) 104.6 84.0 61.8 41.5 27.9 t 1 cSt = 1 x l O - 6 m2/s Table 6.2: Measured kinematic viscosity coefficients of the test oil. showing measured viscosity against temperature is shown in Figure 6.35. Viscosities at other temperatures are determined from the interpolated line, as shown. 73 Chapter 6. Computational Results and Analysis 74 Temperature (°C) Figure 6.35: Kinematic viscosity of the test oil versus temperature. 6.2 Flow in the simulated vertical bearing assembly A computational study was performed using the numerical model of the simulated vertical bearing assembly as shown in Figure 5.33. A sample of typical computed velocity distributions is shown in Figure 6.36. As the runner rotates, oil primarily swirls about the axis of the rotation as shown in Figure 6.36(b). As in a moving boundary layer, the tangential velocity of the oil adjacent to the moving surface is equal to the speed of the runner, but this velocity gradually decreases as the distance from the moving surface increases. At the oilwell tube wall and at other stationary surfaces, the tangential velocity of the oil is zero due to the no-slip boundary condition. In the earlier experimental studies [5], a perforated baffle was used to reduce the angular momentum of the oil as it swirled in the oil reservoir. Chapter 6. Computational Results and Analysis 75 The computational results shown in Figure 6.36(b) indicate that the tangential velocity of the oil is effectively confined to that part of the reservoir enclosed by the perforated baffle. Close to the right wall of the reservoir, the angular momentum of the oil is very small. E N 0.08 0.10 0.12 r(m) (a) (b) Figure 6.36: Velocity distributions in the simulated bearing assembly (WJJ = 15.7 rad/s, v = 1.82 x 10 - 4 m2/s): (a) secondary velocity vectors; (b) primary velocity contours. As a result of this swirling motion, the free surfaces in both the reservoir and in the Chapter 6. Computational Results and Analysis 76 oilwell tube deform from their flat still level and take on a distinctive profile. The shape of the free surface is mainly determined both by the centrifugal acceleration at the free surface and the gravitational acceleration. For a steady axisymmetric flow, the shape of a free surface in the r-z plane should satisfy the classical free surface condition (zero shear stress at the free surface in the r-z plane). In other words, at the free surface, the tangent of the free surface should be normal to the direction of the sum of the body forces (centrifugal and gravitational forces in this case) if flow is not moving in the direction of the body forces. For example, the slope of free surface at a free surface point (rj,Zf) may be expressed as dz dr = w2(rf,zf) for purely rotating flow. The shape of free surface may be modified by secondary flow since the secondary flow may bring in momentum in the body force direction at the free surface. Apart from the deformed free surface, there are secondary flows (in r-z plane) which are also associated with the primary flow, as shown in Figure 6.36(a). For such a bearing assembly, an important aspect of the secondary flow is the flow recirculation and in the region underneath the bottom of the runner, a positive radial pressure gradient is established. Figure 6.37 shows computed static pressure in this region, corresponding to the same condition shown in Figure 6.36. The pressurized oil in this region squeezes into the reservoir through the bearing clearance space and recirculates back into the oilwell tube annulus through the oil return channel. This recirculation of oil prevents the build up of heat between the bearing surfaces. The heat taken away by the recirculating flow is transferred to a cooling system, usually located in the reservoir. The temperature dependent viscosity of the oil is a key factor in effective fluid film lubrication of the bearing. Effective flow recirculation and an efficient cooling system are very important Chapter 6. Computational Results and Analysis 77 Units: Pa Figure 6.37: Static pressure distribution beneath the runner. factors for such a bearing design. In addition to the general secondary flow recirculation, some local secondary vortex flows are formed. From the computational results, it can be seen that a strong vortex flow is formed just beneath the bottom of the runner. This is similar to the well known rotating flow formed in a cylinder with the top and side walls rotating and with a stationary bottom surface [4, 11]. As a result of the rotation of the bottom surface of the runner, oil swirls and a Couette type flow is established in the circumferential direction. Oil particles with high circumferential velocity close to the bottom surface of the runner experience a high centrifugal force, which drives the oil outwards in the radial direction. At the same time, a pressure gradient is established in the radial direction as a result of the restriction at the outside wall of the cavity. This pressure gradient drives oil away from the bottom surface where the centrifugal force is low, causing it to move radially inward. In the annulus between the outer surface of the runner and the perforated baffle, a Chapter 6. Computational Results and Analysis 78 secondary vortex pair is formed. The lower vortex, just above the top surface of the pseudo-guide bearing, is caused by a similar reason described above. The particles, close to the top surface of the stationary bearing body, have relatively low tangential velocity (low centrifugal acceleration) and tend to move radially inward driven by the positive radial pressure gradient of the oil. Therefore, one can conjecture that the presence of the stationary boundary (bearing body) at the bottom of this annulus generates such a secondary vortex in the annulus in the presence of the swirling flow. The upper counter-rotating vortex could be a combined result of the lower vortex and the presence of the deformed free surface. It is worth mentioning that in such a configuration, Taylor type vortices may set in at certain speed of the runner even without the "endwall" effect from the prescribed stationary bearing body. From the critical Taylor number given in Equation 1.3, the critical runner speed for the Taylor vortices to set in is |critical = 20.3 rad/s. The runner speed used for the computation is U>R = 15.7 rad/s, which is close to this critical value. Aeration of the oil is a result of air-oil mixing, which is due to the break-up of the air-oil interface. Experimental observations [5] have shown that as the runner speed increases to a certain limit, oil in the reservoir suddenly goes unstable and becomes aerated, with the oil taking on a milky appearance. From Fig. 6.36, it may be postulated that this aeration takes place at the intersection of the free surface and the runner's outer surface, where the tangent of the free surface and the outer surface (in r-z plane) is at an angle. From the computational results, the angle becomes smaller as this runner speed increases. At the same time, a counter-clockwise vortex ring is formed just below the free surface. As the runner speed increases and exceeds a certain value, the oil's free surface may not be able to sustain the steep profile. Oil beneath the free surface may push back towards the runner surface and air then becomes wrapped into the oil. As soon as this oil contacts the runner surface, it accelerates and becomes separated from the oil behind Chapter 6. Computational Results and Analysis 79 it. This may explain how the oil aeration occurs in real bearings. The modeling of the simulated bearing assembly gave an insight into the flow field in the bearing assembly and provided a general explanation of the physical phenomena observed in the earlier studies [5]. In the following sections, numerical modeling of the test rig's oilwell tube annulus is described and the computational results are compared with data from the flow visualization tests. 6.3 Compu ted flow field in the oi lwell tube annulus In the numerical simulation, the flow field is computed with the A3D-CFD code. In order to compare the numerical results with the experimental data, the following data are used in the computations: The oil density: p = 876 kg/m3 < The oil temperatue: T = 27 °C v = 7.50 x 10 - 5 m2/s The rotating runner speed: WR = 3.14 rad/s Reynolds number is defined as V where R0 is the outer cylinder radius of the oilwell tube annulus and R0 = 0.13035 m. Figure 6.38 shows the computed flow field in the oilwell tube annulus. Secondary flow is mainly occurring in the region underneath the bottom of the rotating runner. At low speeds, the free surface is not severely deformed, but the deformation close to the rotating runner is still visible. The primary tangential velocity contours show a well stratified flow, as shown in Figure 6.38(b). Contour lines above the free surface are out of the fluid domain and should be ignored since it is produced by the plotting procedure and has no real physical meaning. Chapter 6. Computational Results and Analysis 80 Free surface (b) Figure 6.38: Computed flow field (Refl=711): (a) secondary flow vectors; (b) tangential velocity contours. Flow visualization, using micro air bubbles as tracers, was the principal means of ex-perimental observation of the secondary flow. As shown in Figure 6.38(a), the secondary flow is strongest in the region beneath the runner, which is where the flow visualization experiments were conducted. The computed secondary flow streamlines in this region are shown in Figure 6.39. From the secondary streamlines, it can be seen that the strong Chapter 6. Computational Results and Analysis 81 Figure 6.39: Computed secondary flow streamlines (Re#=711) vortex flow beneath the runner transfers momentum further into the oilwell tube annulus. Also, oil recirculation can be observed. 6.4 Path l ine of an air bubble An ideal tracer will follow a streamline of the flow in a steady flow. Figure 6.38 shows a typical computed flow field for the oil, from which the path of an oil particle (ideal tracer) may be computed. Parallel to the experiments, the computed flow field of Figure 6.38 in the flow visualization region is of particular interest. Figure 6.40 shows a close view of the computed flow field (Figure 6.38) in the flow visualization region. Assume that an oil particle is marked and is used as a tracer for the oil flow. As an ideal tracer, it will follow a streamline of the steady oil flow without any deviation. The path followed by the oil particle is identical to the streamline. Let's first examine Chapter 6. Computational Results and Analysis 82 Figure 6.40: Close-up views of computed flow in flow visualization region (Refl=711): (a) secondary flow vectors; (b) tangential velocity contours. Chapter 6. Computational Results and Analysis 83 the pathline of the oil particle. In the computed flow field, an oil particle at initial location (r0, ZQ, 0Q) will travel with the oil flow. After a certain time interval, it will have completed one complete circle (0 = 360°) about the axis of the rotating runner, but will probably end up at a different radial and/or axial locations because of the influence of the secondary flow. This is the moment that the particle reappears in the initial (r-z) plane as set up in the experiment. The particle continues traveling about the axis and reappears in the plane again at a new location 0 = 720°. The particle will keep reappearing in the {r-z) plane in the same manner, but usually at different locations. Based on the computational results (Figure 6.40), the pathline of an oil particle is estimated using the method described in Appendix H. Figure 6.41 shows the estimated location of the oil particle as time changes. In the computation, the (r-z) plane of the first appearance of the particle is chosen as a reference (t = 0, 0 = 0), which corresponds to the illuminated plane in the flow visualization test. Parallel to an experimental case, the initial location of the oil particle in the CFD modeling coordinate system is at (r0 = 131.3 mm, zn = 26.9 mm, 0Q = 0). The second appearance will be at 0 = 360° and the third at 0 = 720° and so on. Figure 6.41 can be re-plotted in terms of the secondary flow pathline in the (r-z) plane and this is shown in Figure 6.42. From the first appearance (t = 0), it takes 3.92 seconds for the particle to reappear in the reference plane for the second appearance (t = 0, 0 = 0). The third appearance of the particle in the plane occurs at t = 7.64 seconds. However, because of the nonhomogeneity of the air bubbles in the oil, their pathline tends to deviate from the prescribed streamline of the oil flow when the bubble is in an accelerating field. An immediate example is the upward movement of an air bubble in a still fluid as a result of the buoyancy force in the gravitational field. According to Chapter 6. Computational Results and Analysis 0.15 N .1 0.10 co o o J D n 0.05 T — 1 — 1 — 1 — 1 — i — 1 — 1 — • — r " " i ' 1 1 ^ -i , , , , i , , , ,: i , , , , 4 6 Time t (sec.) 800 CD 600 2. CD a> o) c co 2 400 . Q) 1 200 O 10 Figure 6.41: Computed location and time history of an oil particle (Reij=711). 0.035 r ~ 0.030 N 0.025 h 0.020 Second appearance (t=3.92s) Third appearance (t=7.64s) First appearance (t=0) _L 0.120 0.130 0.140 0.150 r(m) 0.160 Figure 6.42: Computed pathline of an oil particle (Rejt=711). Chapter 6. Computational Results and Analysis 85 Eq. B.48, a small air bubble at rest will accelerate and reach a terminal velocity given by d2 12 v where da is the air bubble diameter and u is the kinematic viscosity of the fluid. The terminal velocity is in the opposite direction of the gravity g. For an air bubble of 1.0 x 1 0 - 4 m in diameter in the test oil with a kinematic viscosity v = 7.50 x 1 0 - 5 m 2 / s , the resulting terminal velocity is 1.1 x 1 0 - 4 m/s. The motion of an air bubble in a rotating flow experiences a more complicated process compared to that in still oil. Detailed analysis is given in Appendix H . Based on this analysis, the significant body forces on the bubble are gravitational and centrifugal. By analogy with the buoyancy due to gravity, the relative velocity of the air bubble to the oil particle at the same location may be written as: d2 VaT = - t ^ - O b 12 v where a& is the vector form of the body forces. Based on this formulation, the deviation of the air bubble from the oil particle (stream-line) may be evaluated. It is worth noting that in the calculations the oil particle used as a marker is considered to be an imaginary object. As the path of the air bubble con-stantly deviates from the path of the oil particle, the oil particle is always referred to as an imaginary one at the location of the air bubble during the step-by-step calculations. Details are shown in Appendix H . Considering the buoyancy effects due to the body forces as described above, the pathline of an air bubble in the flow field is computed. Figure 6.43 shows how the location of a 0.1 mm diameter air bubble in the oil flow field (as shown in Figure 6.40) varies with time. In this figure, the oil particle location is plotted as a dotted line for the purpose of comparison. It can be seen that the buoyancy effects on the air bubble location are visible, but also very small. Chapter 6. Computational Results and Analysis 86 0.15 E, N . i 0.10 co cj o a) JD - O 3 0.05 i r "* T i | i i i i , , . 1 1 , 1 ! . . 1 ' ' ' l/--• / , , . , 1 . , , , 1 1 1 1 1 1 . ! , 1 . . . . 800 CD 600 2, CD 03 O ) c CO To 400 CD . 0) E y b 200 4 6 Time t (sec.) 10 Figure 6.43: Computed location of the air bubble, and its time history (Refl=711; da = 0.1 mm). Figure 6.44 compares the pathlines of the air bubble with the ideal tracer - oil particle (dotted line). In terms of time history, the air bubble reappears in the starting (r-z) plane at t — 3.92, 7.66 seconds, which deviates slightly from the oil particle (t = 3.92, 7.64 seconds) due to the buoyancy effect. In terms of the locations of these reappearances of the air bubble, both reappearances of the bubble seem to occur ahead of the oil particle. In other words, in this situation, the buoyancy force slows down the secondary motion of the air bubble. Based on the computed pathline of the air bubble, it is noted that, after the first appearance, the air bubble has moved upward, taking a slightly different path in com-parison to the oil particle. This is mainly attributed to the buoyancy effect due to the gravitational acceleration. As the bubble moves, it is exposed to a new velocity of oil Chapter 6. Computational Results and Analysis 87 Second appearance (t=3.92s) Third appearance (t=7.66s) 0.035 r-— 0.030 -N 0.025 -0.020 V- 1 1 _ i i I • • • i I • . • i I • • • • I i • . • L 0.120 0.130 0.140 0.150 0.160 r(m) Figure 6.44: Computed paths of the air bubble and an oil particle (Ren=711). at its new location, since the oil flow field is very location-dependent. As the bubble approaches the outer radial location, it turns and moves radially inward to return back on the pathline of the oil particle. From the earlier analysis, it has been shown that the buoyancy effect on the air bubble due to acceleration of the oil flow is dependent on the size of air bubble and the kinematic viscosity of oil. The terminal velocity of the bubble relative to the bulk oil flow is proportional to the surface area of the bubble. As da —> 0, the buoyancy effects tend to vanish and the pathline of the air bubble approaches the pathline of the oil particle. The computed pathline of the air bubble is compared to the experimental observation. In the experiment, appearances of the air bubble in the illuminated plane are accurately timed by counting the recorded video images frame by frame. The image is recorded at a rate of 30 frames per second, so the accuracy of the timing sequence is ±1 /60 seconds. The location of the air bubble in the illuminated plane is obtained by measuring the location of the bubble in the optically distorted image and then correcting this mea-surement with a restoration function. Figure 6.45 shows the locations of the air bubble, Chapter 6. Computational Results and Analysis 88 as they appeared in the illuminated plane during a flow visualization test. In the same figure, computed results are shown as dotted lines and symbols. From fluid dynamics, 0.035 h 0.030 0.025 h 0.020 Second appearance (t=3.87s) Third appearance (t=7.67s) First appearance (t=0) _L _ l I I I L—1_ J I I I L-J—I—I—I—I—I—I—I—I—I—I—L. 0.120 0.130 0.140 0.150 0.160 r(m) Figure 6.45: Appearances of air bubbles in a flow visualization test (Refl=711). it is clear that the location and timing of the appearance of an air bubble are greatly dependent on the flow field of the oil. Therefore, precise prediction of the primary and secondary velocity components of the oil flow is essential in computing, with reasonable accuracy, the location and timing of the appearance of an air bubble. With regard to this study, the appearance time predicted by the computation agrees very well with the experimental observation. However, the computed locations of the appearance of the air bubble seem to be somewhat ahead of the experimental ones. It should be kept in mind that, from the starting time t — 0, the air bubble has completed one revolution along the secondary pathline before its first reappearance (the second appearance) and has completed two revolutions along the secondary pathline before its second reappearance (the third appearance). Therefore, a reasonable method for evaluating the difference in location between the computation and the experiment is to compare the difference of the two secondary path lengths. From the computation, the air bubble has traveled 63.7 mm Chapter 6. Computational Results and Analysis 89 and 131.6 mm in the (r-z) plane from the first appearance of the bubble to the second and third appearances respectively. From Figure 6.45, it can be seen that the second and third measured appearances were about 2.2 mm and 3.4 mm ahead of the computed values, which results in only 3.5% and 2.6% differences, respectively. There are many factors that may attribute to this difference between the experi-mental and computed results. Experimental errors resulting from measurement of oil temperature and therefore oil viscosity, runner speed and air bubble location recorded at the time of the flow visualization tests will all be important contributors. Their effects on the study will be described in later sections. Numerical assumptions relating to the experimental setup and the size of the air bubble may also be important contributors to the differences between theory and practise. It is of interest to note that, even though there are noticeable differences between the computed and experimentally measured locations of the air bubble (in the illuminated plane), agreement between theory and practise was close in terms of the path taken by the air bubble. As discussed earlier, images taken from the flow visualization tests were optically distorted as a result of the curvature of the rig's outer surface. It is relatively easy to restore some points of the image using the optical correction function (Eq. A.44), while restoring the whole image is somewhat more difficult and time consuming. However, the coordinates of the computed results can be altered easily if the transform function is known. Therefore, in order to compare the computed and experimental results, com-puted "images" are intentionally distorted to simulate the optical distortion. The image distortion was done using the following transform functions. rd = -1.138 + 0.9706 r - 0.01233 r 2 + 0.0001079 r 3 (6.42) zd = z (6.43) Chapter 6. Computational Results and Analysis 90 where r<j, r, zd and z are in millimeters. The image coordinates originate from the lower-left corner of the flow region in the im-age, which is different to the coordinate system used in the CFD modeling. Equation 6.42 is obtained from a least squares regression of data obtained from Equation A.44. Sub-sequently, using the transform functions, the pathline of the air bubble (Figure 6.44) is "distorted" as shown in Figure 6.46. Figure 6.46: Computed pathline of the air bubble after the "distortion" (Re.R=711; da = 0.1 mm). Comparing Figure 6.44 with Figure 6.46, the distortion becomes very obvious. By intentionally "distorting" the computed pathline, it can be seen that the computed path-line is very similar to the measured pathline (Figure 3.15). Hence, from here on forth, computed pathlines will be presented as a distorted image as described above if not specifically mentioned. In accordance with the experimental results presented in Figure 3.16, pathlines are computed for different starting locations. These correspond to the different locations of the air bubble probe in the experiments. Figure 6.47 shows the computed pathlines of the air bubble for different starting locations. In the figure, small circles are starting Chapter 6. Computational Results and Analysis points as interpreted from the experiments (Figure 3.16). 91 (e) (f) Figure 6.47: Computed pathlines ("distorted") of the air bubble at different starting points (Reji=949; da = 0.1 mm). In order to compare pathlines at different runner speeds, computations for the same experimental conditions as in Figure 3.18 have been made. In the experiment, the air bubble probe was placed at a certain fixed location. However, for the different runner speeds, the starting points (that is the first appearances) were not at a fixed location in the images. This is due to the fact that the air bubbles have traveled for a while (traveled 22° in the ^-direction) before their first appearance in the image plane, and the location Chapter 6. Computational Results and Analysis 92 of the appearance is determined by the flow field. Similarly, using the starting points interpreted from the experiments, air bubble pathlines were computed. Figure 6.48 shows the computed pathlines of a 0.1 mm diameter bubble for different runner speeds. (c) Figure 6.48: Computed pathline ("distorted") of the air bubble (da = 0.1 mm) at different runner speeds: (a) ReR - 711; (b) ReR = 949; (c) ReR = 1185. 6.4.1 Effect of bubble size In a still fluid, the terminal velocity of a micro air bubble due to buoyancy is propor-tional to the surface area of the bubble. Considering the same flow conditions as shown in Figure 6.43, it can be shown that different sizes of air bubbles will experience different buoyancy forces and therefore different pathlines. To evaluate the effect of air bubble size, different size air bubbles (da = 0.05, 0.1 and 0.2 mm) are examined. Figure 6.49 shows the computed locations of the air bubbles and how these locations are affected by bubble size. Chapter 6. Computational Results and Analysis 93 0.00 0 2 4 6 8 Time t (sec.) Figure 6.49: Location of air bubble and the effect of bubble size (Ren = 711). The computed locations are re-plotted in terms of the secondary flow pathlines. Fig-ure 6.50 shows the computed (undistorted) pathlines for different sizes of air bubble. The pathlines of the 0.05 mm diameter and 0.1 mm diameter air bubbles are very close to each other, but the pathline for the 0.2 mm diameter air bubble experiences is quite different. Any deviation of the pathline caused by a difference in the size of the air bubble will lead to a different time history of the air bubble in the flow visualization plane. In terms of appearances in the (r-z) plane as described earlier, the time t of reappearance of the air bubble is listed in Table 6.3. As before, the first appearance of the air bubble is set at In Table 6.3, data for da = 0 is from a calculation that pertains to the oil particle, since deviation of an air bubble because of buoyancy tends to vanish as da —> 0. The de-viation in time t for reappearances of the air bubble increases as air bubble size increases. t - 0. Chapter 6. Computational Results and Analysis 94 0.035 ^ 0.030 E. N 0.025 0.020 0.120 0.130 0.140 0.150 0.160 r(m) Figure 6.50: Pathlines affected by air bubble sizes (ReR = 711). Air bubble diameter da (mm) 0 0.05 0.1 0.2 The second appearance 3.92 (s) 3.92 (s) 3.93 (s) 3.94 (s) The third appearance 7.64 (s) 7.66 (s) 7.66 (s) 7.72(B) Table 6.3: Time t of reappearances of air bubbles in the (r-z) plane. Therefore, small air bubbles make good tracers in oil flow. Actually, the terminal veloci-ties for the three bubble sizes (da = 0.05, 0.1, 0.2 mm) are Va, 4Va and 16Va respectively, if Va is the terminal velocity for da = 0.05 mm, as the terminal velocities are proportional to the bubble surface area. 6.4.2 Effect of viscosity As pointed out earlier, air bubble size and the oil's viscosity are important factors in respect to the movement of the air bubble as a result of the buoyancy effect. However, the oil's viscosity also has a large influence on the flow field and therefore on the pathline of the air bubble. da(mm) J • • A Second appearances — • — 0.05 | • o A Third appearances 0.1 — • — o.2 )) First appearances (t=0) \ > ' 1 ' I. • I 1 1 1 1 L Chapter 6. Computational Results and Analysis 95 Figure 6.51 shows, for a constant runner speed of UJR = 3.14 rad/s, computed sec-ondary flow fields for different oil viscosities. At the same runner speed, secondary flow underneath the rotating runner was greatly affected by the viscosity of the oil. It is shown that, as the viscosity decreases, the velocity of the secondary flow increases. By paying close attention to the velocity profiles close to the core of the vortex, the effect of viscosity can be clearly observed. The velocity profiles at r — 150.2 mm, along the z direction, are shown in Figure 6.52. From Figure 6.52, it is clear that viscosity has a large effect on the magnitude of the secondary flow (u, v). As viscosity decreases, the magnitude of both u and v increases. Moreover, flow patterns are also altered, as seen from the velocity profiles of v. The tangential velocity basically has a linear profile, much like a typical parallel Couette flow, although, as the viscosity decreases, the velocity profile does become nonlinear. For example, for v = 0.40 x 10 - 4 m2/s, the upper portion of the velocity w has decreased, while the lower part of w has increased. This is caused by secondary flow convection. As shown in Figure 6.52, the velocity component u is much larger than the velocity component v. If only u is considered, the secondary flow convective interaction with the primary flow (u^j may be described as follows. Basically, the tangential velocity gradient ( § 7 ) is positive (refer to Figure 6.40(b)) in the flow visualization region, while the velocity u is positive in the upper part and negative in the lower part of the region. In the positive (fj^) flow field, positive u tends to slow down the tangential velocity w and vice versa. As the secondary velocity u increases (as viscosity decreases), the interactive effect increases and the non-linearity of the tangential velocity w becomes more and more noticeable as shown in Figure 6.52(c). The strong effect of viscosity on the secondary flow beneath the bottom of the run-ner may even be transferred to the flow in the oilwell tube annulus. As an example, the velocity profiles along the mid-span of the annulus (r = 124.0 mm) are shown in Chapter 6. Computational Results and Analysis 96 0.05 m/s 0.05 m/s (a) 0>) 0.05 m/s (c) (d) Figure 6.51: Computed secondary flow fields at a constant runner speed U>R — 3.14 rad/s: (a) i/=1.235 xlO- 4 m2/s; (b) i/=1.09 xlO" 4 m2/s; (c) i/=0.75 xlO" 4 m2/s; (d) i/=0.40 xlO" 4 m2/s. Chapter 6. Computational Results and Analysis 97 (b) l , , , , l , , , y - 0 . 0 0 8 - 0 . 0 0 4 0 . 0 0 0 v (m/s) Figure 6.52: Computed velocity profiles at r=150.2 mm = 3.14 rad/s; da = 0.1 mm). Chapter 6. Computational Results and Analysis 98 Figure 6.53. At the mid-span location of the oilwell tube annulus, the effect of viscosity on sec-ondary flow is similar to that observed in the region beneath the runner. Secondary flow is strongest at the bottom of the runner region, which is basically the source of secondary flow in the oilwell tube annulus. Further up the oilwell tube annulus, in the z direction, secondary flow eases up but closer to the free surface, it again gets strong. The tangential velocity profiles show no strong effect from viscosity. However, in the upper region of the oilwell tube annulus (above the bottom of the runner), the end effects from the bottom of the runner have a stronger effect on the tangential velocity profile at lower viscosity. In other words, end effects from the region underneath the runner penetrate deeper into the upper part of the oilwell tube annulus as the oil viscosity is lowered. Based on the velocity distributions from the CFD computations, locations of an air bubble can be traced. Figure 6.54 shows computed locations of the air bubble, taking changes in oil viscosity into account. The computed locations are re-plotted in terms of secondary flow pathlines. Fig-ure 6.55 shows the "undistorted" pathlines of the air bubble, as affected by changes to the viscosity of the oil. In terms of appearances in the (r-z) plane (as described earlier), the times (t) are listed in Table 6.4. As before, the first appearance of the air bubble is set at t = 0. Oil viscosity v (m2/s) 0.40 xlO" 4 0.75xl0~4 1.09xl0"4 The second appearance 3.84 (s) 3.93 (s) 3.76 (s) The third appearance 8.08 (s) 7.66 (s) 7.33 (s) Table 6.4: Time t for reappearance of the air bubble in the (r-z) plane. For a constant runner speed, the magnitude of the secondary flow velocities is very Chapter 6. Computational Results and Analysis 99 • ' O . O Q O ' -0.008 -0.004 ffi>00 0.004 0.008 v (m/s) z(m) 0.060 0.040 0.020 O0Q01 1 1 ' 1 1 1 ' 1 1 1 1 1 ' 1 1)000 0.100 0.200 w (m/s) Figure 6.53: Computed velocity profiles at r=124.0 mm (U>R = 3.14 rad/s; da — 0.1 mm). Chapter 6. Computational Results and Analysis 100 0 2 4 6 8 Time t (sec.) Figure 6.54: Locations of the air bubble as affected by oil viscosity (UIR = 3.14 rad/s; da = 0.1 mm). 0.035 ? 0.030 N 0.025 0.020 I • • A D o A Second appearances Third appearances i i 0.120 I i 0.130 v (m2/s) • 0.40x10"4 " 0.75x10"4 "1.09x10^ 0.140 i I i 0.150 0.160 r(m) Figure 6.55: Pathlines of the air bubble as affected by viscosity = 3.14 rad/s; da — 0.1 mm). Chapter 6. Computational Results and Analysis 101 dependent on oil viscosity, which leads to the conclusion that the locations of the air bubble reappearances in the (r-z) plane are strongly affected by the oil viscosity. However, tangential velocity profiles are not severely affected by oil viscosity, as are the secondary flow pathlines, with the result that the times of the reappearances are not greatly affected by changes in the viscosity of the oil. 6.4.3 Effect of runner speed In the experiments, the location of an air bubble at its first appearance varied with runner speed, even though the location of the air bubble probe was fixed. This, in part, was due to the setup of the experimental apparatus, as described earlier. It is of interest to know the effect of runner speed on the oil's flow field, and the pathline of the air bubble. In this section, pathlines of air bubbles, released at the same starting location, are examined. Similarly, the effect of runner speed on the oil's velocity profiles will also be examined. Figure 6.56 shows velocity profiles in the z direction, across the core of the secondary vortex flow. As expected, the magnitude of the secondary flow increases as the runner speed is increased. The non-linearity of the secondary flow can be observed from Fig-ure 6.56(b), and the similarity of the secondary flow patterns observed for increasing the runner speed and decreasing the viscosity is also evident. The primary tangential velocity profiles exhibit a similarity to the cases where the viscosity is changed. Linear tangential velocity profile is altered as runner speed is increased. Figure 6.57 shows the computed velocity profiles along the mid-span of the annulus (r = 124.0 mm). At the mid-span location of the annulus, the magnitude of tangential velocity of the oil increases linearly with runner speed WR. However, the magnitude of the secondary flow velocity increases non-linearly with runner speed U>R. More precisely, it increases at a faster rate than the increasing rate of runner speed. Chapter 6. Computational Results and Analysis -0.080 -0.040 0.000 0.040 u (m/s) coR=5.23 rad/s coR=4.19 rad/s coR=3.14 rad/s co„=2.09 rad/s z(m) • 0.032 -0.016 -0.012 -0.008 -0.004 0.000 v (m/s) 0.020 toB=5.23 rad/s coR=4.19 rad/s coR=3.14 rad/s coR=2.09 rad/s 0.000 0.200 0.400 0.600 0.800 w (m/s) Figure 6.56: Computed velocity profiles at r=150.2 mm [y mm). Chapter 6. Computational Results and Analysis 103 coR=5.23 rad/s wR=4.19 rad/s coR=3.14 rad/s coR=2.09 rad/s z(m) 0.060 -0.010 -0.005 I.000 u (m/s) 0.005 (a) z(m) 0.060 0.040 0.020 °'0<f# coR=5.23 rad/s coR=4.19 rad/s coR=3.14 rad/s C0o=2.09 rad/s 000 0.005 0.010 0.015 0.020 v (m/s) (b) z(m) 0.060 h 0.040 0.020 O.OOj). coR=5.23 rad/s wR=4.19 rad/s coR=3.14 rad/s coR=2.09 rad/s 000 0.100 0.200 0.300 0.400 w (m/s) (c) Figure 6.57: Computed velocity profiles at r=124.0 mm (i/=7.5 x I O - 5 m2/s; «i a=0.1 mm). Chapter 6. Computational Results and Analysis 104 As shown earlier, as the runner speed changes so do all three velocity components of the oil. Figure 6.48 shows computed pathlines for different runner speeds, which can be compared with the results from the flow visualization tests. However, due to practical considerations in relation to the experiments, the location of the first appearance of the air bubble in the (r-z) plane could not be fixed at the same starting point for different runner speeds. This would certainly lead to different pathlines. Now, the effect of the runner speed on the pathline of the air bubble will be numerically examined for the same starting point with other conditions the same as the experiment (refer to Figure 3.18). Based on the velocity distributions from the CFD computations, locations of an air bubble can be traced. Figure 6.58 shows computed locations of the air bubble, as they are affected by runner speed. 4 6 8 Time t(sec.) 1400 Figure 6.58: Locations of the air bubble as affected by runner speed (v = 7.5 x 10 5 m2/s; da — 0.1 mm). Chapter 6. Computational Results and Analysis 105 The computed locations are re-plotted in terms of secondary flow pathlines. Fig-ure 6.59 shows "undistorted" computed pathlines for different runner speeds (O;.R = 2.09, 3.14 and 4.19 rad/s). In terms of appearances of the air bubble in the (r-z) plane as 0.035 0.030 0.025 0.020 • • A Second appearances • oA Third appearances coR (rad/s) - - • - - 2.09 —•— 3.14 —A- - 4.19 _ l_ 0.120 0.130 0.140 0.150 r(m) 0.160 Figure 6.59: Pathlines of the air bubble as affected by runner speed (v = 7.5 x 10 5 m2/s; da = 0.1 mm). described earlier, time t of reappearances is listed in Table 6.5. As before, the first ap-pearance of the air bubble is set to t = 0. When the runner speed changes, all three Runner speed UJR (rad/s) 2.09 3.14 4.29 The second appearance 5.49 (s) 3.93 (s) 2.83 (s) The third appearance 10.67 (s) 7.66 (s) 6.17 (s) Table 6.5: Time t for reappearance of the air bubble in the (r-z) plane. velocity components change accordingly. Therefore, both the location and time of the air bubble reappearances are affected by runner speed. However, air bubble pathlines which have the same starting point are not greatly affected by runner speed. Chapter 6. Computational Results and Analysis 106 6.5 Secondary vorticity at the vortex core The strength of the secondary vortex flow can be presented in terms of the local secondary vorticity u}gc at the core of the vortex given by _ (du dv\ 6 c \dz d r / c o r e o f v o r t e x ' The secondary vorticity is calculated from the computed flow fields for the cases presented above. The results of computed uec (rad/s) are listed in Table 6.6. UJR = 3.14 (rad/s) v (m2/s) 0.40 xlO" 4 0.75 xlO" 4 1.09 xlO" 4 1.24 xlO" 4 <JJ6C (rad/s) 3.74 2.71 1.97 1.75 v = 0.75 x 10"4 (m2/s) U>R (rad/s) 2.09 3.14 4.19 5.23 wee (rad/s) 1.28 2.71 4.37 6.02 Table 6.6: Calculated secondary vorticity u)$c (rad/s) at core of the secondary vortex. As indicated earlier, velocity profiles exhibit similarities for both changes in oil vis-cosity and runner speed. This similarity may be presented in terms of a non-dimensional vorticity (U>0C/U>R). Based on this result, the vorticity at the core of the secondary vortex, in terms of a non-dimensional vorticity, may be correlated to the Reynolds number Re^. From Table 6.6, the non-dimensional vorticity vs. Reynolds number Ren are plotted in Figure 6.60. One can say that this Reynolds number defines the major properties (vortex strength, pathlines) of the secondary flow field. All the other features of the flow not easily quantified are also expected to depend primarily on this Reynolds number. Chapter 6. Computational Results and Analysis 107 Chapter 7 Concluding Remarks and Future Work 7.1 Concluding remarks An oilwell tube annulus test facility was designed for the purpose of performing flow visualization in the oilwell tube annulus region of a vertical bearing assembly. Using a light-sheet visualization technique with micro air bubbles as the tracer, secondary vortex flow patterns that exist underneath the bottom of the runner in the oilwell tube annulus test rig were observed. Appearances of air bubbles in the radial-axial plane have been timed and analyzed to quantitatively examine the oil flow field and compare the experimental results with the numerical results. A three dimensional CFD code for examining laminar axisymmetric flow with free surface was developed for this project. To examine the different aspects of the A3D-CFD code, the code was verified for the following three cases: • Shape of the free surface (oil-air interface) in a rotating annulus. Computed free sur-face shapes have shown very good agreement with the experimental measurements. The effect of the non-dimensional radial acceleration (gR) on the non-dimensional shape of the free surface (AH*) was well predicted numerically in comparison with the experimental results. • Velocity distributions close to a rotating disk in a fluid at rest. Good agreement be-tween the computational velocity profiles and the analytical solution was achieved, except that there were some discrepancies in the magnitude of the velocities (refer 108 Chapter 7. Concluding Remarks and Future Work 109 to Figure 4.27). These discrepancies may be attributed to the fact that the disk radius was modeled as of finite value in the computational approach while infinite disk radius was assumed in the analytical approach. • Taylor vortex flow patterns in a counter-rotating annulus. Taylor type vortices were predicted numerically. The computed vortex flow patterns were similar to the Taylor's descriptions. In terms of the onset of the Taylor vortices (flow instability), the instability in the CFD computation occurred at u>i/v = 4.0 ~ 4.5 x 107 rad/m2, while Taylor experimentally observed the instability at UJI/U = 4.3 x 107 rad/m2. It has been shown that the A3D-CFD code is capable of correctly predicting different features of fluid flow. Using the A3D-CFD code, oil flow in the oilwell tube annulus was examined. In the oilwell tube annulus, the primary swirling flow was well stratified in the oilwell tube annulus as shown in Figure 6.36. A strong secondary vortex flow underneath the bottom of the rotating runner is calculated. This vortex flow expands into the bottom of the oilwell tube annulus. The presence of this vortex flow, together with the free surface, makes the flow unique in comparison to the classical rotating cylinder problem. In the region of the vortex core, the magnitude of the secondary velocity increases, as the Reynolds number Re# increases. As a result of the convective interactions between the secondary and primary flows, the primary flow (i.e. the tangential velocity profile of the oil) was also affected by a change of the oil viscosity at constant runner speed. The secondary vortex strength (in terms of the local vorticity cogc at the core of the secondary vortex) was strongly affected by the Reynolds number (Ren) of the flow. The non-dimensional vorticity (wedwR) was correlated with the Reynolds number (Rejj). Pathlines of air bubbles were computed using the numerical technique developed for tracing an air bubble. The computed pathlines of an air bubble were compared with Chapter 7. Concluding Remarks and Future Work 110 the experimental pathlines and reasonably good agreement was obtained. Although the secondary vortex strength is correlated to the Reynolds number Rejj, the secondary flow pattern, in terms of the secondary flow pathline, is not substantially affected by the Reynolds number. Effect of air bubble size was evaluated and small size bubble makes a better tracer of the flow. 7.2 Future work The work reported in this study is limited to relatively low Reynolds numbers. How-ever, using the A3D-CFD code, an extensive numerical parametric study can be carried out, paying particular attention to the velocity at the bottom boundary region of the annulus. In this numerical parametric study, variable parameters may include oilwell tube geometry, still oil level, oil viscosity and runner speed. These computational results may improve the experimental design of the oilwell tube geometry, which may lead to the final goal - a better industrial design. The A3D-CFD code may be improved in terms of the convergence of the numerical solution at high Reynolds numbers. Special attention will be given to proper presentation of the free surface conditions. At a more extensive level, modifications of the oilwell tube geometry can be made based on the computational studies. Experiments can be further conducted with the modified test rig to determine what triggers the onset of unstable flow. Bibliography [1] Ahlborn, F., "On the Mechanism of Hydrodynamic Drag", abh. a.d. Geb.d. Natur-wiss, Vol.17, 1902, pp.8-37. [2] Andereck, C. D. and Hayot, F., "Ordered and Turbulent Patterns in Taylor-Couette Flow", Edited, Nato ASI Series, Plenum Press, New York, 1992. [3] Berggren, R. and Engelhart, R., "Large Vertical Motor Upper Bearing Reservoir Oil Leakage", General Electric report 75MPL231, Schenectady, USA, 1975. [4] Bertela, M. and Goli, F., "Laminar flow in a cylindrical container with a rotating cover", J. Fluid Eng., Vol.104, 1982, pp.31-39 [5] Brockwell, K. and Kleinbub, D., "Flow Instability in the Oilwell Tube of Vertical Bearing Assemblies - A Parametric Study", National Research. Council of Canada report IMR-T&R-CTR-002, 1993. [6] Burkhalter, J. E. and Koschmieder, E. L., "Steady supercritical flow", J. Fluid Mech., Vol.53, Part 3, 1973, pp.547-560. [7] Cole, J. A., "Taylor vortex instability and annulus length effects", J. Fluid Mech., Vol.75, Part 1, 1976, pp.1-15. [8] Coles, D., "Transition in Circular Couette Flow", J. Fluid Mech., Vol.21, Part 3, 1965, pp.385-425. [9] Couette, M. , "Etudes sur le frottement des liquides", Ann. Chim. Phys., ser.6, vol.21, 1890, pp.433-510. [10] Dukler, A.E. , Wicks, M. and Cleveland, R., "Pressure drop and hold-up in two-phase flow: Part A, A comparison of existing correlations; Part B, An approach through similarity analysis", A.I.Ch.E.J., vol.lO(l), 1964, pp.38. [11] Escudier, M. P., "Observations of the flow produced in a cylindrical container by a rotating endwall", Experiments in Fluids, Vol.2, 1984, pp.189-196. [12] Ferguson, J., "Oilwell Tube Test Rig - Results with. Ring and Pilot Tubes", General Electric Canada Technical Report, 91-JHF-l, 1991. [13] Flowle, T. I., "Aeration in Lubricating Oils", Tribology international, Vol.14, June 1981, pp.151-157. I l l Bibliography 112 [14] Fox, R. W. and McDonald, A. T., "Introduction to Fluid Mechanics", John Wiley & Sons, Inc., New York, 1985, 3rd Edition, pp.684. [15] Fujimura, K., Koyama, H. S. and Hyun, J. M. , "Vortex breakdown in a difFerentially rotating cylindrical container", Proc. 6th int. symp. on transport phenomena and dynamics of rotating machinery, Vol.2, 1996, pp.67-76. [16] Haberman, W. L. and Morton, R. K., "An Experimental Study of Bubbles moving in Liquids", ASCE Proceedings, Vol.80, No.387, 1954. [17] Harlow, F. H. and Welch, J. E. , "Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface", The physics of Fluids, Vol.8, 1965, pp.2182-2189. [18] Hirt, C. W. and Nichols, B. D., "Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries", J. Computational Physics, Vol.39, 1981, pp.201-225. [19] Isbin, H. S., Moy, J. E . and Da Cruz., A. J. R., "Two-phase steam-water critical flow", A.I.Ch.E.J., 1957, Vol.3, pp.361. [20] Joseph, D. D., "Stability of Fluid Motions I", Springer-Verlag, Berlin, 1976, ppl32-136. [21] Kaiser, J. , "Free Surface Stability in a Simulated Vertical Bearing Seal Assembly", University of Waterloo report, 1989. [22] Kobus, H. E . , "Characteristics of self aerated free-surface flow" editors: Rao, N.S.L., and Kobus, H. E. , Erich Schmidt Verlag, 1976, pp.47-65. [23] Lamb., H., "Hydrodynamics" Dover Publications, New York, 6th ed., 1945, ppl23-124. [24] Mallock, A., "Determination of the Viscosity of Water", Proc. Royal Soc, A 45, 126, 1888. [25] McGregor, I., "The Vapor-Screen Method of Flow Visualization", J. Fluid Mech., Vol.11, No.4, 1961, pp.481-511. [26] Monk, T., "Summary of Oilwell Tube Tests on Thrust Bearing Test Rig", General Electric Canada Technical Report, 93-TDM-01, 1993. [27] Nichols, B. D., Hirt, C. W. and Hotchkiss, R. S., "SOLA-VOF: A solution algo-rithm for transient fluid flow with multiple free boundaries", Los Alamos Scientific Laboratory Report, Rreprint LA-8355. Bibliography 113 [28] Patankar, S. V., "Numerical heat transfer and fluid flow", McGraw-Hill Inc., 1980, pp.67-68. [29] Peyret, R. and Taylor, T. D., "Computational Methods for Fluid Flow", Spring-Verlag New York Inc., 1983, pp.324-337. [30] Piao, Y. and Brockwell, K., "Modeling of Oil Flow with Free Surface in the Oilwell Tube of Self-Contained Bearing Assemblies", Proc. 6th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Honolulu, U.S.A. 1996, Vol.1, pp.669-678. [31] Piao, Y. and Brockwell, K., "Secondary Flow Patterns in a Vertical Guide Bearing Assembly", Proc. 2nd Annual Conference of the CFD Society of Canada, Toronto, 1994, pp.155-161. [32] Piao, Y., "A3D-CFD code - a solver for 3D Axisymmetric laminar flow with free surface", Naval Architecture Lab. Technical Report 96-03, University of British Columbia, Canada, 1996. [33] Rayleigh, L., "On the Dynamics of Revolving Fluids", Proc. Roy. Soc. London Ser. A, Vol.93, 1916, pp.148-154. [34] Schlichting, H., "Boundary Layer Theory", 7th ed., McGraw-Hill Inc., 1979, pp.289-293. [35] Schraub, F. A., Kline, S. J. , Henry, J. , Runstadler, P. W. and Little, A., "Use of Hydrogen Bubbles for Quantitative Determination of Time-Dependent Velocity Fields in Low-Speed Water Flows", J. Basic Eng., Vol.87, 1965, pp.429-444. [36] Taylor, G., "Stability of a Viscous Fluid Contained Between Two Rotating Cylinders", Trans, of the Royal Society, London, England, series A, 223, 1923, pp.289-293. [37] Vanyo, J. P., "Rotating fluids in engineering and science", Butterworth-Heinemann, 1993, pp.111. [38] Wei, T., Kline, E. M. , Lee, S. H.-K., Woodruff, S., "Gortler Vortex Formation at the Inner Cylinder in Taylor-Couette Flow", J. Fluid Mech., Vol.245, 1992, pp.47-68. [39] Wendt, F., "Turbulente Stromungen Zwishen Zwei Rotierender Konaxialen Zylindern", Ingenieur-Arcbiv, Vol.4, 1933, pp.577-595. [40] White, F. M. , "Fluid Mechanics", 2nd ed., McGraw-Hill Inc., 1986, pp.142-146. Bibliography 114 [41] Wild, P. M. , Djilali, N. and Vickers, G. W., "Experimental and Computational As-sessment of Windage Losses in Rotating Machinery", J. Fluid Engineering, Vol.118, 1996, pp.116-122. Appendix A Image restoration of optical distortion in the experiment In the present experimental setup, the image recorded with the video camera is optically distorted due to refraction of light as shown in Figure 2.7. In order to evaluate the distortion, a true scale is placed in the flow aligned with the sheet of light in the experiment. The scale is then distorted as shown. Figure A.61 illustrates the true scale and its distorted image. The true scale consists of four dash-dotted 10 mm x 10 mm squares, while the shaded area shows the the distorted scale of the four radially "compressed" squares. The coordinates shown originate from the lower-left corner of the fluid domain of the image (refer to Figure 2.7). Keep in mind that this is different from the origin of coordinates in the CFD modeling of the oil flow. It should also be noted that the distortion is only in the r direction. z (mm) Distorted scales True scales 10 0 10 20 30 40 r (mm) Figure A.61: Illustration of the true scales and its distorted image. 115 Appendix A. Image restoration of optical distortion in the experiment 116 In Figure A.61, the thick-lined boxes just show the distortion of a square in the true scale. Also, the straight slope line (1:4) in the true scale is distorted non-linearly as shown. Using the distortion of the slope line, a point -Pd(r*d,Zd) in the distorted image can be converted to its true coordinates P(r, z). To find out the true coordinates P(r, z), draw a vertical line (parallel to z) passing through Pd and find the intersection point a of this line and the distorted slope line in the distorted image. Since the vertical z coordinate is not distorted, we draw a horizontal line (parallel to r). The intersection point b of this horizontal line and the true slope line is the true image of point a. Finally, we draw a vertical line passing through point b and find the point P from z = zd. Now the image point Pd is transformed into P in the true scale. In order to restore the image, the coordinate reading (rd,zd) in the distorted image need to be transformed into the true coordinates (r, z). This transformation may be described by complicated mathematics. But here a simple physical method is used. To describe the axial optical compression of the image, define Cr as the radial optical compression ratio C = — drd where dr is an actual infinitesimal radial length and drd is its image length after the optical distortion. In Figure 2.7, the optical compression ratio at center of each rectangular cell may be approximated as Az Cr = Ard where Az and Ard are vertical and axial lengths of the rectangular cell in the image respectively. Figure A.62 shows the changes of the optical compression ratio Cr with rd. The symbols were obtained from direct measurements of the distorted cell in Figure 2.7. Using the least square method, a second order polynomial expression of the the optical Appendix A. Image restoration of optical distortion in the experiment 117 2.2 2.0 O 1.8 1.6 1.4 1.2 1.0' 10 15 rd (mm) 20 25 Figure A.62: Optical compression ratio Cr vs. rd. compression ratio for the experimental set-up is obtained: Cr = 0.000996 r\ + 0.008228 rd + 1.295 where rd is in millimeters. This polynomial expression of the optical compression ratio is shown as a solid line in Figure A.62. The radial coordinate of the image can be converted to actual radial coordinate, r = fTd Crdrd = 0.000332 r3d + 0.004114 r\ + 1.295 rd, (A.44) Jo where both r and rd are in millimeters. Figure A.63 is the graphic expression of the trans-formation function (Eq. A.44). This formula can be used in computer image restoration. Appendix A. Image restoration of optical distortion in the experiment 118 rd (mm) Figure A.63: Graphic expression of the optical correction function. Appendix B Dynamic analysis of a micro air bubble in a purely rotating flow Assume that the fluid is rotating around a vertical axis and it only has tangential velocity component (no secondary flow). Here, we examine how a micro air bubble behaves in such a fluid flow. In the following analysis, the properties of the air bubble are subscripted with a and the air bubble is assumed to be of spherical shape [22]. Given the diameter of the air bubble da, the Reynolds number is defined as ' V where Ur is the relative velocity of the air bubble to the oil flow. The drag coefficient Cd for an air bubble in various liquids was studied experimentally by Haberman and Morton [16], which confirmed theoretical solutions to the problem. The two theoretical solutions for small Reynolds number (Rea < 1) are given as follows: 24 Cd = — Stokes solution Re0 16 Cd = zz— Hadamard-Rybczynski solution Rea where the drag coefficient is expressed as where D is the drag on the bubble. Stokes solution assumes that the bubble behaves as a solid sphere, which is the limit as da —>• 0. In the Hadamard-Rybczynski solution, a 119 Appendix B. Dynamic analysis of a micro air bubble in a purely rotating flow 120 fluid sphere is assumed on which equality of the tangential velocities at the bubble-liquid interface is assumed. The drag coefficient for the fluid sphere is less than that for the solid sphere due to the internal air flow in the bubble as discussed by Kobus [22]. Using Hadamard-Rybczynski solution, the drag force applied on the bubble is expressed as Fd = Cd (\pU?) (^j =2irpda Ur. (B.45) It can be seen that the drag force depends linearly on the velocity Ur, but is acting in the opposite direction to the velocity Ur. Since both Fd and UT are vectors, we might directly write the drag force component in directional form. Let's first examine forces acting on the air bubble in the vertical z direction. The forces consist of upward buoyancy force, fluid viscous drag force acting in the opposite direction to the motion of the air bubble relative to the oil flow, and the downward weight of the air bubble. Let m be the mass of the oil of the same volume as the air bubble. According to Archimedes' principle, the buoyancy force experienced by the bubble equals to mg. The forces may be summerized as follows: the air bubble weight: W — mag, the buoyancy force: = mg, the drag force: Fd = 2itpdava where ma = pa and m = p (^<%) • In dynamic analysis, the effect of the fluid pressure on the inertia of the bubble is sub-stantial and can not be neglected. The effect of the fluid pressure is equivalent to an addition to the inertia of the bubble, half of the mass of the fluid displaced [23], namely the added mass. Therefore, the z-momentum equation: n , ( 1 \ dva mg - 2Trpdava - mag = \ma + - m l — . Appendix B. Dynamic analysis of a micro air bubble in a purely rotating flow 121 Since m >^ m a , the above equation may be simplified as dva 2Au „ ... ~ i + i r v - = 2 9 - (B-46) A solution to the above equation is given as: ( 24u > e x p "T* , (B.47) Starting from va = 0, the bubble reaches its maximum velocity as time t goes to infinity. This maximum velocity is the so-called terminal velocity, noted as Va, and is given by V. = fg. (B.48) Before going further to other momentum equations, let's examine the transient char-acteristic of the air bubble. As discussed earlier, the bubble at rest graduately gains its speed and reaches the terminal velocity at last. Time period that an air bubble takes from rest to a speed of 99 percent of the terminal velocity Va, noted as ttc, is a measure of the transient characteristics. Therefore, we have / 24i/ \ 0.99 = 1 - exp I——ttc] and t t c = ±M = 0.192^. (B.49) 24i/ v In the r direction, an oil particle experiences centrifugal force which is balanced by a positive pressure gradient. By analogy with the buoyancy effect due to gravity, an air bubble, in rotating oil flow, will experience "buoyancy force" due to the centrifugal acceleration. The "buoyancy force" equals the force of balancing an oil bubble of the same size of the air bubble at the same location, that is, mass of the oil bubble times the centrifugal acceleration (w2/r), acting towards the axis of the rotation. In addition, Appendix B. Dynamic analysis of a micro air bubble in a purely rotating flow 122 the air bubble also experiences centrifugal force and viscous drag force in the direction opposite to ua. The forces acting on the bubble in the r direction are: the centrifugal force: Fc = ma(wl/r), the centrifugal buoyancy force: F0 = m (w2/r), the drag force: Fd = 2ir p da ua. The r-momentum equation may be written as w2 w2 ( 1 \ dua ma m 27r udaua = I ma + -m ——. r r \ 2 / at (B.50) Since m ^ > ma and w > wa, the above equation may be simplified as 2w2 dua 2Av dt + (B.51) In this equation, the RHS term is not constant as in Equation B.46, but a function of r. Therefore, a simple analytical solution seems impossible. However, if the RHS term changes negligibly slowly with time for a given time duration, it may be treated as a constant in this time period. Then Eq. B.51 may be solved easily as: / 24u \2u r Similarly, the terminal velocity, Ua = exp dl w2 d2 -t (B.52) 12v r Now, let's examine the relative change of the RHS term of Eq. B.51 over time period of ttc, and let ac be w2/r. Then the relative change of ac is and Aar dac dac dr _^  _ ^ dt dr dt Appendix B. Dynamic analysis of a micro air bubble in a purely rotating flow 123 Since dr Tt=va<Ua, daCrr d(w2/r)rT 12w dw w2\ r r Aa c < -^Uattc = -L-J-Lu.tu = — Uattc. dr dr \ r dr r* J Therefore Aac _r_ (2wdw w*\ ( dj w2\ U.bd\\ _ 2.3^ (2w_dw w*\ ac ~ w2 \ r dr r2 J \ \2v r J \ 24i/ J 144i/2 \ r dr r2 J For a narrow gap annulus problem with the inner cylinder of radius r; stationary (a>; = 0) and outer cylinder of radius ra rotating at a speed of co0, the tangential velocity of the fluid may well be approximated as linear profile with respect to r. Assign w = w2/r, the RHS of the above inequality may be written as 2.3da f2(uoro-0) \ 2.3da ( 2u„ \ 2.3da ( 2u0 \ —U) I —1 ; id I = —U) I OJ I < —U) I I . 144i/2 \ r 0 - T i j 144i/2 \l-n/r0 j 144i/2 \l-n/r0) At a stationary wall, u) = 0 for no-slip condition. Obviously, the air bubble will not make a good tracer of flow. Let the minimum value of u be 0.05 u>0 (the air bubble is placed in the region where UJ > 5%^). Then we may finally have an estimated value for the relative change of the RHS term of Equation B.51, Aac 2.3^ / 0.1a>o2 \ ac 144i/2 \1 - ri/r0) ' This inequality is valid only for the assumptions applied in the derivation procedure. In the experiment, the air bubble diameter is about 0.1 mm and the kinematic vis-cosity of the oil is about 0.75 x l O - 4 m2/s. In the region of flow visualization, r; and ra are 117.5 mm and 160.0 mm respectively. The resulting value for u>0 = OJR = 5.23 rad/s, Aa, a — - < 7.92 x IO"9, ac which means that the relative change of the RHS term of Eq. B.51 over the time period for the air bubble to reach 99 percent of its terminal velocity is negligibly small. The term (2w2/r) in Eq. B.51 can be treated as a constant in solving Eq. B.51. Appendix B. Dynamic analysis of a micro air bubble in a purely rotating flow 124 In the tangential ^-direction, the bubble experiences drag force due to the relative motion of the bubble to the oil flow. In addition, there exists the Coriolis acceleration (2u/aua) when the bubble gains a radial velocity component ua. The drag force is given by Fd — 2itdap (w — wa). The 0-momentum equation: 2ndap (w - wa) + ma (2^ua^j = (ma + ^rnj dwa IT The Coriolis force term is negligible compared to the drag force during transient period. Therefore, the momentum equation may be simplified as: dwa 24i/ -(W — WN). dt <Pa { a ) A solution to the above equation may be written as: 24i / "~dl (B .53) W„ — W 1 — exp -t (B .54) The terminal velocity of the air bubble Wa is equal to the tangential velocity w of oil flow. Surprisingly, the three velocity components approaches their terminal velocities at the same rate. A p p e n d i x C A i r bubble size affected by the oi l stat ic pressure In the experimental flow visualization, air bubbles are used as tracers in the oil flow. As an air bubble moves across isobaric lines, the static pressure experienced by the bubble changes. The purpose of this section is to evaluate the size change of an air bubble for a given change of the static pressure. Assume that an air bubble has a spherical shape of diameter da and air in the air bubble is an ideal gas and the oil is isothermal (In the experiment, the temperature change is negligibly small due to very low operating speed of the runner). According to the ideal-gas law, — = constant, (C.55) Pa where pa is the air pressure in the bubble. Considering the surface tension acting on the bubble, Pa-P=4f, (C56) da where p is the oil pressure acting on the air bubble. The conservation of mass gives the air bubble mass = pa ^-irda^J = constant, or pa d\ — constant. (C.57) From Eq. C.55 and Eq. C.57, pa d\ = constant. (C.58) 125 Appendix C. Air bubble size affected by the oil static pressure 126 Combining Eq. C.56 and Eq. C.58, ——Vp) dza = constant. (C.59) Therefore, as the air bubble moves from low to high oil pressure regions, the air bubble size decreases according to Eq. C.59. In the experiments, air bubbles are placed in the region underneath the bottom of the runner. For the runner speed at OJR = 5.19 rad/s, the maximum oil pressure change in the flow visualization plane is approximately from 430 Pa (~50 mmOil) to 560 Pa (~65 mmOil) plus the atmospheric pressure. For the atmospheric pressure pa = 101,300 Pa, the maximum oil pressure change in the fluid domain of the image ranges from Pmin = 101,730 Pa to p m a x = 101,860 Pa. Assume that diameter of the air bubble changes from daX to da2 corresponding to the pressure change from Pmin to pmax- Then from Eq. C.59, Surface tension of the lubricant oil contacting with air is given as cr = 25-35 x 10~3 N/m [14]. If the surface tension a = 30 x 10~3 N/m and the air bubble diameter dai = 0.1 mm. Then a meaningful solution to da may be obtained as da2 = 0.099958 mm by solving the above equation. Therefore, the relative change of the bubble size under the given condition is: size. From the calculations, the air bubble experiences only 0.042% change in size under the experimental conditions. If the surface tension is neglected, the air pressure in the bubble is the same as the oil pressure pa = p. Since the pressure change Ap = pmax — Pmin 1 S small, from the given condition. Therefore, the assumption of a constant bubble size is validated for Appendix C. Air bubble size affected by the oil static pressure 127 Equation C.58, we may write, Ada _ lAp da 3 p Taking the extreme case for the pressure p = Pmhn the relative change of the bubble size yields: ^ x 100% = -0.043%. da Therefore, neglecting and including the surface tension effect result in about 2.3% differ-ence for this calculation. Appendix D Subroutines in Figure 4.23 Each subroutine in Figure 4.23 is listed with a brief description of its major functions. Subroutine BC: • Sets the values of appropriate variables at the wall boundaries, i.e. rigid free slip, no slip conditions, etc. • Sets the values of appropriate variables around the boundaries established by free surfaces. Subroutine DELTADJ: • Computes maximum allowable At for stability. • Adjusts A i according to number of iterations and maximum allowed for stability. Subroutine MESHSET: • Generates the computing mesh from the input data given by the user. • Evaluates all of the necessary geometric variables that are used throughout the code. Subroutine OUTPUT: • Generates output data file in a format required by TECPLOT software for plotting velocity components, pressure, stream function etc. 128 Appendix D. Subroutines in Figure 4.23 129 • Interprets computed results and produces data for plotting the shape of free surface. Subroutine PETACAL: • Determines the slope of the surface in the surface cells. • Determines the cell flag NF(ij) to indicate the interpolation neighboring cells of the surface cell. Subroutine PRESSIT: • Iterates the velocity and pressure field such that mass is conserved in each cell of the mesh. • Computes a free surface cell pressure adjustment based on the applied surface pressure; mass conservation in the surface is not iterated, it is set by application of the free surface boundary conditions. Subroutine SETUP: • Initializes constants necessary for the calculation. • Computes the initial hydrostatic pressure distribution to initialize the p(i, j) pres-sure array. • Sets up the initial velocity u(i,j) and v(i,j). • Computes the relaxation factors that are used in the pressure iteration. • Sets up obstacles. Obstacle definition, in general, must be coded by hand for each problem. Appendix D. Subroutines in Figure 4.23 130 Subroutine STREAM: • Calculates stream function based on the solution to u and v. Subroutine TIDLE: • Computes an explicit solution for r-momentum and z-momentum equations, (i.e., new values of velocities are obtained from the time n values of pressure, advective, and diffusive accelerations). These tilde values are advanced to time n + 1 in the pressure iteration. Subroutine VFCONV: • Computes the solution to the F function (Eq. 4.9). Subroutine WTHETA: • Sets the values of w at wall boundary conditions. • Sets the values of w around the boundaries established by the free surface. • Computes an explicit solution for -^momentum equation. Appendix E Free surface cell and its side ratios In the 2D SOLA-VOF code, shape of free surface in a free surface cell is denned by the following parameters: F(i,j) represents the fraction of fluid in a cell, NF(ij) indicates on which side of a cell the fluid is primarily on, TANTH(ij) shows the slope of the free surface with respect to the cell side. In order to accommodate the free surface boundary conditions for a free surface cell in solving -^momentum equation, the side ratio is defined at each side of the four faces of the cell, as such that the side ratio on a side of the four faces of the cell equals to the ratio of the wetted length on the side (by the fluid) and the length of that side of the cell. Figure E.64 illustrates a free surface cell. According to the definition, the side ratios are: j. le r Lj, ln ls / e ~ A l ; / W ~ A7 ; / e _ A 7 ; u ~A7' The procedure of calculating the side ratios from F(i,j), NF(ij) and TANTH(ij) axe as follows. First we determine which side of the cell the fluid is primarily "sitting" on according to NF(ij). There are four possible sides that the fluid can sit on. For example, if IMF(iJ) = 3 (fluid is sitting on the bottom side of the cell) and TANTH(ij) > 0, there will be three possible cases to be considered separately. Figure E.65 illustrates the three different cases. Based on the known values of F(i,j) and TANTH(ij), we determine which case the free surface cell fits in among these three cases. Then the side ratios are calculated. For 131 Appendix E. Free surface cell and its side ratios 132 j Ie I AZ f w (i.j) Ar ~* • Figure E.64: A free surface cell, other situations, the analysis is similar. The detailed procedure can be found in [32]. (a) (b) (c) Figure E.65: Three different cases for a free surface cell (NF(ij) = 3; TANTH(ij) > 0). Appendix F Dimensional analysis of shape of free surface Consider an annulus with inner and outer radii of Ri and R2 respectively as shown in Figure F.66. Assume the inner cylinder is stationary and the outer cylinder is rotating at a constant angular velocity of u>2. AR AH / \ / \ Rj / \ J \ Figure F.66: Illustration of the rotating annulus. First, let's examine the height difference AH of the free surface. Based on basic fluid mechanics, the height difference of the free surface would depend on geometry parameters AJ? and R2, the angular velocity u;2, gravitational acceleration g, fluid density p and fluid 133 Appendix F. Dimensional analysis of shape of free surface 134 viscosity p. This relation may be expressed in a symbolic form: AH = Fx (AR, R2, u>2, g, p, p) where T\ indicates a functional relationship. Using Buckingham's ir theorem, the follow-ing dimensionless group relationship may be obtained, All (AR u;22R2 _u2 AR~ 1 \ #2 ' 9 'gR2 where V\ indicates another functional relationship. The second non-dimensional variable of function V\ is a form of Froude number [37], which compares the ratio of the centrifugal acceleration to the gravitational acceleration, F r = ^ . 9 Second, let's examine the shape of the free surface. Compared to the height difference of the free surface as described earlier, the radial coordinate r will be an additional independent parameter in the functional relationship. Therefore, we may easily write, h = F2{AH,r) where T2 indicates a functional relationship. Similarly, dimensionless group relationship may be given as " = ^ , M , ^ , J ^ (F.61) AH \AR' R2 ' g ' gR\ where V2 indicates another functional relationship Appendix G Analysis on isobaric lines in a rotating annulus An annulus is formed by two vertical cylinders as shown in Figure F.66. Assume that the infinitely long annulus is filled with homogeneous Newtonian fluid. An analytical solution to the tangential velocity may be obtained for a stable laminar flow. The angular velocity u> in the ^ -direction is given by v 1 - = A + f ? - , (JL>2 T1 (G.62) where the constants A and B are obtained by applying boundary conditions: w — 0 at r = Ri and u; = o>2 at r = R2, A= 1 Rn R\ R\) The swirling fluid experiences centrifugal force in the r-direction. This force is balanced by the radial pressure gradient, ldp p dr u>2r. By integrating the above equation along a constant z line (z = zi) from r = R\ to r = r , the static pressure at point ( r , Z i ) may be written as: p(r , z i ) -p(Ruzx) = pu>\ 2 (G.63) In the vertical z-direction, the force balance is given by ldp p dz -g-(G.64) 135 Appendix G. Analysis on isobaric lines in a rotating annulus 136 By integrating the above equation along a constant r line (r = r) from z = z\ to z = z, the static pressure at point (r, z) may be expressed as: p(r, z) - p(r, zi) = -pg(z - Zi). Combining Eq. G.63 and Eq. G.65, P(r,z)-p(R1,z1) = -pgiz-^ + pu* y ( r 2 - R\) + 2ABln (^-) - y (G.65) 1 (G.66) Therefore isobaric line p = p(Ri,Zi) are given by A2 g(z - zi) = u;2 The above equation may be rewritten as z — zi w\R2 1 2 L 2 2 U2 fl2 (G.67) AR g R2AR Therefore the height difference AH of the isobaric line from the inner cylinder to the outer cylinder may be written as AH u\R2 A 2 2 (R22-Rl) + 2AB\n (|) - f AR g R2AR Now let (f> = R2/Ri, then the above equation becomes AH u2R2 AR g where A!" is a function of c6 only. The annular radius ratio coefficient K is given by, 0 3 + < £ - ^ W (G.68) (G.69) K = (G.70) 2(C6-l)(c6 2-l) It can be seen that for a given annulus (given Ri and R2), AH is proportional to u>|. A graphic expression of Eq. G.70 is shown in Figure G.67. From the figure, as the annular Appendix G. Analysis on isobaric lines in a rotating annulus 137 K 0.55 0.50 0.45 Figure G.67: The annular radius ratio coefficient K vs. the annular radius ratio <6. gap AR becomes smaller, <j> —>• 1 and K -» | . When <f> —>• oo, the diameter of the inner cylinder approaches to zero and the fluid approaches to the solid body rotation, and K —> | . Also, the maximum value Kmaji — 0.515 as <$> — 11.9. However, it should be noted that after K reaches the maximum value, AH still increases as <j> is increased for a constant value of gR, since AR is also increased as a result. Appendix H Computation of trace of an air bubble in the computed flow field Flow field in the oilwell tube annulus of the test rig is computed using the A3D-CFD code. An oil particle in the steady oil flow moves along a streamline in the secondary (r-z) plane and the oil pathline is identical to the streamline. However, if an air bubble is placed in such a flow field, the pathline of the air bubble tends to deviate from the streamline. An air bubble in stationary fluid will accelerate and reach a terminal velocity (refer to Appendix B) given as where da is the air bubble diameter and v is kinematic viscosity of the oil. The terminal velocity is in the opposite direction of the gravity g. Based on the analogy between the gravity and other forms of body forces (i.e. centrifugal force), the deviation of the tracer (air bubble) caused by the body forces may be evaluated. In the A3D-CFD code, Eulerian method is used. The computed flow field is relative to a fixed (r-z) plane. To evaluate the dynamics of an oil particle moving in the rotating oil flow, Lagrangian method has to be used. In order to use the existing computed flow field, define a reference coordinate system rotating at an angular speed u>, the same as the speed of the oil particle around the axis of the runner. Due to the axisymmetric assumption of the flow, it is known that computed secondary flow velocity (u,v) is the relative velocity to the rotating reference frame. Let V be the vector form of the relative 138 Appendix H. Computation of trace of an air bubble in the computed flow field 139 velocity. Then the absolute acceleration of the oil particle is expressed as [40] a = ^ + ^ x r + 2w xV + v x (a> x r), (H.71) dt dt or in directional forms The terms, du 2 dv °> = A ' doj a6 = r — + 2am. at civ c L ; eft ' dt1 dt are accelerations of motion of the oil particle in r, z and 6 directions respectively. The term (—u>2r) is the centrifugal acceleration and the term (2 am) is the Coriolis accelera-tion. The later two act as body forces on the oil particle. However, in the experimental condition for a># = 3.14 rad/s and v = 7.5 x 10~5 m 2 /s , the computed results give | - o;2r| < 0.1 g and |2om| < 0.013g. From earlier analysis, the effect of the Coriolis acceleration will be only about 1.3% of that from the gravitation g in terms of the terminal velocity. Compared to the primary tangential velocity, such a small effect can be neglected. Only the centrifugal and gravi-tational body forces are considered in the computation of the pathline of the air bubble. Body forces are directly acting on the mass of the particle and are proportional to the mass of the particle. These body forces are normally balanced by the pressure gradient. Significant body force on the particle is given by ab = ac + g (H.72) Appendix H. Computation of trace of an air bubble in the computed £ow field 140 where ac is centrifugal acceleration acting radially outward and g is gravitational accel-eration acting downward. If the oil particle is substituted with an air bubble of the same size at the same place in the oil, the body forces acting on the air bubble are proportional to the mass of the air bubble while the pressure gradient balancing the oil particle still exists, which results in unbalanced forces on the air bubble. Therefore, the air bubble will experience movement relative to the oil in the direction opposite to the body force since density of the air bubble is much smaller than that of the oil. This is in analogy with the buoyancy effect due to gravity. The body forces in this situation are acting in (r-z) plane only. Therefore, buoyancy effects will only directly affect velocity components in the (r-z) plane. The centrifugal acceleration (w2/r) depends on the location of an air bubble and it varies as the air bubble travels. We assume that the micro air bubble is so small that the air bubble gains its terminal velocity instantly. In other word, the transient time of the air bubble is assumed to be negligibly small. Then, the deviation of the path of an air bubble from the path of an oil particle due to acceleration of motion can be neglected. Now, the air bubble will travel mainly with the oil flow, but also with a relative velocity to the oil flow. This relative velocity is the terminal velocity due to body forces acting on the bubble. By analogy with the buoyancy effect due to gravity, the relative velocity of the air bubble to the oil flow is given by d? Then, the air bubble velocity can be obtained by adding the relative velocity Var to the velocity of the oil particle. Therefore, step by step, the pathline of an air bubble may be traced by obtaining the pathline of an oil particle (streamline) and adding the "buoyancy" velocity effect. In order to calculate the pathline of an oil particle, oil flow velocity at any point Appendix H. Computation of trace of an air bubble in the computed £ow field 141 of the flow domain needs to be known. However, the CFD computed flow field only gives velocity at discrete points. Therefore, velocity at an arbitrary point needs to be interpolated. In the CFD computation, velocity components are averaged at the center of each cell. An arbitrary point P within the rectangle defined by centers of the four finite-difference cells (i, j), (i+1), (t+1, j+1) and is shown in Figure H.68. Velocities (i+1, j+1) - 9 -6 0+1.J) 0.5[AzG)+Az(j+1)] (i. j+1) 9" " Ar P P 1 4 6-(i-j) Arp 0.5[Ar(i)+Ar(i+1)] Figure H.68: Interpolation of velocity at an arbitrary point P. V(i,j),V(i+l),V(i+l,j+l) and V(i,j-\-l) at centers of the cells are known from the CFD computation. Velocity at point P is linearly interpolated as: v * = v ( i ' j ) + 0 . 5 ^ ) +Mi +1)] W + ^ ~ v ^ + 0.25[Ar(«) + Ar(i + l)][A*(j) + A z 0 + 1)] [ V ( * ' i ) + V { ' + 1 J + " V{i + l,})-V{i,} + i)]. Appendix H. Computation of trace of an air bubble in the computed flow field 142 From the continuous flow field, represented by the above equation, pathline of an air bubble can be traced numerically, by following the procedure below. Assume that at time t an air bubble is at a known location A in the flow and its loca-tion after a time interval At is to be determined. Figure H.69 illustrates the calculation procedure of the secondary path of the air bubble during the time interval At. Velocities shown are velocities of the oil flow field. Z - r Figure H.69: New location of the air bubble after the time interval At. With an infinitesimal At, path of the air bubble can be traced by marching the air bubble since velocity of oil flow and the relative velocity of the air bubble to the oil flow are known. However, for a finite time step At, the marching technique may lead to significant errors, especially for a curved path as is in this secondary vortex flow. To reduce the numerical error, a new technique with iterative procedure is developed and used, which is described in the following three steps. Appendix H. Computation of trace of an air bubble in the computed flow field 143 First, find the secondary path of an oil particle (streamline). Assume that in a small time interval At, the oil particle travels along a small circular arc. At first, find the center of the arc and the angular speed of the traveling oil particle along the arc. The first guess of the new location is from direct marching the oil particle, which is shown as point B in Figure H.69. Velocity VB is then known. Find the intersection point B' of a line passing through point A and perpendicular to VA and another line passing through point B and perpendicular to Vg. A finite point B' implies the existence of curving path of the oil particle. Calculate averaged angular speed Then, second guess of the new location C is obtained by traveling the particle from A, at the angular speed of LOO,, around the center point B'. Repeat the above procedure until a convergent solution to the new location of the oil particle is obtained. This convergent location is considered to be the new position of the oil particle after the time interval At. The second step is to evaluate the buoyancy effect on the air bubble and find the new location of the air bubble. Assume that point C is the new location of the oil particle at time t + At. Then centrifugal acceleration of the air bubble may be approximated by where and wc are tangential velocity of the oil flow at points A and C respectively. With known body forces ac and g, the relative velocity of the air bubble to the oil particle where (= ac+g) is the vector form of the body forces. With this relative velocity, the new location of the air bubble D is obtained by moving the bubble a distance (Var At) from C. If D is very close to C compared to point A, D can be considered as the new is given by 12 v Appendix H. Computation of trace of an air bubble in the computed flow field 144 location of the air bubble. Otherwise, further iterative procedure should be performed due to the possible change of the averaged angular velocity for calculating the point C. The third step is to calculate the new circumferential location of the air bubble. Since the new location D in the secondary (r-z) plane is known. The circumferential angle traveled by the air bubble in the time interval At is given by 2 V rA rD J Advancing the time step and repeating the above procedure, pathline of the air bubble is then traced. It is worth noticing that the oil particle is an imaginary object. As the air bubble deviates from the oil particle constantly, the oil particle is always referred to the imaginary one at the location of the air bubble. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080938/manifest

Comment

Related Items