TURBULENT SWIRLING FLOW IN SHORT CYLINDRICAL CHAMBERS By ARDESHIR RIAHI B.Sc, The University of Cardiff, South Wales,U.K. 1983 M.A.Sc., The University of B r i t i s h Columbia, 1985 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (MECHANICAL ENGINEERING DEPARTMENT) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July 1990 Ardeshir Riahi, 1990 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. 1 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. The University of British Columbia Vancouver, Canada Department DE-6 (2/88) ABSTRACT The effects of aspect ratio (L/D) on the rate of decay of swirl in a cylindrical chamber were experimentally studied using the Laser-Doppler-Anemometry technique. PreUrninary measurements revealed that water should be used as working fluid; the results pertaining to air were inferred from Reynolds number similarity. The steady-state measurements revealed that a solid body type of rotation can be generated by a disc whose surface has been uniformly roughened. The effect of aspect ratio on the rate of decay of such flow field was studied in three chambers with aspect ratios in the range of interest to engine combustion. Experimental results showed a faster decay rate in the shorter chamber (i.e. smaller aspect ratio). This was attributed to the stronger swirl driven secondary flow pattern in the shorter chamber. A mathematical model describing axi-symmetric, decaying, turbulent swirling flow inside a short cylindrical chamber was also developed. The model was numerically solved, using the control-volume analysis, to provide insight on swirl decay in engines. The model validation was based on experimental observations. Turbulence parameters were represented by a two-equation turbulence model, modified for streamline curvature effects. The ad-hoc curvature modification of the standard k-e model proposed by Launder et al. and the mixing energy model developed by Saffman-Wilcox-Traci (SWT) were used to account for curvature effects. The analysis of steady flow between two long concentric cylinders, established the superiority of the latter over the former method. The SWT model was also successfully used in reproducing previously published experimental results, pertaining to decaying swirling flow field (mean velocity and turbulence intensity) in a short cylinder. The calculated turbulence intensity profile revealed that swirl promotes anisotropic turbulence. The validated numerical model was used to predict the effect of aspect ratio on the rate of decay of the flow field observed by the experimental measurements in the present study.The overall prediction of decay rate was successful, leading to the conclusion that Wilcox and Chambers model can be used in predicting the behaviour of two-dimensional transient turbulent swirling flows. T A B L E OF CONTENTS ABSTRACT i NOMENCLATURE v LIST OF TABLES viii LIST OF FIGURES ix ACKNOWLEDGEMENTS xiv CHAPTER 1. INTRODUCTION 1.1 General 1 1.2 Physical behaviour of turbulent swirling flows 2 1.3 Literature review 8 1.3.1 Swirling flows in internal combustion engines 10 1.3.2 Flows with streamline curvature 14 1.4 Purpose of present study 18 CHAPTER 2. MATHEMATICAL FORMULATION 2.1 General 20 2.2 Instantaneous conservation equations 20 2.3 Ensemble-averaged conservation equations 22 2.4 Turbulence modeling 24 2.4.1 Mixing length model 25 2.4.2 Curvature-modified mixing length model 27 2.4.3 One-equation k-1 model 28 2.4.4 Curvature-modified one-equation k-1 model 31 2.4.5 Two-equation k- e model 31 2.4.6 Curvature-modified k- e model 34 2.4.7 Two-equation e-Q model 35 2.5 Initial and boundary conditions 38 ii CHAPTER 3 NUMERICAL PROCEDURE 3.1 Introduction 45 3.2 The grid arrangement 45 3.3 Control-volume formulation 47 3.3.1 Discretization of the transport equation 52 3.3.2 False diffusion error 58 3.3.3 The bounded skewed hybrid differencing scheme 59 3.4 Implementation of boundary conditions 63 3.4.1 Inlet boundary 63 3.4.2 Outlet boundary 63 3.4.3 Axes of symmetry 64 3.4.4 Wall 64 3.5 The method of solution 69 CHAPTER 4 RESULTS OF NUMERICAL TESTS 4.1 General 71 4.2 Laminar rotating fluid near a stationary disk 71 4.3 Decay of a free laminar vortex 76 4.4 Turbulent flow between rotating cylinders 79 4.4.1 Inner cylinder rotating, outer cylinder stationary 84 4.4.2 Outer cylinder rotating, inner cylinder stationary 97 4.4.3 Conclusions 99 4.5 Decay of a turbulent swirling flow 99 CHAPTER 5 EXPERIMENTAL SET UP AND MEASUREMENTS TECHNIQUE 5.1 Introduction 114 5.2 Constant-volume chamber 115 5.3 Laser Doppler Anemometry 120 5.3.1 Introduction 120 5.3.2 Differential doppler mode 121 iii 5.3.3 Major components of LDA 127 5.3.4 Scattering particles 134 5.3.5 LDA data collection and analysis 135 5.4 Choice of the working fluid 137 CHAPTER 6 EXPERIMENTAL RESULTS AND DISCUSSION 6.1. General 140 6.2. Reynolds number similarity test 140 6.3. Steady-state measurements with smooth disc 141 6.4. Steady-state measurements with roughened disc 154 6.5. Transient measurements in short chambers 166 CHAPTER 7 CONCLUSIONS AND RECOMMENDATION 182 REFERENCES 185 APPENDIX A DERIVATION OF MEAN TRANSPORT EQUATIONS 190 APPENDIX B DERIVATION OF REYNOLDS STRESS TRANSPORT EQUATIONS 195 APPENDLXC NON-ISOTROPIC EDDY VISCOSITY MODEL 200 APPENDIX D SIMILARITY BETWEEN NORRIS-REYNOLDS AND Van-DRIEST LENGTH SCALE MODIFICATION 202 APPENDIX E EQUIVALENCY OF k-e AND c-Q. MODEL FOR ISOTROPIC, HOMOGENEOUS TURBULENCE 205 APPENDIX F WALL FUNCTIONS 208 APPENDIX G ERROR ANALYSIS 212 APPENDIX H THE EFFECT OF ASPECT RATIO ON THE MEAN VELOCITY 214 iv N O M E N C L A T U R E Ai Total of diffusive and convective flux C Convective flux c p Specific heat capacity J/kgK D Diffusive flux e Mixing energy m2/s2 F Dimensionless radial velocity G Dimensionless tangential velocity g Acceleration due to gravity m/s2 H Dimensionless axial velocity h Convective heat transfer coefficient W/m2K k Turbulence kinetic energy m2/s2 1 Length scale m P Pressure N/m2 Pk Production terms in kinetic energy equation W/m3 Pe Peclet number Pr Prandtl number R Radius of curvature m r Radial distance m Re Reynolds number Ri Richardson number So Source terms in the transport equation T Temperature K t Time s u Axial velocity m/s u + Dimensionless axial velocity u-c Friction velocity m/s V u'2 Axial turbulence intensity m2/s2 v Radial velocity m/s V Numerical cell volume m 3 v'2 Radial turbulence intensity m2/s2 w Tangential velocity m/s w + Dimensionless tangential velocity w'2 Tangential turbulence intensity m2/s2 x Axial distance m y Normal distance to the wall m Y + Dimensionless distance Z Roughness height m Greek Symbols P Bradshaw's constant e Turbulence kinetic energy dissipation rate m2/s3 P Density kg/m3 M- Viscosity N s/m2 § Arbitrary variable v Eddy viscosity m2/s x Shear stress N/m K von Karman constant <t> Diffusivity coefficient N s/m2 £ Dimensionless distance X Wavelength m CI Mixing energy dissipation rate s_1 vi Subscripts E N NE NW S SE SW W Superscripts n o East North North-East North-West South South-East South-West West New Old vii LIST OF TABLES Table 2.1. Listing of turbulence models and the calculated parameters pertaining to each model 26 Table 2.2. Production terms in the k-e and e-Q models 33 Table 2.3 Constants in the k-e and e-Q models 33 Table 3.1 Diffusion coefficients and source terms 46 Table 3.2 Discretized source terms Su and Sp 57 Table 5.1 High and low pass filter settings on the TSI counter 135 viii LIST OF FIGURES Figure 1.1 Stability of a fluid element in a curved path of motion 4 Figure 1.2. Motion of fluid element over a curved surface 6 Figure 1.3. Secondary flow pattern in a closed chamber of three aspect ratios 9 Figure 2.1. Correlation of law of the wall constant C with SR 44 Figure 3.1. A typical discretized solution domain 48 Figure 3.2. A typical scalar cell 49 Figure 3.3. Axial (u) and radial (v) velocity cells 50 Figure 3.4. Control-volume cell over which the general transport equation is integrated 51 Figure 3.5. A typical cell pertaining to skewed upwind differencing scheme 61 Figure 3.6. A typical cell adjacent to the wall 66 Figure 3.7. Scalar cell used to calculate the term ^ OVeffw) 68 Figure 4.1. Schematic of rotating fluid near a stationary disc of finite radius R 73 Figure 4.2. Variation of dimensionless axial (H), radial (F) and tangential (G) velocities with respect to the dimensionless distance £ 74 Figure 4.3. Schematic of the numerical calculation domain and boundaries pertaining to rotating fluid near a stationary disc 75 Figure 4.4. Comparison between analytical solution of rotating fluid near a stationary disc and the corresponding numerical results, using 20 by 20 grids 77 Figure 4.5. Comparison between analytical solution of rotating fluid near a stationary disc and the corresponding numerical results, using 30 by 33 grids 78 Figure 4.6. Comparison between the exact and numerical solution of decaying.laminar, free vortex after 80 ms. 80 Figure 4.7. Stability of annular flow presented by Taylor [ 1923] 82 Figure 4.8. Radial distribution of angular momentum in an annulus of two concentric cylinders with inner cylinder rotating and outer one stationary 83 Figure 4.9. Characteristic dimensions and grid arrangement pertaining to simulation of flow field in an annulus of concentric cylinders 85 Figure 4.10. Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with original form of the source term in the angular momentum equation (inner cylinder speed of 18.86 m/s) 87 ix Figure 4.11. Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with modified source term in the angular momentum equation (inner cylinder speed of 18.86 m/s) 89 Figure 4.12. Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with curvature corrections (inner cylinder speed of 18.86 m/s) 91 Figure 4.13. Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with curvature corrections (inner cylinder speed of 24.4 m/s) 92 Figure 4.14. Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with curvature corrections (inner cylinder speed of 30 m/s) 93 Figure 4.15. Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results using Wilcox and Chambers model with refined grids (inner cylinder speed of 24.4 m/s) 95 Figure 4.16. Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results using Wilcox and Chambers model with refined grids (inner cylinder speed of 30 m/s) 96 Figure 4.17. Characteristic dimensions and grid arrangement pertaining to simulation of flow field in an annulus of concentric cylinders (Taylor flow field) 98 Figure 4.18. Comparison between experimental velocity profile due to Taylor and numerical results with and without curvature corrections 100 Figure 4.19. Comparison between experimentally measured decaying velocity (by Inoue et al.) and the corresponding calculated field using k-e model with large initial dissipation rate 102 Figure 4.20. Comparison between experimentally measured turbulence intensity (by Inoue et al.) and the corresponding calculated field using k-e model with large initial dissipation rate 103 Figure 4.21. Comparison between experimentally measured decaying velocity (by Inoue et al.) and the corresponding calculated field using k-e model with zero initial dissipation rate 105 Figure 4.22. Comparison between experimentally measured turbulence intensity (by Inoue et al.) and the corresponding calculated field using k-e model with zero initial dissipation rate 106 Figure 4.23. Comparison between experimentally measured decaying velocity (by Inoue et al.) and the corresponding calculated field using k-1 model 108 Figure 4.24. Comparison between experimentally measured turbulence intensity (by Inoue et al.) and the corresponding calculated field using k-1 model 109 Figure 4.25. Comparison between experimentally measured decaying velocity (by Inoue et al.) and the corresponding calculated field using Wilcox and Chambers model 111 Figure 4.26. Comparison between experimentally measured turbulence intensity (by Inoue et al.) and the corresponding calculated field using Wilcox and Chambers model 112 Figure 5.1. Schematic of a constant volume cylindrical chamber assembly 116 Figure 5.2. Schematic of the chamber's end plate cross-section 117 Figure 5.3. Schematic of labyrinth seal details 118 Figure 5.4. Schematic of quartz glass holder plate 119 Figure 5.5. A typical arrangement of laser doppler anemometry in reference beam mode of operation 122 Figure 5.6. Characteristic dimensions of laser beams and ellipsoidal measuring volume 124 Figure 5.7. A typical LDA signal representing the combination of pedestal and doppler signal 126 Figure 5.8. A schematic of LDA set up in the differential doppler, forward scatter mode of operation 128 Figure 5.9. A typical photo-multiplier out put signal 131 Figure 5.10. Continuous and total burst mode of operation settings on the TSI counter 133 Figure 5.11. A typical frequency distribution of 5000 velocity data points collected by the LDA technique 138 Figure 6.1. Dimensionless tangential velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/2 142 Figure 6.2. Dimensionless tangential velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/2, water as working fluid 143 Figure 6.3. Dimensionless rms velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/2, water as working fluid 144 Figure 6.4. Dimensionless tangential velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 3/10, water as working fluid 145 Figure 6.5. Dimensionless rms velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 3/10, water as working fluid 147 Figure 6.6. Dimensionless tangential velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/10, water as working fluid 148 Figure 6.7. Dimensionless rms velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/10, water as working fluid 149 Figure 6.8. The comparison between experimental and calculated velocity profiles in a chamber of aspect ratio equal to 1/2, using Wilcox and Chambers model (rpm = 2090) 151 Figure 6.9. The comparison between experimental and calculated velocity profiles in a chamber of aspect ratio equal to 3/10, using Wilcox and Chambers model (rpm = 2090) 152 Figure 6.10. The comparison between experimental and calculated velocity profiles in a chamber of aspect ratio equal to 1/10, using Wilcox and Chambers model (rpm = 2090) 153 xi Figure 6.11. Dimensionless tangential velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 1/2, water as working fluid 155 Figure 6.12. Dimensionless rms velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 1/2, water as working fluid 156 Figure 6.13. Dimensionless tangential velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 3/10, water as working fluid 157 Figure 6.14. Dimensionless rms velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 3/10, water as working fluid 158 Figure 6.15. Dimensionless tangential velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 1/10, water as working fluid 159 Figure 6.16. Dimensionless rms velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 1/10, water as working fluid 160 Figure 6.17. Variation of mid-radius dimensionless tangential velocity profile across the chamber, with a rough disc rotating at 1180 rpm 162 Figure 6.18. Comparison between the calculated and the corresponding experimental dimensionless velocity profile using a rough disc in a chamber of aspect ratio equal to 1/2 163 Figure 6.19. Comparison between the calculated and the corresponding experimental dimensionless turbulence intensity profile using a rough disc in a chamber of aspect ratio equal to 1/2 164 Figure 6.20. Comparison between the calculated and the corresponding experimental dimensionless velocity profile using a rough disc in a chamber of aspect ratio equal to 1/10 165 Figure 6.21. Comparison between the calculated and the corresponding experimental dimensionless turbulence intensity profile using a rough disc in a chamber of aspect ratio equal to 1/10 167 Figure 6.22. The mean tangential velocity profile at dimensionless radial location of 0.9 in a chamber with aspect ratio equal to 1/2 at different time interval after the motor has been stopped 169 Figure 6.23. The radial profile of mean tangential velocity in a chamber with aspect ratio equal to 1/2 at different time intervals after the disc is stopped 171 Figure 6.24. The radial profile of turbulence intensity in a chamber with aspect ratio equal to 1/2 at different time intervals after the disc is stopped 172 Figure 6.25. The radial profile of mean tangential velocity in a chamber with aspect ratio equal to 1/10 at different time intervals after the disc is stopped 173 Figure 6.26. The radial profile of turbulence intensity in a chamber with aspect ratio equal to 1/10 at different time intervals after the disc is stopped 175 Figure 6.27. The comparison between experimentally measured velocity profile and the corresponding calculated profile (using Wilcox and Chambers model) after 50 ms and 100 ms decay (chamber aspect ratio = 1/2) 177 Figure 6.28. The comparison between experimentally measured turbulence intensity profile and the corresponding calculated profile (using Wilcox and Chambers model) after 50 ms and 100 ms decay (chamber aspect ratio = 1/2) 178 xii Figure 6.29. The comparison between experimentally measured velocity profile and the corresponding calculated profile (using Wilcox and Chambers model) after 50 ms and 100 ms decay (chamber aspect ratio =1/10) 179 Figure 6.30. The comparison between experimentally measured turbulence intensity profile and the corresponding calculated profile (using Wilcox and Chambers model) after 50 ms and 100 ms decay (chamber aspect ratio = 1/10) 180 xiii ACKNOWLEDGEMENT The author wishes to express his gratitude towards his supervisor Prof. P.G. Hill for all his help and encouragement throughout the course of this research. His expert advice has been invaluable in achieving the goals of present study. The author is also indebted to his co-supervisor Prof. M. Salcudean for generously offering her technical expertise in the CFD work when it was needed most. I would also like to thank Mr. R. Pierik for designing the experimental model. Special thanks goes to Tony Bessel for making the experimental rig and solving the mechanical problems encountered in the course of experimental study. I am also thankful to all the other technicians at UBC for all their help. The financial support provided by NSERC and the University fellowship provided by UBC in the later stages of the program is acknowledeged. The cooperation of the University of Toronto CRAY computer center is appreciated. I am also indebted to the NSERC grant for CRAY computer time. I would like to thank my family for supporting me throughout these five years. I would also like to thank my wife Afarin for her tireless moral support and encouragement and also her invaluable help in preparing this manuscript. I dedicate this work to her for being so understanding. xiv CHAPTER 1. INTRODUCTION 1 31.1 General Swirling flows are generally divided into three categories depending on the swirl intensity; these are weakly swirling flows, strongly swirling flows and purely swirling flows. Weakly swirling flows are identified as those whose dominant stream-wise velocity component is in the axial direction, superimposed by a weak tangential velocity component. Under such circumstances, the maximum tangential velocity is normally an order of magnitude less than the maximum axial velocity. On the other hand, strongly swirling flows are those flows whose maximum tangential velocity is of the same order of magnitude as the maximum axial velocity. Purely swirling flows can be considered as the limiting case of strongly swirling flows, in which the stream-wise velocity component is in the tangential direction with no, or minimum, axial flow. Swirling flows form a major constituent of many industrial devices; typical examples are turbo-machinery, internal combustion (IC) engines, industrial burners and cyclones. Swirl is normally generated by means of guide vanes, rotating perforated disks or tangential inlet injection. Even-though the presence of swirl in such devices is normally beneficial, it can induce some undesirable adverse effects. A characteristic example of beneficial swirl effect is the large axial pressure gradient associated with the strong swirling flow in industrial furnaces. In such devices, pressure fields are induced in order to balance the centrifugal forces due to swirl. However, as the swirl decays in the axial direction, an adverse pressure gradient is set up in that direction. If the swirl is strong enough, the forces due to this adverse axial pressure gradient exceed the forward momentum and the flow near the center changes direction; thus recirculation zones are formed. Such zones play an important role in flame stabilization. Swirling flow in internal combustion (IC) engines provide both beneficial and adverse effects. Swirl is believed to enhance turbulence which causes higher flame speed. At the same time higher turbulence promotes higher wall shear stress and consequently greater heat losses to the wall. Close to the top-dead-center portion of the engine cycle, the flow can be assumed to be 2 purely swirling, and hence no axial pressure gradient is anticipated. However, due to the nature of the flow, swirl-driven secondary flows1 are set up. Secondary flows can increase the rate of decay of mean velocity and turbulence intensity. Based on the above example, it is clear that engine designers should strive for an optimum swirl at which the benefits outweigh the disadvantages. To do so would require a detailed understanding of fluid mechanics and thermodynamics of turbulent swirling flows. Experimental and numerical investigations are the essential tools for acquiring such an understanding. The purpose of the present study is the experimental observation and numerical investigation of turbulent swirling flow in a short cylindrical chamber. The objective is to simulate conditions near the top-dead-center (TDC) portion of the cycle in internal combustion engines. The goal is to develop a mathematical model which adequately represents the fluid dynamics of swirling motion in short chambers of varying length-to-diameter ratio. The distribution of mean motion and turbulence intensity are both of concern; of special interest are curvature effects on turbulence, secondary flows and transient swirl behaviour. The experimental observations have been used in the development of the model. 1.2. Physical behaviour of turbulent swirling flows Swirling flows have been categorized, by their nature, under the general heading of flows with curved streamlines. In this retrospect, theoreticians treat swirling flows as "complex" turbulent flows2. Based on extensive review of the literature on this topic, Bradshaw [1973] concluded that curvature has substantial effects on the flow behaviour. Rayleigh [1916] was one of the pioneers in the study of the curvature effect; he postulated that curved flows can be sub-divided 1 The mechanism b e h i n d t h e g e n e r a t i o n o f s e c o n d a r y flow i s e x p l a i n e d i n t h e f o l l o w i n g s e c t i o n . 2 B r a d s h a w . P [1973] d e f i n e s c o m p l e x t u r b u l e n t f l o w s as s h e a r l a y e r s w i t h c o m p l i c a t i n g i n f l u e n c e s s u c h as e x t r a s t r a i n r a t e o r i n t e r a c t i o n w i t h o t h e r s h e a r l a y e r s . 3 into two categories, namely stable and unstable curved flows. Rayleigh's flow categorization can be best understood by considering a curved flow in an ideal fluid. Fig. 1.1 shows a fluid element at a radial location r from the center of curvature which is moving at a tangential velocity w in its circular path. At any location the centrifugal force on this element (due to the tangential velocity) is balanced by the radial pressure gradient viz: 3p _ p w^ = p (w rft, Br r r 3 (1.1) If a fluid particle which is originally at location r 0 is displaced to a new radial location r, it would conserve its angular momentum and hence would attain a velocity wp given by : (wr)p = constant (1.2) Rayleigh suggested that two situations might arise. In one case the angular momentum of the flow (wr)fl0W decreases with increasing radius and the other where the angular momentum of the flow increases with increasing radius. In the former case, as the particle is displaced in the positive radial direction, the centrifugal force experienced by it would be: p w2 = p (w rg ( i 3 ) but P (w r)£ o w < P (w rg (1.4) Substitution of equation 1.1 in 1.4 would result in : | < ^ i ,i.5) 1.1. Stability of a fluid element in a curved path of 5 So, as a fluid particle is displaced radially outwards, the centrifugal force in the new location outweighs the radial pressure gradient and the particle is displaced still further. Under such circumstances the flow is considered to be unstable. It can be shown in a similar manner that if the mean angular momentum increases with increasing radius, the force on the displaced particle due to radial pressure gradient is greater than the centrifugal force and as a result the particle is forced back to its original location. Under such circumstances the flow is said to be stable. Hence, the Rayleigh stability criterion can be expressed in terms of the radial profile of the mean flow angular momentum. If the angular momentum decreases in the radial direction, the flow is unstable, while if the angular momentum increases in the radial direction, the flow is stable. von-Karman [1934] suggested that Rayleigh's stability criterion also holds for the more general case of viscous turbulent flows. Based on the Rayleigh stability condition, stability criterion is based on the radial gradient of angular momentum. In the case of viscous flows, stability argument can be based on the rate of augmentation or reduction of turbulence intensity. This mechanism is governed by the magnitude of shearing stress which represents the rate of transfer of momentum. von-Karman suggested that the shearing stress can be expressed as x = pv'i^(wr) It can thus be concluded that the stabihty criterion can be explained in terms of gradient of angular momentum and also in terms of the normal component of velocity fluctuation v (figure 1.2). Bradshaw [1973] draws an analogy between buoyancy and curvature effects to explain the physics of curved flows. By considering a fluid element which is displaced in a stratified fluid in which the density decreases in the direction opposite to the gravity vector, it is shown by Tennekes and Lumley [1972] that the element would oscillate in a simple harmonic motion at a characteristic frequency fB where: 6 Figure 1.2. Motion of fluid element over a curved surface. 7 fB = V :f_dp p dy (1.6) This is the so called Brunt-Vaisala frequency which is the frequency of a gravity wave in a stable atmosphere. In the above equation g is the acceleration due to gravity, p is the mean density dP and — ^ is the mean density gradient in the direction opposite to the gravity vector. For an unstable stratified fluid, the density gradient is positive and the gravity wave breaks up to turbulence. Bradshaw shows that, in analogy with buoyant flows, an element in its curved path of motion would oscillate with a simple harmonic motion at a frequency of fr if it is displaced radially outwards and then released; fr is given by: It can be seen that as suggested by Rayleigh, the radial gradient of angular momentum (wr) is the governing parameter for stabihty. For negative gradient of angular momentum, the amplitude of oscillation increases and hence the flow is unstable. The reverse is true for positive gradient of angular momentum. Bradshaw suggests that, as for the stratified fluid in a gravity field, a curvature parameter can be used to quantify the stream-line curvature effects. This parameter is the so called gradient Richardson number which for buoyant flows is defined as: (1.7) Ri = -g_dp p dy (1.8) For curved flows the analogous parameter would be: Ri = (1.9) for stable flows Ri is positive and for unstable flows Ri is negative. 8 Because swirling flows have the general characteristics of curved flows, the stability criterion is applicable to these flows. The swirling flow field inside a closed chamber, which is the subject of the present study, is comprised of both stable and unstable regions. The flow near the center tends to be of solid-body rotation, which is stable, while the flow near the curved walls of the chamber tends to be unstable. Another interesting feature of the swirling flow inside a closed chamber is the swirl-driven secondary flow pattern which is illustrated in Fig. 1.3. An axial survey of mean tangential velocity within the chamber will show that the profile is very nearly uniform across the chamber, except in thin regions near the wall (shown experimentally by Dyer [1979]). This observation is the consequence of axial velocity being much smaller than the tangential velocity (u«w); the axial pressure gradient is much smaller than the radial one and hence the radial pressure gradient is essentially independent of x. Since at any radial location r, the p w^ centrifugal forceI—^— is balanced by the radial pressure gradient, it can be said that the tangential velocity w is independent of x (as observed by Dyer [1979]). However near the radial end walls of the chamber (see Fig. 1.3), the tangential velocity decreases (due to no-slip condition) and as a result there will be an imbalance between the centrifugal force and the radial pressure gradient. This causes the fluid to accelerate radially inward along the wall toward the swirl axis. For continuity this motion is compensated by a slow upward motion of fluid in the bulk of the chamber. The net result is the formation of swirl-driven secondary flow pattern shown in Fig 1.3. for three values of aspect ratio (chamber length to diameter ratio). The distribution of this secondary flow depends strongly on aspect ratio. The secondary flow in such chambers can have a direct influence on the wall shear stress and hence on the rate of decay of mean tangential velocity. It also affects the distribution of turbulence intensity . 1.3. Literature review Since the present study is motivated by the need to know more about swirl in IC engines, the literature on this topic is reviewed first. The literature on the general effect of streamline curvature on turbulence is reviewed in the second part of the present section. C i r c u m f e r e n t i a l W a l l A x i s o f R o t a t i o n (o) (b) Figure 1.3 Secondary flow partem in a closed chamber of three aspect ratios, (a) L/D = 1/10 (b) L/D = 1/5 (c) L/D = 1/2 10 13.1 Swirling flows in internal combustion engines: Understanding the role of in-cylinder fluid motion, in particular the effect of mean swirl velocity and turbulence intensity, in determining the combustion characteristics has been the focus of numerous investigations. The hostile environment of IC engines and the lack of accessibility however, make any detailed experimental measurement a formidable task. Many researchers have thus resorted to simple experimental models which closely simulate the real engine cycles or parts of a cycle; single-cylinder engines and combustion chambers are two typical examples. Arcoumanis and Whitelaw [1985] give an overview of some of these investigations; the general conclusion of these studies is that swirl has substantial effects on the flow turbulence intensity and hence indirect influence on the combustion rate. Swirl increases the rate of combustion by increasing the flame front area and also increases the rate of heat losses. Dao et al. [1973] studied the heat transfer rate at the interface of the working gas and walls of a motored single-cylinder engine. The swirl level was varied by adjusting the angle of the shrouded inlet valve. Their experimental results revealed that the peak heat flux consistently increased with swirl. An improved cycle-to-cycle repeatability of heat flux traces with increasing swirl was also reported. Davis et al. [1986] investigated the effects of swirl on the rate of combustion in a single-cylinder engine with a disc-shaped combustion chamber. Their experimental results revealed that swirl increased the turbulence intensity by as much as 20 per cent, and also increased the burning rate by up to 50 per cent. The increase in the burning rate was observed by comparing the mass fraction burned deduced from the pressure record at different swirl intensities. To predict the above-mentioned swirl influences, the turbulence quantities were calculated by Davis et al. using the standard k-e model with an extra production term due to swirl introduced in the standard model. This term was implemented in the kinetic energy equation and was accompanied by a constant which was optimized to attain best agreement with experimental results. Their numerical model predicted a 20 per cent increase in turbulence intensity when high swirl was imparted to the cylinder flow as compared with the zero swirl case. 11 Piston motion, intake process and cylinder head geometry are a few parameters which have direct and inter-related influences on the in-cylinder charge motion in an engine. Incorporating all of these parameters in a suitable numerical model would make analysis cumbersome, even with a single-cylinder engine. To isolate the effect of swirl on the turbulence intensity profile and combustion characteristics, constant-volume combustion bombs can be used. Dyer [1979] studied the combustion of methane-air mixture in a constant-volume cylindrical chamber. Swirling flow was generated through a tangential injection of mixture. Laser Doppler Anemometry was used to measure the mean velocity and turbulence intensity. Dyer reported a faster burning rate with increasing swirl. He drew this conclusion from his pressure data as well as by tracking the flame front, using laser refraction technique. Inoue et al. [1980] studied the role of swirl on combustion of propane-air mixtures inside a disc-shaped combustion chamber. Swirl intensity was varied by adjusting the inlet valve orientation. Mean velocity and turbulence intensity were measured using hot-wire anemometry in cold flow. The results revealed slower decay of turbulence as swirl intensity was increased. This observation suggested that the turbulence generated in the induction process of a real engine can be maintained through the combustion period by increasing the swirl intensity. The study of combustion in a single-cylinder engine was also carried out by the same authors. They concluded from their study that cycle-to-cycle repeatability was improved with increasing swirl. Hanson and Thomas [1984] studied the flame development in swirling flows between two long concentric cylinders. Lean propane-air mixture was used as working fluid. Swirl was generated by rotating one or both of the cylinders. Photographic records of flame development showed that rotation induced a cylindrical shape into what would otherwise have been a spherical flame. This increased the flame front area and resulted in a faster burning. Faster burning was also reflected by steeper pressure rise, as rotation was increased. However, increased rotation resulted in a steeper pressure decay which indicated higher heat losses. Wahiduzzaman [1985] studied the effect of swirl on the rate of convective heat transfer to the wall of constant-volume cylindrical chambers. Through dimensional analysis, he arrived at a correlation between the dimensionless heat transfer coefficient and the Reynolds number based on 12 the mid-radius mean velocity. The correlation suggested that the heat transfer increased with increasing swirl. He also calculated the global heat flux for different aspect ratios of the chamber. His results showed an increase in the heat flux with decreasing aspect ratio. However, his attempt to collapse the data into a single correlation was unsuccessful#M5\Afla$-parfietn^^ for large aspect ratios. He speculated that the reason was due to the axial variation of tangential velocity at higher aspect ratios; hence there was some ambiguity in defining the mid-radius velocity as the characteristic velocity for Reynolds number calculation. Recent advances in numerical capabilities have provided the ground for detailed study of in-cylinder motion. Numerical models comprising the conservation equations for momentum, energy and species concentration have been developed. The k-e turbulence closure approximation3 has been the favorite amongst the model developers (see for example Kondoh et al. [1985], Wakisaka et al. [1986], Ramos [1980], and Brandstatter et al. [1985]). In this turbulence model, k represents the total turbulence kinetic energy and e represents the rate at which the energy is dissipated by the action of viscosity. Validation of these numerical models are based on reliable experimental data. Elkotb et al. [1982] developed a numerical model of the turbulent swirling flow in a single-cylinder engine. The k-e model was used to calculate the turbulence quantities. A reasonable agreement between model prediction and experimental results was only attained after an extra term was introduced in the dissipation equation of k-e model. The addition of this empirical term was originally proposed by Launder et al. [1977] to account for curvature effects. Wang and Berry [1985] studied the effect of swirl on the rate of heat transfer between the working fluid and the coolant through cylinder walls of a single-cylinder engine under motoring conditions. The heat transfer was calculated based on the assumption of quasi-steady one-dimensional heat flow. The total heat transfer between the gas and cylinder walls included the convective and radiative heat transfer. The convective heat transfer coefficient was calculated from the following equation: 3 T u r b u l e n c e m o d e l s i n g e n e r a l a r e u s e d t o a p p r o x i m a t e t h e unknown s t a t i s t i c a l c o r r e l a t i o n t e r m s i n t h e mean c o n s e r v a t i o n e q u a t i o n s . A b r i e f i n t r o d u c t i o n o f t h e s e m o d e ls i s g i v e n i n n e x t C h a p t e r 13 h=Re°- 8 Pr where Re is the flow Reynolds number and Pr is the fluid Prandd number. They compared their model prediction of heat transfer with those reported by Woschni. The comparison revealed that their model predicted a peak heat transfer 5 times larger than Woschni's prediction which had not allowed for the presence of swirl. They concluded that swirl had a significant influence on the convective heat transfer coefficient. It should be pointed out that Wang and Berry only compared their model against a previous analytical correlation; no comparison with experimental results were reported by the authors. Wang and Ferguson [1985] utilized the KIVA code, developed by Amsden et al. [1986], to simulate the angular momentum and turbulence intensity decay reported by Dyer. The code used the k-1 turbulence model to solve for turbulence parameters. They observed that reasonable agreement between the predicted and experimental velocity field was attained only after optimizing the constant in the turbulence model. They also calculated the wall heat flux from the following equation: The constant 3 in the above equation was 200 per cent larger than that implemented in the original KIVA code. This constant adjustment was found necessary by the authors to attain good agreement between the predicted and measured wall heat flux. This discrepancy was attributed to the inability of the standard model to account for streamline curvature effects. It was explained in the previous section that the secondary flow pattern is a characteristic of swirling flows in short cylindrical chambers. The flow pattern is associated with slow radial and axial velocities in thin regions close to the chamber boundaries. Such velocities can only be resolved by allocating a large number of calculation grid points near the boundaries. The 14 computation grid used by Wang and Ferguson did not satisfy this criterion; hence it can only be assumed that the presence of secondary flow was neglected by those authors. In summary, the review of literature revealed that swirl has substantial effects on the flow turbulence intensity and hence indirect influence on the combustion rate. Swirl increases the rate of combustion by increasing the flame front area and also increases the rate of heat losses. The literature review revealed that the standard k-e turbulence closure approximation is inadequate for the numerical simulation study of swirling flows in IC engines. It was pointed out previously that swirling flows belong to the more general topic of flows with streamline curvature. The literature on this general topic is reviewed next. 132 Flows With Streamline Curvature: Swirling flows and flows over curved surfaces are examples of flows with curved streamline. There is now abundant evidence of the sensitivity of turbulence to curvature, excellent review of which have been given by Sloan et al. [1986] and Lakshminarayama [1986]. The general conclusion of these studies is that streamline curvature has substantial influence on the wall shear stress and wall heat fluxes. Some of the more relevant studies will be reviewed here. Thomman [1968] studied the heat transfer in turbulent boundary layer on curved surfaces. His experimental results suggested that the heat flux to a concave wall was 30 per cent larger than that to a flat wall. He also observed a 15 per cent reduction in the heat flux through a convex wall as compared to a flat wall. Mayle et al. [1979] observed a similar effect on the heat transfer rate. However, the curvature effect in the case of concave wall was not as drastic as that reported by Thomann. Yowakim and Kind [1988] studied the annular turbulent swirling flows. The measured wall shear stress on the concave wall of the duct was as much as 15 per cent higher than that at the convex wall. The difference between the two increased with increasing swirl. These experimental findings can be explained qualitatively in terms of the Rayleigh stability criterion which was reviewed in the previous section. It might be expected that for stable flows, mixing and turbulence would be inhibited, while for unstable flows mixing and turbulence would be augmented. The stability argument suggests that flows over concave walls are unstable; this 15 would promote mixing, increase turbulence and thus wall shear stresses. If the Reynolds analogy holds, one would also expect higher wall heat fluxes. The flows on convex surfaces on the other hand are stable so one would expect a reduction in the wall shear stresses and heat fluxes. Reasonable prediction of curvature effects would require adequate modeling of the turbulence field. The standard eddy-diffusivity-based models, such as the Prandd mixing length or the standard k-e model, have been shown to be inadequate for predicting the curvature effects (see for example Ramos [1984], Gupta et al. [1984] and Rodi [1979]). One reason for such inadequacy could be the assumption of isotropic eddy viscosity; this was shown by Lilley and Chiger [1971] in their study of swirling flows. Koosinlin and Lockwood [1974] used Reynolds stress modeling to predict the ratio of stream-wise to the cross-stream eddy diffusivity coefficients. They concluded that the ratio is near unity in the fully turbulent part of the swirling boundary layer but departs significantly from unity for free swirling jets. Ramos [1984] studied incompressible, confined swirling flows using the standard k-e model and compared his numerical results against the experimental axial and tangential velocity profiles in a model combustor reported by Vu and Gouldin [1980]. The calculated results underpredicted the maximum and minimum axial velocity peaks. Ramos attributed these discrepancies to the assumption of isotropic eddy diffusivity used in the standard k-e model. The experimental results of annular swirling flows reported by Yowakim [1985] also support this argument, though he concluded that the anisotropy was not pronounced. The second reason for the inability of standard models to predict adequately the curvature effects could be due to the extra strain rate4 caused by curvature. The extra strain rate appears naturally in the turbulence kinetic energy equation for curved flows. However, as was pointed out by Bradshaw, these extra terms do not fully represent the curvature effects. Bradshaw showed that these terms amount to about 1 per cent of the total production of energy for curvature that causes 10 per cent change in shear stress. Various modifications to the standard models to account for the effect of extra strain-rate have been proposed. Most of these modifications have been based on the idea originally proposed by Bradshaw [1973]. In his study of flows with curved streamline, 4 T e r m i n o l o g y u s e d b y Bradshaw [1973] t o e x p l a i n c u r v a t u r e e f f e c t s 16 Bradshaw used the analogy between buoyancy and curvature effects on turbulence. In so doing, he defined a Richardson number and heuristically modified the Monion-Oboukov type length scale by the following relation: l = l 0( l - P Ri) (1.10) where 1 is the length scale for curved flows and (3 is the accompanying constant whose value changes for different flows. So [1976] showed that Bradshaw's hypothesis could be derived from the Reynolds stress equations. In arriving at this conclusion, he had to invoke the assumption of local equilibrium5 and also to assume that the Richardson number was small. His calculation showed that P in Bradshaw's length scale equation took on the value of 4.5 which was of the same order of magnitude as that originally proposed by Bradshaw. The Richardson number is positive for stable flows; if its magnitude is greater than — , then P the modified length scale attains an unrealistic negative value. To avoid such a problem Bradshaw suggested the following general form of modification ( l + PRi) Castro and Bradshaw [1976] studied the turbulence structure of highly curved mixing layers. Their results showed that so long as the Richardson number was small, the linear modification of length scale was reasonable; however at larger values of Richardson number, any modification to the mixing length in terms of local Richardson number would be inadequate. 5 The c o n v e c t i o n a n d d i f f u s i o n o f t h e t u r b u l e n c e k i n e t i c e n e r g y a r e n e g l i g i b l e ; p r o d u c t i o n o f k i n e t i c e n e r g y i s b a l a n c e d by i t s r a t e o f d i s s i p a t i o n 17 Many successful applications of Bradshaw's length scale model have been reported in the literature (see for example Sharma et al. [1976]). Nevertheless there are shortcomings of this model. Firsdy, all the disadvantages associated with the Prandd mixing length model are equally applicable to any extension of it. Secondly, the value of the empirical constant may differ from flow to flow. Thirdly, a hnear modification is only at best adequate for small curvatures. Launder et al. [1977] extended Bradshaw's idea to incorporate the effect of curvature in the more general two-equation k-e model. This was achieved by introducing an extra term (a function of Richardson number) in the dissipation equation of the k-e model. Their study of flow on spinning cones and curved surfaces showed that the introduced modification improved the performance of the k-e model. However, the extra term was ad-hoc, with an empirical constant which was tuned to attain best agreement with specific experimental results. Other researchers have accounted for curvature effects by making the k-e model coefficients a function of local turbulence properties (see for example Pourahmadi and Humphrey [1983]). In the most recent studies, the trend has been to switch to higher order turbulence closure approximations for prediction of flows with curved streamline. Fu et al. [1988] resorted to Reynolds stress modeling in their prediction of strongly swirling recirculating jet flows. Jones and Pascau [1989] used the same type of modeling in their study of confined turbulent swirling flows. It has been argued that the curvature effects appear naturally in the Reynolds stress equations and thus no ad-hoc modification should be necessary. The improvement observed in these recent studies seemed to have established the superiority of Reynolds stress modeling. The so called e-Q model which was originally proposed by Saffman [1970], and later extended by Wilcox and Traci [1976], was used by Wilcox and Chambers [1977] in their study of flow over curved surfaces. Wilcox and Chambers claim that this model takes account of curvature in a natural way. They argue that, as was pointed out by von Karman, the flow stability criterion can be based on the behaviour of particles fluctuating in the direction normal to the stream-wise flow direction. Based on this observation, they favour the use of the equation which describes the transport of normal component of turbulence intensity (*) as opposed to the total kinetic energy. The extra term due to streamline curvature appears naturally in this equation and Wilcox suggests 18 that this term increases the magnitude of v'2 for unstable flows and decreases it for stable flows. He also argues that curvature gives rise to the redistribution of total kinetic energy between the normal and stream-wise intensity components, leaving the total energy relatively unaffected. It does not appear that the Saffman-Wilcox-Traci model has been used extensively in the study of flows with streamline curvature. In summary, this literature review provided evidence that streamline curvature has substantial influence on wall shear stress and wall heat fluxes. It also revealed that the standard length scale and k-e models can not take account of curvature effects. Ad-hoc modification of these standard models have proven to be satisfactory in predicting curvature effects for special cases. Though the Wilcox and Chambers model seems to be free of any empirical modification to account for curvature, it has not been extensively used by previous investigators. 1.4. Purpose of Present Study One objective of the present study was to examine the effect of swirl on turbulence and its rate of decay using the k-1 model, k-e model with curvature accounted for by Launder et al. modification, and also the Wilcox and Chambers model. The numerical model of axisymmetric, swirling flow was developed. The model was first tested by simulating both steady and transient laminar flows for which analytical results exist. These were the laminar rotating fluid near a stationary disk, and the decay of a free laminar vortex, respectively. The k-e model with the Launder et al. modification and also the Wilcox and Chambers model were then used to simulate the flow in the annuli of two concentric cylinders, either of which was rotating. The two models were lastly employed to simulate the decay of turbulent swirling flow in a cylindrical chamber for which experimental results of Inoue et al. exist. The second objective of the present study was a detailed experimental investigation of cold swirling flow in short chambers of varying length-to-diameter ratio. The aim of this work was two-fold. Firstly, the effect of aspect ratio on the rate of decay of swirl and turbulence intensity was studied. Secondly, the measurements under steady-state conditions were used to validate the numerical simulation of this flow. The validated solution was then used to completely define the 19 initial conditions for the subsequent transient calculation. The transient calculation were validated against experimental measurements of decaying swirling flow in short chambers. In this experimental study, swirl was generated by a rotating impeller disc. It was observed that swirl generated by a smooth disc was in a transitional region (between laminar and turbulent flow); the numerical model used in the present study was not able to predict transition. This difficulty was solved by using a disc whose surface was uniformly roughened to promote early transition to turbulence. The mean velocity and turbulence intensity measurements were carried out with the aid of the non-intrusive Laser-Doppler-Anemometry (LDA) technique. 20 CHAPTER 2 M A T H E M A T I C A L F O R M U L A T I O N 2.1 General The aim of the present investigation is the study of swirling flows in short chambers. Under such a condition all three components of velocity, namely axial (u), radial (v) and tangential (w) velocity are present. In this Chapter the governing equations describing unsteady, incompressible, isothermal, turbulent axisymmetric swirling flows are derived. The Navier-Stokes equations describe laminar as well as turbulent flows. A numerical solution of these equations will therefore provide the results which can fully describe the fluid flow under laminar conditions. However, turbulent flows consist of random velocity fluctuations which should be treated with statistical methods (Tennekes and Lumely [1972]). The simplest statistical method is the decomposition of all instantaneous quantities (velocity, pressure etc.) into mean and superimposed fluctuating parts (originally proposed by Reynolds in [1895]). It is assumed that the instantaneous variables satisfy the Navier-Stokes equations. The instantaneous quantities in the Navier-Stokes equations are decomposed into mean and fluctuating parts; the equations are then suitably averaged. The averaging procedure gives rise to unknown statistical correlations due to the nonlinearity of the equations. These unknown correlations are approximated in terms of quantities which are either known or can be calculated. This process of approximation is usually known as "turbulence modeling". The ideology behind the derivation of such approximations will be explained fully in the next sections. In the following sections the instantaneous conservation equations are presented first, the ensemble-averaged equations follow. Then the variety of turbulence approximations, and their extensions to incorporate the streamline curvature effects, will be reviewed. Finally the appropriate boundary and initial conditions, used in the course of numerical simulation will be described. 2.2. Instantaneous Conservation equations The governing equations for incompressible, isothermal fluid flow are written in cylindrical {x,r,o) co-ordinate systems. The flow is quasi two-dimensional because the variations in the 21 tangential direction can be neglected, but all three components of velocity are present. The differential equations for mass and momentum are (Schlichting [1979]): fa) Continuity: | ( P » ) 4 S ^ - 0 (b) Axial momentum: du du du —+ v —+ u — dt dr dx dr dx (c) Radial momentum: (2.1) (2.2) dv dv dv w 2 + v-=^ + u dt dr dx r (d) Tangential momentum: (2.3) 3w , 3w , 9w , v w ^T- + V — + U^r—+ dt dr dx r (2.4) where T is the viscous stress which for Newtonian fluid can be expressed in terms of strain-rates viz: Xrr = 2 u ^ (2.5.a) ^6Q = 2]xf ~ du x x x = 2 u ^ (2.5.b) (2.5.c) (2.5.d) 3u 3v dr dx (2.5.c) 22 Xxe = H ^ (2.5.f) The above equations describe the instantaneous conservation of dependent variables6. The instantaneous variables are decomposed into mean and fluctuating parts (as was originally proposed by Reynolds) namely: 0 = 0 + 0 ' (2.6) where 0 can be any of the dependent variables; 0 represents the corresponding ensemble-averaged value which is defined by: 0 = l im N _ > o o i=l (2.7) The use of ensemble-averaging corresponds to a laboratory situation where a large number of data points are averaged at a fixed location and a fixed time span (time-window) over a reasonably large number of cycles. The cycles in this context refer to the number of repetitions of the same measurement. The ensemble-averaging and the time-averaging are equivalent for steady mean motion. 2.3. Ensemble-averaged conservation equations The instantaneous variables in the conservation equations(2.1-2.4) are replaced by their mean and fluctuating parts. The equations are then averaged, noting that by definition 0 =0. The resulting equations will be of the form (refer to Appendix A for detail) 6 Dependent v a r i a b l e s a r e a x i a l , r a d i a l and t a n g e n t i a l v e l o c i t y c o m p o n e n t s . 23 (a) Continuity: | - ( p u ) + l 3 { i p v ) a s 0 (2-8) dx ^ ' r 3r (b) Axial momentum: | (pu) + f f ( p r u v ) + |(p«>) = - f t + 1 | [ r k - p uVJ) + | k - p ^ (2.9) (c) Radial momentum: I t ( p v ) + r aT ^ p r v ^ + Jx~ ^ p u v ) " p : T" = " ^ " r ( T e 0 ' P w ' 2 ) + ^ [ r k - P v^] + ^ -[xrx - p u'v] (2.10) (d) Tangential momentum: ir ( p w ) + T ih ( p r v w^+Ir ( p u w ^ + p y r = ^ ir t r 2 l T r f ) ~ p w ' v ' ) l (2.11) For convenience the bars over mean variables have been dropped. The averaging process has given rise to unknown parameters in the mean conservation equations. Physically, these terms represent the momentum fluxes generated by turbulent velocity fluctuations; a momentum flux can generally be thought of as a stress and hence the terminology Reynolds stress is attributed to these unknown quantities. The conservation equations in their present form are not closed. The unknown terms should be related to quantities that are either known or can be calculated, in order to close the equations. This procedure constitutes the theme of turbulence closure approximations (or turbulence modeling). These models will be reviewed next 24 2.4. Turbulence modeling The essence of turbulence modeling is to represent the unknown Reynolds stresses in terms of known parameters. There has been a variety of proposed models, the merits and shortcomings of which have been reviewed by Rodi [1982]. There are two main categories of modeling approach; one category is based on the suggestion of Boussinesq [1877] that Reynolds stresses can be represented in terms of mean strain-rates (in analogy with laminar Newtonian flows). The second category is based on the development of differential equations describing the transport of individual stresses (Launder and Spalding [1972]). Hereafter, the first category is referred to as turbulent-viscosity modeling and the second category as stress modeling. The turbulent viscosity modeling can be represented typically by the following equation: Stress = p:t (strain-rate) (2.12) where fa is the so called turbulent viscosity. Boussinesq's approximation by itself does not constitute a turbulence closure model; this is because the turbulent viscosity in the above equation, unlike its counter part laminar viscosity, is not a fluid property and is generally unknown. The task of relating fa to variables which are either known or can be calculated, is known as turbulence modeling. The major drawbacks associated with turbulent viscosity models are two-fold. Firstly such models are based on the assumption of isotropic turbulent viscosity i.e. the turbulent viscosity is the same for all components of stress. This might not be true for swirling flows as was pointed out by Lilly [1973]. The second shortcoming of turbulent viscosity modeling arises in the cases where the surfaces of maximum velocity and zero shear stress do not coincide. This phenomenon was reported by Wattendorf [1935] in his study of flow in curved channels; such a phenomenon can not be incorporated in the Boussinesq type of turbulent viscosity modeling. An alternative approach to the turbulent viscosity concept is to calculate the individual Reynolds stresses from their respective transport equations which can be obtained directly from 25 mean momentum equations (refer to Appendix B). These transport equations in turn contain further unknown higher order statistical correlations which have to be modeled in terms of known mean parameters (see Hanjalic and Launder [1972], Launder et al. [1975]). Rodi's [1982] review of these higher order turbulence modeling suggests that there is still much work to be done to attain a complete and accurate representation of statistical correlations. A complete solution of the mean momentum equations is obtained after the Reynolds stresses are calculated. This necessitates the solution of six partial differential equations which describe the transport of individual Reynolds stresses. The turbulence model incorporated in the present study is based on turbulent viscosity approach. The curvature effects are incorporated in two ways; one way is through implementation of empirical modifications in the standard models and the second way is by the use of a model which incorporates the effects of curvature in a natural way (i.e. based on physical argument). These methods will be reviewed in the following sections. Table 2.1 summarizes different turbulent models and the quantities that each model calculates. 2.4.1. Mixing length model The simplest turbulence model, based on turbulent viscosity concept, is that proposed by Prandd [1925]. This is the so called mixing length model. In analogy with the kinetic theory of gases7, Prandtl defined the turbulent viscosity in terms of fluid density p a length scale l m and a characteristic velocity vt viz Ut = plmv t (2.13) The calculation of the turbulent viscosity requires prior knowledge of l m and v t . The turbulent viscosity is a local variable whose magnitude vary from one location to the other. Prandd suggested that the characteristic velocity vt should be represented by: ^ M- — 1 P lfree path vmolecular 26 Turbulence Model uV wv wu v'2 W'2 k Mixing length model D k-1 model I I I I I I D k-e model I I I I I I D e-Q model I I I I D I Reynolds stress D D D D D D D D : Direct calculation I : Indirect calculation Table 2.1 Listing of turbulence models and the calculated parameters pertaining to each model 27 V t = p l n 9u dy (2.14) The mixing length is normally represented in terms of some characteristic length. In the case of wall bounded flows, Prandd suggested that l m is proportional to the distance from the wall, viz l m = K y (2.15) where K is the Von-Karman constant. Near the boundary, where the laminar viscosity effects can not be neglected, the Prandd mixing length may be modified by Van Driest [1956] modification, viz: l m = K y 1 - exp • p y u x (2.16) where ux is the wall friction velocity and is equal to 'y . I t can be seen that at large distances8 away from the wall, equation (2.16) yields Prandd's mixing length. Mixing length modeling has been used extensively for simple flows (e.g. boundary layer flows). One major drawback of mixing length modeling is its inability to account for convection and diffusion of turbulence. The other major drawback is that in its standard form, the mixing length can not account for curvature effects. Several modifications to the standard mixing length model have been proposed to accommodate curvature effects. These will be reviewed next. 2.4.2. Curvature-modified mixing-length models The literature review in the previous Chapter, confirms that streamline curvature has pronounced influence on turbulence. Based on the Rayleigh stability criterion, it was shown that for stable flows turbulence is inhibited, while for unstable flows turbulence is augmented. Stable flows are those in which the radial gradient of angular momentum is positive while in the case of unstable flows the gradient is negative. 8 A t d i m e n s i o n l e s s d i s t a n c e Y + = 100 t h e e x t r a e x p o n e n t i a l t e r m c a u s e s 2 p e r c e n t change i n m i x i n g l e n g t h . 28 It was explained in the first Chapter that Bradshaw used the analogy between curvature and buoyancy effects to describe curvature effects. Two major drawbacks are associated with Bradshaw's modification. Firstly the linear modification is only applicable for small Richardson numbers. Secondly the empirical constant P has been found to be flow dependent (refer to Bradshaw [1973], Koosinlin et al. [1974], Mayle et al. [1979]). The mixing length model assumes an isotropic length scale (and hence isotropic turbulent viscosity). It was suggested by Lilley [1973] and Yowakim and Kind [1985].that isotropic eddy viscosity assumption might not be justified for swirling flows. To avoid the use of the isotropic eddy viscosity assumption for swirling/rotating flows, some investigators have suggested other approaches. Lilly [1973], Koosinlin and Lockwood [1974] favour the use of two different turbulent viscosities for such flows. The r-x viscosity is defined (analogous to Prandtl's suggestion) in terms of the mean resultant strain-rate, while the r-9 viscosity is defined through the use of viscosity ratio fy-e , viz: #r6 The viscosity ratio is usually represented in terms of Richardson number which brings Bradshaw's suggestion into consideration (see Koosinlin and Lockwood [1974]). Kind et al [1989], in their study of annular swirling flows, used two different length scales (and hence two eddy viscosities) based on the suggestion of Galbraith et al [1977] (refer to Appendix C). Kind et al. study suggested that the ratio of the two viscosities is near unity close to the boundary. In the present study, the merits of Bradshaw's modification against more advanced models was investigated; no attempt was made to introduce non-isotropic viscosity. It was found that Bradshaw's suggestion produced reasonable results, bearing in mind its simplicity and restrictions. 2.4.3 One-equation k-1 model As pointed out by Prandd, one of the shortcomings of the mixing length concept is its inability to take account of convective and diffusive nature of turbulence. By defining the turbulence viscosity in terms of velocity fluctuation (vt) and then relating v t to the local velocity 29 gradient, the influences of downstream and upstream happenings are neglected. In the case of recirculatory flow this can be a serious deficiency of the model. Kolmogorov [1942] and later Prandd [1945] suggested that the characteristic velocity (vt) should be ideally defined in terms of a turbulence property rather than a mean flow property. The square root of the total turbulence kinetic energy k was suggested as the characteristic velocity, k being defined as: in which u'2 , v'2 and w'2 are the energy components (or Reynolds stress components) in the axial, radial and tangential directions respectively; with this as characteristic velocity the turbulent viscosity becomes: in which 1 is a length scale taken to be proportional to Prandtl mixing length and is obtained through some algebraic relation; k is obtained by solving the transport equation for turbulent kinetic energy.. The transport equation is obtained by adding the transport equations for u'2 , v'2 and w'2 which in turn are obtained from the Navier-Stokes equations (see Appendix B); the general form of this equation is : (2.17) pt = P 1VF (2.18) Dk = Dt a) =- u=--lpu k + up, / 9k -r —\" y ^ - P vk-vpj ( I D ,'2 Y. r (in) -pV (u,v,w) (2.19) 30 Term I in the above equation accounts for convection of energy by the mean flow. The second group of terms (group LT) are diffusive terms because their integral for the whole flow are zero. Hence, these terms do not contribute to the total production or dissipation of energy; they merely act to redistribute the energy from one location to the other. The third group of terms are referred to as production terms; the reasons are two-fold. Firsdy they are of the form stress times strain-rate which can be interpreted as rate of work done and hence rate of production of energy. Secondly these terms appear in an identical form, except of opposite sign, in the mean kinetic energy equation; thus they represent the rate at which energy is extracted from the mean flow to enhance turbulence. This energy transfer from the mean flow is cascaded down to smaller scales of turbulence and is dissipated by viscosity at the smallest scale. The dissipation rate is represented by the last term in the turbulence kinetic energy equation. There are a number of unknown terms in equation (2.19) which can be approximated (based on their physical interpretations) in terms of kinetic energy and turbulence length scale, (refer to Launder and Spalding [1972]). The final form of the equation is : where [Leff = p t + M-ian -hi the above equation, represents the summation of all production terms and is listed in table 2.2. The constant c^ in the above equation is equal to 0.09 (refer to Appendix B). The length scale in the above equation is normally prescribed. Norris and Reynolds [1975] suggested the following definition of length scale: d [ dkl i d r -\. 1 _ I (2.20) 1 = 6.41 min (icy , X d) (2.21) in which A, = 0.085 and 5 is the boundary layer thickness. They also suggested that the following relation should be used to define the turbulent viscosity 31 ut = p V k " ^ (2.22) where D is defined as D = 1 + UL2 7 Pl (2.23) M-lam It can be shown (refer to Appendix D) that the term D effectively modifies the length scale near the wall in much the same way as the Van-Driest modification of the mixing length. In the present study, the one equation k-1 model with Norris and Reynolds definition of length scale was used. 2.4.4 Curvature-modified one-equation k-1 model By comparing the turbulence kinetic energy transport equation for curved flows and the corresponding equation for linear flows, it can be seen that in the case of former there is an extra production term, namely P * v y • However, this term does not fully account for large effects of curvature, as was pointed out by Bradshaw. Extra empirical terms can be added to the standard turbulence kinetic energy equation to account for such large effects. However, since the equation is exact (it is derived from Navier-Stokes equation), it might be more feasible to account for curvature through modification of the length scale 1. Bradshaw's modification of the mixing length can be equally applicable to the dissipation length scale of the k-1 model. Indeed this modification was implemented in the present investigation, purely for comparison purposes. 2.4.5. Two-equation k- r model The one-equation k-1 model can be successfully applied to predict flows in which prescribing a feasible length scale distribution is not difficult. However, situations might arise when the nature of the flow makes any prior knowledge of length scale unlikely; recirculatory flows and strongly accelerated flows are typical examples of such instances. 32 An alternative approach is the calculation of length scale distribution (or a combination of k and 1) through the use of an appropriate transport equation. This task forms the basis of the so called two-equation turbulence modeling. The so called k-e model, developed by Launder and Spalding [1972], is one example of two-equation models. In this model the transport equation for k i .5 the combination — j — is used as the second equation; this combination is the so called dissipation of kinetic energy which is represented by : e = C n Y - ( 2 - 2 4 ) here c^ is a constant of proportionality whose magnitude is equal to 0.09 (Launder and Spalding [1972]). Based on the Prandtl-Kolmogrov suggestion the following definition of turbulent viscosity is adopted: p t = p Vk~ 1 = c ^ p fe2- (2.25) Equation (2.25) is obtained after substituting for the length scale from equation (2.24) in equation (2.18). The transport equation of dissipation rate ( e ) is direcdy obtained from Navier-Stokes equations (refer to Harlow and Nakayama [1968]). However the equation contains a number of unknown correlation terms which should be approximated in terms of known quantities. One set of modeling of these terms have been given by Launder and Spalding [1972] which would result in the following equation: K Dt dx M-eff de o> dx r dr \ a e / dr + CiC^ p k P k - c2 p e2 . (2.26) In the above equation, Pk represents the production of kinetic energy and is given in table 2.2. The constants ci , C2 , o~e are listed in table 2.3; these constants have been obtained based on experimental observation of decay of grid generated turbulence and near-wall turbulent flows (refer to Launder and Spalding [1972]). The constants E and K are used in the implementation of the boundary conditions (refer to page 40). 33 Model Production terms P k-e model [9u dvl2 e-Q model 3w pu 9v| j 8x |3r dxj J Table 2.2 Production terms in the k-e and e-Q models. k-e model C l C2 Ok E K 1.44 1.92 1 K 2 (c2-ci)VcfT 9 0.4187 0.09 e-Q model P Ce a K a* 3 20 0.09 2 2 1 3 0.4187 0.3 Table 2.3 Constants in the k-e and e-Q. models. 34 Solution of the two transport equations for k and £, completely define the turbulence parameters which can subsequendy be used to close the Navier-Stokes equations. 2.4.6 Curvahire-moHified k-emodel In analogy with Bradshaw's modification of length scale, Launder et al. [1977] suggested the implementation of an extra Richardson number dependent term in the dissipation equation of k-e model. The reason for implementing this term in the e equation is that it seems to be an obvious place to incorporate the curvature effects on turbulence length scale distribution9. The proposed term is of the form extra term = cc Ri t (2.27) where cc is an empirical constant whose recommended value by Launder et al. is 0.2; Ri t is the turbulent Richardson number defined as : (2.28) By comparing the above definition of Richardson number with that given by equation (1.9), it can be seen that the difference in the two definition is the turbulence time scale. This time scale is represented by — in the Launder et al. definition; they feel that this time scale is superior to the time scale given in terms of mean flow since it is defined in terms of turbulence parameters. The extra term is implemented in the dissipation equation in the following manner P ° f = c 2 p £ ( l - c c R i t ) (2.29) 9 This i s because e = ^ y— . So any effect on the length scale can be incorporated in the e equation. 35 For positive Richardson number (stable flows), the extra term tends to increase the rate of dissipation of turbulence kinetic energy and hence turbulence is inhibited; this agrees with the Rayleigh stability criterion. The reverse is true for negative Richardson number (unstable flows). Launder et al [1977] successfully applied their modification of the k-e model to predict boundary layer flows on spinning and curved surfaces. However, one major disadvantage of this type of modification is the need for the case-adjusted empirical constant. Priddin [1975] had to reduce the above cited value of cc by as much as 50 per cent in his study of flow over curved wall with strong blowing through the wall. 2.4.7 Two-eqiintinn e-Q morid The momentum transfer process in a thin turbulent shear layer is governed by the behaviour of motion perpendicular to shear layer. Motivated by this observation, Saffman [1970] developed a two-equation turbulence model with the characteristic turbulence energy represented by the normal component of intensity rather than total kinetic energy (k). The Saffman model was later modified by Wilcox and Traci [1976] to make the model more universal (in the absence of curvature). In this model the mixing energy e is defined as : The proportionality constant ^ , has been chosen so that the mixing energy and turbulence kinetic energy reduce to the same value in the near wall region of a flat-plate boundary layer. In this region The rate of dissipation of mixing energy is represented by Q. Wilcox suggested the following definition of Cl: (2.30) the total kinetic energy is equal to j v , as shown by the experimental results of Klebanoff [1955]. 3u U_y (2.31) r'2 36 where P is a constant whose suggested value by Wilcox and Traci is equal to 0.09. Wilcox arrived at this definition of dissipation rate by comparing the dissipation term in the mixing energy equation with that appearing in the exact equation for the normal component of the Reynolds stress ( * ) . Wilcox and Chambers [1977] went one step further and suggested that the same model can be applied to curved flow analysis without any further empirical modification. This is obviously superior to any ad-hoc modification of the standard models to account for curvature effects. The physical reasoning for such a claim has been given by Wilcox and Chambers [1977]; this is briefly reviewed here. The transport equations for u'2 , v'2 and w'2 components of Reynolds stresses for flows over a curved surface with radius of curvature R can be written in the following general form : p D ^ + 2 p ^ 7 7 f f i . = . 2 p ^ v | ^ + ( 2 . 3 2 . a ) P ¥ " 2 P W ' V ' R = 2 P W ' V ' R + ( 2 - 3 2 , B ) = ( 2 . 3 2 . C ) The curvature related terms are of two origins. Those appearing on the left hand side of equations stem from body forces and those appearing on the right hand side are part of production terms. For flows over convex surfaces, the radius of curvature R is positive. Hence, it can be seen that the body forces tend to increase the w'2 component, decrease the v'2 component and leave u'2 unaffected. However, as was suggested by Von-Karman, the normal velocity fluctuations are responsible for mixing and hence it can be concluded that in flows over convex surfaces, mixing (and hence turbulence) is inhibited. This is in accord with the general conclusion drawn from experimental observations. The situation is exactly reversed in the case of flows over concave 37 surfaces where R is negative and body forces tend to reduce w'2 , increase v'2 but still leave u'2 unaffected. The interesting feature is that upon summation of all three equations, which gives the turbulence kinetic energy equation, the body forces cancel out and the only curvature related term is components, leaving the total energy minutely affected. This is in accord with Bradshaw's conclusion that the extra production term (due to curvature) which appears in the turbulence kinetic energy equation does not account for substantial effects of curvature. Because the curvature effects appear naturally in the v'2 equation and also because normal fluctuating velocities are responsible for mixing, Wilcox and Chambers conclude that curvature effects can best be represented by e-Q, model. The following two transport equations have been used in the present study (based on Wilcox and Chambers equation): The universal constants appearing in these equations have been given in table 2.3. These constants have been obtained based on study of wall equilibrium layer and other observed properties of turbulent flows. In this model the turbulent viscosity u.t is defined as: 2 p wV Y . This suggests that curvature tends to redistribute the energy between different energy (2.33) (2.34) p.t = p- e- (2.35) It can be shown that e-Q and k-e models reduce to the same form for homogeneous, isotropic turbulence field (refer to Appendix E). 38 2.5. Tnitial and houndarv conditions The solution of the governing partial differential conservation equations is obtained based on specified initial conditions and also conditions applied at the boundaries of the calculation domain. In the case of steady flow, the solution within the domain is only influenced by the specified boundary conditions. In the case of transient flow, the solution depends on both initial and boundary conditions. However, in some cases it might not be possible to specify the initial profiles of all dependent variables. In the present investigation, for example, the situation arose where the initial profiles for axial and radial velocity components were not available; one way to deal with this problem is to start with zero initial profile for these velocity components and check that the solution at later times is not significantly affected by this lack of initial information. The same concern is also valid for initial kinetic energy and its dissipation rate profiles. The effect of arbitrary specification of these profiles will be discussed in Chapter 4. Four types of boundary conditions were encountered in the present investigation; these are explained next. (a) Inlet boundary condition: This type of boundary condition specification is only necessary when there is a flow across the inlet boundary and into the solution domain. This situation was encountered in the numerical simulation of laminar rotating fluid near a stationary disk. This case, whose exact analytical solution is available, was chosen to test the numerical code performance. The profiles of all 3 velocity components, at the inlet boundary were specified based on analytical solution. (b) Outlet boundary condition: This type of boundary condition is necessary when there is a net flow across the exit plane and out of the solution domain . The standard implementation of the outlet condition is to set the normal gradient of all dependent variables equal to zero except for the velocity component normal 39 to the boundary. This component of velocity is adjusted until the overall mass conservation is satisfied. (c) Axis of symmetry: At the axes of symmetry, the normal gradient of all dependent variables, except the radial velocity component is set to zero. Since there is no flow across the symmetry line, the radial velocity itself would be zero. (d) Wall boundary condition: The walls are treated as impermeable surfaces and the usual no-slip condition should be applied. Near the wall, the velocity gradient is steep. An accurate representation of these profiles would require substantial number of grid points near the surface. Also the turbulence models (k-1 and the two-equation models) are generally restricted to high Reynolds number conditions, where the effects of laminar viscosity can be neglected. Very close to the surface, in the laminar sublayer region, this assumption is no longer valid. The turbulence models should then incorporate the effects of laminar viscosity; an example of this is the Van-Driest modification of mixing length. Jones and Launder [1972] have devised an extended k-e model for low-Reynolds number flows which is equally applicable in the laminar-sublayer region. However, these modifications do not alleviate the problem of large number of grid points required to resolve steep near-wall variations. An alternative approach is to use the,so called "wall functions" (described fully by Launder and Spalding [1974] ). The essence of this method is to bridge the laminar sub-layer region by locating the first calculation grid point at a location sufficiendy remote from the wall where the flow is completely turbulent. This location is determined in terms of dimensionless distance from the wall defined by: Y - P - ^ l (2.36) Ulan ' According to Hinze [1959], the near wall flow region is divided into 3 distinct layers, the thickness of each of which is given in terms of Y + . The laminar sublayer region extends from the 40 wall up to Y + of about 5. The fully turbulent region, where the turbulent stresses completely overwhelm their laminar counter parts, extends from 30 < Y + < 200. In this region the stress can be assumed to be constant and equal to that prevailing at the wall (the so called one dimensional, steady, Couette flow). The region in between these two layers is called transition or buffer region. In the fully turbulent region, the velocity is characterized by the dimensionless distance Y + in the following manner (Schlichting [1979]): u + = ^ l n Y + + B (B = 5.5) (2.37) where u + = .In the laminar sub-layer region the following simple relation holds (Schlichting Ux [1979]): u + = Y + (2.38) In the case of swirling flows and flow over curved surfaces, different authors have proposed modified logarithmic profile (see for example Wilcox and Chambers[1977], Koyama et al [1979], Lawn and Elliott [1972], Kind et al [1988]). However, the experimental results of Kind et al. revealed that logarithmic profile is valid for tangential velocity up to Y + of 200. In the present study, the logarithmic profile was used for both tangential and axial velocities. These profiles are used to calculate the wall shear stresses. The axial and tangential components of shear stress can be calculated from (refer to Appendix F): (TwaU)x = p c / 2 5 V F n g ^ (2.39) (Twaii)e = p c ^ 2 5 V ^ i - ^ (2.40) The constants E and K are listed in table 2.3.In the laminar sublayer region, the following relations are used: (Twall)* = H £ 41 (2.41.a) Ccwaii)e = H y (2.41.b) These stresses are then used to calculate the generation terms in the kinetic energy equation. The dissipation term in the kinetic energy equation is the integerated value of e over a volume extending from the wall up to the point under consideration viz (refer to Appendix F): The boundary condition used for dissipation equation is based on equilibrium condition which is achieved by forcing the dissipation rate to take on the following value (refer to Appendix F): (e) Rough wall It was explained in Chapter 1 that in order to increase the turbulence intensity of the flow generated by a rotating-disc in a cylindrical chamber, it was found necessary to roughen the disc surface. Saffman and Wilcox [1974] suggested that the law of the wall can be written in the following genral form: (2.42) (2.43) = ^ ln Y + + C (2.44) K where C is a function of wall roughness height. The function SR is defined as: 42 (2.45) where Z + is the dimensionless roughness height, defined by the following equation: Z + = ^ (2.46) v where Z is the roughness height. The variation of C with SR cited by Saffman and Wilcox is shown in Fig. 2.1. Depending on the roughness height, a value of SR and hence C can be obtained. In the present study, the value of C was found by trial and error to be -8.5. The procedure was to vary the constant until adequate agreement between the experimental and numerical results was attained. The roughness height associated with this constant was found to be 1.54 mm. The roughness height on the disc used in the experimental study was also measured to be 1.5 mm. Schlichting [1979] suggested that for a completely rough pipe, the law of the wall can be written as: w+ = 5.75 log Y + + Di (2.47) where is calculated from: Di = 8.5 - 5.75 log ^ (2.48) The constant D : associated with the roughness height of 1.54 mm was calculted from the above equation to be -8.75. It can be seen that the two constants Dj and C obtained from two different sources agree to with in 3 per cent Hence it was concluded that in the case of rough wall numerical calculation, the value of the universal constant used in the law of the wall should be equal to -8.5. 43 Figure 2.1 Correlation of law of the wall constant C with SR 44 The value of constant C lies outside the range of the curve shown in Figure 2.1. The value was obtained by extrapolating the curve up to the SR value which corresponded to the roughness height in the present experimental set up. 45 CHAPTER 3 N U M E R I C A L PROCEDURE 3.1 Introduction The governing equations describing the conservation of mass, momentum and turbulence quantities for an axi-symmetric, isothermal, incompressible flow were presented in the previous Chapter. These equations (except for continuity equation) can be written in the following general form: where (]) is any of the dependent variables, T$ is the diffusivity constant and S<j> is the source term pertaining to each equation. The diffusivity coefficient and the source terms for each variable are listed in table 3.1 .It should be pointed out that because the flow is axi-symmetric, the tangential component of momentum equation can be rearranged so that the dependent variable 4> is the angular momentum (wr) rather than tangential velocity w; this is necessary to transform the tangential momentum equation to the form given by equation 3.1. The numerical procedure adopted in the present study to solve the general transport equation is based on the control-volume analysis, described by Patankar [1980]. The large scale numerical code TEACH II, developed by Gosman group at Imperial college, incorporates the control-volume discretization technique. This code was used to solve the transport equations described in the previous Chapter. In this Chapter the grid arrangement is introduced first, followed by an overview of control-volume analysis. The difficulties associated with the discretization method and the suggested remedies will also be described. The numerical representation of the boundary conditions will be reviewed in the last section. 3.2 The grid arrangement The domain of calculation is discretized into a finite number of four-sided control-volumes (or cells). The size of these cells might vary from one location to the other. Normally, a large number of these cells are concentrated at locations where steep gradients of velocity and/or 46 Conserved property <t> r * Axial momentum u Heff 3p 9 / 3u\ i 3/ 9v\ Radial momentum V Heff 3p 3/ 3u\ i 3/ 3v\ pw2 Angular momentum wr M-eff -fr^lrpeffw) Turbulence K.E. k Me£f 0\ Pk-pe Dissipation rate e Heff OP c i ^ P k - c 2 p ^ Mixing energy e Heff Oe. P e -pP*Qe Dissipation rate Q 2 M-eff Table 3.1 Diffusion coefficients and source terms 47 turbulence parameters are anticipated (for example near walls). Fig. 3.1 shows a typical solution domain which has been discretized in this manner. The grid points which lie within these four-sided cells are the locations at which each of the variables under consideration are calculated. Following Patankar's suggestion, the grid points at which the scalar quantities (such as pressure, angular momentum etc.) are calculated are located at the center of the so-called scalar cells. Fig. 3.2 represents a typical scalar cell with grid point P (or node P) located at its center and the surrounding nodes designated by E, W, N and S located in a compass fashion. Velocity components u and v are calculated at the grid points located at the center of the faces of the scalar cells; these points are represented by small letters e,w,n and s. The control-volumes which enclose the velocity grid points are called velocity cells; a typical velocity cell is shown in Fig. 3.3. This type of arrangement is normally referred to as staggered grid arrangement. The reasons for adopting such an arrangement are two-fold. Firstly the convective fluxes through the faces of a scalar cell are based on the velocities across the cell boundaries. With a staggered grid arrangement, the velocities at these locations are readily available. Secondly, the pressure difference between the two adjacent scalar grid points becomes the natural driving force for the velocity component between these points. If on the other hand, the pressure and the velocity were to be calculated at the same grid point, a physically unrealistic solution might arise (Patankar [1980]). Additional fictitious cells are located at the boundaries to impose the boundary conditions. The faces of these cells coincide with the actual boundary to facilitate an accurate implementation of fluxes. 3.3. Control-volume formulation The finite volume form of the general transport equation (equation 3.1) is derived by integration of this equation over a finite-volume such as that shown in Fig. 3.4. The integration can be expressed as: Figure 3.1 A typical discretized solution domain. 49 Figure 3.2 A typical scalar cell. t N 1 ^ 1 i 4 ( b) 1 Figure 3.3 Axial (u) and radial (v) velocity cells . (a) u velocity cell (b) v velocity cell 51 Figure 3.4 Control-volume cell over which the general transport equation is integrated. 52 S*dV = 0(3.2) The second and third volume integrals in the above equation can be transformed into surface integrals using Gauss theorem. The equation after some re-arrangement would be: The first integral represents the transient term. The second, third, forth and fifth integrals represent the total convective and diffusive fluxes across the cell boundaries. The last integral represents the summation of all the other terms in each of the individual transport equation. In the following section the discretization of the above equation is explained. 3.3.1 Discretization of the transport equation: In the present study, a fully implicit time discretization technique was used to represent the transient terms. The first integral in equation 3.3 is thus represented as: (3.3) (3.4) 53 In the above equation, the subscript 'P' denotes the nodal point at which the variable is calculated. The superscripts 'n' and 'o' denote the new and old values of the variable under consideration. For convenience the superscript 'n' is omitted hereafter. The second, third, forth and the fifth integrals represent the total convective and diffusive fluxes across the boundaries of the control-volume. As an example, the discretization of the second integral is considered in detail. The integral can be expressed as: The subscripts 'e' and 'P' in the above expression denote the locations at which the variables are specified. It should be pointed out that the magnitude of a variable at location 'e' (i.e. the cell interface) is evaluated by assuming a linear variation of that variable between the closest neighbouring points; for example : in the case of a uniform grid spacing, fe = 1/2; for simplicity a uniform grid spacing arrangement is assumed in the following derivation. The simplest discretization scheme that can be adopted to represent the convective and diffusive fluxes, is the so-called 'Central Differencing (CD)' scheme; in this scheme a linear variation of variable <}> between the nodes E, P is assumed. The CD scheme would result in the following discretized form of convective and diffusive fluxes : Tc>e = fe T(|)E + (1 " fe) T<t>p (3.5) where fe is a weighting factor defined by : fe = 0.5 Ax p 6xn (3.6) 54 (pu<)))e = convective flux = C e = (pu)e ^ E * (3.7) r * ^ | = diffusive flux = D e = r e ^ ^ (3.8) I 3 x / e (5x)e Patankar [1980] shows that the CD scheme can be satisfactorily applied for diffusive dominated flows (low Reynolds number flows). However, in situations where convective fluxes are dominant (high Reynolds number flows), CD scheme gives rise to numerical instabilities. To remedy the difficulty, Runchal and Wolfshtein [1969] suggested the use of the so-called 'Upwind Differencing (UD)' scheme in discretizing the convective fluxes. The UD scheme amounts to representing the value of <j) at any cell boundary by the value of <|> at the upstream grid point, depending on the direction of the flow, viz: ()>e = §p for u>0 <t>e = (J)E for u<0 (3.9) hence, the convective flux can be represented by: (pu<J>)e = max[(pu)e , 0] <j>p - max [-(pu)e , 0] <|>E (3.10) The symbol 'max' in the above equation, denotes the greater of the two enclosed parameters. It can be easily seen that for u>0 in the positive x direction, the above expression reduces to : (pu<t>)e = (pu)e <}>p (3.11) It should be emphasized that CD scheme is second order accurate10 in x, while UD scheme is first order accurate (this point will be clarified in the next section); hence, it is advantageous to 1 0 The truncation error in Taylor series expansion is O (fix)2 55 use the CD scheme whenever possible. Parankar [1980] shows that from the exact solution of the steady, one-dimensional, convection-diffusion problem, the dividing line between the use of CD and UD scheme can be defined in terms of the cell Peclet number, the cell Peclet number is defined as: Pe = P ^ (3.12) The so-called 'Hybrid Differencing (HD)' scheme, developed by Spalding [1972], uses the combination of UD and CD discretization schemes. The method amounts to the use of CD scheme for diffusion dominated flows (-2 < Pe < 2) and the use of UD scheme for convection dominated flows (|Pej > 2) for discretization of the convective fluxes. Spalding suggests the following implementation of HD scheme: (P^)e - (i* | | ) e = (Pu)e [(1-Ye) <t>E + Ye Op] (3.13) where Ye =[o.5 + ^max (^ , =E*, 1 j] (3.14) It can be seen that three regions of interest can be defined, based on the above expression: (i) Pe < -2 (convection dominated flow) (PU<J))e- =(Pu)e(t)E (3.15) which is the UD representation of convective flux, (ii) -2 < Pe < 2 (diffusion dominated flow) M - (r, | [=(PU), [ • s i t ] + v te] (3.16) e Sxg 56 which is the CD representation of convective and diffusive fluxes. (iii) Pe > 2 (convection dominated flow) which is the UD representation of convective flux. The last integral in equation 3.3 is represented as : S^dV=[sp*())p + Su*](rArAx)p (3.18) Sp1*" and Svfi pertaining to different variables are listed in table 3.2. The linearization of the source term S<j> as shown above, has been suggested by Patankar [1980] to enhance numerical stability. The final form of the control-volume version of the general transport equation is p * P _ * E ( R A R A X ) + ( P U ) E [ ( 1 . Y E ) ( F E . ^ ) ] ( R A R ) St + (pu)w [YW (<t>P - 0w)] (r Ar) + (pvV.[(l-Yn) (<t>N - Op)] Ax (3.19) + (pv)s [Ys (<t>P - <|>s)] Ax - [Sp* 0 P + Su-"] (r Ar Ax) = 0 The equation can be written in a more compact form as suggested by Patankar [1980]: A P < l ) P = Z(Ai<t>i) + AP,<l>p + Su<l, (3.20) i where A p = ^ A i + A P ' S P* »• X A i = A N + A S + A w + A E i i A N = (pv)n(l-Yn)Ax , A s =-(pv)s Ys Ax , A E = (pu)e(l-Ye)rAr , A w =-(pu)w Yw r Ar 57 Variable Sp* Sn«f u 0 ( p ^ r A r + L f f ^ r A r + j r p e f f ^ A x ' 8x / w ' 5x's V 0 { p l A A x+jpeff ^ j 6 r A r + U e f f ^ f A x ' 8r ' w \ 8rJ s rw -2(rpeff w^Ax 2(rpeff w)s Ax k -Cap 2 k r A r A x Pk r Ar Ax e —f-—r Ar Ax k ci ^-Pk r A r A x e -p |3* ft r Ar Ax P e r A r A x a 2 -p ft p+2an(|J-]2 r A r A x P ^ Pe r Ar Ax Table 3.2 Discretized source terms Su and Sp 58 3.3.2 False diffusion error: It was explained in the previous section that if the cell Peclet number exceeds ±2, the convective fluxes should be discretized by the Upwind Differencing (UD) scheme. Raithby [1976.a] shows that such a scheme produces a numerical error which is commonly called "False Diffusion". The reason for adoption of such a terminology can be explained by considering the Taylor series expansion of the variables § v and <()£ which appear in the discretized convective flux representation. The convective flux across the east boundary of a typical cell was shown to be (refer to previous section): u w x -i J (pu)e(()ErAr forUD Fe=[(pu<|>)e(rAr)pJ = <t>E + <t>p (Pu)e F V ^ r A r f 0 r C D (3.21) The Taylor series expansion of $pand (|>E in the above equation are: (3.22.a) (3.22.b) In the UD scheme, it is assumed that <t>e = <t>E • From Taylor series expansion of 0E , it can be seen that the error associated with such approximation is of the order ^r-=r~ ; also the 2 dx A Y "(J) truncation error in UD representation of convective flux would be (puu^^—r Ar. This 2 dx truncation error has the form of a diffusion flux with a diffusion coefficient of: rfalse = (p U ) e ^ 59 (3.23) The Central Differencing (CD) approximation of convective flux is free of this false diffusion error; M this is because the gradient term I cancels out (adding the two series would give the CD approximation). Patankar [1980] argues that this interpretation of false diffusion error is misleading; this is because the Taylor series expansion is only meaningful for very small values of 8x , in which case the cell Peclet number would be small and UD scheme need not be applied. Raithby and Patankar suggest that the basic cause of false diffusion can be attributed to the practice of treating the flow across each control-volume cell as locally one-dimensional. They suggest that false diffusion occurs when there is a strong cross-stream gradient of <(> or when the flow is oblique to the direction of grid lines. Vahl-Davis and Mallinson [1972] show that an approximate magnitude of the false diffusion coefficient can be obtained from: where U is the resultant velocity vector and 9 is the inclination of this vector with respect to the x direction. False diffusion error can be significant in the weak swirling flows, due to the presence of a recirculation bubble near the center (explained in Chapter 1). It can also be significant in the case of swirling flows in short chambers due to high cross-stream gradient of tangential velocity (for high-turbulence-intensity flow). 33.3 The Bounded Skewed Hybrid Differencing (BSHD) scheme: One of the schemes that reduces the false diffusion error is the Skewed Upwind Differencing (SUD) scheme, proposed by Raithby [1976.b]. This differencing scheme takes the skewness of the velocity vector with respect to the grid lines into account. The full explanation of p U Ax Ay sin (20) (3-24) 60 this scheme has been given by Raithby; the method is briefly explained here with reference to the discretization of the convective flux through the west boundary of a typical cell shown in Fig 3.5. Assuming the cell Peclet number is greater than +2, then the Upwind Differencing scheme prevails. The skewness of the velocity vector passing through the west boundary is taken into account by extending the vector upstream until it intersects one of the grid lines shown in Fig 3.5 If the vector intersects the SW-W grid line, <|>w is calculated based on a linear interpolation of the two grid point values <J>sw and Owviz: <pw = (1 - kw) <pw -»- kw <psw (3.25) here k w is a suitable interpolation function which is given by (assuming u>0, v>0): (3.26) k w = min 1 ilYw '2 Ax AyJ If on the other hand, the projected velocity vector intersects SW-S grid line, (j)w is equated to the closer of the two 0sw or <j>s . Benodekar et al. [1983] showed that even though the SUD scheme reduces the false diffusion error, it might give rise to numerical instability in the form of non-physical oscillation in the solution. The cause of these oscillations has been attributed by Benodekar et al. to the negative coefficients appearing in the difference equation (equation 3.20) n» Lai and Gosman [1982] have developed a flux-blending scheme which eliminates negative coefficients by blending the skewed differencing and upwind differencing; this results in the so-called Bounded Skewed Upwind Differencing (BSUD) scheme. Also to take advantage of the smaller truncation error associated with CD scheme, a hybrid of BSUD and CD scheme is used; this is the so-called Bounded Skewed Hybrid Differencing (BSHD) scheme. The full account of this scheme can be found in the work of Benodekar et al. [1983]; for clarity, this scheme is briefly 1 1 Patankar [1980] shows that if any of the coefficients in the difference equation is negative, unrealistic solution may be obtained. 61 Figure 3.5 A typical cell pertainig to skewed upwind differencing scheme. 62 explained with reference to the calculation of the total flux through the west boundary of the cell shown in Fig. 3.5. The total flux is given by : F w = Xw F wcD +(l • A-w)F wSUD (3.27) where F W C D = total flux through the west boundary using CD. FwCD = (pu}w <l>W + <t>p 6xw (3.28) and FwSUD = total flux through the west boundary using SUD. FwSUD = ( p u ) w <t>w ( p u ) w - 2-<(>w 8xv k w (<|>w - <t>sw) (3.29) and X = 1 for -2<Pe< 2 0 for Pe>2 (3.30) It can be seen that if the velocity vector is parallel to the x direction (i.e. no skewness), from equation 3.26, k w = 0 and the total flux FWSTJD reduces to : FwSUD = ( p u ) w <t>w = F W U D (3.31) where F W U D is the UD representation of the total flux through the west boundary of the cell. The fluxes across the other boundaries are calculated based on similar expressions (given by Benodekar et al. [1983]); the final difference equation is then written in the following form: 63 (Ap-Sp*)<|>p = A N (t>N + A s <(>s + A E <t>E + Aw <l>w + A N E <t>NE + A N W $ N W + Asw 0sw + A S E <1>SE + Su* (3.32) The expressions for all the coefficients have been given by Benodekar et al. [1983]. 3.4 Implementation of boundary conditions It was explained in the previous Chapter that 4 types of boundary conditions were encountered in the present study. In this section, the manner in which these boundary conditions are implemented into the numerical procedure is briefly explained. 3.4.1 Inlet boundary: The inlet boundary condition amounts to the specification of values of all variables at the inlet boundary. 3.4.2 Outlet boundary: At the oudet boundary, the mass flow rate is fixed by the overall mass balance. This is achieved by correcting the velocity component normal to the exit plane with an incremental value at each iteration step, viz: where u p is the velocity at node p, u w corresponds to the velocity at the upstream node and U i n C is calculated from: Up = U w + Uinc (3.33) 'inc — X (puA)jn - X (P"A)bi pA0ut lUt (3.34) 64 This method of incremental correction should theoretically be applied only when the flow at the boundary is everywhere outwards. Once the solution has converged ujnc becomes equal to zero (or a very small number). 3.4.3 Axes of symmetry: For the cells adjacent to the axis of symmetry, the normal gradients of all variables are set to zero. This is implemented by setting the fluxes across the axis to zero; for example if the east boundary of the solution domain coincides with the line of symmetry then : A E =0 . Also the values on either side of the axis are set to be equal. The velocity normal to the boundary is also set to zero; in this case ue=0. 3.4.4 Walls: The usual no-slip condition prevails at the walls. The implementation of the wall boundary condition is carried out in two steps. First, the fluxes at the node adjacent to the wall are set to zero for all variables. For example, in the case of a wall coinciding with the north boundary of the solution domain, A N for all variables is set to zero. Next, the diffusive flux which is calculated based on the wall functions (explained in Chapter 2) is implemented in the source term of the discretized equation. The incorporation of the wall functions for each variable is different; this warrants some explanation which is given next (a) u, v velocities: The diffusive flux for these velocity components is due to the wall shear stress. This is incorporated in the source term of the discretized equation as: (volume)p = (volume)p + xw a U A w a U (3.35) where A w a u is the area over which the shear stress acts. The wall shear stress is calculated from (for u velocity as an example): 65 Twall = A, Up (3.36) where A, is the coefficient which can take on two values, depending on the location of the grid point with respect to the wall (refer to Fig. 3.6), viz M-lam y P piccj}-25 Vk" In (EY+) for Y + < 11.63 for Y + > 11.63 (3.37) here Y + is the dimensionless distance between the nodal point P and the wall which is defined by: Y + = ypUx v (3.38) where uT = Y ^ = cJ-25Vk" (3.39) (b) Turbulence kinetic energy: The production terms in the turbulence kinetic energy equation are incorporated in the source term of this equation (refer to Chapter 2). At the nodal point adjacent to the wall, the production terms are expressed in terms of the wall shear stresses; for example for the north wall (refer to Fig. 3.6): Pk = |it (3.40) this is replaced by: Pk = (Xwau)u ^ + (TwallXv TT (3-41) 66 A= 0 N — ... i W a | | Figure 3.6 A typical cell adjacent to the wall. 67 where (TWaii)u »fawaiijw represent the wall shear stresses associated with the axial and tangential velocity components. (c) Dissipation rate (e): The dissipation rate of kinetic energy near the wall is forced to take the following value, based on the equilibrium condition (explained in Chapter 2): (d) Angular momentum: It was explained previously that due to axi-symmetric nature of the flow, it is advantageous to solve for angular momentum rather than tangential velocity. This formulation gives rise to a source term in the angular momentum equation which is of the following form : S ^ f ^(rueffw) (3.43) The source term is discretized as (refer to Fig. 3,7): Y5T(r M-eff w) = -=2_ {rn u e f f „ wn - r s p e f f s ws) (3.44) o r r p Ar At the wall one of the two terms is set to zero; for example in the case of a north wall, the first term is set to zero. The diffusion term should also be replaced by a wall shear stress which is incorporated in the source term of the discretized equation, viz: XwaU = X Wp (3.45) where Xwas defined by equation 3.37. Figure 3.7 Scalar cell used to calculate the term Y ^ (rH*ff w) 69 3.5 The method of solution The general finite volume equation for each variable is solved using the PISO algorithm, developed by Issa [1982]. The main feature of this algorithm is the splitting of the solution procedure into a series of steps whereby the pressure and velocity fields are decoupled at each step; this amounts to a two-stage, predictor-corrector procedure. In the first stage, the pressure field pertaining to the n* iteration step is used to obtain a velocity and a new pressure field. These intermediate values are then corrected in two steps in the second stage in such a way that mass continuity is satisfied. The corrected velocity and pressure field after the second correction are taken to stand for (n+1)* iteration values. The sequence of iterative operation which is followed in the PISO algorithm can be summarized as follows (refer to Benodekar et al [1983] for further detail): (1) Guess a pressure field. (2) A n intermediate set of axial and radial velocity fields are obtained, based on the guessed pressure field. The pressure and velocity fields are then corrected using the two corrector stages of PISO algorithm. (3) The general discretized equation for angular momentum, kinetic energy and its dissipation rate are then solved. (4) The total residual for each variable is calculated from : R«t, = (A p<l) p)-X(Ai(t>i)-Su1» (3.46) i where A i are A N , As etc. (equation 3.20). If the total residual is less than some prescribed value, the solution is assumed to have converged. Otherwise the values prevailing at the end of each iteration step are used as guessed field for the next iteration step. In the case of transient problems, the solution procedure consists of two iteration loops, namely outer time iteration loop and inner loop. The inner loop is similar to that outlined above. The converged solution of the inner loop is used as the initial field for the outer iteration loop. The 70 procedure is initiated by either a specified or a guess field or even a combination of the two for all variables of interest. Depending on the number of grid points within the solution domain, a set of algebraic fmite-difference equations for each variable are obtained. These simultaneous set of equations are solved using the alternate direction iterative procedure (refer to Benodekar et al. [1983]), which utilizes the tri-diagonal matrix algorithm. 71 CHAPTER 4 RESULTS OF N U M E R I C A L TESTS 4.1 General The TEACH II code, described by Benodekar et al [1983], was used in the present study to investigate the decay of swirling flow in short chambers. This code was originally intended for numerical simulation of axi-symmetric, steady swirling flows, with the k-e closure approximation used to model the turbulence parameters. The suitability of the code for swirling flows was tested by simulating a simple steady, laminar flow, for which exact solutions were available. Since the focus of the present investigation was on transient swirling flows, the code had to be modified to incorporate the transient terms. The modified code was tested against a simple transient swirling flow, namely decaying free laminar vortex, for which exact analytical solutions were available. The curvature modification to turbulence closure approximations were also implemented in the code. Three sets of modifications were implemented; namely Bradshaw's modification of the k-1 model, the Launder et al. modification of the k-e model and the Wilcox-Chambers e-Q model. The merits of each model were investigated against previously published experimental results. The present chapter includes a summary of the above-mentioned numerical tests and the conclusions drawn from them. 4.2. Laminar rotating fluid near a stationary disc Figure 4.1 shows the motion generated by a rotating fluid near a stationary disk. The fluid at large distances away from the disk, rotates at a constant angular speed . The main feature of this flow is the secondary flow pattern set up near the wall owing to the centripetal forces. At large dp distances away from the wall, the radial pressure gradient is balanced by the centripetal force P-^ —. Near the wall, owing to the no-slip condition, the fluid velocity decreases and as a result there is an imbalance between the forces. The fluid is accelerated radially inward and, for continuity, causes an upward axial flow at the symmetry axis. The prediction of this swirl-driven 72 secondary flow pattern is of direct interest to the present study; the same flow pattern is generated in the case of decaying, swirling flow in short chambers. The exact solution of two-dimensional, axi-symmetric Navier-Stokes equations is given by Nydahl. J.E. (cited by Schlichting [1979]) in terms of dimensionless velocities defined as : F = — (v = radial velocity) (4.1 .a) rco ^ = no ^ w = ^ S 6 1 1 0 ^ velocity) (4.1 .b) H = -JLr (U = axial velocity) (4.1 .c) Vvco where F, G and H are dimensionless radial, tangential and axial velocities respectively; these are functions of dimensionless distance % defined by % = * ^ (4.2) Figure 4.2 shows the variation of F, G and H with respect to c\ given by the exact solution. The TEACH code was used to simulate this flow field. Figure 4.3 shows a schematic of the geometry used. The solution domain was sub-divided into 20-20 non-uniform grids (as shown). The boundary designated by A-B represents the surface of the disk. The usual no-slip condition was applied along this boundary. The boundary designated by B-D represents the axis of symmetry at which the normal gradients of axial and tangential velocity were set to zero. From the exact solution (figure 4.2) it can be seen that the radial velocity approaches zero at t\ equal to 12.5. This was used to set the boundary designated by C-D. Boundary A-C is along the tip of the disk. The exact solution was used to specify all three components of velocities at A-C. The specified radial velocity at A-C defined the total mass flux into the solution domain. The axial velocity, normal to the exit boundary C-D, was corrected at each iteration step in such a way to Figure 4.1 Schematic of rotating fluid near a stationary disc of finite radius R (from Schlichting [1979]). 74 -1 -r 2 - r 4 " T B 10 12 14 Figure 4.2 Variation of dimensionless axial (H), radial (F) and tangential (G) velocities with respect to dimensionless distance % ~~ ~* Axial velocity H Radial velocity F *~~ Tangential velocity G 75 B D Figure 4.3 Schematic of the numerical calculation domain and boundaries pertainig to rotating fluid near a stationary disc. 76 satisfy the overall continuity. The solution was assumed to have converged when there was a 0.1 per cent difference between the overall mass influx and efflux. Figure 4.4 shows the comparison between the numerically calculated dimensionless velocities and the corresponding analytical solutions. It can be seen that the agreement for the axial and tangential velocities deteriorated at large distances away from the disk. This was attributed to the coarse grid size used in that region. The effect was not reflected in the radial velocity component since its magnitude (lirninished at large distances away from the disk. It was decided to increase the number of grid points to 30-33 in the axial and radial directions respectively. Figure 4.5 shows the comparison between the new numerical results and the corresponding analytical results. The agreement between the two results were excellent. It should be pointed out that the same agreement was observed with 60-60 equally spaced grid points. This proved the importance of grid size and distribution on the converged solution. 4.3 Decay of a free laminar vortex The transient terms, based on the fully implicit time differencing scheme, were implemented in the TEACH II code. To examine the validity of the implementation, the code was used to simulate the decay of a free laminar vortex due to the action of viscosity. The exact solution of the one-dimensional, non-steady Navier-Stokes equation has been given by Ossen. W.C. (cited by Schlichting), viz: ) - « P ( £ 1 (4-3) where To is the vortex strength at time t equal to zero (before the viscosity is set into action) which was calculated by: To = 2 7t r 2 coo (4.4) here COQ is the initial rotational speed. 77 Figure 4.4 Comparison between analytical solution of rotating fluid near a stationary disc and the corresponding numerical results, using 20 by 20 grids. — — — Numerical axial velocity profile O Analytical axial velocity profile •' 1 1 •• Numerical radial velocity profile X Analytical radial velocity profile — — Numerical tangential velocity profile A Analytical tangential velocity profile 78 • <o <Sk H 1 i $ ri K..." \ Ii ! ¥> ! ( 1 j F \ K, A'A»< •? 0 2 4 6 8 10 1 2 U Figure 4.5 Comparison between analytical solution of rotating fluid near a stationary disc and the corresponding numerical results, using 30 by 33 grids. — — Numerical axial velocity profile O Analytical axial velocity profile — Numerical radial velocity profile X Analytical radial velocity profile Numerical tangential velocity profile A Analytical tangential velocity profile 79 The tangential velocity w was set initially to and the calculation proceeded at 0.5 ms time interval up to dimensionless time |^rj equal to 0.1. Figure 4.6 shows the comparison between the numerical and the exact analytical solution after 80ms decay; the observed agreement was within 1 per cent It was thus concluded that the transient terms had been correctly implemented. 4.4 Turbulent flow between rotating cylinders To exarnine the validity of the turbulence modelling for simulation of swirling flows, it was decided to simulate the turbulent flow field between two long concentric cylinders, either of which can be rotating. Such a decision was made for three reasons. Firstly, because the flow was generated by the action of wall shear stresses, accurate prediction of the flow field was only possible if the wall shear stresses were calculated correctly. Hence, this was a good test of the numerical procedure and also of the turbulence modelling, both of which directly influence the wall shear stress calculation. Secondly, in the case of a flow between long concentric cylinders, the radial distribution of shear stress can be approximated directly from the tangential component of momentum equation, which for axi-symmetric steady flow is: In the above equation, the axial variations can be neglected; the radial velocity is small and can also be neglected. It can be seen, from the above equation, that the radial distribution of shear stress is: Tjfl r 2 = constant (4.6) which states that the torque is constant. This relation can be used to examine the performance of the code. 80 Figure 4.6 Comparison between the exact and numerical solution of decaying, laminar, free vortex (coo = 300 rpm). Initial velocity distribution at t=o Analytical velocity profile after = 0. 1 Numerically calculated velocity profile after | ^ 2 j =0-1 V L 81 Thirdly, the effects of curvature clearly manifest themselves in such a flow field. The turbulent flow near the cylinder walls can be either stable or unstable, depending on which one of the cylinders is rotating. The flow is considered to be unstable if any type of disturbance in the flow causes the flow to become turbulent and remain turbulent. The reverse is true for the case of stable flow. The stability condition for such flow fields was first presented by Taylor [1923] who studied a laminar flow between two concentric cylinders of radii Ri and R 2. The flow being generated by either the inner cylinder rotation (fti) or outer cylinder rotation (Q2). Figure 4.7 shows a summary of some of Taylor's results. It can be seen that for rotating outer cylinder and stationary inner cylinder, the flow is stable, i.e. any disturbance in the flow quickly dies out and the flow remains laminar. In the case of rotating inner cylinder and stationary outer cylinder, the flow becomes unstable, when a certain Reynolds number (or Taylor number12) has been exceeded. At such a Taylor number, disturbances (for example in the form of vortices) appear in the flow; as the inner cylinder rotational speed is increased further, the flow becomes unstable and eventually turbulent. The same phenomenon can be observed in the case of fully turbulent flow between concentric cylinders. According to the Rayleigh stability (explained in chapter 1), the flow in the case of the inner cylinder rotating, the outer one being stationary, is unstable. The reverse is the case for rotating outer cylinder and stationary inner cylinder. This can be explained by examining the radial distribution of angular momentum. Figure 4.8 shows a typical angular momentum distribution reported by Wattendorf [1935] for the case of inner cylinder rotation. It can be seen that the angular momentum decreases radially outwards away from the inner cylinder, is nearly constant close to the center and decreases near the outer stationary cylinder. This would mean that the flow is unstable and turbulence is augmented, since the angular momentum decreases with increasing radius. 12Taylor number = ^ £ = (ObE-f = 1700 where h = R 2 - Ri R v 2 R * v ' 82 Figure 4.7 Stability of annular flow presented by Taylor [1923]. 83 1*6 1.2 wr/(ur) l •8 .4 L 0 1 2 3 4 5 Distance from inner wall (cm) Figure 4.8 Radial distribution of angular momentum in an annulus of two concentric cylinders with inner cylinder rotating and outer one stationary (by Wattendorf [1934]). (ur)j = core angular momentum 84 Both inner and outer cylinder rotations were investigated in the present study; the accuracy of the predictions were tested against published experimental results. 4.4.1 Inner cylinder rotating, outer cylinder stationary: Experimental velocity and heat transfer measurements within an annular channel, having a rotating inner cylinder and stationary outer cylinder were reported by Zmeikov and Ustimenko [1966]. Characteristic dimensions of their apparatus have been given in Fig. 4.9. Because the length of the chamber was almost eight times larger than its height, no axial variations were anticipated; hence only a few grid points were located along the axial direction. Due to the symmetry condition, only i- of the channel was considered. 16-5 grid points in the radial and axial directions were used respectively (shown in Fig. 4.9). Three sets of rotational speeds were considered, namely 18.86 m/s, 24.4 m/s and 30 m/s. Zmeikov and Ustimenko have reported wall shear stresses which were measured using two different methods, namely hot wire and Preston tube. They have also reported that their measured tangential velocity profiles can be approximated by the following equations: w+ = Y + (Laminar sublayer) 5 > Y + w+ = - L - In (0.5433 Y +) (Transition) 50 > Y + > 5 ° - 2 (4.7) w + = 1.2 In Y + + 10.3 (Fully turbulent) Tjo > Y + > 50 where rjo = ^ ^ y ^ 1 u<tt f ° r ^ e ^ n Q T w a ^ ; U x l * s & c w a ^ friction velocity, Ri and R 2 are the inner and outer radii respectively. It should be pointed out that Zmeikov and Ustimenko reported a good agreement between their results and the empirical relation used in the transition region which was originally proposed by Spalding [1961]. However, according to Spalding, the transition region extends from Y + of 5 to Y + of 30 (as opposed to the value 50 which was cited by Zmeikov and Ustimenko). It was thus assumed that there might have been an error made by Zmeikov and Ustimenko in reporting the limits; this point is clarified later. 85 E S o "3< Axis of symmetry i •= 200 m m Figure 4.9 Characteristic dimensions and grid arrangement pertaining to simulation of flow field in an annulus of concentric cylinders (Zmeikov and Ustimenko). 86 The TEACH code with the standard k-e turbulence model was used to predict the flow field reported by Zmeikov and Ustimenko. Fig. 4.10 shows the comparison between the calculated dimensionless velocity W against the experimental profile near the inner cylinder, rotating at 18.86m/s. The results revealed that the experimental wall shear stresses were underpredicted by almost 30 per cent The constant torque condition was not satisfied either. The maximum error in the calculated velocity profile, as compared with the corresponding experimental profile, was 25 per cent. Such large discrepancies between the calculated and experimental values were attributed to two sources of error; these were discretization of the source term in the angular momentum equation and the effect of curvature on turbulence modelling. As explained in the previous chapter, in order to use a common form for all transport equations, it was necessary to solve for angular momentum rather than the tangential velocity. For this reason, part of the source term in the angular momentum equation becomes (equation 3.43): S * = f ^(rwpeff) (4.8) The finite volume representation of this term is : S<!>=— (rnWnPeffn-rsWsUeffs) (4.9) r Ar where the subscripts n and s refer to north and south surfaces of the control volume. The term rwpeff at any nodal point is thus represented by a linear interpolation of two neighbouring points. Close to the boundary, where there is a large gradient of peff, a linear interpolation might not be adequate, particularly in the case of a moving boundary. This is because if the wall is stationary, one of the two terms in the above equation drops out; for example if the north face of the control volume coincides with the wall, then wn is zero and the whole term r n wn peff n cancels out. However, if the wall is moving, then wn has a finite value and hence peff n has to be approximated arbitrarily. 87 54 -i 20-y y 16- • • W+ 12-8-4-A _ 10 100 1000 Y+ Figure 4.10 Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with original form of the source term in the angular momentum equation (inner cylinder speed of 18.86 m/s) • Experimental results — — Numerical results 88 This problem can be alleviated to some extent by re-arranging the source term into a more suitable form. By taking derivatives, the differential term in the source term can be expanded as : S* = -2 ^ M l i N + w ar (4.10) The second term in the above equation is the shear stress. The first term is negligible compared to the remaining terms, hence: S* = -2 au eff ar (4.11) Townsend [1956] has suggested the following expression for ueff in the vicinity of the wall: Ueff 5 p K y u x (4.12) where K is the von Karman constant and u x is the wall friction velocity. The substitution of the above expression in the original equation would yield S ^ - l ^ i p w K u , ] (4.13) The positive sign is used for the north wall and the negative sign is used for the south wall. The above expression is of the same form, whether the wall is moving or stationary. The corrected source term was implemented in the code and the simulation was repeated for inner cylinder rotational speed of 18.86 m/s. Fig. 4.11 shows the comparison between the calculated results, using the alternative form of the source term, and the corresponding experimental values. A considerable improvement was observed. The constant torque condition was now satisfied. However, the calculated wall shear stresses under-estimated the experimental values by as much as 10 per cent There was also a'maximum of 6 per cent difference between the calculated velocity profile and the experimental values. It was anticipated that this discrepancy could be due to curvature effects on turbulence modelling. 8 9 24 20-Figure 4.11 Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with modified source term in the angular momentum equation (inner cylinder speed of 18.86 m/s) • Experimental results Numerical results 90 It was decided to implement the Launder et al.[1977] proposed modification of standard k-e model. The modification caused a new coupling between the mean flow and turbulence parameters due to Richardson number. This led to divergence during the iterative solution procedure, whenever the sign of Richardson number was in error. To prevent divergence, a curve fit of angular momentum between four adjacent nodal values was used to obtain the correct sign of the Richardson number. The sign was obtained from the gradient of the fitted curve during each iteration step. If the sign was not correct, the Richardson number was set to zero for that iteration i.e. effectively no curvature correction was applied for that iteration. Fig. 4.12 compares the experimental velocity profde and the profile calculated with the modified k-e model. The discrepancy between the calculated and measured wall shear stresses was only 1.3 per cent. The maximum error on the calculated velocity field was 5 per cent and occurred near the wall. This might have been partially due to the logarithmic profile used by Zmeikov and Ustimenko to represent their experimental velocity in the range of 5<Y+<50. If the upper limit of Y+equal to 30 were used to represent the transition region, as was originally proposed by Spalding, the discrepancy near the wall would have reduced to less than 1 per cent. Other than this difference, the comparison was excellent. It was also decided to simulate the same flow field using the Wilcox and Chambers model. The model was implemented in the code and was then used to simulate the flow field. Fig. 4.12 shows the comparison between the calculated and experimental profiles. Excellent agreement was observed, except in the near wall region. This discrepancy was again attributed to the upper limit of Y + used by Zmeikov and Ustimenko in defining the transition region. It should be pointed out that the boundary condition used in the e-ft model (logarithmic law of the wall) was the same as that used in the modified k-e model; hence the same near wall discrepancy was observed with the two different models. These models were also used to simulate the flow field for two more inner cylinder rotational speeds of 24.4 m/s and 30 m/s. The comparisons are shown in Fig. 4.13 and 4.14 respectively. It was seen that the disagreement between the calculated and experimental values increased with increasing rotational speed, away from the inner wall. In both cases a closer 91 20 12-10 10 100 1000 Y+ Figure 4.12 Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with curvature corrections (inner cylinder speed of 18.86 m/s) • Experimental results Numerical results with Launder et al. modification ——— Numerical results with Wilcox and Chambers model 92 i o 4 -10 11 100 1000 Figure 4.13 Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with and without curvature corrections (inner cylinder speed of 24.4 m/s) • Experimental results Numerical results with no curvature correction — —. — Numerical results with Launder et al. modification Numerical results with Wilcox and Chambers modifications 93 12-10 ! i 10 100 1000 Figure 4.14 Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results with and without curvature modification (inner cylinder speed of 30 m/s) • Experimental results — — Numerical results with no curvature correction Numerical results with Launder et al. modification • Numerical results with Wilcox and Chambers modification 94 observed with the Wilcox and Chambers model. The maximum error observed was within 5 per cent. It was also observed that with inner cylinder rotational speed of 30 m/s, the disagreement between the numerical and experimental values near the wall disappeared completely. In this case, the first calculation nodal point was located at Y + of 51. According to Zmeikov and Ustimenko, the flow at this Y + is fully turbulent; the expression used by them to represent their experimental values in this region gave a similar value to that obtained from the law of the wall used in the present numerical model. Hence if the upper limit of Y + equal to 30 were used to separate the fully turbulent region from the transition region in the previous two cases (18.86 m/s and 24.4m/s), the same good agreement would have been observed. In order to obtain grid independent results, and to reduce the discrepancy between the calculated and the corresponding experimental results, the number of grid points in the radial direction was increased from 16 to 36. However, the nodal points adjacent to the walls were kept at Y + greater than 30. This was based on the observation that the logarithmic law of the wall is only applicable at these Y + . Also at locations close to the wall, where the turbulent Reynolds number is low, the turbulence models do not perform well; extra terms have to be added to account for low Reynolds number effects (see for example Jones and Launder [1972]). To avoid the use of low Reynolds number modelling and to be certain that the logarithmic law of the wall holds, the first nodal point was always kept at Y + greater than 30. Fig. 4.15 and 4.16 show the comparison between the numerical results, obtained with refined grid points, with the corresponding experimental results for inner cylinder rotational speeds of 24.4 m/s and 30 m/s respectively. The results in both cases were obtained using the Wilcox and Chambers model. The agreement between the numerical and experimental results did not improve as the number of grid points were increased. It was thus concluded that the results with 16-5 grid points were grid independent. 95 Figure 4.15 Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results using Wilcox and Chambers model with refined grids (inner cylinder speed of 24.4 m/s) • Experimental results —— Numerical results 96 Figure 4.16 Comparison between experimental dimensionless velocity due to Zmeikov and Ustimenko and numerical results using Wilcox and Chambers model with refined grids (inner cylinder speed of 30 m/s) • Experimental results Numerical results 97 4.22 Outer cylinder rotating, inner one stationary: As proposed by Taylor (Fig. 4.7), the flow under such conditions is stable; also based on Rayleigh stability criterion (refer to Chapter 1), because the angular momentum increases radially outwards from the inner to outer wall, the flow is stable and this confirms Taylor's suggestion. To test the modification proposed by Launder et al. and the model proposed by Wilcox and Chambers for stable conditions, it was decided to simulate the flow in an annulus between two concentric cylinders with the outer cylinder rotating. The experimental results of Taylor [1936] were used to validate the simulation. Fig. 4.17 shows the schematic of the configuration used by Taylor in his experimental study with the outer cylinder rotational speed of 33.3 rps. Fig. 4.18 shows the comparison between the experimental and the corresponding numerical results using the standard k-e model. It was observed that the standard k-e model over-predicted the wall shear stresses by as much as 45 per cent. The curvature in this case tends to reduce turbulence and hence the wall shear stress. This suggested that the discrepancy could be due to curvature effects. The predicted velocity profile was also in large disagreement with the experimental profile. The Launder et al. modification of the k-e model was used to account for curvature effects. Fig. 4.18 shows the comparison between the experimental and the numerical results pertaining to the modified k-e model. Excellent agreement was observed. The constant torque condition was also satisfied up to 1 per cent. However, these agreements were only obtained after optimizing the empirical constant associated with Launder et al modification term (Q in equation 2.28). It was found that the constant had to be reduced by as much as 50 per cent (Cc = 0.1), as compared to the value originally suggested by Launder et al. (C c = 0.2). It was mentioned earlier that Priddin arrived at the same conclusion in his study. Launder's explanation of this discrepancy was that the modification should be incorporated in the production terms of dissipation equation (with an opposite sign) as opposed to the decaying parts. The interpretation of this suggestion by the present author is that the extra term should increase the dissipation rate for stable flows and decrease the dissipation rate for unstable flows; this can only be achieved by implementing the modification as part of the production term in the dissipation equation. No attempt was made to test the validity of Launder's suggestion in the present study. 98 £ E co E E CN C 5 ~ = 40.2 mm Si Figure 4.17 Characteristic dimensions and grid arrangement pertaining to simulation of flow field in an annulus of concentric cylinders (Taylor flow field) 9 9 The same flow field was simulated using the Wilcox and Chambers turbulence model. Fig. 4.18 shows the comparison between the numerical and experimental results. Identical results were obtained with both modified k-e and Wilcox-Chambers models, with the difference that no constant optimization was necessary in the case of the Wilcox and Chambers model. 4.4.3 Conclusions: The general conclusion drawn from these test cases was that curvature has substantial effects on turbulence, velocity field and wall shear stresses. The standard k-e model is not well equipped to reproduce these effects. The test cases have proven that the Launder et al. modification of the k-e model and also the Wilcox and Chambers mixing energy model can both be effectively used to account for curvature effects. It was also shown that the Wilcox and Chambers model was superior to the Launder et al. modification in that the former has no need for additional case-tuned empirical constant 4.5 Decay of a turbulent swirling flow As it was explained in Chapter 1, at near top-dead-center, the IC engine flow field can be simulated in constant-volume chambers. The swirling flow in such chambers is normally generated by injection of fluid, through specially designed valves, though other methods are possible. Once the valve is closed, the flow starts to decay until a quiescent state is attained. The rate of decay is governed by the wall shear stresses and is influenced by the chamber geometry; the rate of decay is expected to increase as the aspect ratio (length/diameter) is decreased. An accurate prediction of wall shear stresses is thus an essential factor in determining the rate of flow decay. This would require an accurate prediction of the turbulent field. To examine the merits of turbulence models in predicting such flow fields, it was decided to simulate the decay of an isothermal, turbulent swirling flow inside a closed chamber. Experimental results of Inoue et al [1980], who studied the decay of swirling flow in a chamber of 120 mm diameter and 27 mm length were used to validate the simulation. 100 Figure 4.18 Comparison between experimental velocity profile due to Taylor and numerical results with and without curvature corrections. • Experimental results — — Numerical results with no curvature correction Numerical results with Launder et al. modification Numerical results with Wilcox and Chambers model 101 The experimental results reported by Inoue et al. included mean tangential velocity and turbulence intensity profiles at several time delays after completion of injection process. Turbulence was suspected to have been generated by two mechanisms. One by injection process and the other through turbulence generation in the shear layer close to the wall of the chamber. The experimental profiles at 20 ms after the decay had commenced, were used as initial conditions to start the simulation. However, no initial data for axial and radial velocity profiles were available. It was decided to start the simulation with zero values for both radial and axial velocities. The physical domain was divided into 29-19 grid points in the r and x directions, respectively. The fully implicit time discretization technique, with 1 ms time steps, was used to solve for the transient terms. The effect of assuming zero initial axial and radial velocity field on the final results was tested by comparing the results after 30 ms decay with those when the velocity distribution after the first time step (after 1 ms decay) were used as an initial profile for axial and radial velocities. The maximum difference of 1.5 per cent in the calculated tangential velocity profile between the two cases was observed. The difference was anticipated to increase as the aspect ratio was reduced. This was because with a lower aspect ratio, the secondary flow pattern plays a more important role on the rate of decay of mean tangential velocity. The k-e model, which included the Launder et al curvature correction term, was used for this simulation. The initial dissipation rate profile was obtained from the following relation : e = Cn k 1 5 (4.14) 0.09 S in which 8 was arbitrary taken to be equal to 9 mm. The initial turbulence kinetic energy profile was obtained from the measured stream-wise intensity, assuming isotropic turbulence k=1.5w'2 Fig. 4.19 shows the comparison between the experimental tangential velocity and the calculated profile at the mid-plane after 30 ms decay. The maximum error of the order 2 per cent close to the center was observed. Fig. 4.20 shows the measured and calculated tangential component of turbulence intensity after 30 ms decay. To obtain the stream-wise intensity, the turbulence was again assumed to be isotropic, so intensity was directly calculated from the total 102 30 Figure 4.19 Comparison between experimentally measured decaying velocity (by Inoue et al.) and the corresponding calculated field using k-e model with large initial dissipation rate. Initial velocity profile fitted to experimental values • Experimental velocity profile after 30 ms ~~ — Numerical velocity profile after 30 ms 103 1.8-1.6-1.4-1.2-1-» ' (? ) 0.6-0.6-0.4-0.2-/. m J 0 U 1 • 1 • • 1 A \ A 1 A IA • ! A 0-c > 10 20 SO 40 SO 6 R (mm) Figure 4.20 Comparison between experimentally measured turbulence intensity (by Inoue et al.) and the corresponding calculated field using k-e model with large initial dissipation rate. • Measured intensity profile at time zero A Measured intensity profile after 30 ms Numerical intensity profile after 30 ms 104 kinetic energy. The computation indicated an extremely fast decay of turbulence up to about 45mm radius. This substantial discrepancy in the turbulence field was not reflected severely in the velocity profile for the following reason. The turbulence effect on velocity field manifests itself through either turbulent viscosity or gradient of turbulence kinetic energy. The velocity near the center exhibited a nearly solid-body type of rotation, for which the strain-rate is zero. Hence, the turbulent viscosity effect is negligible. Also, as is apparent from the numerical results, the turbulence kinetic energy is very uniform near the center and hence the gradient effects are negligible. Thus the error in the calculated turbulence intensity was not reflected as much in the velocity field. The discrepancy between the measured and calculated intensity, near the center, was attributed to the over-estimation of the initial dissipation rate. Any turbulence energy generated at the wall, and consequently transported towards the center, was quickly dissipated. It was thus decided to start with a very small initial dissipation rate profile and also proceed with a smaller time step, namely 0.5 ms. This allowed the diffusion process to yield a reasonable distribution of dissipation rate after the first few time steps. By the same token, a smaller time step was expected to alleviate the problem of unknown initial radial and axial velocities to some extent. Fig. 4.21 shows the comparison between the experimental and calculated velocity profiles after 30 ms decay with zero initial dissipation rate. The agreement improved, which resulted in a smaller maximum error of the order 1 per cent. Fig. 4.22 shows the corresponding turbulence intensity profiles after 30 ms decay. A remarkable improvement, compared to the previous case was obtained, particularly close to the center. The large discrepancy between the measured and calculated intensity values near the wall was attributed to the isotropy assumption which was necessary for calculating the stream-wise intensity from the total calculated kinetic energy. The decay of velocity and turbulence intensity after 80 ms is also shown in Fig. 4.21 and 4.22 respectively. An interesting feature pertaining to the flow under investigation was the existence of a combination of stable flow, near the center, and unstable flow close to the wall of the chamber. It 105 SO 0 10 20 30 40 50 60 R (mm) Figure 4.21 Comparison between experimentally measured decaying velocity (by Inoue et al.) and the corresponding calculated field using k-e model with zero initial dissipation rate. Initial velocity profile fitted to experimental values • Experimental velocity profile after 30 ms Numerical velocity profile after 30 ms O Experimental velocity profile after 80 ms Numerical velocity profile after 80 ms 106 1.6 -1 1 1 I I 10 20 30 40 50 R (mm) Figure 4.22 Comparison between experimentally measured turbulence intensity (by Inoue et al.) and the corresponding calculated field using k-e model with zero initial dissipation rate. • Measured intensity profile at time zero A Measured intensity profile after 30 ms Numerical intensity profile after 30 ms • Measured intensity profile after 80 ms Numerical intensity profile after 80 ms 0 + o 107 was noticed that the Launder et al. modification of k-e model in the stable region, gave rise to a faster decay of turbulence than that exhibited by experimental results. The profiles shown in Fig. 4.21 have been obtained by setting the correction term to zero in the stable region. The feasibility of Bradshaw's curvature correction of length scale was also studied in the context of decay of swirling flow in short chambers. The Inoue et al. experimental results were used to validate the numerical prediction. The one equation k-1 model, explained in chapter 2 was used. The experimental profiles of velocity and tangential intensity were used to define the initial conditions for numerical simulation (as for the previous case). The length scale proposed by Norris and Reynolds [1936] (refer to chapter 2) and modified by Bradshaw's curvature correction was used to calculate the initial length scale distribution within the boundary layer on the curved wall of the chamber. The constant (3 used in Bradshaw's modification was set equal to 10. Previous experimental results on flows in short chambers (e.g. in cylinders of IC engines) have suggested that turbulence field in such chambers is isotropic and homogeneous with the length scale equal to chamber width. The length scale outside the boundary layer was thus assumed to be constant and equal to the chamber width, in the present simulation. Fig. 4.23 shows the comparison between the calculated and experimental velocity profiles after 30 ms decay. Good agreement was observed near the center. However the model predicted a thicker boundary layer than that exhibited by experimental profile. This discrepancy was attributed to under-prediction of the length scale near the wall, which was caused by large value of P used in Bradshaw's modification of length scale. This under prediction of length scale gave rise to a lower turbulent viscosity value and hence lower wall stresses. The comparison between the corresponding calculated and experimental turbulence intensity is shown in Fig. 4.24. The agreement between the experimental and calculated profiles after 30 ms decay deteriorated as compared with the calculation using the modified k-e model. The intensity was under-predicted near the center and over-predicted near the wall. The kink observable in the calculated intensity, around 2. radius, was caused by the sudden change in the length scale distribution i.e. from constant value outside the boundary layer to that proportional to the distance 108 30 - i R (mm) Figure 4.23 Comparison between experimentally measured decaying velocity (by Inoue et al.) and the corresponding calculated field using curvature modified k-1 model. —— Initial velocity profile fitted to experimental values • Experimental velocity profile after 30 ms Numerical velocity profile after 30 ms 109 1 . 8 -1 . 6 -1 . 4 -1 . 2 -» ' ( ? ) " 0 . 8 -0 . 6 -0 . 4 -0 . 2 -fs -7 • • 1 • 1 • • A A A a . . A w " ( > 1 0 2 0 3 0 4 0 S O 6 0 R (mm) Figure 4.24 Comparison between experimentally measured intensity profile (by Inoue et al.) and the corresponding calculated field using curvature modified k-1 model • Measured intensity profile at time zero A Measured intensity profile after 30 ms Numerical intensity profile after 30 ms 110 from the wall inside the boundary layer. However, in sum, the exercise suggested that with a good choice of value, reasonable prediction can be obtained with Bradshaw's modified k-1 model. The Wilcox and Chambers model was also used to simulate the same flow field. This was helpful in determining the validity of the model for decaying swirling flow prediction and also in investigating the non-isotropic nature of turbulence for swirling flows. As explained in Chapter 2, the effect of swirl is the redistribution of total kinetic energy between the normal and stream-wise components of intensity. Hence it was expected that the normal intensity calculated by e-Q model would under-predict the experimental stream-wise intensity near the center, and would over-predict it near the wall, where swirl is expected to increase the normal intensity and decrease the stream-wise intensity. Fig. 4.25 shows the comparison between the experimental and calculated velocity profde, using the Wilcox and Chambers model, after 30 ms and 80 ms decay. Good agreement was observed, suggesting that the wall shear stresses had been predicted correcdy since these were responsible for decay of swirl. Fig. 4.26 shows the comparison between the calculated normal intensity and the experimental stream-wise intensity after 30 ms and 80 ms decay. The flow near the center can be considered to be stable, while the flow near the chamber's wall is unstable. Hence, as suggested by Wilcox and Chambers turbulence energy redistribution criterion, the calculated normal intensity was smaller than the stream-wise intensity near the center and was larger than the stream-wise intensity near the wall. It should be pointed out that because mixing energy (e=^ v^ J and kinetic energy k are equal near the wall, the wall value shown in Fig. 4.26 was smaller than the The general conclusion drawn from the simulation of decaying swirling flow in short chambers was that the decay rate can be reasonably predicted using the modified k-1 model, the modified k-e model and e-Q model. However, the Wilcox and Chambers e-Q model is superior to the other two models because it does not need an extra adjustable empirical constant. In the case of the k-e model, modified by the Launder et al. proposal, the modification term had to be set to zero in the stable part of the flow. In the case of Bradshaw's modified k-1 model, the accompanying corresponding value shown in Fig. 4.22, which was equal to U l 3 0 1 i i l l l l T 0 10 20 30 4 0 50 60 R (mm) I Figure 4.25 Comparison between experimentally measured decaying velocity (by Inoue et al.) and the corresponding calculated field using Wilcox and Chambers model. • Initial velocity profile fitted to experimental values • Experimental velocity profile after 30 ms -—— Numerical velocity profile after 30 ms A Experimental velocity profile after 80 ms — — Numerical velocity profile after 80 ms 112 2.5 2 -cH 1 1 1 1 1 0 10 20 30 40 50 60 R (mm) Figure 4.26 Comparison between experimentally measured turbulence intensity (by Inoue et al.) and the corresponding calculated field using Wilcox and Chambers model. • Measured intensity profile at time zero A Measured intensity profile after 30 ms Numerical intensity profile after 30 ms • Measured intensity profile after 80 ms Numerical intensity profile after 80 ms 113 empirical constant needed to be adjusted to attain good agreement; also the k-1 model suffered from the need to prescribe the length scale. It was also concluded that the initial distribution of dissipation rate could have substantial effects on the overall prediction of the flow field and turbulence intensity. It was concluded that with sufficiendy small values for dissipation rate and correspondingly small time step, the problem of unknown initial profile can be alleviated. 114 Chapter 5. EXPERIMENTAL SET UP AND MEASUREMENTS TECHNIQUE 5.1 Introduction Swirl can be generated in a number of different ways; tangential injection of fluid, specially designed inlet valves and rotating impeller inside a chamber are a few examples of the methods by which swirl is generated. The third method of swirl generation was chosen for the present experimental investigation. This method has an advantage and a disadvantage over the swirl generation by injection process. The advantage is that a truly steady, two-dimensional flow field can be generated in the case of former whereas the swirling flow generated by steady injection of fluid is not two-dimensional. The disadvantage associated with the present method of swirl generation is the low turbulence intensity level. The turbulence is generated through the shear action at the walls of the chamber and if the disk is smooth the overall intensity level is low compared to the case where turbulence is generated by the injection process13. To promote turbulence generation in the present study, the disc surface was roughened. The experimental study of steady, two-dimensional swirling flow is beneficial, because the accuracy of the experimental results can be tested against an already validated numerical model which has been developed in the present study. If the agreement between the experimental results and the corresponding numerical results is satisfactory, the latter can be used to define the initial conditions for transient analysis. The mean velocity and the turbulence intensity measurements in the present study were carried out using the Laser-Doppler-Anemometry (LDA) technique. This technique has several advantages over the more conventional hot-wire anemometry technique. The LDA method does not require any calibration. It is not intrusive i.e. there is no probe to alter the flow characteristics. It can also be used for hot flow measurements e.g. in the presence of combustion. 1 3 I n t h i s c a s e t h e t u r b u l e n t l e n g t h s c a l e i s c h a r a c t e r i z e d by t h e s i z e o f t h e v a l v e o p e n i n g . 115 In the present Chapter the experimental setup is explained first, followed by a brief review of the LDA technique. The decision on the choice of the working fluid is also discussed. 5.2 Constant volume chamber A cylindrical chamber14 of 4 in diameter and 3.5 in length was bored inside a 6 in by 5.4 in by 4.5 in aluminum block (Fig 5.1). The chamber length was reduced to 2 in by adding a removable end plate which closed the chamber (Fig 5.2). The length to diameter ratio of the chamber was varied by adding extra blocks of aluminum inserts to the end plate. In overall three length-to-diameter ratios were considered, namely ^ and . The original length to diameter ratio was ^ ; inserts of 0.8 in and 1.6 in lengths were added to the end plate to obtain the required length to diameter ratios. The impeller was made of a 0.03 in thick aluminum disk whose diameter was 3.868 in. The disk was mounted on the end plate with a clearance of 0.05 in by screwing its hub to the plate. The disk was also connected to a ^ in diameter stain-less steel shaft through its hub (Fig. 5.1). To accommodate the shaft, a shaft housing of 0.2503 in I.D. and 0.75 in O.D. was drilled inside the end plate. The sealing between the shaft and the housing was provided by labyrinth seal (Pierik [1987]). The details of the seal design are given in Fig. 5.3. The other end of the shaft was connected to a D.C. motor via a universal joint which made the alignment procedure less critical. The motor speed was varied from 20 to 25000 rev/min by adjusting its input power supply through a voltage regulator. The LDA technique requires optical access to the measuring point To accommodate such a need, two windows of 2.7 in by 1.7 in were cut out on the curved walls, on either side of the chamber. Quartz glasses were inserted into these openings and were kept in place with holder plates (Fig 5.4). The chamber and the motor were mounted on a traversing table which was able to move in all three x-y-z directions. Because the flow was believed to be axi-symmetric, the measurements The chamber was d e s i g n e d by P i e r i k . R. Figure 5.1 Schematic of a constant volume cylindrical chamber assembly Figure 5.2 Schematic of the chamber's end plate cross-section A l l dimensions i n inches. Figure 5.3 Schematic of Labyrinth seal details 119 A l l dimensions in inches. Figure 5.4 Schematic of quartz glass holder plate 120 were only made in the x-r plane, along the mid-plane. The details of the measurements are explained in the following sections. 5.3 Laser Doppler Anemometry (LDA) 53.1 Introduction In 1842, Doppler, an Austrian physicist considered a phenomenon of the change in the frequency of acoustic and electromagnetic waves owing to the movement of the source, receiver or intervening reflectors and scatterers. The principle behind the operation of Laser-Doppler-Anemometry technique is based on this doppler frequency shift phenomenon. The frequency of a given light source can be altered by the presence of scattering or reflecting particles in its path. The change in frequency of the light beam depends on the velocity of intervening particles and the direction of scattering. However, the change in the frequency is very small15 and is not direcdy detectable by optical spectrometry. The principle of optical hetrodyning or mixing is used to measure this small frequency difference. This principle is based on the addition of two light waves which are then passed through a detector with a non-linear response16. The detector output contains the difference in the frequency of the two scattered waves which is a measure of particle velocity. The two light sources are generated from a single beam. The LDA technique is generally divided into two modes of operation, based on the methods by which these two beams are used to generate a doppler frequency. These are the reference beam mode and the differential doppler mode. In the case of the reference beam mode, the frequency of one beam is shifted (due to intervening scatterers) while the frequency of the other remains unchanged and is used as a reference frequency. The two beams are then passed through a detector whose output signal 1 5 7 parts in 108 for velocity of 10 m/s (Drain [1972]). 1 6 This w i l l be discussed in more d e t a i l in the later parts of the present Chapter 121 includes the difference between the frequencies; a typical arrangement of this mode is shown in Fig 5.5. In the differential doppler mode, the frequency of both beams is shifted by the intervening scatterers. The change in the frequency is associated with the direction of scattering and the velocity of scattering particles; since the two beams are scattered at different angles, their frequency changes by different amounts. The output of the detector includes the difference between these two frequencies. The differential doppler mode is superior to the reference beam mode for two reasons. Firstly, the frequency of the detector output signal in the case of former is independent of the direction in which the detector receives the light (see Drain [1972]); hence the scattered light can be collected in any direction. Secondly in the reference beam mode, the detector aperture size is restricted by the coherence requirement, as explained by Drain [1972]. The smaller aperture size would result in a lower signal intensity. The coherence requirement is automatically satisfied in the differential doppler mode (Drain [1972]). If the scattered light is collected in the forward direction, the LDA arrangement is said to be in the forward scatter mode. If on the other hand the scattered light is collected in the backward direction, the arrangement is said to be in the back scatter mode. In the case of former, the detector has to be positioned opposite to the light source, while in the case of latter, the source and the detector are along side each other. The back scatter mode is beneficial in situations where access to the measurement location from opposite sides is not possible. However, as was demonstrated by Adrian and Earley [1975], a higher portion of light is scattered in the forward direction; hence a higher laser power is required for back scatter arrangement. In the present study, the differential doppler mode of LDA, in the forward scatter arrangement was used. 532 Differential doppler mode: In this mode of operation, a single laser beam is split into two beams. The diameter of each beam is defined as the location at which the intensity is e 2 of maximum intensity which occurs at the beam axis. This diameter is referenced by the symbol De-z • The intensity profile across the Figure 5.5 A typical arrangement of laser doppler anemometry in reference beam mode of operation 123 beam diameter is a Gaussian function of the radial distance from the beam axis. The beam converges and then diverges in its path of propagation, with its minimum diameter occurring at the so-called "beam waist"; at this point the light waves are very nearly planar. If the two beams are focussed in such a way that the cross-over occurs at the beam waist, a series of uniformly spaced bright and dark fringes are formed; the fringe spacing df is defined by: (refer to Fig. 5.6.a) d f = — ^ — (5.1) 2 sin 0 where 0 is half angle between the two beams. The region of cross-over is normally called the measuring or scattering volume. This is an ellepsoidal surface with the major and minor axes defined by : (refer to Fig. 5.6.b) l m = - ^ - (5.2.a) sin 0 0^ = - ^ - (5.2.b) cos 9 where de-2 is the beam diameter at the cross over point and is equal to here X. is the wavelength and F is the focal length of the focusing lens. The total number of fringes in the control volume is then : N f r = d f =1-2]^-2 (5.4) df De-2 v ' Figure 5.6 Characteristic dimensions of laser beams and ellepsoidal measuring volume (a) Measuring volume and fringe pattern cross-section (b) Typical beam dimensions 125 The physical interpretation of the LDA technique can be based on the fringe concept. It takes a particle, moving perpendicular to the fringes, t seconds to traverse one set of fringes which are located at distance df apart. The velocity of the particle is then : w = df = f d f = _f_A_ (5.5) 1 2 sin 9 The frequency in the above equation is the doppler frequency which is measured by some suitable means; the fringe spacing is obtained from equation 5.1 and hence the velocity can be calculated from the above equation. A typical LDA signal is shown in Fig. 5.7. The shape of the signal can be explained in terms of variation in intensity of scattered light, as the particle traverses the fringes. Let's consider a particle entering the measuring volume; as it traverses the measuring volume, it goes through a series of light and dark regions. The amplitude of the output signal from the detector is directly proportional to the intensity of the scattered light. When the particle is in the light region, the amplitude of the output signal increases; as the particle moves towards the dark region, the signal amplitude decreases. The amplitude of the whole signal attains a Gaussian envelope which is due to the Gaussian distribution of light intensity across the measuring volume. The above description does not involve the concept of doppler frequency, which is the central theme of the LDA technique; nevertheless it serves as a simple explanation of the LDA technique. The rigorous analysis of the LDA technique which is based on the doppler frequency phenomenon will be discussed in the following section. It is noted that the relation between the frequency and the speed of the particle is the same irrespective of the direction in which the particles traverse the fringes. In some instances it is desirable to know the direction of the motion; for example in the case of separated flows, the point of reattachment can be obtained by locating the position at which the velocity changes direction. The frequency shifting technique is the most common method used for identifying the flow direction. If the frequency of the two beams are identical, the fringes at the point of crossing are stationary. If on the other hand, the frequency of one of the beam is different to that of the other, 126 0.2 _ Figure 5.7 A typical LDA signal representing the combination of pedestal and doppler signals. 127 the fringes move at a speed proportional to the difference between the frequencies (fs). The direction of fringe movement is from the higher frequency to the lower frequency beam. If the particle which enters the measuring volume, moves in the direction opposite to the fringe movement, the output signal at the detector would be at the frequency fd+fs , where fd is the doppler frequency. If the particle is moving in the same direction as the fringes, the resulting frequency would be fd-fs In this way the direction of the flow can be clearly determined. Bragg cells are usually used to attain the required frequency shift in one or both of the beams, as explained by Drain [1972]. Based on the foregoing brief introduction, it can be seen that the velocity measurements can be carried out using a single laser beam, appropriate optics, light detector and a device which measures the frequency of the signal from the detector. The knowledge of frequency, the characteristics of the light source and those of the optics would yield the velocity magnitude by the use of the following equation: w = f d —^- (5.6) 2 sine The essential components required to accomplish such measurements are explained in the following section. 5.3.3 Major components of LDA: All the following components have been manufactured and supplied by TSI. Fig. 5.8 shows the overall arrangement of the system (a) Laser source : The He-Ne laser was used as the light source for the LDV system. The supplied light was at 632.8 nm wavelength and had the characteristics of red light, emitted by Neon gas. The rated power under operating condition was 15 mw. The laser beam was emitted at its TEM mode (Adrian et al. [1981]); at this mode, the intensity profile across the 1 mm diameter beam was very nearly Gaussian. LASER RECEIVING OPTICS PHOTODETECTOR TRANSMITTING OPTICS SIGNAL PROCESSOR Figure 5.8 A schematic of LDA set up in the differential doppler, forward scatter mode of operation 129 (b) Transmitting optics : The laser beam exiting the source was passed through a collimator. The function of the collimator was to position the beam waist (the location of minimum diameter) at one focal point ahead of the focusing optic. This assured the two beams were focussed at their waist which was a necessary condition to attain a uniform fringe spacing at the cross-over point The collimated beam was divided into two beams of equal intensity and 50 mm apart, with a beam splitter. The beam splitter could be rotated at 2° intervals to accommodate velocity measurements at different angles. The two beams were focussed into a single point with a lens of 250 mm focal length and 60 mm clear aperture. With the standard 50 mm beam spacing and the lens focal length of 250 mm, the angle between the two beams at the point of crossing was 11.4°. The measuring volume dimensions were (refer to Fig. 5.6.b): d m = 0.18 mm l m = 1.9 mm The fringe spacing pertaining to this optical arrangement was : d f=3.19um and hence the typical number of fringes across the measuring volume was: N f r = 57 (c) Receiving optics and photo-multiplier: The tight scattered by the particles (traversing the measuring volume) were collected by the receiving optics. These comprised of two lenses of 250 mm focal length; these focussed the collected light onto the tight detector's aperture of 0.2mm size. The receiving lens position was adjustable to insure that it was focussed exacdy on the measuring volume. This was a critical part of the LDV setup; other lights scattered by surroundings and/or reflected by intervening surfaces can also reach the detector which give rise to erroneous 130 signals. These erroneous signals, can be omitted by focussing the receiving lens exactly on the measuring point. The photo-multiplier (light detector) was an RCA vacuum tube with 13 per cent quantum efficiency (see Adrian et al. [1981]). The principle of LDA operation lies in the theory of photo-multiplier non-linear response; the terminology non-linear response stems from the fact that the output of this device is proportional to the intensity of the incident light which itself is proportional to the square of the amplitude of the light wave. It is shown by Adrian et al. [1981] that the output signal from the photo-multiplier is of the form : I = C ((E? + El) + Ei E 2 cos [(f! - f2) t + (q>i - 92)]} (5.7) It can be seen, from the above equation, that the detector output signal contains a constant term (E2 + E2) and a term which includes the difference between the frequency of the two scattered light waves. The first part of the signal is called the pedestal and the second part is called the doppler term. The doppler part of the signal is a measure of the particle velocity. Fig. 5.9.a shows a typical output signal from the photo-multiplier. The pedestal part of the signal is of low frequency and is a measure of the time it takes the particle to traverse the whole measuring volume. The omission of this part of the signal results in the doppler signal shown in Fig. 5.9.b. (d) Frequency shifter: The doppler part of the photo-multiplier signal is proportional to the difference between the frequencies of the two scattered light waves. If one of the original beams is at a higher frequency than the other, the difference in the frequency is superimposed on the doppler frequency. As explained previously, frequency shifting of one of the original beams has a number of advantages; to make use of these advantages, a Bragg cell was used to shift the frequency of one of the beams by 40 MHz. The frequency of the photo-multiplier output signal was then equal to fs+fd. The signal was passed through a down-mixer, which electronically reduced the frequency to a measurable range. Figure 5.9 A typical photo-multiplier out put signal (a) Before filtering (b) After high pass filtering 132 (e) Counter: this unit was used to measure the frequency of the photo-multiplier output signal with a 2 ns resolution clock. This was achieved by counting the period of a pre-selected number of cycles within each signal. The unit consisted of three major modules. These were the input conditioner, timer and data transfer module. The function of input conditioner was to reduce the photo-multiplier signal into a form which could be detected by the timer module. It consisted of a series of filter settings, with the high pass filters used to remove the pedestals and low pass filters used to remove the inherent noise in the signal. The signal was then amplified above a threshold level which fired the Schmitt trigger. This initiated the timer module whose function was to measure the period of the signal. The counter could be operated in two modes. In the continuous mode, the number of cycles were specified and the corresponding period was measured; however if the total number of cycles in each signal exceeded this prescribed number, a new set of measurement was done on the same signal. In this mode, a few measurements per signal could be obtained. In the total burst mode, the counter measured the period of the signal by counting the total number of cycles in the signal; however, the measurement was only validated, if the prescribed number of cycles were exceeded, (refer to Fig. 5.10). In each mode of operation, the counter carried out two sets of time measurements per data point. The period for N cycles and also the period for some multiple of N cycles were measured. If the difference between the two measured periods were within some prescribed tolerance, in this case 1 per cent, the result was validated and a data ready signal was issued. The data was then latched onto a buffer where it was temporarily stored until it was read by an external device (in this case an on line computer). During the process of data transfer, the counter was internally inhibited from updating its buffer. This was achieved by a data hold signal which informed the counter that the data was being read by an external device. The data hold signal was useful for the purpose of hand-shaking between the counter and any external device. For example in the case of transient measurements, the data transfer was inhibited for some desired time, by externally invoking the data hold signal; as a result, though the computer was ready in the data transfer mode, it was not able to collect any data since no data was Figure 5.10 Continuous and total burst mode of operation setting on the TSI counter 134 available at the counter buffer. Once the pre-selected time was exceeded, the data transfer started instantaneously. The data were latched into the counter buffer in a series of three 16 bit words. Word A included the information on the number of cycles, word B included the information on the period of the signal and word C included the information on the time between each data point 5.3.4 Scattering particles The scattering particles were the integral part of the LDA system. The system, in principle, measured the velocity of scattering particles; based on the assumption that these particles closely follow the flow, the fluid velocity was inferred. To satisfy such a criterion, the particles must be light and be of suitable sizes. The particle size were restricted by the size of fringe spacing; if the particles were larger than a typical fringe spacing the associated signal was not validated by the counter. At the same time, as noted by Adrian and Early [1975], larger particles scatter more light; so there should be an optimum particle size at which all requirements are satisfied. The signal intensity was directly influenced by the intensity of scattered light, so particles had to be good light scatterers. Particles concentration was another influential factor. Low data rates were associated with low particle concentration. At the same time a large concentration of particles at any location, resulted in a low data rate; this was because of the overlap of different signals, due to the presence of multiple particles in the measuring volume. In the case of liquid, the naturally occuring particles in the fluid produce adequate signals. For air measurements however, the fluid has to be seeded in order to get any measurable doppler signal. Particle densities are normally 3 orders of magnitude larger than air density; under such a condition the particles might not follow the flow. Particle sizes of 1pm or less have been shown to follow the flow at moderate Reynolds number (Re = 2<105) quite accurately (see Dring [1982] and also Feller and Meyers [1976]). In the case of swirling flows, where there is an extra body force (due to centrifugal effects) on the particles, the problem of particle lag might be more pronounced. Dring and So [1978] studied the particle trajectories in swirling flows. Their study suggested that particles of diameter 135 less than 4 um were not affected by centrifugal forces and hence were expected to closely follow the flow. In the present study, velocity measurements with air as working fluid, were carried out using Zirconium dioxide particles, with specific density of approximately 5. The particle size distribution was measured to be in the range of 0.5 um to 10 um . The large particles were centrifuged out towards the walls of the chamber, thereby reducing the data rate in those locations; smaller particles could also easily escape through the shaft housing and thereby again reduce the data rate. The particle consideration had a significant influence on the choice of the fluid used in the present study. This point is considered in more detail in the following section. 53.5 LDA data collection and analysis It was explained previously that the photo-multiplier output signal, comporised of a low frequency pedestal and a high frequency doppler part. The pedestals were removed by the high pass filter and the inherent noise in the signal was removed by the low pass filter. The importance of the filter setting can thus be appreciated. The filter settings on the TSI counter are hsted in table 5.1. High pass filters 1 3 10 30 100 300 1 3 10 Khz Khz Khz Khz Khz Khz Mhz Mhz Mhz Low pass filters 10 100 1 3 10 30 Khz Khz Mhz Mhz Mhz Mhz Table 5.1 High and low pass filter settings on the TSI counter 136 The large filter bandwidths, made the signal cleaning more difficult. The frequency shifter was used to remove the pedestals completely and also to bring the expected frequency into the range where it could be easily cleaned (e.g. in the 1-3 MHz range). In the case of transient measurements, the velocity magnitudes occupied a large range of frequencies and hence frequency shifter was essential to cover such a large range. Frequency shifting was also used in locating the geometric center of the chamber; this was the location at which the tangential velocity changed sign. Without the shifter, the center could be located by traversing the beam crossing point from one end of the chamber to the other and thereby locating the center. With the frequency shifter, this method of locating the center was double checked by noting the change in the flow direction. Biasing is one source of error associated with the LDA measurements. Velocity biasing is attributed to the higher probability of faster particles entering the measuring volume and thereby producing more data points. The mean velocity obtained by a simple arithmetical averaging would be biased towards the higher velocity points. The biasing error can be reduced by means of appropriate weighting of the data points. Maclaughlin and Tiederman [1973] suggested that for one-dimensional flows, the reciprocal of the velocity magnitude is a good weighting function. Based on this suggestion, the average velocity can be obtained by the following equation: I-(WT) w = M _ ( 5 . 8 ) y Wj i=i The more general form of velocity biasing correction can be expressed as [Dimotakis [1976]): N X w i te w = i ^ (5.9) S te i=l 137 where ts is the total time taken by the particle to traverse the measuring volume; this is the time measured by the counter in its total burst mode of operation which was previously explained; both methods of biasing correction reduce to the same form in the continuous mode of operation. The mean velocity with and without biasing corrections were obtained for two sample cases in the present study. The results revealed a 0.5 per cent difference between the two methods of averaging (refer to the next Chapter for further details). Due to the random nature of LDA signals, a large collections of data points17 are expected to have a normal distribution. The FIND software (provided by TSI) was used to generate the frequency plot for each set of data points. The skewness and flatness factors for each set were calculated and compared against the corresponding values for normal distribution18 (see Fig. 5.11). If the agreement were close, The data points were accepted and recorded for further processing; otherwise the set was rejected and the measurement was repeated until satisfactory conditions were attained. 5.4 Choice of the working fluid Water and air were considered as two choices of working fluid in the present investigation. In deciding on the right choice, three major parameters were taken into considerations; these were the need for scattering particles, the problem of sealing and the problem of motor stoppage for transient measurements. It was previously explained that water, in the context of LDA measurements, was considered to be naturally seeded. In the case of air however, seeding particles had to be introduced into the flow for the purpose of LDA measurements. Because the fluid was enclosed within the chamber, the particles had to be deposited in the chamber prior to any measurements. Also with air as working fluid, high motor rpm was necessary to attain a reasonable flow Reynolds number (20000 rpm was the maximum speed considered); this would necessitate the use 1 7 5000 d a t a p o i n t s were c o l l e c t e d a t e a c h m e a s u r i n g l o c a t i o n i n t h e p r e s e n t s t u d y . 1 8 F o r n o r m a l d i s t r i b u t i o n skewness f a c t o r i s z e r o a n d f l a t n e s s f a c t o r i s 3 . 138 Filenue: c:\flrfcsfcir\Mter7B.sfS I Processor*: 1 Mode: fUadon , ^ j I i 1 f 1.868 1,688 2.286 2.868 3.488 4.888 Velocity (n/s) for Processor 1 Position (cul/in) (8, 488, ) Velocity ri cursor = 1.888 Velocity wan - 2.386 Percent at cursor - 8.88 <F1> Help screei for bistogrM Figure 5.11 A typical frequency distribution of 5000 velocity data points collected by the LDA technique. 139 of labrynth seal. Small seeding particles could easily leak through this type of a seal. It was also found that, after large number of test cases, the particles damaged the seal and had to be eventually replaced. However, with water as a working fluid, the same Reynolds number could be attained with a much smaller motor rpm19. At low motor speeds, the use of a normal water seal proved to be satisfactory. This provided a complete seal between the chamber and the surroundings. In the case of transient measurements, it was ideal to stop the disk rotation instantaneously. The technique used in the present study to achieve such a goal was to stop the motor (by short-circuting) rather than stopping the disk directly. With water as working fluid, the short-circuting of motor, stopped the disk rotation in 100 ms; a much higher delay time was anticipated for air flow measurements. This is because to attain a representable Reynolds number, with air as a working fluid, the motor has to run at a much higher rotational speed; hence because of the higher motor inertia, it takes a longer time for it to come to complete rest Based on the above three considerations, it was decided to use water as working fluid and the corresponding flow characteristics for air were inferred based on Reynolds number similarity. The Reynolds number similarity was experimentally confirmed, as will be shown in the following Chapter. The total error associated with the measurements presented in the following chapter were calculated to be within 12 per cent. It is shown in Appendix G that 2 per cent of this error is due to the LDA technique and the rest is attributed to experimental uncertainty (eg. repeatability). maximum s p e e d o f 2190 rpm was u s e d i n t h e p r e s e n t s t u d y 140 CHAPTER 6 Experimental Results and Discussions 6.1. General It was pointed out in the previous chapter that water was found to be a better choice for working fluid as compared with air. This choice necessitated the establishment of Reynolds number similarity between the air and the corresponding water measurements. In the present Chapter the experimental results establishing the Reynolds number similarity criterion are presented. This is followed by presenting the experimental results on the effect of aspect ratio on the mean tangential velocity and turbulence intensity for steady-state rotation. The accuracy of the measurements is then examined against the numerical calculation of the flow field. The effect of disc surface roughness on the velocity and turbulence intensity profiles for steady-state rotation is also presented in this chapter. Lastly the experimental results and the corresponding numerical calculations pertaining to decaying flow field in short chambers of different aspect ratios are presented. 6.2 Reynolds number similarity test Owing to the difficulties encountered in measurements of air velocity and turbulence intensity field, generated by a rotating disc inside a chamber, it was decided to use water as working fluid and use Reynolds number similarity hypothesis to infer the corresponding air measurements. The Reynolds number similarity was established by comparing the air measurement with the corresponding water measurement at the same Reynolds number. The air measurement was carried out at one flow Reynolds number in a chamber with aspect ratio equal to 1/2. The results presented in this section, pertain to the disc rotational speed of 10840 rpm; this rotational speed corresponds to Reynolds number 0^-) equal to 224 x 103 . The difficulties encountered in the air measurements were mostly due to the lack of seeding particles in the chamber. Originally it was thought that prior deposition of seeding particles (Zi 02) in the chamber would produce a reasonable data rate. However, it was found that the particles escaped the channel, through the shaft housing, before a reasonable number of data points could be 141 collected. The alternative was to introduce seeding particle through a radial port while the disk is rotating; because the particles were small and light, they followed the flow once they were introduced into the chamber. Seeding particles were added whenever the data rate slowed down. The disc rotational speed at which the same Reynolds number with water as working fluid was attained corresponded to 610 rpm. Figure 6.1 shows the comparison between the radial profiles of mean dimensionless tangential velocity at the mid-plane of the chamber with water and air as working fluids. The agreement between the two profiles was within 15 per cent which was within experimental error20 ; this established the Reynolds number similarity between the two flows. One interesting feature exhibited by both profiles is the sudden change in the slope of the velocity profile at dimensionless radius equal to 0.4. The solid line drawn through the experimental data points clearly shows this change in the slope. This feature will be examined in more detail in the following sections. 6.3 Steady-state measurements with a smooth disc The effects of rotational speed and the chamber aspect ratio on the flow field generated by a smooth disc were investigated. Three rotational speeds of 1080, 1700 and 2090 rpm in three chambers with aspect ratios equal to 1/2, 3/10 and 1/10 were considered. Figure 6.2 shows the radial profile of dimensionless tangential velocity at the mid-plane of chamber with aspect ratio equal to 1/2, for different disc rotational speeds. The maximum speed attained by the bulk of fluid was 24 per cent of the tip speed, occuring at the maximum disc rotational speed. The sudden change in the slope of velocity profile is again apparent at the lowest disc rotational speed. Figure 6.3 shows the corresponding dimensionless rms velocity profile. It can be seen that the profile was reasonably uniform near the center and increased towards the wall of the chamber. Figure 6.4 shows the dimensionless velocity profile at the mid-plane of the chamber with aspect ratio equal to 3/10 for different disc rotational speeds; the maximum velocity was 30 per cent 2 0 R e f e r t o A p p e n d i x G 142 0.20 0.18-R Figure 6.1 Dimensionless tangential velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/2 • Air as working fluid (disc rpm = 10840) A Water as working fluid (disc rpm = 610) 143 • A • • • A • • • ft • a • i • i ft • 5 0 0.2 0.4 0.6 0.8 1 r_ R Figure 6.2 Dimensionless tangential velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/2, water as working fluid. • disc rotational speed = 1080 rpm A disc rotational speed = 1700 rpm • disc rotational speed = 2090 rpm 144 O.OIOi 0.009-0 . 0 0 8 -0.007 0 .006 -w W 0.005 0 .004-0 . 0 0 3 -0 .002-0 .001-0 .000 -• 0... A A , , MM . . . M » km • • 1 • • • • : A A • • • • • A g A A -0.2 0.4 0.6 0.8 Figure 6.3 IJimensionless rms velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/2, water as working fluid. • disc rotational speed = 1080 rpm A disc rotational speed = 1700 rpm • disc rotational speed = 2090 rpm 145 0.35 0.30-0.25 0.20 W W 0.15 0.10-0.05-0.00-• A • • A • • • A J • • • • • 1 • 1 1 0.2 0.4 0.6 0.8 1 r R Figure 6.4 Dimensionless tangential velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 3/10, water as working fluid. • disc rotational speed = 1080 rpm A disc rotational speed = 1700 rpm • disc rotational speed = 2090 rpm 146 of the tip speed which again occured at the highest rotational speed. The change in the slope of velocity profile was again present; however, the location at which the change occured had shifted towards the wall |^ -= 0.8). Figure 6.5 shows the corresponding dimensionless rms velocity profiles. Again the turbulence intensity tended to be uniform near the center and increased towards the wall. Figure 6.6 shows the dimensionless velocity at the mid-plane of the chamber with aspect ratio equal to 1/10 for different disc rotational speeds. The maximum fluid velocity in this case was 40 per cent of the tip speed. The change in slope of velocity profiles appeared at all three rotational speeds for this aspect ratio. Figure 6.7 shows the corresponding dimensionless rms velocity profiles; same observation as with the previous two rms velocity profiles can be made. Comparison of these dimensionless velocity and turbulence intensity profiles pertaining to different disc rotational speeds, in chambers of different aspect ratios reveals the following features: (a) The maximum measured mean tangential velocity for the cases considered, always occured at the highest disc rotational speed. The maximum velocity increased with the decrease in the chamber length (smaller aspect ratio). This was because at the smallest aspect ratio, the slowing down effect of the chamber curved walls on the bulk flow was least (refer to Appendix H). (b) The velocity profiles tended to be almost Reynolds number independent at higher rotational speeds. The Reynolds number independency was most apparent at the smallest aspect ratio. (c) The turbulence intensity appeared quite uniform near the center of the chamber and increased towards the wall. This was because the mean flow was stable (solid body type of rotation) near the central region of the chamber, and was unstable near the walls of the chamber. (d) The turbulence intensity level was small; the maximum intensity was 0.85 per cent of the disc tip speed . (e) There was a change in the slope of velocity profile at all three investigated aspect ratios. The sudden change in the slope occured at lower speeds for aspect ratios of 1/2 and 3/10; the change in the slope was apparent at all three rotational speeds for aspect ratio of 1/10. 147 0.010 0.009-0.008 0.007 0.006 W 0.005 0.004-0.003-0.002-0.001-0.000-A A 0.2 0.4 0.6 0.8 r R Figure 6.5 Dimensionless rms velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 3/10, water as working fluid. A disc rotational speed = 1080 rpm • disc rotational speed = 2090 rpm 148 Figure 6 . 6 Dimensionless tangential velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/10, water as working fluid. • disc rotational speed = 1080 rpm A disc rotational speed = 1700 rpm # disc rotational speed = 2090 rpm 149 0.010-1 0 . 0 0 9 -0 . 0 0 8 0 . 0 0 7 -0 . 0 0 6 -w W 0 . 0 0 5 -0 . 0 0 4 -0 . 0 0 3 0 . 0 0 2 -0 . 0 0 1 -0 . 0 0 0 -1 1 • A : • A A A • •;• « 1 0.2 0.4 0.6 r R 0.8 Figure 6.7 Dimensionless rms velocity profile generated by a smooth disc at mid-plane of a chamber with aspect ratio equal to 1/10, water as working fluid. • disc rotational speed = 1080 rpm A disc rotational speed = 1700 rpm • disc rotational speed = 2090 rpm 150 This change in the slope was attributed to the local transition of the boundary layer from laminar to turbulent. It appeared that the transition was very much aspect-ratio dependent. It was observed that as the aspect ratio was increased, the transition occured at a lower rotational speed. The same conclusion was made by Daily and Neece [1960] in their study of chamber dimension effects on induced flow and frictional resistance of enclosed rotating discs. This phenomenon will be explored in more detail in the following section. At this stage it was decided to numerically calculate the flow-field generated by a smooth rotating disc inside a chamber at three different aspect ratios. It was decided to limit the calculation to the highest rotational speed of 2090 rpm. The Wilcox and Chambers model was used to calculate the turbulent quantities. The solution domain was divided into 56 by 25 grid points in the x and r directions, respectively. Figure 6.8 shows the comparison between the experimental and the corresponding calculated dimensionless velocity profile in a chamber of aspect ratio equal to 1/2. Excellent agreement was observed with the maximum error less than 5 per cent. In the shorter chamber (aspect ratio = 3/10) it was necessary to change the number of grid points, in order to make sure that the nodal point adjacent to the boundary was placed at a dimensionless distance (X + or Y + ) greater than 30. Figure 6.9 shows the comparison between the experimental and the corresponding calculated dimensionless velocity profiles in this chamber. The number of grids used for this calculation were 36 by 30 in the x and r directions, respectively. Good agreement was observed with the maximum error of 3 per cent cnxurring near the wall. The same calculation was carried out for the chamber of aspect ratio equal to 1/10. The number of grid points for this calculation were changed to 12 by 35 in the axial and radial directions, respectively. It should be pointed out that the experimentally measured velocity profiles in this chamber exhibited the sudden change in slope at all three rotational speeds considered. Figure 6.10 shows the comparison between the calculated and measured velocity profile. The disagreement between the two profiles was significant. The disagreement was attributed to the speculation that the flow in this chamber undergoes a local laminar to turbulent transition; the numerical model used in this study, at its present form, 151 0.30 0.25-0.20-W W 0.15-0.10-0.05-0.00 r R Figure 6.8 The comparison between experimental and calculated velocity profiles in a chamber of aspect ratio equal to 1/2 using Wilcox and Chambers model (rpm = 2090) • Experimental results — Numerical results 152 0.40 0.35-Figure 6.9 The comparison between experimental and calculated velocity profiles in a chamber of aspect ratio equal to 3/10 using Wilcox and Chambers model (rpm = 2090) • Experimental results Numerical results 153 Figure 6.10 The comparison between experimental and calculated velocity profiles in a chamber of aspect ratio equal to 1/10 using Wilcox and Chambers model (rpm = 2090) • Experimental results Numerical results 154 was not equipped to account for laminar to turbulent transition and hence the disagreement between the experimental and calculated velocity profiles. In order to investigate the validity of such speculation, it was decide to promote turbulence in the flow by roughening the disc surface. The following section summarizes the results pertaining to the flow generated by this roughened disc. 6.4 Steadv-state measurements with roughened disc In order to keep the flow generated by the disc axisymmetric, the disc surface was uniformly roughened by fastening a thin disc (cut out of carrot grating) to the original smooth disc. The roughness element height on the grating was measured to be approximately 1.5 mm. The measurement was limited to two rotational speeds of 2090 and 1180 rpm. The measurements were carried out inside chambers of aspect ratio equal to 1/2, 3/10 and 1/10. Figure 6.11 shows the dimensionless tangential velocity generated by roughened disc, rotating inside a chamber with aspect ratio equal to 1/2. The change in the slope of velocity which occurred with the smooth disc rotating at 1180 rpm did not appear in the present case. The maximum speed attained by the fluid was 50 per cent of the tip speed which was two times that with the smooth disc. The velocity profiles tended to be Reynolds number independent. Figure 6.12 shows the corresponding dimensionless turbulence intensity profiles. These tended to be uniform near the center and increased towards the wall. The maximum intensity in this case was three times larger than that obtained with the smooth disc. Figures 6.13 and 6.14 show the dimensionless tangential velocity and turbulence intensity at two rotational speeds, generated by roughened disc inside a chamber of aspect ratio equal to 3/10. The same observations as in the previous case (aspect ratio = 1/2) on the velocity and turbulence intensity field could be made. The maximum speed was almost 2 times that generated by a smooth disc; no change in the slope of the velocity profile could be seen. Figure 6.15 shows the dimensionless velocity profile generated by a roughened disc surface in a chamber with aspect ratio equal to 1/10. The velocity profile exhibited a solid body rotation with no change in slope. The profiles were Reynolds number independent and attained a maximum value 1.5 times of that with a smooth disc. It was thus concluded that the sudden change 155 Figure 6.11 Dimensionless tangential velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 1/2, water as working fluid. • disc rotational speed =1180 rpm A disc rotational speed = 2090 rpm 156 0.030-1 0.025-0.020-w W 0.015-0.010 0.005-0.000-• • i 1 1 r— 0 0.2 0.4 0.6 0.8 r R Figure 6.12 Dimensionless rms velocity profile generated by a rough disc at mid-plane chamber with aspect ratio equal to 1/2, water as working fluid. A disc rotational speed = 1180 rpm • disc rotational speed = 2090 rpm 157 Figure 6.13 Dimensionless tangential velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 3/10, water as working fluid. • disc rotational speed = 1180 rpm A disc rotational speed = 2090 rpm 158 o.oi8n 0.016-0.014 0.012-0.010-w W 0.008-0.006-0 . 0 0 4 -0 . 0 0 2 -0.000-• A * 0.2 I I 0.4 0.6 r R i 0.8 Figure 6.14 Dimensionless rms velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 3/10, water as working fluid. • disc rotational speed = 1180 rpm A disc rotational speed = 2090 rpm 159 Figure 6.15 Dimensionless tangential velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 1/10, water as working fluid. • disc rotational speed = 1180 rpm A disc rotational speed = 2090 rpm 160 0.025-0.020-0.015-w W 0.010-0.005-0.000 + 0 . 2 — I 1— 0.4 0.6 r R o.e Figure 6.16 Dimensionless rms velocity profile generated by a rough disc at mid-plane of a chamber with aspect ratio equal to 1/10, water as working fluid. • disc rotational speed = 1180 rpm A disc rotational speed = 2090 rpm 161 in the slope of velocity profile was due to local transition and could be overcome by promoting turbulence. It was explained in Chapter 1 that, because the axial pressure gradient is negligible, the tangential velocity is expected to be uniform across the chamber in the axial direction, except in thin regions near the walls. To examine the validity of this point, the tangential velocity across the chamber at the mid-radius was measured. Figure 6.17 shows the axial profile of the tangential velocity generated by a rough disc rotating at 1180 rpm inside two chambers with aspect ratios equal to 1/2 and 1/10. The profile was very nearly uniform in the shorter chamber. In the longer chamber, the velocity tended to decrease from the stator (stationary wall on the left) to the rotor (rotating disc on the right). This decrease was attributed to the secondary flow pattern which transferred the high momentum fluid close to the curved wall towards the stationary wall and at the same time transferred the low momentum fluid near the center towards the rotating disc; hence giving rise to the profile shown in Fig 6.17. The experimental results were compared with the flow field calculated by the numerical model. It was explained in Chapter 2 that in the case of rough walls, the logarithmic profile used in the law of the wall incorporates a different constant. The value of this constant was found to be equal to -8.5 which agreed with the value proposed by Schlichting (as explained in Chapter 2). The comparison was limited to only two aspect ratios, namely 1/2 and 1/10. Figure 6.18 shows the comparison between the experimentally measured dimensionless velocity in a chamber with aspect ratio equal to 1/2 and the corresponding calculated profile. Good agreement was observed near the center, and deteriorated towards the wall with the maximum difference of 8 per cent. Figure 6.19 shows the corresponding comparison between the calculated and measured turbulence intensity profiles. It should be pointed out that the Wilcox and Chambers model calculates the normal intensity component; so from Fig.. 6.19 it can be concluded the turbulence field was isotropic near the center of the chamber and became non-isotropic in the vicinity of the wall. Figure 6.20 represents the comparison between the experimentally measured velocity and the corresponding calculated profile in a chamber with aspect ratio equal to 1/10. The agreement 162 W 0.35 • • • • • ! B - _ A :A A A " " 1 . A A A A A • • • • • • • a i • . u - l 1 1 1 1 I I -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 X Figure 6.17 Variation of mid-radius dimensionless tangential velocity profile across the chamber, with a rough disc rotating at 1180 rpm. • Aspect ratio L/D = 1/2 A Aspect ratio L/D = 1/10 163 Figure 6.18 Comparison between the calculated and the corresponding exr>erimental dimensionless velocity profile using a rough disc in a chamber of aspect ratio equal to 1/2. — Calculated profile • Measured profile 164 0.03 0.02 w'/W 0.01 0.00 < i i E • • | » • • * • I ! ! i 1 ! i t i 1 1 ..... ........... 0.0 0.2 0.4 0.6 0.8 1.0 r / R Figure 6.19 Comparison between the calculated and the experimental dimensionless turbulence intensity profile using a rough disc in a chamber of aspect ratio equal to 1/2. Calculated profile • Measured profile 165 0.8-0.7- • • / 0.6-• / 0.5-w 0 4 -w • 0.3-• 0.2-/m 0.1- /m 0.0-( ) I 0.2 0.4 r R 0.6 0.8 I Figure 6.20 Comparison between the calculated and the corresponding experimental dimensionless velocity profile using a rough disc in a chamber of aspect ratio equal to 1/10. Calculated profile • Measured profile 166 between the two profiles is not as good as the previous aspect ratio. However, a substantial improvement as compared to the case of smooth wall calculation has been attained, providing further evidence for the speculation that the disagreement between the calculation and measurement was due to the laminar to turbulent transition phenomenon. Figure 6.21 shows the corresponding comparison between the calculated and measured intensity profiles. The disagreement near the wall was attributed to the non-isotropic nature of turbulence in that vicinity. 6.5 Transient measurements in short chambers In this section the velocity and turbulence intensity measurements pertaining to decaying flow fields in cylindrical chambers of different aspect ratios are presented. The aim of this study was to investigate the effect of aspect ratio on the rate of decay of flow. Two aspect ratios, namely 1/2 and 1/10 were considered. In each case, the steady-state flow field was generated by a rotating rough disc. The on-line data aquisition program (purchased from TSI) was used to monitor the mean tangential velocity at the mid-plane in real time. Once the steady-state velocity at each measuring radial location reached the value reported in the previous section, the rotation was stopped. This was achieved by short-circuting the motor which drove the shaft. Instantaneous motor stoppage could not be achieved, owing to the inertia of the motor armature and the disc. Two choices of data collection schemes were available. In one scheme the collection of data commenced immediately after the motor was short-circuted. The second choice was to externally inhibit the TSI counter from collecting data for a predetermined adjustable time delay. In the former case, the initial measured decay of flow could be used to approximate the time taken by the disc to completely come to rest, once the motor is switched off. In the latter case, this portion of data which pertains to initial decay (disc speed is decelerating to zero) could be omitted by inhibiting the counter during this time period. It was decided to adopt the first scheme of data collection. The post processing of measured velocity field data showed a very slow decay in the first 100 ms, followed by a much faster decay. 167 r/K Figure 6.21 Comparison between the calculated and the experimental dimensionless turbulence intensity profiles using a rough disc in a chamber of aspect ratio equal to 1/10. Calculated profile • Measured profile 168 It was thus concluded that the disc gradually comes to rest in the first 100 ms after the motor has been switched off. Data were collected at seven radial locations in each of the two chambers (L/D=l/2 and L/D=l/10). At each location, ten batches of 5000 data points were collected to ensure an accurate averaging process. The data reduction process consisted of extracting the decay of mean velocity and turbulence intensity at different radial locations from the raw collected data. This procedure can be summarized by the following steps: (a) The best mean velocity curve for each of the ten data files, at each of the seven locations, was determined over a time span of 100 to 600 ms after the motor had been switched off. The exponential curve of the form We(t)=exp ( a + bt + ct2 + dt3 ) was found to be a good fit to the raw data (see Fig. 6.22) (b) The total time span was then divided into 20 equal 25 ms time windows and the mean velocity at the mid-point of each window was obtained (Wej). The rms velocity (WnnS) at these mid-points (for each of the 10 data files) were then calculated by the following equation: where N corresponds to the total number of data points in the k* time window interval. (d) The mid-point mean and rms velocities for each window were then ensemble-averaged over the 10 data files, using the following equations: (6.1) m (W)en s k n T l W ^ (6.2) m (6.3) j=i 169 Veloc i ty (m/s) • A l A A2 O A3 A A4 • A5 + A6 K A7 a A8 O A9 O A10 MEAN Time (s) Fieure 6 22 The mean tangential velocity profile at dimensionless radial location of 0.9 in a chamber with aspect ratio equal to 1/2 at different time interval after the motor has been stopped. 170 These are the ensemble averaged mean and rms velocities at the mid-point of the k* window; m is equal to the total number of data files per location. As an example, figure 6.22 shows 10 mean velocity curves at dimensionless radius equal to 0.9 in a chamber of aspect ratio equal to 1/2 over the time span of 600 ms after the motor has been stopped. The corresponding ensemble averaged mean velocity curve at mid-point of 25 ms time windows has also been shown in this figure. Following the averaging procedure outlined above, the radial profiles of mean tangential velocity and turbulence intensity at different time intervals after the motor has been stopped, have been obtained. Figure 6.23 shows the dimensionless velocity profiles pertaining to a decaying flow field inside a chamber with aspect ratio equal to 1/2. It was pointed out that the flow decays slowly in the first 100 ms, due to the fact that the disc stoppage is not instantaneous. The figure shows the decay from 100ms to 565 ms after the disc has completely stopped. By comparing the steady state measured profile and the profile after 112.5ms decay, it can be seen that the peak velocity has decayed by only 2 per cent. It can also be seen that the maximum velocity decays by 12.5 per cent in the first 50 ms after the disc has come to complete rest. The overall decay in the 560 ms is as much as 64 per cent. Figure 6.24 shows the corresponding decay of turbulence intensity profile in the tangential direction. It appears that the maximum drop in the intensity occures in the 300 to 400 ms interval. The intensity profile tended to be uniform near the center and increased towards the wall in the first 350 ms. The profile flattened out from 400 ms onwards. Figure 6.25 shows the decay of dimensionless tangential velocity profile in a chamber of aspect ratio equal to 1/10 at different time intervals. It is interesting to note that the peak velocity has dropped by as much as 33 per cent in the first 50 ms time interval. Also it should be pointed out that the peak velocity has decayed by 40 per cent from its steady-state value in the first 100 ms. Overall it can be seen that the decay is much faster in the shorter chamber (L/D = 1/10). The faster decay rate associated with the smaller aspect ratio chamber was attributed to the stronger secondary flow pattern. The secondary flow tend to transfer the slower moving fluid near 171 0.6 0.5' 0.4 03 • 0.2 0.1 - -• A 0.0 + x • — t - 0 ms • t-112.5 ms A t-162.5 ms • t = 212.5 ms • t •» 262.5 ms A t- 312.5 ms + t - 362.5 ms x t-412.5 ms o t = 462.5 ms • t = 562.5 ms 0.0 0.2 0.4 0.6 0.8 1.0 r/R Figure 6.23 The radial profile of mean tangential velocity in a chamber with aspect ratio equal to 1/2 at different time intervals after the disc is stopped. 172 Figure 6.24 The radial profde of turbulence intensity in a chamber with aspect ratio equal to 1/2 at different time intervals after the disc is stopped. 173 0.8' 0.6 • 0.4 - •• 0.2 • A 0.0 0.0 0.2 0.4 i X o 1 •H J 0.6 t - 0 ms • t-112.5ms A t = 162.5 ms Q t-212.5ms + t - 262.5 ms x t = 312.5 ms o t - 562.5 ms 0.8 1.0 r/R Figure 6.25 The radial profile of mean tangential velocity in a chamber with aspect ratio equal to 1/10 at different time intervals after the disc is stopped. 174 the wall to the center and this gives rise to a thinner boundary layer at the wall. As a result, the wall shear stress would be larger and consequendy the decay rate will be faster. Figure 6.26 shows the decay of turbulence intensity in a chamber with aspect ratio equal to 1/10. It can be seen that the decay of turbulence intensity is much faster in the shorted chamber.The peak intensity near the wall decayed by as much as 47 per cent in the first 100 ms after the disc has completely come to rest (i.e. between 112.5 to 212.5 ms time interval). The faster decay of turbulence intensity in a shorter chamber is consistent with the suggestion that the dissipation length scale is proportional to the length of the chamber. Consequently the dissipation length scale is smaller in a shorter chamber which gives rise to a higher dissipation rate of mean turbulence kinetic energy. At this stage, it was decided to use the numerical model to calculate the experimentally measured decay in the two chambers considered (L/D = 1/2 and L/D = 1/10). It should be pointed out that one of the objectives of the steady-state calculation (which was reported in the previous section) was to be able to define a complete initial condition for all the variables prior to transient calculation. The steady-state calculation produced satisfactory results and hence it could be concluded that the initial conditions are completely known. However, these calculated initial conditions are only useful if the disc can be stopped instantaneously; this condition could not be achieved experimentally. For the purpose of comparison however, it was decided to start the transient calculation with the measured velocity and turbulence intensity profiles after the disc has completely come to rest (i.e. after the first 112.5 ms). This is hereafter referred to as time zero. The effect of incomplete knowledge of initial conditions on the calculation of rate of decay of mean velocity and turbulence intensity was investigated by adopting two different approaches and then comparing the corresponding results. The first approach consisted of initiating the calculation using the experimentally measured velocity and turbulence intensity profiles after the first 112.5 ms decay; the axial and radial velocities were initially set to zero. In the second approach, the axial and radial velocities calculated after 5 ms,by the first approach, were used as initial profiles in conjunction with the measured velocity and turbulence intensity profiles after the first 112.5 ms decay. The difference between the calculated velocity and turbulence intensity 175 Figure 6.26 The radial profile of mean tangential velocity in a chamber with aspect ratio equal to 1/10 at different time intervals after the disc is stopped. 176 profiles after 50 ms decay was less than 0.5 per cent using the two approaches; this observation was true for both aspect ratio calculations (i.e. 1/2 and 1/10). It was thus concluded that the incomplete knowledge of initial conditions does not influence the subsequent calculation of decay rate. Figure 6.27 shows the comparison between the measured and calculated velocity profiles, in the chamber with aspect ratio equal to 1/2, after 50 and 100 ms decay. Good agreement was observed near the center of the chamber. The maximum difference of 8 per cent was observed near the wall of the chamber. Figure 6.28 shows the corresponding comparison between the measured and calculated turbulence intensity profiles after 50 and 100 ms decay. It can be seen that the experimentally observed quasi-steady nature of the turbulence field was reasonably predicted by the numerical model. The difference between the calculated and the corresponding measured intensity profiles near the wall was again attributed to the non-isotropic turbulence in that vicinity. The transient calculation was also carried out in the chamber with aspect ratio equal to 1/10. The comparison between the predicted and measured velocity profiles after 50 and 100 ms decay is shown in Figure 6.29. The difference between the predicted and measured decay rate increased to 15 per cent in the present case. It is speculated that this source of error could be due to the inaccurate representation of initial velocity and turbulence intensity profiles near the walls of the chamber. These profiles were obtained by smooth curve fit to the experimental results. However, the experimental velocity profile did not reveal any information about the thickness of boundary layer on the curved of the chamber; the existence of boudary layer on the walls of the chamber was thus ignored in the specification of initial conditions for transient calculation. The other source of error could be the under-estimation of roughness height which results in lower wall shear stress value and hence underprediction of decay rate. Figure 6.30 shows the corresponding comparison between the calculated and measured turbulence intensity profiles after 50 and 100 ms decay. The overall decay rate near the center of the chamber was reasonably predicted. Near the wall however, the direct comparison between the calculated and measured profiles is not adequate due to non-isotropic turbulence. 177 Figure 6.27 The comparison between the experimentally measured velocity profile and the corresponding calculated profile (using Wilcox and Chambers model) after 50 ms and 100 ms decay (chamber aspect ratio = 1/2). • Initial velocity profile • Experimental results after 50 ms Numerical results after 50 ms • Experimental results after 100 ms Numerical results after 100 ms 178 0.04 r/R Figure 6.28 The comparison between die experimentally measured turbulence intensity profile and the corresponding calculated profile (using Wilcox and Chambers model) after 50 ms and 100 ms decay (chamber aspect ratio = 1/2). • Initial turbulence intensity profile A Experimental results after 50 ms Numerical results after 50 ms o Experimental results after 100 ms Numerical results after 100 ms 179 0.6 0.5 0.4 WW 0 3 0.2 1 • i • 1 • / A t \ • -» { ^f.. * o m - o • y 0'' y 0.0 0.2 0.4 0.6 0.8 1.0 Figure 6.29 The comparison between the experimentally measured velocity profile and the corresponding calculated profile (using Wilcox and Chambers model) after 50 ms and 100 ms decay (chamber aspect ratio = 1/10). • Initial velocity profile A Experimental results after 50 ms Numerical results after 50 ms o Experimental results after 100 ms Numerical results after 100 ms 180 Figure 6.30 The comparison between the experimentally measured turbulence intensity profile and the corresponding calculated profile (using Wilcox and Chambers model) after 50 ms and 100 ms decay (chamber aspect ratio = 1/10). • Initial turbulence intensity profile A Experimental results after 50 ms — Numerical results after 50 ms o Experimental results after 100 ms Numerical results after 100 ms 181 The difference between the calculated and measured flow field decay rate can be minimized by initiating the transient calculation with a more accurate representation of flow field. This can only be achieved by using the calculated steady-state results, the accuracy of which has been established by direct comparison with the measurements. Unfortunately this condition was not met in the present calculation, due to the difficulty in achieving an instantaneous disc stoppage in the course of transient measurements. 182 C H A P T E R 7. C O N C L U S I O N S A N D R E C O M M E N D A T I O N S The main objective of the present study was to provide some insight on the behaviour of decaying, turbulent swirling flow in short cylindrical chambers. Of particular interest were the effects of chamber aspect ratio on the rate of decay of mean swirl velocity and turbulence intensity; also of interest were the effects of streamline curvature on turbulence. The effect of chamber aspect ratio was investigated through experimental study of swirling flow in chambers of different length to diameter ratio. A numerical model was also developed to predict the transient behaviour of such flow field. The following general conclusions have been drawn in the course of this study. The experimental study of steady-state swirling flow field, generated by a rotating disc inside a closed chamber, revealed that a solid body type of rotation can be generated by a uniformly roughened disc. The maximum speed attained by the core fluid was a function of chamber aspect ratio; the maximum speed increased as the aspect ratio was reduced. It was also observed that the radial distribution of the mean tangential velocity was axially invariant in the shortest chamber. The tangential component of turbulence intensity was relatively low and uniform near the center of the chamber and increased towards the chamber wall. This was because turbulence was mainly generated by the shear action at the walls. The study of decay of swirling flow in these chambers revealed that the rate of decay was directly influenced by the chamber aspect ratio. It was observed that the mean flow decayed much faster in the shorter chamber. This behaviour was attributed to the stronger secondary flows in the shorter chamber which gave rise to a thinner boundary layer on the curved wall and consequently larger wall shear stress and hence faster decay rate. It was thus concluded that for time scales typical of engine combustion, the decay of swirl and turbulence is substantial. A numerical model of decaying, turbulent swirling flow was developed to provide a tool for predicting the flow behaviour. Based on the calculation of steady flow in annuli of two concentric cylinders, either of which could be rotating, it was concluded that streamline curvature substantially affects the turbulent shear stress and velocity field. The turbulence model should account for such effects. It was found that Wilcox and Chambers mixing energy model can account 183 for streamline curvature effects. It was also found that ad-hoc modification of the standard two-equation turbulence model can also account for streamline curvature, with the disadvantage that the extra curvature related term is accompanied by an empirical "case-tuned" constant. It was thus concluded that Wilcox and Chambers model should be used to predict two-dimensional turbulent swirling flows. Based on such a conclusion, the model was used to predict the decay of swirling velocity and turbulence intensity in a short cylinder. The decay of flow field reported by Inoue et al. was successfully predicted, with the provision that the initial dissipation should not be over-predicted. Based on the comparison between the turbulence intensity profile calculated by the Wilcox and Chambers model and that reported by Inoue et al., it was concluded that swirl promoted non-isotropic turbulence. The same model was used to predict the steady-state flow field generated by a roughened disc. It was concluded that the model can successfully predict this flow field, provided that the surface roughness effect is adequately accounted for. The transient calculation of this flow field was only possible by using the experimentally measured profiles after the first 100 ms decay as an initial condition. The overall agreement between the measured and predicted decay rate was fair. It was thus concluded that the numerical model, developed in the present study, can be successfully applied in predicting the effects of aspect ratio on the rate of decay of mean flow and turbulence intensity field in short cylindrical chambers. It should be pointed out that the model has not been tested for fully three dimensional flows with curved streamline. Also the model at its present form is not suitable for predicting the relaminarization of stable flows; the relarninarization prediction would require the implementation of the so called low Reynolds number modelling (refer to Wilcox and Chambers). Based on the observations in the course of present study, the following areas of further research are recommended: (i) To provide a direct comparison between the experimentally measured decay rate of swirling flow field reported in the present study and that predicted, it is recommended that a mechanism be developed to stop the rotating disc instantaneously. This procedure would be more 184 advantageous than using the experimentally measured velocity field and turbulence intensity after the first 100 ms decay as initial conditions for subsequent calculations. An alternative and probably easier approach is the measurement of the disc deceleration rate wich can subsequently be incorporated in the numerical model. In this way the transient calculation can be initiated with the fully known steady-state calculated profiles. (ii) It is believed that turbulence field generated by swirling flows is non-isotropic. A two-components LDA system can be used to measure both tangential and radial components of turbulence intensity simultaneously. These measurements can provide physical evidence of such anticipated behaviour. Wilcox and Chambers model predicts the normal component of turbulence intensity which can be validated against the above mentioned measurements. (iii) It is also believed that turbulent shear stress is increased in the unstable region of swirling flow. The two-component LDA system can be used to examine the validity of such phenomenon. REFERENCES 185 Adrian, R J . , Fingerson, L.M., "Laser Anemometry Theory, Application and Techniques", TSI. Lecture Notes, 1981. Adrian, R.F. and Earley, W.L., "Evaluation of LDV Performance Using Mie Scattering Theory", Proc. of the Minnesota Symp. on Laser Anemometry, PP.426-454, 1975. Amsden, A.A., Ramshaw, J.D., O'Rourke, P J . and Dukowicz, J.K., "KIVA: A Computer Program for Two-Dimensional and Three-Dimensional Fluid Flows with Chemical Reaction and Fuel Sprays", Los-Alamos Report, 1984. Arcoumanis, C. and Whitelaw, J.H., "Fluid Mechanics of Internal Combustion Engines: A Review", Int. Symp. on Flows in Internal Combustion Engines, Vol 28, pp.1-17, 1985. Benodekar, R.W., Gosman, A.D. and Issa, R.I., "The TEACH-U Code for the Detailed Analysis of Two-Dimensional Turbulent Recirculatory Flow", Report No. FS/83/3, Imp. College of Sci. and Tech. Eng. Dept. London 1983. Boussinesq, J. Cited by Schlichting, 1877. Bradshaw, P. "Effect of Streamline Curvature on Turbulent Flow", AGARD-ograph, 1973. Brandstatter, W., Johns, R. and Wigley, G., "Calculation of the Flow Produced by a Tangential Inlet Port", S.A.E. Symp. on Flows in Internal combustion Engines, Nov.1985. Castro, LP. and Bradshaw, P., "The Turbulence Structure of a Highly Curved Mixing Layer", J.Fluid.Mech, Vol 73, pp.265-304, 1976. Daily, J.W. and Nece, R.E., "Chamber Dimension Effects on Induced Flow and Frictional Resistance of Enclosed Rotating Disc", J. Basic Eng. March 1960, pp 217-229. Dao, K., Uyehara, O.A. and Myers, P.S., "Heat Transfer Rates at Gass-Wall Interface in Motorid Piston Engine", S.A.E. Trans. Vol 82.3, pp.2237-2258, 1973. Davis, G.C., Mikulec, A., Kent, J.C. and Tabaczynski, R J . , "Modelling The Effect of Swirl on Turbulence Intensity and Burn Rate in S.L Engines and Comparison with Experiment", S.A.E. Technical Paper, 860325, 1986. Dimotakis, P.E., "Single Scattering Particle Laser Doppler Measurements of Turbulence", AGARD Conf. No. 193, 1976. Drain, L.E., "The Laser Doppler Technique", John Wiley and Sons, 1972. Dring, R.P. and Suo, M., "Particle Trajectories in Swirling Flows", J.Energy, Vol.2, No.4, pp.232-237, 1978. Dring, R.P., "Sizing Criteria for Laser Anemometry Particles", J.Fluids Eng., Vol.104, pp.15-17, 1982. Dyer, T.M., "Characterization of One and Two-Dimensional Homogeneous Combustion Phenomena in a Constant Volume Bomb", SAE Paper 790353, Vol.88, pp. 1196-1216, 1979. 186 Elkotb, M.M., Abou-Ellali, M.M.M. and Salem, I.S., "Rotating Flow and Non-isotropic Turbulence in Reciprocating Engines", A.S.M.E. Symp. on Flows in Internal Combustion Engines, pp.81-91, Nov. 1982. Feller, W.V., and Meyers, J.F., "Development of a Controllable Particle Generator for LV Seeding in Hypersonic Wind Tunnels", Minnesota Symposium on Laser Anemometry, pp.345-357, 1975. Fu, S., Launder, B.E. and Leschziner, M.A., "Modelling Strongly Swirling Recirculatory Jet Flow With Reynolds-Stress Transport Closures", Proc.6th Symp. on Turbulent Shear Flows, Toulouse, France, 1987. Galbraith, R.A., Sjolander, S.A. and Head, M.R., "Mixing Length in the Wall Region of Turbulent Boundary Layers", Aeronautical Quarterly, VoLXXVII, pp.97-110,1977. Gupta, A.K.,LiIley, D.G. and Syred, N., "Swirl Rows", Abacus Press, 1984. Hanjalic, K. and Launder, B.E., "A Reynolds Stress Model of Turbulence and Its Application to Thin Shear Flows", J. Fluid Mech. Vol 52, pp.609-638, 1972. Hanson, R J . and Thomass, A., "Flame Development in Swirling Flows in Closed Vessels", Combustion and Flames, Vol. 55, pp.255-277, 1984. Harlow, F.H. and Nakayama, P.I., "Transport of Turbulence Energy Decay Rate", Los Alamos. Sci. Lab. of University of California, pp.3-7, 1968. Hinze, J.O., "Turbulence", 2nd Edn., McGraw and Hill Book Company,1975. Inoue, T., Nakanishi, K., Noguchi, H., and Iguchi, S., "The Role of Swirl and Squish in Combustion of the SI Engine", VDI, No 370, pp.181-188, 1980. Issa, R.I., "Solution of Implicitly Discretized Fluid Flow Equations by Operator-Splitting", Dept. Mech. Engng., Imperial College, Fluid Section Rept. FS/82/15, 1982. Jones, W.P. and Launder, B.E., "The Prediction of Laminarization with a Two-Equation Model of Turbulence", IntJ.Heat and Mass Transf. Vol 15, pp.301-314, 1972. Jones, W.P. and Pascau, A., "Calculation of Confined Swirling Rows with a Second Moment Closure", J. of Fluids. Eng. Vol.111, pp.248-255, 1989. Kind, R J . , Yowakin, F.M. and Sjolander, S.A., "The Law of the Wall for Swirling Flow in Annular Ducts", Journal of Fluids Eng., Vol.111, pp. 160-164, 1989. Klebanoff, P.S., "Characteristics of Turbulence in a Boundary Layer with Zero Pressure Gradient", Report 1247, NACA, 1955. Kolomogorov, A.N., "Equations of Turbulent Motion of an Incomressible Ruid", 1942, Cited by Launder and Spalding 1972. Kondoh. T., Fukumoto, A., Ohsawa and Ohkubo, Y., "An Assessment of a Multi-Dimensional Numerical Method to Predict the Row in Internal Combustion Engines", S.A.E. Paper 850500, 1985. Koosinlin, M.L. Launder, B.E. and Sharma, B.L, "Prediction of Momentum, Heat and Mass Transfer in Swirling, Turbulent Boundary Layers", J. Heat Transf. Vol.96, pp.204-209, 1974. 187 Koosinlin, M.L. and Lockwood, F.C., "The Predication of Axi-Symmetric Turbulent Swirling Boundary Layers", AIAA, Vol.12, No.4, pp.547-554, 1974. Koyama, H., Masuda, S., Ariga, L and Watanabe, L, "Stabilizing and Destabilizing Effects of Coriolis Force on Two-Dimensional Laminar and Turbulent Boundary Layers", J. Eng. for Power, Vol.101, pp.23-29, 1979. Lai, K.Y., and Gosman, A.D., "Finite Difference and Other Approximations for the Transport and Navier-Stokes Solutions", Dept. Mech. Engng., Inperial College Rept. FS/82/16,1982. Lakshminarayana, B. "Turbulence Modeling for Complex Shear Rows", AIAA, Vol 24, pp.1900-1917, 1986. Launder, B.E., Reece, G J . and Rodi, W., "Progress in the Development of a Reynolds-Stress Turbulence Closure", J. Fluid Mech. Vol.68, pp.537-566, 1975. Launder, B.E., Priddin, C.H., Sharma, B.I., "The Calculation of Turbulent Boundary Layers on Spinning and Curved Surfaces", J. of Fluids Eng, Vol.99, No.l, pp.231-239, 1977. Launder, B.E. and Spalding, D.B., "Mathematical Models of Turbulence", Academic Press, London, 1972. Lawn, C J . and Elliott, C.J., "Fully Developed Turbulent Flow Through Concentric Annuli", J. Mech. Eng. Sci, Vol.14, No.3, pp. 195-204, 1972. Lilley, D.G. and Chigier, N.A., "Nonisotropic Turbulent Stress Distribution in Swirling Flows from Mean Value Distribution", IntJ. Heat and Mass Transf. Vol.14, pp.573-585, 1971. Lilley, D.G., "Prediction of Inert Turbulent Swirling mows", AIAA, Vol.11, No.7, pp.955-960, 1973. Mayle, R.E., Blair, M.F. and Kooper, F.C. "Turbulent Boundary layer Heat Transfer on Curved Surfaces", J. Heat. Transf. Vol.101, pp.521-525, 1979. McLaughlin, D.K. and Tiederman, W.K., "Biasing Correction for Individual Realization of Lazer Anemometer Measurements in Turbulent Flows", Phys. of Fluids, Vol.16, pp.2082-2088, 1973. Norris, L.H., and Reynolds, W.C., "Turbulent Channel Flow With a Moving Wavy Boundary", DepL Mech. Eng. Stanford University, Rept. FM-10, 1975. Patankar, S.V., "Numerical Heat Transfer and Fluid Flow", McGraw-Hill Publication, 1980. Pierik RJ . , "Swirling Combustion of Premixed Gaseous Reactants in a Short Cylindrical Chamber", M.A.Sc. Thesis, University of B.C. 1987. Prandtl, L. 1925, Cited by Launder and Spalding 1972. Prandtl, L. 1945, Cited by Launder and Spalding, 1972. Priddin, C.H., "The Behaviour of the Turbulent Boundary Layer on Curved, Porous Walls", Ph.D. Thesis, University of London, 1975. Pourahmadi, F. and Humphrey, A.C. "Prediction of Curved Channel Flow With an Extended K-Model of Turbulence", AIAA, Vol.21, No. 10, pp.1365-1373, 1983. Raithby, G.D., "A Critical Evaluation of Upstream Differencing Applied to Problems Involving Fluid Row", Comp. Meth. Appl. Mech. Engng., Vol. 9, pp. 75-103, 1976.a. 188 Raithby, G.D., "Skew Upwind Differencing Schemes for Problems Involving Fluid Row", Comp. Meth. Appl. Mech. Engng., vol. 9, pp. 153-164, 1976.b. Ramos, J.L, "Turbulent Non-reacting Swirling Flows", AIAA, Vol.22, No.6, pp.846-848, 1984. Ramos, J.L and Sirignano, W.A., "Axi-Symmetric Flow Model With and Without Swirl in a Piston-Cylinder Arrangement With Idealized Value Operation", S.A.E. Paper 800284,1980. Rayleigh, J.W.S., "On The Dynamics of Revolving Fluids", Proceedings of the Royal Society, V0I.6A, pp.148-154, 1916. Reynolds, O. "On the Dynamical Theory of Incompressible Viscous Fluids and the Determination of the Criterion", Phil. Trans. Roy. Soc. A, 186, pp.123-164, 1895. Rodi, W., "Examples of Turbulence Models for Incomressible Flows", AIAA, Vol.20, No. 7, pp.872-879, 1982. Rodi, W., "Influence of Buoyancy and Rotation on Equations for the Turbulent Length Scale", Second Int. Symp. on Turbulent Shear Flows, London, 1979. Saffman, P.G. "A Model for Inhomogeneous Turbulent Flow", Proc. Royal Soc. Lon. A, 317, pp.417-433, 1970. Schlichting, H., "Boundary Layer Theory", McGraw and Hill Company, 1979. Sharma, B.I., Launder, B.E. and Scott, C J . "Computation of Annular, Turbulent Flow with Rotating Core Tube", J.Fluids. Eng. Vol.98, No.4, pp.753-758, 1976. Sloan, G.D., Smith, P J . and Smoot, L.D., "Modeling of Swirl in Turbulent Flow Systems", Prog. Energy Combust. Sci, pp. 163-250, 1986. So, R.M.C., "Turbulence Velocity Scales for Swirling Flows", Turbulence in Internal Flows, Edited by Murthy, pp.347-369, 1976. Spalding, D.B., "A Novel Finite-Difference Formulation for Differential Expressions Involving Both First and Second Derivatives", Int. J. Num. Methods Eng., Vol. 4, pp. 551, 1972. Spalding, D.B., "A Single Formula for the Law of the Wall", Trans. ASME, Series E, pp.455-458, 1961. Taylor, G.L, "Fluid Between Rotating Cylinders", Proc. Roy. Soc, Series A, CLVIJ, pp.565-578, 1936. Tennekes, H. and Lumley, J.L. "A First Course in Turbulence", MIT. Press, 1972. Thomann, H., "Effect of Stream-Wise Wall Curvature on Heat Transfer in a Turbulent Boundary Layer", J.Fluid Mech., Vol.33, pp.283-292, 1968. Townsend, A.A., "Structure of Turbulent Shear Flow", Academic Press, New York, 1956. Vahl-Davis, G. and Mallinson, G.D., "False Diffusion in Numerical Fluid Mechanics", Dept. Mech. Engng., New South Wales Univ. Rept. 1972/FMT/l, 1972. Van Driest, E.R., "On Turbulent Flow Near a Wall", J.Aero. Sci, Vol.23, 1956. 189 Von Karman, T., "Some Aspects of the Turbulence Problem", Proc. of the 4th Int. Congress of App. Mech. pp.54-91, 1934. Vu, B.T. and Gouldin, F.C., "Flow Measurement in a Model Swirl Flow", AIAA Paper 80-0076, 1980. Wahiduzzaman, S., "A Study of Heat Transfer Due to a Decaying Swirling Flow in a Cylinder with Closed Ends", Ph.D. Thesis, Purdue University, 1985. Wakisaka, T., Shimamoto,Y. and Isshiki, Y., "Three-Dimensional Numerical Analysis of In-Cylinder Flows in Reciprocating Engines", S.A.E. Paper 860464,1986. Wang, CS. and Berry, G.F., "Heat Transfer in Internal Combustion Engines", A.S.M.E., Technical Paper, 85-WA/HT-23,1985. Wang, O.C. and Ferguson, C.R., "Multi-Dimensional Modeling of the Decay of Angular Momentum and Internal Energy in a Constant Volume Cylindrical Vessel", A.S.M.E. Symp. on Flows in Internal Combustion Engines, FED, Vol.28, pp.79-86, 1985. Wattendorf, F.L., "A Study of Effect of Curvature on Fully Developed Turbulent Flow", Proc. Roy. Soc. of London, A, Vol. CXLVIII, pp.565-598, 1935. Wilcox, D.C. and Chambers, T.L., "Streamline Curvature Effects on Turbulent Boundary Layers", AIAA, Vol.15, No.4, pp.574-580, 1977. Wilcox, D.C. and Traci, R.M., "A Complete Model of Turbulence", AIAA Fluid and Plasma Dynamics Conf. Paper 76-351, 1976. Yowakin, F.M., "Experimental Investigation of Turbulent Swirling Flow in an Annulus", Ph.D. Thesis, Carleton University, Canada, 1985. Yowakin, F.M. and Kind. R.J., "Mean Row and Turbulence Measurements of Annular Swirling Flows", Journal of Fluids Eng. Vol.110, pp.257-263, 1988. Zmeikov, V.N. and Ustimenko, B.P., "Hydrodynamics and Heat Transfer of Rotating Row Between Two Coaxial Cylinders", Proc. All Soviet Union Conf. on Heat Transf., Vol. 1, pp. 164-174, 1966. APPENDIX A. Derivation of Mean Transport equations: 190 In the following sections, the differential equations describing the conservation of mass and momentum pertaining to isothermal, incompressible, transient, axisymmetric turbulent swirling flows are derived. The conservation equations are first written in terms of the instantaneous dependent variables (u,v,w,p). Following Reynolds suggestion, the variables are decomposed into a mean and a superimposed fluctuating part. The equations are then ensemble-averaged to obtain the mean conservation equations. A . l . Conservation of mass The conservation of mass equation in terms of instantaneous velocities is ^u- + §Y_ + v dx dr r = 0 (A.l) now u = u + u' and v = v + v' so d_(-dx (u + u') + |.(v + v') + (^ )" = 0 (A.2) Taking the average of the above equation and noting that by definition u' = v' = 0 then the mean equation describing the conservation of mass would be du dv y dx + dr + r" = 0 (A.3) 191 A.2. Conservation of x momentum The differential equation describing the conservation of x component of momentum, in terms of instantaneous velocities and pressure is: 3u 3u 3u -=r- + V-r- + U=-3t dr dx 3x + r d t ( T X n ) + 3xxx (A.4) , [3u 3v where t r a = p + — and x x x = 2u du dx the left hand side of equation A.4 can be written as 3u 3u du ^ r - + VTT- + UTT-3t 3r 3x = P 3u 1 3 , x uv dv 3u 3 u 2 + f — (ruv) - up.. u— + u^ — - K r — 3t r 9 r v ' r dr dx 3x (A.5) using equation A. l the above equation can be simplified and substituted in the original momentum equation viz: 3u 1 3 / \ 3 u 2 ^ - + -L-(ruv) + 3— 3t r 3 r l ' dx dx + r d v ( T X T x ) + dx (A.6) writing all the variables in above equation in terms of their mean and fluctuating parts 3 ( -(u + u') +11; [r (u + u) (v + v')] + J- (u + u f 3 t v ~ ' ' r 3 r L M " ' ' ' / J ' 3x - ^ (p + P') + T I: Ir + 4)1 + Ir ta + T ' J 3x (A.7) , 3u 3v where T„ = u + 3r dx and x x x = 2u 3u 3x 192 taking the mean of equation A.7, noting that the mean of products of fluctuating components (uV or u'2) are not zero, then dx dx r 9r dx or =- H+F ^ H -^PuV)]+Ir te - p ^ < A - 9 > This is the mean x component of momentum equation. A 3 Conservation of radial momentum The differential equation describing the conservation of radial component of momentum in terms of instantaneous velocities and pressure is: 3v 3v dv w 2 dt dr 9x r r (A. 10) the steps followed to simplify the x component of momentum equation can be followed for the present case; i.e. the left hand side of the equation is simplified first, using the conservation of mass. The instantaneous variables are then replaced in terms of their mean and fluctuating parts. Finally the equation is averaged which gives the following equation ( A . l l ) 193 A.4 Conservation of tangential momentum The differential equation describing the conservation of the tangential component of momentum equation in terms of instantaneous velocities and pressure is : dt dr dx r (A. 12) The same procedure as the previous two momentum equations can be followed. The final form the mean tangential component of momentum equation is: ^ + i | [ r (vw)] + A(nw) + m ] = -L|[r2 fc-pvV)]+ | - f e - p u V ) (A.13) It was explained in Chapter 3 that the tangential component of momentum equation is replaced by the conservation equation for angular momentum. The latter equation is obtained direcdy from the former, by the following simple manipulations. The left hand side of equation A.13 can be written as: -J-p 4(wr) + 1- ^ -(rvwr) + ^-(uwr) dv ' r dr ' dxy ' (A. 14) The mean and the Reynolds stresses on the right hand side of equation A.13 are combined and are expressed in terms of strain-rates viz: 3w T i e - P V W = U e f f l - ^ - f and x x 9 - pu w = peff 3w 3x~ also 194 3w r 2 3r rUeff^(rw) - ^ ( r u e f f w ) also Substituting all of the above in equation A. 14, the following equation is obtained 3(wr) i 3(rvwr) 3(uwr) " 3 T + r dv 3x r 2 3r = T T Veff 3(wr) 3r r 3x Ueff 3(wr)' 3x ^|{ru e f fw) (A.15) multiplying the above equation by r, the angular momentum transport equation which is similar to the other transport equations is obtained viz: 3(wr) i 3(rvwr) 9(uwr) at r ^ ax r 3r 9 H ' 3r + 3x^ e f f 3(wr) ax 2-|{rueffw) (A. 16) 195 A P P E N D I X B . D E R I V A T I O N O F R E Y N O L D S STRESS TRANSPORT E Q U A T I O N S B .l. u' 2 transport equation: The axial component of instantaneous momentum equation is written in terms of mean and fluctuating velocities as | (u + u) + (v + v) J^u + u) + (u + u) | ^ u + u) d (- ,\ i d Icflu+u') ow+v'M d L d(u+u' (B.l) Multiply the above equation by u' and taking the mean, 3tV 2 / diu'2\ . _ ,3u' . , idu . 'du' .— du ,03u ^du 3x + V U -rr— + U V r - + U V hU U h U + U dr dr dx dx i dp u' 3 I du'i „'d{ dv\ ,dL du = - — u^f- + ^- — |rv—1 + Vr^ -lrv— + u^ -l2v p 3x r dr \ 9r / r dr\ dx / 3x\ 9x (B.2) D _ d d d rearranging the above equation and defining £ x ~9t + v a r ~ + u a ^ " m e n m e a D O v e equation can be simplified to D k 2 . ) = . u v §H. . u ' 2 ^ - A (ul) - 1 A lu&A JL U^P_ DtA 2 i dr 3x 3x I 2 / r dA 2 ) o dx 2„'2 |3_u^ dr du_ dx, (B.3) u vr2 3 2 1 ^ 5 2 where V = — + — — + 8r2 r 3 r 8x2 196 B.2. v'2 transport equation: The radial component of momentum equation in terms of mean and fluctuating velocity components is: ^(v + v ) + (v + v ) ^ v + v ) + (u + u ) ^ ( v + v ) - ^ W " r W ^ = - | ( P + P ) + J -3 3r| rp \ 3 r (B.4) 3x multiply the above equation by v and taking the mean, the following is obtained *h/ _ dv' ,r, dv ,03v' 8t + v v 3r + v ar + v ,3v' / <3v 8r r r dx 3 x 1 -dp = —v~ + v p 8r fav'j2 v'2 12/ i a r ) r 2 a v 3x A2 (B.5) D _ a a a rearranging the above equation and defining fx - 37 + V 3 r + U 3 x ^ t n t n e a D o v e equation can be simplified to D J v l L Dt\ 2 / / * dv ^ dv _ ' ' w o i n v 8x 8r r dx r 3r \ 2 ' w z v p 3r V 2 M -3v' law dv_ dx >\2 (B.6) B.3. w'2 transport equation The tangential component of momentum equation in terms of the mean and fluctuating velocity components is: 197 J^ (w + w') + (v + v') ^ ( w + w') + (u + u') + w ) " I-. . \ 9 u (w + w') (v + v') r dr 9(w+y/)\ d I d(w+w') dr dxf dx (B.7) w + w Mul t ip ly the above equation by w' and taking the mean, the fol lowing equation is obtained ,dw — ,3w - ,3w , ,3w W ' 2 V ' - w ' 2 ' ' w - -9w / /3w W - = r — + V W - T — + W V ^ r — + W V - T — + w v + V _ + W V + U W — + U W ^r— 3t dr dr dr r r r dx dx , /dw + u w — = v dx w , d 2 w ' ' dw' w ' 2 , d 2 w ' dr2 r ^r d x 2 (B.8) rearranging the above equation and defining D t " ^ ^ ^ m e n t n e above equation can be simplified to D k 2 Dt\ 2 = - u w "3 w dx" -w' 2^- - w'v r dr r . - v dw' "37 A(lLw'2 3x r 3r r w'2y' "7 w z . 3w' dx" w z v (B.9) B.4. Turbulence kinetic energy transport equation T h e turbulence kinetic energy k is defined as : k = l ( u ' 2 + v 2 + w'2) (B.10) T h e transport equation for the turbulence kinetic energy is obtained by adding the transport '2 '2 '2 equations for ^ » ^ and (i.e. equations B.3, B.6 and B.9). A d d i t i o n o f al l these three equations results in the following: 198 Dk = Dt Tdv - v - S r — + 2 wV — - u'v ^ — w'2 -^ - w'v dr r dx r 3w dr , -3w ^3u du - U W^r U - U V T -dx dx dr r dr r yl + r w|yl + R H Z V ! + j . D v 2 p 3_ 3x u l + v / V + v2uL4-lp'u' 2 2 2 p F _ pV ,oV du + v (viscous terms) The fifth group of terms on the right hand side of the above equation (V) cancels out due to conservation of mass. The third group of terms (TH) can be written as follows: r [v'k + jkp'v Also the fourth group of terms (IV) can be written as follows uk + — p u P The second group of terms (H) are the so-called production terms and are hereafter referred to by the symbol P. The viscous terms are expressed as: v (viscous terms) =v V 2 (k) - fcf . I^LF . ( l ^ L f . fexlF . l^!L .F - f o w ' l 2 . v l _ w'2 d r / I d r / \ d r / \ d x | \ d x j \ dx r2 r2 199 . w 2 , M 9*k i 9k 92k j 3 / 9k\ 9 \ „ u . . L . r u . where V (k) = + -7-5- + = -TKT r ^ H + . Combining this part of the viscous terms 9r2 r * 9x2 r 9r\ 9r/ 3 x 2 with the corresponding groups (IH) and (IV) and substituting the other terms in the original transport equation, the following equation is obtained: 200 APPENDIX C. Non-isotropic eddy viscosity model Kind et al. suggested that the mixing length can be defined by the following relation : (Cl) where x is the shear stress and xw is the corresponding value at the wall. It can be seen that the assumption of constant shear layer near a wall will reduce equation C l to the well known Prandd mixing length relation. Kind et al. also suggested that for a more general case of swirling annular flow, the length scale in the r-0 and r-x directions should be defined in analogy with equation C l as: (C.2) (C.3) where x^ and ^rx are the shear stresses in the r-0 and r-x direction, Te,w and xx,w are the corresponding wall values. The shear stresses are defined by the following relations: (C.4) 201 r a_ (C.5) Substitution of equations C.2 and C.3 into equations C.4 and C.5 results in the following relations for the wall shear stresses: (C.6) (C.7) where xw = V x >^w + w . From equations C.4 and C.5 it can be seen that M-rx and are d r ^ (C.8) M-rx = P & A/ ( f ^ + r|-(^) dr *r/ (C.9) Hence the two eddy viscosities are generally not equal, because according to equation C.2 and C.3 the two length scales are not equal. APPENDIX D. Similarity between Norris-Reynolds and van Driest length-scale modification According to Norris and Reynolds [1975], the dissipation rate e is defined as 202 1,1.5 k 0 5 l V (D.l) where 1 = CD min [K y , X, S] (D.2) the constant ce is equal to 13.2 and the constant CD is equal to 6.41. The turbulent viscosity v t also been defined by the same authors as: v t = k 0 5 1 CD c^ (D.3) Substituting equation D.3 and D.2 in equation D.l, the following expression for the dissipation rate is obtained: 1,1.5 e _ _ K — c DKy / + c e c D c ^ Yi v (D.4) Spalding and Launder [1974] definition of dissipation rate is based on the following relation 203 e = c ^ (D.5) 1.5 The dissipation length scale in the above equation is related to the mixing length by the following equation (refer to Spalding and Launder) l^o - 2 5 1, m where the mixing length is defined as: l m = Ky 1 - exp 26 (D.6) the bracketed term in the above definition of length scale is the modification proposed by van Driest which takes account of molecular viscosity effect in the viscous dominated part of the flow (near wall region). Substituting for length scale in equation D.5 from equation D.6, the following expression for dissipation rate is obtained: e = c 0.75 k 1 . 5 Ky 1 - exp 26 11 (D.7) Comparing equations D.7 and D.4, it can be seen that in the region where viscous effects are dominant (eg Y + = 4) the two relation are equal to within 8 per cent. This verifies that the length 204 scale modification proposed by Norris and Reynolds are of the same order of magnitude of that proposed by van Driest in the region of interest 205 APPENDIX E. Equivalency of k-e and e-Q model for isotropic, homogeneous turbulence The transport equation for turbulent kinetic energy for homogeneous, isotropic turbulent flow is: ! - « <*•» where by definition e = ^ y^ -;substitute for e in equation E. l , then, the following equation is obtained: (E.2) Now for isotropic turbulence, k = ; so after some simple arrangement, the above equation can be written as: The dissipation equation of e-Q model for isotropic, homogeneous turbulent flow is : 8 Q 2 = - p n 3 (E.4) 1 2 where by definition Q =-&-. Substitute for Q in equation E.4, then However, e = ^ v'2 , so equation E.5 can be written in terms of the normal component of intensity as: 1 For convenience the dissipation length scale in this definition has been given a different symbol namely 1^ 2 0 6 For isotropic turbulent flows, the two equations E.3 and E.6 should be identical, so it can be seen that the right hand side of these two equations should be equal, viz: 1 « = V ? T " 1 ( E ' 7 ) For homogeneous, isotropic turbulent flow, the transport equation for the mixing energy e is but by definition e = 2- v'2 so the above equation can be written as: ^ = -P*"v ' 2 (E.8) By comparing Equation E.8 with the exact transport equation for the normal component of intensity, it can be observed that: -P*Qv'2 = - l e or 0 = — ^ T Z T (E.9) 3 B* v'2 2 1 5 But by definition, Q = -f- and also e = ^ —, so substitute these back into equation E.9, and 1 1 replacing k and e by v'2, then the following equation is obtained: (E.10) 207 Comparing equation E. 10 and E.7 it can be deduced that P * = | (E.11) Wilcox and Traci suggest that P = 3-, So from the above equation it can be seen that p = 0.075 as apposed to the original value of 0.09 suggested by Wilcox and Traci. So equation E.l 1 which was obtained based on the assumption of isotropic and homogeneous turbulence is satisfied to within 16 percent 208 APPENDIX F. Wall Functions Close to a solid boundary, the molecular viscosity effects predominates. To take account of such effects, the turbulence model (e.g. the k-e model) need to be modified; one example of such modifications is that proposed by Jones and Launder [1972]. Low Reynolds number modelling requires large number of grid points in the vicinity of the wall. An alternative approach is to bridge the region where the molecular viscosity effects predominate by the use of the so-called wall functions. These functions are briefly explained in this section. F. l . Velocity functions Following Launder and Spalding [1974], close to the wall, the shear stress can be assumed to be constant and equal to the value prevailing at the wall; viz (Piam + Pt)^ u r = x s x w (F.l) Defining the dimensionless parameters u+=-S-u x (F.2) where y is the distance from the point of calculation to the wall and u^ is the friction velocity. Equation F.l can be written in terms of the dimensionless parameters defined by equation F.2 as: (F.3) 209 Two region of interest prevail. Close to the wall where the laminar viscosity is dominant (Pt « P l a n ) the above equation can be integrated to give: u + = Y + for Y + < 11.63 (F.4) Away from the wall, where the flow can be assumed to be fully turbulent (Pt »M-iam) equation F.3 can be written as (F-5> and following the suggestion by Townsend p t = P K Ux y substituting for l^ t i n equation F.5 a n d t n e n integrating, u + =^ln(EY + ) (F.6) rC where E is a universal constant whose value is equal to 9. F.2. Wall shear stress In the laminar sublayer where Pt « M - t m , the wall shear stress can be obtained directly from equation F.4 viz: u + = Y + U _ P y u x UT p. 210 but uT = <yj — so the wall shear stress can be obtained from the above equation to be *w = U « r (F.7) In the fully turbulent region, Launder and Spalding [1974] show that U X = C<P5VK" (F.8) Also from equation F.6 it can be seen that „ - K u x"ln(EY+) or n2 = _ U L K _ U _ ^ ln(EY+) or pu* K u In (EY+) Substituting for u x on the right hand side of the above equation from equation F.8 , the following equation for wall shear stress is obtained: Turbulence kinetic energy and dissipation rate The production terms in the kinetic energy equation are of the form stress times the strain-rate. Near the wall the stresses in these terms are replaced by equations F.7 or F.9 depending on the value of Y + . 211 Also near the wall, the equilibrium condition prevails. Under such a condition, the production of turbulence kinetic energy is balanced by dissipation rate, viz: " • i f f - i * < f i o> However Pt = P K U r y as suggested by Townsend. Also the velocity gradient in the above equation can be obtained by differentiating equation F.6; after substitution equation F.10 can be written as: p K U x y ( ^ f = p e or e = ^ (F.ll) The dissipation rate in the kinetic energy equation is replaced by the above expression near the wall. APPENDIX G. Error Analysis 212 The experimental data presented in the present study were mainly those using LDA technique. According to Wahiduzzaman [1985], there are a number of uncertainties associated with this technique. Velocity biasing is one source of error. It was explained in Chapter 5 that Dimotakis [1976] averaging process is believed to reduce this error. This method of averaging can only be used in the total burst mode of LDA operation. The comparison between the measurements using continuous and total burst mode of operations revealed that biasing error was minimal in the present study. The second source of error is due to velocity lag of seeding particles. In the present study, this was not of concern, since water which was used as working fluid is naturally seeded and hence particles follow the flow without any lag. The velocity is calculated from the following equation in the LDA analysis: W = | f d -f s |d f (G.l) where fa = doppler frequency fs = shift frequency fd = fringe spacing The error in the velocity measurements can be obtained by the method proposed by Kline and McClintoock [1953] which suggests that: (G.2) where °~w = error in w = independent variables = error in the independent variables From equation G.l , it can be seen that fd, fs and df are the independent variables. The error in the shift frequency was measured to be 0.2 per cent. The error in the doppler frequency as estimated 213 by Wahiduzzaman was of the order 0.25 per cent. The error in the specification of fringe spacing was estimated to be 2 per cent. For a typical velocity of 5 m/s, the doppler frequency is 1.56 Mhz, the shift frequency used to detect this doppler signal is 1 Mhz, and the fringe spacing in the present set up is 3.19 M-m. Using equation G.2 then: rjw = V|df (Ofc - (Jf,f + lo* (fd - Uf Substituting the appropriate numbers in the above equation, it can be seen that the error in the velocity is approximately 0.1 m/s or the error in the measured velocity is approximately 2 per cent This is the total error associated with the velocity measurements, using LDA technique. The other source of error is the random error associated with experimental procedure, such as motor rotational speed measurement, repeatability etc. The random error was calculated by repeating the same measurements up to 10 times and calculating the standard deviation of the mean velocity under steady-state condition. This was calculated to be of the order of 12 per cent. APPENDIX H . 214 The Effects of Aspect Ratio on the Mean Velocity It was pointed out in chapter 6 that the maximum speed attained by the bulk of the flow in a chamber was a function of the chamber aspect (L/D) ratio. To clarify this point, Fig. H.1 shows the experimentally measured dimensionless tangential velocity versus aspect ratio at dimensionless radius equal to 0.9 under steady state condition. It can be seen that the velocity decreases with an increase in the aspect ratio. It was also pointed out that the mean velocity decays faster in the shorter chamber. Fig. H.2. shows the decay of dimensionless tangential velocity at dimensionless radius equal to 0.9 for two aspect ratios, namely 0.5 and 0.1. It can be clearly observed that the mean velocity decays much faster in the shorter chamber. 215 1.0 0.8 H w/W 0.6 H • I ' 1 ' 1 ' 1 - 1 1 I 1 1 0 .0 0.1 0 . 2 0 .3 0.4 0 . 5 0 . 6 L/D Figure R l . Variation of dimensionless tangential velocity with chamber aspect ratio at dimensionless radius equal to 0.9. 216 0.8 t (s) Figure H.2. The decay of mean dimensionless tangential velocity at dimensionless radius equal to 0.9 in chambers of different aspect ratios, -e- L/D= 1/2 - * - L/D = 1/10
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Turbulent swirling flow in short cylindrical chambers
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Turbulent swirling flow in short cylindrical chambers Riahi, Ardeshir 1990
pdf
Page Metadata
Item Metadata
Title | Turbulent swirling flow in short cylindrical chambers |
Creator |
Riahi, Ardeshir |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | The effects of aspect ratio (L/D) on the rate of decay of swirl in a cylindrical chamber were experimentally studied using the Laser-Doppler-Anemometry technique. Preliminary measurements revealed that water should be used as working fluid; the results pertaining to air were inferred from Reynolds number similarity. The steady-state measurements revealed that a solid body type of rotation can be generated by a disc whose surface has been uniformly roughened. The effect of aspect ratio on the rate of decay of such flow field was studied in three chambers with aspect ratios in the range of interest to engine combustion. Experimental results showed a faster decay rate in the shorter chamber (i.e. smaller aspect ratio). This was attributed to the stronger swirl driven secondary flow pattern in the shorter chamber. A mathematical model describing axi-symmetric, decaying, turbulent swirling flow inside a short cylindrical chamber was also developed. The model was numerically solved, using the control-volume analysis, to provide insight on swirl decay in engines. The model validation was based on experimental observations. Turbulence parameters were represented by a two-equation turbulence model, modified for streamline curvature effects. The ad-hoc curvature modification of the standard k-e model proposed by Launder et al. and the mixing energy model developed by Saffman-Wilcox-Traci (SWT) were used to account for curvature effects. The analysis of steady flow between two long concentric cylinders, established the superiority of the latter over the former method. The SWT model was also successfully used in reproducing previously published experimental results, pertaining to decaying swirling flow field (mean velocity and turbulence intensity) in a short cylinder. The calculated turbulence intensity profile revealed that swirl promotes anisotropic turbulence. The validated numerical model was used to predict the effect of aspect ratio on the rate of decay of the flow field observed by the experimental measurements in the present study. The overall prediction of decay rate was successful, leading to the conclusion that Wilcox and Chambers model can be used in predicting the behaviour of two-dimensional transient turbulent swirling flows. |
Subject |
Turbulence |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-01-25 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0100446 |
URI | http://hdl.handle.net/2429/30810 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1990_A1 R52.pdf [ 9.07MB ]
- Metadata
- JSON: 831-1.0100446.json
- JSON-LD: 831-1.0100446-ld.json
- RDF/XML (Pretty): 831-1.0100446-rdf.xml
- RDF/JSON: 831-1.0100446-rdf.json
- Turtle: 831-1.0100446-turtle.txt
- N-Triples: 831-1.0100446-rdf-ntriples.txt
- Original Record: 831-1.0100446-source.json
- Full Text
- 831-1.0100446-fulltext.txt
- Citation
- 831-1.0100446.ris
Full Text
Cite
Citation Scheme:
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>
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-0100446/manifest