Numerical simulations of swirling flows

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Request accessible transcript

Transcript (if available)

Content INFORMATION TO USERS This manuscript has been reproduced from the microfilm master. UMI films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter face, while others may be from any type of computer printer. The quality of this reproduction is dependent upon the quality of the copy subm itted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely affect reproduction. In the unlikely event that the author did not send UMI a complete manuscript and there are missing pages, these will be noted. Also, if unauthorized copyright material had to be removed, a note will indicate the deletion. Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, beginning at the upper left-hand comer and continuing from left to right in equal sections with small overlaps. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6" x 9” black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact UMI directly to order. ProQuest Information and Learning 300 North Zeeb Road, Ann Arbor, Ml 48106-1346 USA 800-521-0600 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. NUMERICAL SIMULATIONS OF SWIRLING FLOWS by Peiji Chen A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Aerospace Engineering) May 2000 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. UMI Number: 3017990 _ ( ft UMI UMI Microform 3017990 Copyright 2001 by Bell & Howell Information and Learning Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. Bell & Howell Information and Learning Company 300 North Zeeb Road P.O. Box 1346 Ann Arbor, Ml 48106-1346 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. UNIVERSITY OF SOUTHERN CALIFORNIA THE GRADUATE SCHOOL UNIVERSITY PARK LOS ANGELES, CALIFORNIA 90007 This dissertation, written by Peij i Chen under the direction of h. ia„... Dissertation Committee, and approved by aU its members* has been presented to and accepted by The Graduate School, in partial fulfillment of re quirements for the degree of DOCTOR OF PHILOSOPHY Dean of Graduate Studies April 17. 2000 DISSERTATION COMMITTEE COMMITTEE S j U R eproduced with permission of the copyright owner. Further reproduction prohibited without permission Contents List Of Figures v Acknowledgm ents xxi A bstract xxii 1 Introduction 1 1.1 M o tiv a tio n s.................................................................................................... 1 1.2 Experimental O b se rv a tio n s....................................................................... 2 1.3 Theoretical Works ...................................................................................... 3 1.3.1 Quasi-Cvlindrical T heory............................................................... 3 1.3.2 Stability T h e o r y ............................................................................ 4 1.3.3 Wave T h e o r y ................................................................................... 5 1.3.4 Other Recent Theoretical A pproaches.......................................... S 1.4 Criteria for Vortex B re a k d o w n ................................................................ S 1.5 Previous Numerical Sim ulations................................................................ 9 1.6 Features of Present Numerical C o d e ........................................................ 10 1.7 Outline of the S t u d y ................................................................................... 11 2 M athem atical M odel 13 2.1 Governing E q u a tio n s ................................................................................... 13 2.2 Boundary C o n d itio n s.................................................................................... 14 3 Num erical M ethods 18 3.1 Spatial Differencing...................................................................................... IS 3.2 Temporal D isc re tiz a tio n ............................................................................. 19 3.3 Coordinate Transform in the Radial D ire c tio n ..................................... 19 4 Validation of Com puter Code 21 4.1 Check Computational Domain S i z e ......................................................... 21 4.2 Check Spatial R esolution............................................................................. 23 ii R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4.3 Check time s t e p ............................................................................................. 23 4.4 Check Solution C onvergence..................................................................... 25 4.5 Check the Flow Divergence ..................................................................... 29 4.6 Check Downstream Boundary C o n d itio n .............................................. 29 4.7 Check the Total Head on the Symmetric Line .................................... 29 4.8 Comparisons with Previous R esults......................................................... 31 5 A xisym m etric Solutions o f Re = 200 35 5.1 Basic Features of Solutions ..................................................................... 35 5.2 The Effects of S w irl...................................................................................... 63 5.3 Uniqueness of Solution Re = 200 ............................................................ 66 6 Som e Comparisons with Previous Studies 74 6.1 Comparisons with Experimental O b se rv a tio n s.................................... 74 6.1.1 Breakdown Response to Increase in Swirl L e v e l ...................... 74 6.1.2 Swirl A n g le ........................................................................................ 74 6.2 Comparisons with Theoretical A p p ro ach es........................................... 75 6.2.1 Quasi-Cylindrical A p p ro x im atio n ................................................ 75 6.2.2 Wave T h e o r y ..................................................................................... 76 6.3 Absolute/Convective Instability of Solutions Re = 200 79 7 A M echanism for A xisym m etric Vortex Breakdown 85 8 Verifications of Other Breakdown Criteria 93 9 Other A xisym m etric Solutions with High Re 95 9.1 A dynamical systems representation of the flow ................................ 95 9.2 Topologies of the steady flows for Re = 500 ........................................... 97 9.3 Topologies of the steady flows for Re = 1000 ........................................ 105 9.4 Topologies of temporally periodic flows for Re = 1000 ......................... 112 10 Parallel 3D Sim ulations o f Vortex Breakdown 127 10.1 Transition from bubble to spiral breakdown and its mechanism . . . 128 10.2 3D case of Re = 200. a = 1.0, V = 1 .2 .......................................................146 11 Sim ulations of Jet-like Profiles 150 11.1 Axisymmetric cases of Re = 200 .............................................................. 152 11.2 Axisymmetric Solutions of High R e .......................................................... 155 11.3 3D simulation of case Re = 200. V = 1 . 5 ................................................ 159 12 Conclusion 162 12.1 Future W o rk s ...................................................................................................164 iii R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 13 A p p e n d ice s 13.1 Appendix A: Spatial Finite Differencing ........................ 13.1.1 Treatment of the a x is .............................................. 13.2 Appendix B: Temporal Finite D ifferen ce........................ 13.3 Appendix C: Solving Poisson Equ. about $ using FFT 13.4 Appendix D: The Transport Equation for Passive Scalar R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. List Of Figures 2.1 Boundary c o n d itio n s ................................................................................ 2.2 Axial velocity at inlet. Dashed line: a = 1.0; solid line: a = 2.0; dash-dotted line: a = 0.5........................................................................... 2.3 Azimuthal velocity at inlet. Dashed line: V = 1.0; solid line: V = 0.S944; dash-dotted line: V = 0.63.......................................................... 2.4 Azimuthal vorticity at inlet. Dashed line: a = 1.0: solid line: a = 2.0; dash-dotted line: a = 0.5................................................................... 2.5 Axial vorticity at inlet. Dashed line: V = 1.0: solid line: V = 0.S944; dash-dotted line: V ’ = 0.63.......................................................... 3.1 A typical grid in r-z p l a n e ...................................................................... 4.1 Stream function contours of two domains for case Re = 200. a = 1.0, V = 1.0. Grid size for upper figure: (257, 51); for lower figure: (194,41) ....................................................................................................... 4.2 Stream function contours in the same region of two domains; case parameters Re = 200, a = 1.0, V = 1.0. upper figure: portion of domain [25,12.5] and grid [257,51]; lower figure: portion of domain [20,10] and grid [194,41]. 0-levels in the figures are not identical. 4.3 Axial velocity at r = 0.031 vs. z. Solid line: axial velocity in domain [20,10], dashed line: axial velocity in domain [25.12.5]. Case parameters Re = 200, a = 1.0, V = 1.0................................................. 4.4 Axial velocity at r = 0.031 vs. z in a small region. Solid line: axial velocity in domain [20,10], dashed line: axial velocity in domain [25,12.5]. Case parameters Re = 200, a = 1.0, V '’ = 1.0..................... 4.5 Stream function contours of two spatial resolutions for case Re = 200. a = 1.0. V = 1.0. Upper figure: grid size [257,51]: lower figure: grid size [193.41]. Domain sizes are the same....................................... 14 16 16 16 16 20 22 22 23 23 24 v R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4.6 Stream function contours in the same region of two spatial reso lutions; case parameters Re = 200, a = 1.0. V = 1.0. Upper figure: portion of grid size [257,51]; lower figure: portion of grid size [193,41]. Domain sizes are the same....................................................... 4.7 Axial velocity at r = 0.031 vs. z. Solid line: axial velocity for grid size 193*41, dashed line: axial velocity for grid size 257*51. Case parameters Re = 200, a = 1.0, V = 1.0.................................................. 4.S Axial velocity at r = 0.031 vs. z in a small region. Solid line: axial velocity for grid size 193*41, dashed line: axial velocity for grid size 257*51. Case parameters Re = 200. a = 1.0, V = 1.0....................... 4.9 Axial velocity at r = 0.031 vs. z. Symbol o: axial velocity for dt = 0.05, symbol +: axial velocity for dt = 0.02,0.01. Case parameters Re = 200, a = 1.0, V = 1.0. ' .................................................................. 4.10 Axial velocity at r = 0.031 vs. z in a small region of [0,20]. Symbol o: axial velocity for dt = 0.05, symbol +: axial velocity for dt = 0.02,0.01. Case parameters Re = 200, a = 1.0. V ’ = 1.0................... 4.11 The tim e history of axial velocity at point (1.77,0.139). Symbol o: for dt = 0.05, symbol +: for dt = 0.01. Case parameters Re = 200, Q = 1.0, V = 1.0.......................................................................................... 4.12 Axial velocity at r = 0.031 vs. z. Symbol o: axial velocity for error tolerance = 10-4 checked at every 5 nondimensional time units, symbol +: axial velocity for error tolerance = 10“ 1 0 checked at every 5 nondimensional time units. Case parameters Re = 200, Q = 1.0, V = 1.0.......................................................................................... 4.13 Axial velocity at r = 0.031 vs. z in a small region. Symbol o: axial velocity for error tolerance = 10-4 checked at every 5 nondi mensional time units, symbol +: axial velocity for error tolerance = 10-10 checked at every 5 nondimensional time units. Case pa rameters Re = 200, a = 1.0, V = 1.0. .............................................. 4.14 Axial velocity convergence history. Symbol +: error tolerance = 10-4 ; symbol o: error tolerance = 10~6; solid line : error tolerance = 10- u . All errors are checked every 5 nondimensional time units. Case parameters Re = 200, a = 1.0, V = 1.0....................................... 4.15 Time history of the flow divergence. Case parameters Re = 200. Q = 1.0, V = 1.0.......................................................................................... 4.16 Axial velocity at r = 0.031 vs. z. Symbol o: axial velocity for C = 0.6. symbol +: axial velocity for C = 0.5.0.7 and O.S. Case parameters Re = 200. ft = 1.0. V = 1.0................................................. 24 25 25 26 26 26 27 27 28 29 30 vi R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4.17 Axial velocity at r = 0.031 vs. z in a small region. Symbol o: axial velocity for C = 0.6. symbol +: axial velocity for C = 0.5.0.7 and O.S. Case parameters Re = 200, a = 1.0, V = 1.0.................................. 4.IS Variations of total head H = Pc/p 4- V ’2/2 on the symmetric line with Re for the fixed swirl rate V = 0.75. For comparison purpose the three lines are aligned into one point at entrance........................... 4.19 Axial velocity at r = 0.031 vs. z for three different domains. Dashed line: domain size [20,10], grid size [193,41]; Solid line: domain size [35,10], grid size [361,41]; dash-dotted line: domain size [40.10], grid size [385,41]. Case parameters Re = 500, a = 1.0. V = 0.8944. . . . 4.20 Axial velocity at r = 0.031 vs. z for two different spatial resolu tions. Dashed line: domain size [20,10], grid size [193.41]: Solid line: domain size [20,10], grid size [257,51]. Case parameters Re = 500, a = 1.0. V = 0.8944........................................................................................ 4.21 Stream function contour with super-imposed axial velocity at r — 0.031 vs. z. Domain size: [35,10]; grid size: [361.41]. Case parame ters Re = 500. a = 1.0. V = 0.8944............................................................ 4.22 Stream function contour with super-imposed axial velocity at r = 0.031 vs. z. Domain size: [50,10]; grid size: [501.41], Case parame ters Re = 1000, a = 1.0, V ’ = 0.8944......................................................... 4.23 Comparisons of axial velocities at r = 0.031 vs. z for various Re. Case parameters a = 1.0. V = 0.63............................................................ 4.24 Comparisons of axial velocity at r = 0.031 vs. z for various Re. Case parameters a = 1.0, V = 0.8944........................................................ 5.1 Axial velocity at r = 0.031 Vs. Z. Case parameters Re = 200, a = 1.0. V = 0.63. Domain and grid sizes are the same as in case V = 1.0. .......................................................................................................... 5.2 Stream function contour. Case parameters Re = 200. a = 1.0, V = 0.63. Domain and grid sizes are the same as in case V = 1 .0 . . 5.3 Axial velocity profiles. Re = 200, a = 1.0. V '’ = 0.63. Domain and grid sizes are the same as in case V = 1.0................................................ 5.4 Swirl velocity profiles. Re = 200, a = 1.0, V = 0.63. Domain and grid sizes are the same as in case V ' = 1.0................................................ 5.5 Axial velocity variation with z at fixed radial positions. Case Re = 200. q = 1.0. I ’ = 0.63. Domain and grid sizes are the same as in case I ’ = 1.0. (a) r=0: (b) r=0.463: (c) r=0.697 (d) r=l.O 30 31 33 33 34 34 34 34 37 38 38 38 39 vii R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5.6 Swirl velocity variation with z at fixed radial positions. Case Re = 200, q = 1.0, V ’ = 0.63. Domain and grid sizes are the same as in case V ' = 1.0. (b) r=0.463; (c) r=0.697 (d) r=1.0 39 5.7 Core pressure at r = 0.031. core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V = 1.0. Case Re = 200, a = 1.0, V = 0.63........................................................................ 39 5.S Azimuthal vorticity contour. Re = 200. a = 1.0. V = 0.63. Domain and grid sizes are the same as in case V ’ = 1.0............................................ 39 5.9 Axial vorticity contour. Re = 200, a — 1.0. V * = 0.63. Domain and grid sizes are the same as in case V ' = 1 .0 ................................................ 40 5.10 Radial vorticity contour. Re = 200, a = 1.0. V = 0.63. Domain and grid sizes are the same as in case V = 1.0............................................ 40 5.11 Z-Y side view of vortex line. Re = 200. a = 1.0. V ' = 0.63. Domain and grid sizes are the same as in case V ' =1.0................................................................... 40 5.12 X-Y side view of vortex line. Re = 200. a = 1.0. V = 0.63. Domain and grid sizes are the same as in case V ' =1.0....................................................................... 40 5.13 Vortex line (normalized view). Re = 200. a = 1.0, V ’ = 0.63. Domain and grid sizes are the same as in case V = 1.0................................................................... 41 5.14 Vortex line (real size view). Re = 200. a = 1.0. I' = 0.63. Domain and grid sizes are the same as in case V ' =1.0....................................................................... 41 5.15 Axial velocity at r = 0.031 Vs. Z. Case parameters Re = 200. a = 1.0. V ' = 0.S5. Domain and grid sizes are the same as in case V ’ = 1.0. Domain and grid sizes are the same as in case V = 1.0. . . 43 5.16 Stream function contour. Case parameters Re = 200. ci = 1.0. V ' = 0.S5. Domain and grid sizes are the same as in case V ’ = 1.0. . 44 5.17 Axial velocity profiles. Re = 200, a = 1.0. V * = 0.S5. Domain and grid sizes are the same as in case V = 1.0................................................ 44 5.IS Swirl velocity profiles. Re = 200. a = 1.0. V = 0.S5. Domain and grid sizes are the same as in case V ' = 1.0................................................ 44 5.19 Axial velocity variation with z at fixed radial positions. Case Re = 200, q = 1.0, V ' = 0.S5. Domain and grid sizes are the same as in case V ’ = 1.0. (a) r=0: (b) r=0.463: (c) r=0.697 (d) r= 1.0 45 5.20 Swirl velocity variation with z at fixed radial positions. Case Re = 200. o = 1.0. V = 0.S5. Domain and grid sizes are the same as in case V = 1.0. (b) r=0.463: (c) r=0.697 (d) r=l.O ................. 45 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5.21 Core pressure at r = 0.031. core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V — 1.0. Case Re = 200, a — 1.0, V = 0.85..................................................................... 5.22 Azimuthal vorticity contour. Re = 200, a = 1.0, V ' = 0.85. Domain and grid sizes are the same as in case V ’ = 1.0..................................... 5.23 Axial vorticity contour. Re = 200, a = 1.0. V = 0.S5. Domain and grid sizes are the same as in case V ' = 1 .0 ............................................. 5.24 Radial vorticity contour. Re = 200, a = 1.0, V = 0.85. Domain and grid sizes are the same as in case V = 1.0..................................... 5.25 Z-Y side view of vortex line. Re = 200. a = 1.0, V = 0.85. Domain and grid sizes are the same as in case V = 1.0................................................................ 5.26 X-Y side view of vortex line. Re = 200. a = 1.0. V = 0.85. Domain and grid sizes are the same as in case V' = 1.0................................................................ 5.27 Vortex line (normalized view). Re = 200. a = 1.0. V = 0.85. Domain and grid sizes are the same as in case V = 1.0................................................................ 5.28 Vortex line (real size view). Re = 200. a = 1.0. V = 0.85. Domain and grid sizes are the same as in case V = 1.0................................................................ 5.29 Axial velocity at r = 0.031 Vs. Z. Case parameters Re = 200. a = 1.0. V ' = 0.8944. Domain and grid sizes are the same as in case V ' = 1.0........................................................................................................... 5.30 Stream function contour. Case parameters Re = 200, a = 1.0. V = 0.8944. Domain and grid sizes are the same as in case V = 1.0 5.31 Axial velocity profiles. Re = 200. a = 1.0. V' = 0.8944. Domain and grid sizes are the same as in case V = 1.0..................................... 5.32 Swirl velocity profiles. Re = 200. a = 1.0, V ' = 0.8944. Domain and grid sizes are the same as in case V = 1.0..................................... 5.33 Axial velocity variation with z at fixed radial positions. Case Re = 200, q = 1.0, V — 0.8944. Domain and grid sizes are the same as in case V = 1.0. (a) r=0; (b) r=0.463; (c) r=0.697 (d) r=1.0 5.34 Swirl velocity variation with z at fixed radial positions. Case Re = 200, q = 1.0. V ’ = 0.8944. Domain and grid sizes are the same as in case V = 1.0. (b) r=0.463; (c) r=0.697 (d) r=1.0 ............. 5.35 Core pressure at r = 0.031. core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V * = 1.0. Case Re = 200. o = 1.0. V = 0.8944................................................................. 45 45 46 46 46 46 47 47 4S 49 49 49 50 50 50 ix R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. lO * 0 » 0 » Q l O » Q iO » Q » C .36 Azimuthal v r orticity contour. Re = 200, a = 1.0. V = 0.8944. Domain and grid sizes are the same as in case V ’ = 1.0............................. 50 .37 Axial vorticity contour. Re = 200, a = 1.0, V — 0.8944. Domain and grid sizes are the same as in case V = 1.0.............................................51 .38 Radial vorticity contour. Re = 200, a = 1.0, V = 0.8944. Domain and grid sizes are the same as in case V = 1.0............................................. 51 .39 Z-Y side view of vortex line. Re = '200, a — 1.0, V = 0.8944. Domain and grid sizes are the same as in case V — 1.0.................................................................... 51 .40 X-Y side view of vortex line. Re = 200, a = 1.0, V = 0.8944. Domain and grid sizes are the same as in case V = 1.0................................................................... 51 .41 Vortex line (normalized view). Re = ‘ 200, a = 1.0, V ’ = 0.8944. Domain and grid sizes are the same as in case V = 1.0............................................................ 52 .42 Vortex line (real size view). Re = ‘ 200, a = 1.0, V = 0.8944. Domain and grid sizes are the same as in case V = 1.0................................................................... 52 .43 Axial velocity at r = 0.031 Vs. Z. Case parameters Re = 200, a = 1.0, V = 1.0. Domain size: [20.10]; grid size: [193,41]........................53 .44 Stream function contour. Case parameters Re = 200, a = 1.0, V ’ = 1.0. Domain size: [20.10]; grid size: [193,41]................................... 54 5.45 Axial velocity profiles. Re = 200. a = 1.0, V = 1.0. Domain and grid sizes are the same as in case V = 1.0................................................ 54 5.46 Swirl velocity profiles. Re = 200, a = 1.0. V ' = 1.0. Domain and grid sizes are the same as in case V = 1.0................................................ 54 5.47 Axial velocity variation with z at fixed radial positions. Case Re = 200, a = 1.0, V = 1.0. Domain and grid sizes are the same as in case V = 1.0. (a) r=0; (b) r=0.463; (c) r=0.697 (d) r=1.0 55 5.48 Swirl velocity variation with z at fixed radial positions. Case Re = 200, a = 1.0, V = 1.0. Domain and grid sizes are the same as in case V = 1.0. (b) r=0.463; (c) r=0.697 (d) r=l.O ................. 55 5.49 Core pressure at r = 0.031, core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V ’ = 1.0. Case Re = 200, a — 1.0, V — 1.0.......................................................................... 55 5.50 Azimuthal vorticity contour. Re = 200. a — 1.0, V ’ = 1.0. Domain and grid sizes are the same as in case Y = 1.0.............................................55 5.51 Axial vorticity contour. Re = 200. a = 1.0. V = 1.0. Domain and grid sizes are the same as in case V = 1.0................................................ 56 x R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. i O > C 1 0 1 0 1 0 1 C 5.52 Radial vorticity contour. Re = 200, a = 1.0, V = 1.0. Domain and grid sizes are the same as in case V = 1.0................................................ 56 5.53 Z-Y side view of vortex line. Re = " 2 0 0 , a = 1.0, V ’ = 1.0. Domain and grid sizes are the same as in case V = 1.0................................................................... 56 5.54 X-Y side view of vortex line. Re = ‘ 200. a = 1.0, V — 1.0. Domain and grid sizes are the same as in case V = 1.0................................................................... 56 5.55 Vortex line (normalized view). Re = ‘ 200, a = 1.0, V = 1.0. Domain and grid sizes are the same as in case V = 1.0................................................................... 57 5.56 Vortex line (real size view). Re = ‘ 200. a = 1.0, V — 1.0. Domain and grid sizes are the same as in case V = 1.0........................................................................... 57 5.57 Upper figure: axisymmetric bubble followed by spiral breakdown (from Sarpkaya 1971). Lower figure: Comparison of stream function contours (from Grabowski & c Berger 1976). The lower are uniformly scaled................................................................................................................. 58 5.58 Axial velocity at r = 0.031 Vs. Z. Case parameters Re — 200, q = 1.0, V ’ = 1.095. Domain and grid sizes are the same as in case V = 1.0.............................................................................................................. 59 5.59 Stream function contour. Case parameters Re = 200, a = 1.0. V = 1.095. Domain and grid sizes are the same as in case V ' = 1.0. 60 5.60 Axial velocity profiles. Re = 200, a = 1.0. V = 1.095. Domain and grid sizes are the same as in case V = 1.0................................................. 60 .61 Swirl velocity profiles. Re = 200, a = 1.0, V ’ = 1.095. Domain and grid sizes are the same as in case V = 1.0.................................................. 60 .62 Axial velocity variation with z at fixed radial positions. Case Re = 200, a = 1.0, V * = 1.095. Domain and grid sizes are the same as in case V = 1.0. (a) r=0: (b) r=0.463; (c) r=0.697 (d) r=1.0 61 .63 Swirl velocity variation with z at fixed radial positions. Case Re = 200, a = 1.0. V = 1.095. Domain and grid sizes are the same as in case V = 1.0. (b) r=0.463; (c) r=0.697 (d) r=1.0....................... 61 .64 Core pressure at r = 0.031, core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V ' = 1.0. Case Re = 200, a = 1.0, V = 1.095...................................................................... 61 .65 Azimuthal vorticity contour. Re = 200. a = 1.0, V ' = 1.095. Do main and grid sizes are the same as in case I = 1.0 61 .66 Axial vorticity contour. Re = 200. a = 1.0. I = 1.095. Domain and grid sizes are the same as in case V ' = 1.0....................................... 62 xi R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. •5.67 Radial vorticity contour. Re = 200, o = 1.0. V = L .095. Domain and grid sizes are the same as in case V = 1.0........................................ 62 5.68 Z-Y side view of vortex line. Re = ‘ 200, a = 1.0, V = 1.095. Domain and grid sizes are the same as in case V ' =1.0................................................................... 62 5.69 X-Y side view of vortex line. Re = 200, a = 1.0, V = 1.095. Domain and grid sizes are the same as in case V — 1.0................................................................... 62 5.70 Vortex line (normalized view). Re = 200, a = 1.0. V = 1.095. Domain and grid sizes are the same as in case V = l.Q........................................................... 63 5.71 Vortex line (real size view). Re = 200, a = 1.0. V = 1.095. Domain and grid sizes are the same as in case V — 1.0................................................................... 63 5.72 Streamlines contour. Case parameters Re = 200, a = 1.0, V = 1.20. 64 5.73 Convergence history of maximum velocity change. Case parameters Re = 200, a = 1.0.* V = 1.20........................................................................ 64 5.74 Streamlines contour. Case parameters Re = 200, a = 1.0, V ' = 1.25. 64 5.75 Convergence history of maximum velocity change. Case parameters Re = 200, a = l.0 ,V = 1.25........................................................................ 65 5.76 The first minimum axial velocity on the axis vs. V .......................... 66 5.77 Comparison of results using different initial conditions. Case pa rameters Re = 200, a = 1.0, V = 0.8944.................................................. 67 5.78 Time history of axial velocity change at r = 0. Case parameters Re = 200. a = 1.0. V = 0.8944.................................................................... 68 5.79 Snapshots of axial velocity change at several different positions. Case parameters Re = 200, a = 1.0, V = 0.8944.................................... 69 5.80 Comparison of results using different initial conditions. Case pa rameters Re = 200, q = 1.0, V = 1.095.................................................... 70 5.SI Time history of axial velocity change at r = 0. Case parameters Re = 200, q = 1.0, V = 1.095...................................................................... 71 5.82 Snapshots of axial velocity change at several different positions. Case parameters Re = 200, a = 1.0. V = 1.095...................................... 72 5.83 Comparison of steady state results using different initial conditions. Case parameters Re = 3000, a = 1.0, V = 0.80...................................... 73 6.1 The breakdown parameters vs. R e ........................................................... 76 6.2 Streamfunction contours for t = 500.700,720 and 740. Case pa rameters Re = 1000. V = 0.S2. Domain size: [50.10]: grid size: [501.61] 77 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 6.3 Streamfunction contours for t = 750, 760,780 and 800. Case param eters Re = 1000, V = 0.82. Domain size: [50,10]; grid size: [501,61], .................................................................................................................................... 77 6.4 Streamfunction contours for t = 820,S40,S50 and 900. Case pa rameters Re = 1000. V = 0.82. Domain size: [50,10]: grid size: [501.61 ]............................................................................................................... 78 6.5 Streamfunction contours for t = 1000,1300,1600 and 2500. Case parameters Re = 1000, V = 0.S2. Domain size: [50,10]; grid size: [501.61 ]............................................................................................................... 78 6.6 Streamlines of solitary-wave solutions and the emergence of the re circulation region (From Leibovich & Kribus 1990), (a) X/^oo = 0.80, (b)A/^00 = 0.70, a small recirculation bubble appears.(c)A//i0 o = 0.001, a large recirculation bubble, (d) Detail of the bubble in (c). . S O 6.7 a-q param eter plane. Case parameters Re = 200. a = 1.0, V = 0.S5. 82 6.S a-q param eter plane. Case parameters Re = 200. ct = 1.0. V ' = 0.8944................................................................................................................. 83 6.9 a-q param eter plane. Case parameters Re = 200. a = 1.0. V ’ = 1.095 84 7.1 Time history of axial velocity change at r = 0. Case parameters Re = 200, a = 1.0. V = 0.85.’ ................................................................... 87 7.2 Snapshots of axial velocity change at several different positions. Case parameters Re = 200, a = 1.0, V = 0.85.............................................S8 7.3 Time history of axial velocity change at r = 0. Case parameters Re = 200, a = 1.0. V = 0.S944.................................................................... 89 7.4 Snapshots of axial velocity change at several different positions. Case parameters Re = 200, a = 1.0, V = 0.S944.................................... 90 7.5 Time history of axial velocity change at r = 0. Case parameters Re = 200. q = 1.0, V = 1.095...................................................................... 91 7.6 Snapshots of axial velocity change at several different positions. Case parameters Re = 200, a = 1.0, V ’ = 1.095...................................... 92 9.1 The history of the axial velocity at point (2.5,0.2297). Re = 500, a = 1.0, V = 0.83. Domain size: [50,10]; grid size: [501,61]................ 95 9.2 Streamfunction contour Re = 500, a = 1.0. V = 0.83. Domain size: [50.10]: grid size: [501.61].................................................................... 98 9.3 Streamfunction contour Re = 500. a = 1.0. V = 0.835. Domain size: [50.10]: grid size: [501.61].................................................................... 9S R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9.4 Streamfunction contour Re = 500. o = 1.0. V ' = 0.85. Domain size: [50,10]; grid size: [501,61] 9.5 Streamfunction contour Re size: [50,10]; grid size: [501,61 9.6 Streamfunction contour Rc size: [50,10]; grid size: [501,61 9.7 Streamfunction contour Re size: [50,10]; grid size: [501.61 = 500, q = 1.0, V ' = 0.S944. Domain = 500, a = 1.0, V ’ = 0.95. Domain 9.8 Streamfunction contour Re size: [50,10]; grid size: [501.61 9.9 Streamfunction contour Re size: [50,10]; grid size: [501.61 9.10 Streamfunction contour Re size: [50.10]; grid size: [501,61 9.11 Streamfunction contour Re size: [50,10]; grid size: [501.61 9.12 Streamfunction contour Re size: [50,10]; grid size: [501,61 9.13 Streamfunction contour Re size: [50,10]; grid size: [501,61 9.14 Schematics Re = 500. a = grid size: [501,61]................... 9.15 Schematics Re = 500. a = grid size: [501,61]................... 9.16 Schematics Re — 500, a = grid size: [501,61]................... 9.17 Schematics Re = 500, a = 1 grid size: [501,61]................... 9.18 Schematics Re = 500, a = grid size: [501.61]................... 9.19 Schematics Re = 500. o = grid size: [501.61]................ = 500, a = 1.0, V = 0.87. Domain = 500. o = 1.0. V ’ = 1.0. Domain = 500. a = 1.0. V ’ = 1.05. Domain = 500. a = 1.0. V — l.l. Domain = 500. ct = 1.0. V ’ = 1.15. Domain = 500. a = 1.0. V = 1.2. Domain = 500, a = 1.0. V ’ = 1.3. Domain .0, V ' = 0.835. Domain size: [50,10]; 1.0, V = 0.85. Domain size: [50,10]; 1.0, V ' = 0.87. Domain size: [50.10]; .0, V ’ = 0.8944. Domain size: [50.10]; 1.0. V = 0.95. Domain size: [50.10]; 1.0. I = 1.0. Domain size: [50.10]: 99 99 99 99 100 100 100 100 101 101 101 101 102 102 103 103 xiv R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9.20 Schematics Re = 500. a = 1.0. V ’ = 1.05. Domain size: [50.10]: 9.21 grid size: [501.61]. size: Schematics Re = 500, a grid size: [501,61]. 9.22 9.23 Schematics Re = 500. a = grid size: [501,61]................... Schematics Re — 500, a = grid size: [501.61]................ grid size: [501.61]. 9.24 Schematics Re = 500. a = 9.25 9.26 Streamfunction contour Re size: [50,10]; grid size: [501,61 Streamfunction contour Re size: [50,10]; grid size: [501.61 9.27 9.2S 9.29 9.30 9.31 9.32 9.33 9.34 9.35 Streamfunction contour Re size: [50,10]; grid size: [501.61 Streamfunction contour Re size: [50,10]; grid size: [501.61 Streamfunction contour Re size: [50,10]: grid size: [501,61 Streamfunction contour Re size: [50,10]; grid size: [501.61 Streamfunction contour Re size: [50,10]; grid size: [501,61 Streamfunction contour Re size: [50,10]; grid size: [501,61 Streamfunction contour Re size: [50,10]; grid size: [501,61 1.0, V ’ = 1.1. Domain size: [50.10]; 1.0, V ' = 1.15. Domain size: [50,10]; 1.0. V = 1.2. Domain size: [50.10]: 1.0, V = 1.3. Domain size: [50.10]; = 1000, q = 1.0. V ' = 0.S1. Domain = 1000. q = 1.0. V ’ = 0.S2. Domain = 1000. o = 1.0. V = 0.S35. Domain = 1000. q = 1.0, V ' = 0.85. Domain = 1000. q = 1.0. V ' = 0.875. Domain = 1000, q = 1.0. V = 0.95. Domain = 1000. a = 1.0. V ’ = 1.0. Domain = 1000. a = 1.0, V ' = 1.03. Domain = 1000, a = 1.0, V ' = 1.05. Domain 103 103 104 104 104 106 106 106 106 107 107 107 107 10S Power spectra of the axial velocity at (2.5,0.2297) for V as indicated. The spectra have been taken over 20 < t < 4000 using 21 7 samples. 108 Schematics Re = 1000. a = 1.0. 1 = 0.S1. Domain size: [50.10]: grid size: [501.61]............................................................................................... 109 xv R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9.36 Schematics Re = 1000. a — 1.0, V ' = 0.82. Domain size: [50.10]; grid size: [501,61]............................................................................................... 109 9.37 Schematics Re = 1000. a = 1.0, V = 0.835. Domain size: [50,10]: grid size: [501,61]............................................................................................... 110 9.3S Schematics Re = 1000, a = 1.0, I ' = 0.85. Domain size: [50,10]: grid size: [501,61]............................................................................................... 110 9.39 Schematics Re = 1000, a = 1.0, V' = 0.875. Domain size: [50,10]: grid size: [501,61]............................................................................................... 110 9.40 Schematics Re = 1000, o = 1.0, V ’ = 0.8944. Domain size: [50.10]; grid size: [501,61].................................................................................110 9.41 Schematics Re = 1000, a = 1.0, V ’ = 1.0. Domain size: [50.10]; grid size: [501,61]............................................................................................... I l l 9.42 Schematics Re = 1000. a = 1.0, V ’ = 1.03. Domain size: [50,10]; grid size: [501,61]............................................................................................... I l l 9.43 Schematics Re = 1000, a = 1.0, V = 1.05. Domain size: [50.10]: grid size: [501,61] I l l 9.44 The convergency history of the axial and radial velocity at point (z.r)=(2.5,0.2297). Case parameters Re = 1000, a = 1.0. V = 1.05. Domain size: [50.10]; grid size: [501,61]......................................................... 114 9.45 The convergency history of maximum velocity change and maxi mum local divergency checked every 5 nondimensional time. Case parameters Re = 1000, a = 1.0, V = 1.05. Domain size: [50,10]; grid size: [501.61]............................................................................................... 115 9.46 Stream function contour at t = 4000. Case parameters Re = 1000, q = 1.0, V = 1.05. Domain size: [50,10]; grid size: [501.61]..................... 115 9.47 The convergency history of the axial and radial velocity at point (z,r)=(2.5,0.2297). Case parameters Re = 1000. a = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]......................................................... 116 9.48 The period of the axial and radial velocity at point (z,r)=(2.5,0.2297). . Case parameters Re = 1000, a = 1.0, V ' = 1.11. Domain size: [50.10]; grid size: [501,61]................................................................................ 116 9.49 Stream function contour at t = 3880. Case parameters Re = 1000. a = 1.0. I ’ = 1.11. Domain size: [50.10]: grid size: [501.61].................... 117 9.50 Stream function contour at t = 3SS4. Case parameters Re = 1000. a = 1.0. V = 1.11. Domain size: [50.10]: grid size: [501.61].....................117 xvi R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9.51 Stream function contour at t = 388S. Case parameters Re = 1000, a = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61].................. 118 9.52 Stream function contour at t = 3892. Case parameters Re = 1000, a = 1.0. V ' = 1.11. Domain size: [50,10]; grid size: [501.61]................... 118 9.53 Stream function contour at t = 3896. Case parameters Re = 1000, a = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501.61]...................118 9.54 Stream function contour at t = 3900. Case parameters Re = 1000. a = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61].................... 118 9.55 Stream function contour at t = 3901.25. Case parameters Re = 1000, a = 1.0. V = 1.11. Domain size: [50.10]; grid size: [501,61]. . 119 9.56 Stream function contour at t = 3904. Case parameters Re = 1000, q = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]................... 119 9.57 Stream function contour at t = 390S. Case parameters Re = 1000, a = 1.0. V = 1.11. Domain size: [50,10]; grid size: [501.61]...................119 9.58 Stream function contour at t = 3912. Case parameters Re = 1000. ct = 1.0, V ’ = 1.11. Domain size: [50.10]: grid size: [501,61]...................119 9.59 Stream function contour at I = 3916. Case parameters Re — 1000. a = 1.0, V = 1.11. Domain size: [50,10]; grid size:...[501.61]................... 120 9.60 Stream function contour at t = 3920. Case parameters Re = 1000. a = 1.0, V = 1.11. Domain size: [50.10]; grid size: [501,61].................... 120 9.61 Stream function contour at t = 3922.5(same as t = 3SS0). Case parameters Re = 1000, a = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]...............................................................................................120 9.62 Schematic at t = 3880. Case parameters Re = 1000, a = 1.0, V ' = 1.11. Domain size: [50,10]; grid size: [501,61]............................. 121 9.63 schematic at t = 3884. Case parameters Re = 1000. a = 1.0. V ' = 1.11. Domain size: [50,10]; grid size: [501,61]............................. 121 9.64 schematic at t = 3888. Case parameters Re = 1000, a = 1.0. V = 1.11. Domain size: [50,10]; grid size: [501,61]....................................122 9.65 schematic at t = 3892. Case parameters Re = 1000, o = 1.0, V ’ = 1.11. Domain size: [50,10]; grid size: [501,61]....................................122 9.66 schematic at / = 3896. Case parameters Re = 1000. a = 1.0. 1 = 1.11. Domain size: [50.10]: grid size: [501.61]....................................122 x v ii R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9.67 schematic at t = 3900. Case parameters Re = 1000. a = 1.0. V = 1.11. Domain size: [50,10]; grid size: [501,61].....................................122 9.68 schematic at t = 3904. Case parameters Re = 1000, a = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]...................................123 9.69 schematic at t = 3908. Case parameters Re = 1000. a = 1.0. V = 1.11. Domain size: [50,10]; grid size: [501.61]...................................123 9.70 schematic at t = 3912. Case parameters Re = 1000. a = 1.0, V * = 1.11. Domain size: [50.10]; grid size: [501.61].....................................123 9.71 schematic at t = 3916. Case parameters Re = 1000, a = 1.0. V '’ = 1.11. Domain size: [50,10]; grid size: [501,61].....................................123 9.72 schematic at t = 3920. Case parameters Re = 1000. a = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61].....................................124 9.73 schematic at t = 3922.5. Case parameters Re = 1000. a = 1.0. V = 1.11. Domain size: [50,10]; grid size: [501.61]....................................124 9.74 The period of the axial and radial velocity at point (z.r)=(2.5,0.2297). . Case parameters Re = 1000, a = 1.0, V = 1.20. Domain size: [50.10]; grid size: [501.61]................................................................................. 124 9.75 The period of the axial and radial velocity at point (z.r)=(2.5,0.2297). . Case parameters Re = 1000, a = 1.0, V = 1.30. Domain size: [50.10]; grid size: [501.61].................................................................................124 9.76 The period of the axial and radial velocity at point (z,r)=(2.5,0.2297). . Case parameters Re — 1000, a = 1.0. V = 1.50. Domain size: [50.10]; grid size: [501.61]................................................................................. 125 9.77 The period of the axial and radial velocity at point (z,r)=(2.5,0.2297). . Case parameters Re = 1000. a = 1.0. V ’ = 2.0. Domain size: [50.10]; grid size: [501,61].................................................................................125 9.78 Power spectra of the axial velocity at (2.5,0.2297) for V as indicated. The spectra have been taken over 20 < t < 4000 using 21' samples. 126 10.1 Code Scalability................................................................................................ 127 10.2 Temporal fluctuation of the maximum velocity change. Case pa rameters Re = 200. c v = 1.0. V ’ = 1.095..................................................... 130 10.3 Temporal fluctuation of the three velocity component changes at r = I. Case parameters Re = 200. o = 1.0. I = 1.095......................... 131 xviii R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10.4 Time history of streaklines, T=5, 50 and 400. Case parameters Re = 200. a = 1.0, V = 1.095. The gaps in the lowest figure are due to the attractions of the rear stagnation point.....................................133 10.5 Time history of streaklines, T=500, 600 and 650. Case parameters Re = 200, q = 1.0, V = 1.095...................................................................... 134 10.6 Time history of streaklines. T=750, 850 and 1000. Case parameters Re = 200, a = 1.0, V = 1.095...................................................................... 135 10.7 Sequence of axial vel. changes (a), Case parameters Re = 200. Q = 1.0, V ' = 1.095.......................................................................................... 136 10.8 Sequence of axial vel. changes (b), Case parameters Re = 200, o = 1.0, V = 1.095.......................................................................................... 137 10.9 Sequence of axial vel. changes (c), Case parameters Re = 200, Q = 1.0. V = 1.095.......................................................................................... 13S lO.lOPower spectra density of the axial velocity at ~ = 5. It's sampled by 21 8 samples over 550 < T < 1000.......................................................... 139 10.11 Power spectra density of the axial velocity at : = 5. It s sampled by 21 8 samples over 1500 < T < 2000........................................................ 139 10.12Power spectra density of the axial velocity at z = 10. It's sampled by 21 8 samples over 550 < T < 1000............................................................ 139 10.13Power spectra density of the axial velocity at c = 10. It’s sampled by 21 8 samples over 1500 < T < 2000........................................................ 139 10.14Upper figure: Streaklines at T — 1000. lower part: Vortex lines at T = 1000. Case parameters Re = 200. a = 1.0, V = 1.095.................. 141 10.15Upper figure: pressure contour of plane 0 = 0° and ISO1 3 at T = 1000, lower part: streaklines at T = 1000. Case parameters Re = 200, a = 1.0, V = 1.095.......................................................................................... 142 10.16Upper figure: vorticity magnitude contour of plane 0 = 0° and 180° at T = 1000, lower part: streaklines at T = 1000. Case parameters Re = 200, a = 1.0, V = 1.095...................................................................... 143 10.17Passive scalar concentration at T = 5 0 .....................................................144 lO.ISPassive scalar concentration at T = 300 .................................................. 144 10.19Passive scalar concentration at T = 700 .................................................. 145 10.20Passive scalar concentration at T = 1500 ............................................... 145 10.2lTime history of streaklines, T=5, 100 and 350. Case parameters Re = 200, q = 1.0. V = 1.2.............................................................................147 10.22Time history of streaklines. T=400. 450 and 475. Case parameters Re = 200. a = 1.0. I ' = 1.2.............................................................................14S 10.23Time history of streaklines. T=500, 550 and 600. Case parameters Re = 200. o = 1.0. V = 1.2.............................................................................149 xix R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 11.1 The initial velocity profiles at inlet ...............................................151 11.2 The initial vorticity profiles at i n l e t ....................................................151 11.3 The streamfunction contour for case Re = 200, V * = 1 .2 .............152 11.4 The streamfunction contour for case Re = 200, V ’ = 1 .4 .............153 11.5 The time history of velocity convergence, case Re = 200, V = 1.4 . 153 11.6 The streamfunction contour for case Re = 200, V ’ = 1 .5 .........................153 11.7 The time history of velocity convergence, case Re = 200, V ’ = 1.5 . 154 11.8 The time history of velocity convergence. Case Re = 500. V ’ = 1.2 . 155 11.9 The Streamfunction contour. Case Re = 500. V ’ = 1 .2 ................... 156 ll.lO T he time history of velocity convergence, case Re = 500, V ’ = 1.5 . 156 ll.llT h e Streamfunction contour. Case Re = 500, V = 1 .5 ................... 157 11.12The time history of velocity convergence, case Re = 1000, V = 1 .5 . 157 11.13The Streamfunction contour. Case Re = 1000. V = 1 . 5 ................158 11.14Temporal fluctuation of the maximum velocity change. Case Re = 200, V = 1 .5 ...................................................................................................... 160 U.153D simulation. Time sequence of streaklines. T=30. 75, 200 and S00. Case Re = 200, V = 1 .5 .........................................................................161 13.1 A typical computational c e l l .........................................................................173 13.2 Accuracy checking for Poisson Eqn. about $ ......................................... 179 xx R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Acknowledgments I would like to take this opportunity to thank my advisor Prof. Eckart Meiburg for his encouraging support in my difficulties and the helpful advice. It’s joyful experience to study in this departm ent with so many nice and famous professors. Ell always carry the warm feelings of every Friday's group meeting in my future life. It’s in this wonderful group that I have learnt a lot and had a lot of funs. Somehow. I wonder I was spoiled too much in our group and may have difficulties to face the real outside world when I go out from here. Also, it’s the constant encouragement from my family, especially my wife Linda, that brings me at this stage. Thanks for her understanding about the difficulties of this nonlinear problem which took me so long to make some conclusions. Fm really not sure about the achievements of my research vet, but 1 do have a great feeling of our new born son, Allen Chen, who was born at almost S full pounds and gave us a big surprise. That makes me feel great about this life project. The author thanks Dr. Robert Verzicco for the helpful discussions on the numerical method and the permit to use part of the code from his research team. Also, discussions with Prof. Maxworthy and Prof. Redekopp are always quite helpful and encouraging. This is a joint project in collaboration with Prof. Maxworthy, and is funded by NSF. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Abstract The primary goal of this dissertation is to try and improve the fundamental understanding of the mechanism about vortex breakdowns. Vortex breakdown in spatially growing flows are investigated numerically by simulating the three dimensional incompressible, time-dependent Navier-Stokes equations in cylindrical coordinates. At relatively small values of Reynolds and swirl numbers, axisymmetric vortex breakdown solutions of bubble type are ob served and are stable under axisymmetrv assumption. This kind of breakdown is found to be strongly related to the upstream propagation of disturbance waves inside the vortex core, thereby lending support to recent attem pts to relate vortex breakdown to the emergence of absolute instabilities! Delbende et al. 199S. Rusak et al. 1998). Carrying the simulations to the 3D simulations, the axisymmetric bubble-type solution becomes unstable, and the transition to a temporally periodic, helical vortex breakdown flow is observed. The finite length region of absolute instability observed around the wake of bubble breakdown, together with the temporal peri odicity. indicate the possible existence of a global mode. The transition process of bubble type to spiral type breakdown observed in our simulations also agrees well with Rusak et al.'s theoretical arguments that axisymmetric breakdown represents an im portant step in the evolution of three-dimensional finite-amplitude distur bances and unsteady three-dimensional breakdown flows. At higher swirl number, a double-helix vortex breakdown is also observed. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 1 Introduction 1.1 M otivations Swirling flows with various degrees of swirl are of considerable physical interest and great practical significance. They feature in a variety of propulsion systems involving turbo-machinery, such as jet engines and turbo-pumps. In combustion applications, their ability to create reverse flow regions near the jet nozzle is ex ploited for the purpose of swirl- stabilizing the flame. A lot of chemical reactors and mixing devices are enhanced by making use of the increased spreading rates of swirling jets compared to non-swirling ones. High-Reynolds-number swirling flows are characterized by abrupt changes and regions of flow reversals and instabilities known as vortex breakdown phenomena, an important generic phenomenon for which a universally accepted explanation is still elusive. Therefore, the ability to understand the complicated structures of swirling flows and to predict the flow conditions that lead to these phenomena is essential for future utilization of swirling flows in the design of advanced aerody namic configurations, combustion chambers, nozzles and other flow devices where swirl has a dominant influence. The vortex breakdown phenomena have been widely studied and several review papers on this subject have been presented, including the reports by Hall (1972). Leibovich ( 197S. 1984) and Escuclier (19SS). Although there has been extensive research, the fundamental nature of these phenomena remains largely unexplained. I R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1.2 Experimental Observations Experimental efforts have focused on the studies of wing tip flows and vortex flows in various tubes (Harvey 1962; Sarpkaya 1971, 1974; Bruecker Althaus 1992. 1995). The two most common breakdown types observed in all experiments are spiral and nearly axisymmetric bubble breakdowns. The various types of break down can develop in flows with the same Reynolds number, with only a small change in the swirl ratio which is the one of azimuthal velocity to axial velocity upstream. All the experimental results show the following vortex breakdown behaviors. One or several stagnation points in the swirling flows, followed by regions of flow reversals and/or turbulence behind it; the axial flow decelerating along the vortex core centerline rapidly to a stagnation; and the deviation of swirl from solid-body like rotation near the stagnation point. Vortex breakdown also exhibits widespread hysteresis and bistability behav ior. The classic photograph of Lambourne & Bryer (1961) which simultaneously displays on each side of a delta wing, for the same imposed flow conditions, the spiral and the bubble breakdowns, constitutes a striking example of bistability. In confined tubes, Sarpkaya (1971) also identified flow regimes where the spiral and the bubble were highly unstable and the flow configuration changed over time from one to the other. Although the regime was not observed by Faler < V Lei- bovich (1977) with the same apparatus, the latter authors did reported random and unpredictable transitions between different kinds of bubbles. Vortex breakdown is also a temporally evolving phenomenon. Unsteady and asymmetric motions occur within the bubble region. Observed by Sarpkaya (1971), fluid is injected into and ejected out of the internal bubble region by an emptying and filling process taking place in the downstream part of the bubble. Furthermore, breakdown states display secondary temporal dynamics and wander back and forth along the axis in random fashion (Faler Leibovich 1977). The experiments about the role of the helical disturbances on the dynamics of breakdown is an open issue. As mentioned by Billant. C'homaz ic Huerre (1997). 0 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. two points of view have existed: either breakdown itself results from the devel opment of helical disturbances (Leibovich 197S) or breakdown is an independent phenomenon over which the helical disturbances may grow (Escudier 19SS). Lam- bourne (1965), Escudier Sc Zehnder (1982) and Althaus Sc Krause (1990) reported that the bubble is the basic form but it becomes unstable against helical dis turbances and develops into the spiral form, while Sarpkaya (1971) and Faler Sc Leibovich (1977) have observed that the bubble is a secondary development of the prim ary spiral form. In another approach, based on the observations of kinked solitary waves by Hopfinger, Browand & Gagne (1982) in their oscillating grid in a rotating fluid, Maxworthv, Hopfinger Sc Redekopp (1985) and Maxworthy. Mory Sc Hopfinger (1983) studied traveling kinked solitary wave motions and axisymmetric waves to determ ine their properties and their roles in vortex breakdown of subcritical core flow. They found that axisymmetric traveling waves cause breakdown when their am plitude exceeds a certain magnitude. 1.3 Theoretical Works The theoretical analyses of the vortex breakdown phenomenon have suggested three basic classes: quasi-cylindrical theory (analogy to boundary-layer separa tion), stabilities and wave theory. 1.3.1 Quasi-Cylindrical Theory The analogy to flow separation (Hall 1972) is based on the failure of the quasi- cylindrical approximations to the N-S equations to describe a swirling flow with large streamwise gradients. The topology of the flow near the breakdown point changes drastically and looks similar to the separation of a boundary layer over a rigid wall. Before the failure, assume the gradients in the axial direction are much smaller than those in the radial direction and Re is 0 (e -1 ). Then the X-S equations can be simplified as (Mager 1972) 3 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. (ur)r + (wr)~ = 0. 2/ 1 l’ l r = oPr. (uvr2)r + (l'ilt2)- = -i- [r2(ur - - ] r. tit r (uii7r)r + (w2r), = -^ rp .- + ( L//?e)[r(tt>)r]r. where (u, u, w) are the velocity components in cylindrical coordinates (r, 0. z). Sim ilar to the calculation of the separations of the boundary layers, the calculation of the vortex flow will be stopped when it meets a large axial gradient which in validates the initial assumptions. This stopped point is considered as a vortex breakdown position. Although this theory demonstrates the effects of swirling velocity, the rising of pronounced axial retardation, and the occurrence and posi tion of abrupt change, it still leaves too much unexplained; for example, it cannot provide any details near or after vortex breakdown. 1.3.2 Stability Theory The stability analyses study the tendency of small disturbances imposed on a vor tex flow to grow or decay in time and space. The analyses of Lessen, Singh k Paillet (1974) and Leibovich k Stewartson (1983) define several criteria for the stability of a vortex flow to generate axisymmetric and non-axisymmetric pertur bations. It is found that a columnar vortex with a large rotational core is stable to axisymmetric perturbations when the swirl ratio q (O.G39</ is the ratio of the maximum swirl velocity to the maximum axial velocity) is greater than 0.403, but it is unstable to non-axisymmetric helical perturbations when the swirl ratio is less than about 1.6. A review of vortex stability criteria is given in Leibovich (1984). Recent works done by Delbende. Chomaz c A Huerre (1996) and Lim .A Redekopp (1996) all demonstrate that the introduction of a moderate amount of 4 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. swirl promote absolute instability. Tsai & : Widnall (1980) indicate that incoming swirling flows exhibiting vortex breakdown are supercritical with respect to helical mode m = ±1. Furthermore, in the case of axisymmetric breakdown, the flow displays a transition from supercritical to subcritical with respect to helical modes m = +1. Finally, in the case of spiral breakdown, the supercritical-subcritical tran sition only occurs for the helical mode m = + 1, the axisymmetric mode m = 0 remaining supercritical throughout. Based on their works. Delbende. Chomaz k Huerre (1996) and Loiseleux. Chomaz k Huerre (1998) connect the now widely used concepts of absolute/convective instabilities to the supercriticality and sub- criticality by examining the linear impulse response of a given basic flow for large time. A property of all supercritical flow is that it cannot bear infinitesimal sta tionary waves according to Benjamin (1962), whereas such waves are possible for subcritical flows. If the response wave packet is convected away from the source location, the flow is said to be convectively unstable, i.e. supercritical. If it moves both upstream and downstream from the source location, the flow is said to be absolutely unstable, i.e. subcritical. However, as mentioned by Rusak (1997), the relation between vortex breakdown and vortex stability is still unclear. As was also pointed out by Leibovich (1984), breakdown can occur in a vortex flow with just slight signs of instability; a vortex flow can also become unstable without any breakdown phenomena. 1.3.3 Wave Theory Squire is the first who connected the vortex breakdown with wave theory. He thought the vortex breakdown phenomenon is the result of the propagation of small amplitude waves inside the vortex core. So he tried to find the stationary wave by considering inviscid, axisymmetric flow which is governed by the equation d2a - d2ib 1 du.' _ 2dH dC dx2 dr2 ^ r dr dv dv R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. where v is streamfunction, H is the total head defined as p/p + V 2/2 and 2t tC is the circulation. By assuming the perturbation streamfunction t/’ as j/’(x .r) = f(r) cos(ax), the above equ. becomes # S I df 2C dC , dr‘ rdr rH l'1 dr ’ where W is the axial velocity, a is the axial wave number. Squire considered the upstream velocity profiles in two basic flows: Rankine vortex and Burgers vortex. By solving the above characteristic equation numerically, he found that under certain conditions, there exists a stationary wave whose wave length is infinite, a = 0. Similar results were obtained for other upstream profiles. Based on these findings he concluded that there exists a critical state in which the downstream perturbations can propagate upstream that will bring up the vortex breakdown eventually. The criticism about Squire's theory is that it cannot describe the abrupt tran sition in the vortex core during the vortex breakdown. Moreover, according to his calculations, although the phase speeds will propagate upstream , the group veloc ity still moves downstream. So it's impossible to bring the perturbation energy upstream to cause the happening of vortex breakdown. Benjamin (1962) considered the ideas of the “conjugate-flow” theory, which proposed that vortex breakdown is fundamentally a structure transition from a uniform state of swirling flow to one featuring stationary waves of finite ampli tude. The original flow is assumed to be supercritical (i.e. incapable of bearing infinitesimal stationary waves), and the mechanism of the transition is explained on the basis of physical principles that are analogous to supercritical-flow phe nomenon of the hydraulic jum p or bore. The basic equations are the same as used by Squire. Benjamin proved that if there is a supercritical solution for a given set of H and C values, there must exist a conjugate solution which is subcritical. 6 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Benjamin's (1962) finite transition theory of vortex breakdown leaves open the im portant question of the structure of the transition region joining two conjugate flow states. This region is of very great interest, particularly since, in some cases at least, it remains fairly well ordered and involves backflows. It is likely that dissipation must be included to properly account for the transition region. There fore nonlinear wave theory is applied to vortex breakdown by Benjamin (1962, 1967) , Leibovich (1970, 1973, 1980) and Leibovich Sc Ma (1983). etc. Randall Sc Leibovich (1973) proposed that the bubble type vortex breakdown is the result of the propagation of an axisymmetric, finite am plitude perturbation wave inside the vortex core. To simulate the nearly axisymmetrical bubble type breakdown. Leibovich con sidered the long wave inside the axisymmetrical backflow region. The governing equations are D 2ii't + 2il'yD 2ii'z + -IT ; - 2yu'z{y~x D 2tl')y = 0 y r t - 2 v :r y + 2inyr z = o where ip = p (r.z,t), y = r2, D 2ip = 4yip.v + ti'yy + ipzz, T(i/) = rv. Perturbing the basic flow l'V'(y), T(y). after complicated calculations, the amplitude of perturbation wave .4 (r,t) satisfies the equation 4 t = c0.4. - f * eci.4.4- -I- k ‘c2.Azzz The above equation is the famous KdV equation in nonlinear wave theory. It has an analytical solution .4(r,f) = a se c /i(0 .5 (^ -)I/,2((r — c0t ) + ^actet)) c2 3 When (r — c0t) — > ±oc. .4 — > 0. So the perturbation is concentrated in a finite region, which is called a soliton. The existence of a soliton in the vortex flow strongly supports the wave theory for vortex breakdown which states that as time develops, the small perturbation R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. will develop into the finite am plitude wave which will eventually change the vortex core structure and bring up the vortex breakdown. Leibovich(19S4) summarized a theoretical scenario to explain the onset of vortex breakdown in a vortex flow in a tube. It is suggested that several interacting mechanisms are involved: the axisymmetric, bubble-type breakdown phenomenon is associated with the devel opment of large-amplitude axially symmetric waves in a basic supercritical vortex flow. This hypothesis is based on the above weakly nonlinear 'trapped wave’ the ory of Randall & : Leibovich (1973) and the recent solutions for the development in space and time of large-amplitude axially symmetric solitary waves on columnar vortices (Leibovich Kribus 1990; Ivribus & Leibovich 1994). 1.3.4 Other Recent Theoretical Approaches Recently, trying to give more details about breakdown point and the transition pro cedures, Rusak (1996) and Wang & : Rusak (1997) focus more on the region near vortex breakdown point and the dynamics of transition to axisymmetric break down. Wang & Rusak (1997) presented a theoretical approach which describes the ax isymmetric vortex breakdown process. This approach employs the unsteady Euler formulation, which is an inviscid model. Their results show that the axisymmetric vortex breakdown's evolution is a result of the interaction between disturbances propagating upstream and the inlet flow to the downstream. Later, our numerical results from the full N-S equations reveal in details about this view. 1.4 Criteria for Vortex Breakdown Several criteria have recently been proposed to predict the onset of breakdown in a vortex flow. Spall, Gatski & : Grosch (19S7) suggested a criterion for the breakdown of vortex flows in tubes in terms of the local Rossbv number (the inverse of swirl ratio) of the flow, that is defined at the radial distance of maximum swirl velocity. Based on the previous experimental, numerical and theoretical studies, they found that for any Reynolds number greater than 100. vortex breakdown occurs only 8 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. when the local Rossby number of the upstream flow is less than a critical value of about 0.65. Brown k Lopez (1990) analyzed the Euler equations that describe an incom pressible and inviscid axisymmetric swirling flow and proposed a theoretical crite rion for the axisymmetric breakdown that is based on the generation of negative azimuthal vorticity on some stream surfaces in the flow. It relates the tangents of the helix angles a 0 and 30 for the velocity and vorticity vectors at an upstream sta tion. A stagnation point will appear in the flow only when the condition a 0 > 30 is satisfied. This idea showed good agreement with numerical results of vortex breakdown in tubes. Marshall (1991) studied area-varying waves on curved vortex tubes by using a Polhausen type of integral momentum approach. The vortex breakdown is de scribed as a ’shock jum p’ in the vortex flow across which the core radius and the circulation of the vortex can be discontinuous. Theoretical criteria for both sta tionary axisymmetric and helical breakdown were presented in terms of local swirl ratio of the flow. Billant, Chomaz k Huerre (1997) derived a vortex breakdown criterion in the spirit of Brown & : Lopez(1990)’s breakdown criterion. Mahesh (1996) in his award winning paper proposed a very simple criterion for vortex breakdown of both incompressible and compressible flows. His breakdown criterion is that the breakdown is assumed to occur if the axial pressure rise exceeds the upstream streamwise momentum flux. It demonstrates the dependence of critical swirl number on upstream parameters. More experimental and numerical results are needed to test this simple criterion. 1.5 Previous Numerical Simulations A large number of numerical studies have been devoted to the simulation of vor tex breakdown phenomenon: the steady or unsteady axisymmetric incompressible Navier-Stokes equations have been numerically integrated by kopecky Torrence (1973). Grabowski k Berger (1976) and Krause. Shi .L ' Hartwich (19S3). etc. All 9 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. the results can describe flow fields that resemble vortex breakdowns which depend strongly on the initial disturbances, the upstream (or inlet) conditions and also the boundary conditions along the outer lateral boundaries of the computational domain. However, there are two limits for most of the previous studies. First, most of the simulations are steady cases; Secondly, most of the simulations are for pipe flows or closed containers which have more fixed boundaries than unbounded swirling jets and hence will meet less difficulties in numerical implementations. Recently, some 3D simulations are reported by Althaus k Weimer (199S), Lucas (1997). etc. Those results are time-dependent three-dimensional simulations that reproduced not only the unsteadiness and asymmetry of the internal flow within the bubble, but also the spiral breakdown; Lucas ( L997) also simulates the bubble- spiral breakdown transitions and revealed that the near-periodic oscillations in the spiral breakdown is the mechanism for the bubble-spiral breakdown transition. 1.6 Features of Present Numerical Code The same tip vortex model and inflow conditions as in Grabowski k Berger (1976) is used here. As mentioned by Grabowski. while the energy and angular momen tum considerations require that an unconfined vortex appears as a member of an oppositely rotating vortex pair, a single vortex may still be studied with the assumption that its opposite is far enough away to have little effect. The present numerical study differs from most of previous numerical investi gations in several respects. The first feature is that the solutions are of high res olution, 2nd-order accuracy everywhere, and time-dependent. Calculations using N-S equations by Grabowski k Berger (1976) gave numerical solutions depicting a closed recirculation zone which is very similar to the ’axisymmetric breakdown' bubbles found by Harvey (1960) and Sarpkaya (1971). Its numerical method is through solving the time-independent equations. It has not been shown whether the solutions are steady solutions of the time-dependent equations. Also, by the limitations of numerical method, the maximum Rt based on vortex core radius can be achieved is around 200. which is obviously too small, moreover, the numerical 10 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. accuracy decreases on the symmetrical axis where vortex breakdown phenomenon happens, because of the only first order approximation limited by its numerical method. The second feature is its capability of simulating 3D unbounded swirling flows. Recent numerical studies (Lopez 1990. etc. ) have solutions with higher Re and more accurate resolutions. However, as mentioned before, these kinds of studies are focused more on the vortex flows in pipes or in rotating closed cylindrical containers which are relatively easier to simulate than unbounded vortex flows. Based on the present code, bubble to spiral breakdown and the double helix breakdown are simulated successfully in the time dependent approach, which will reveal more dynamical behaviors of vortex breakdowns and shed some light on the explanations about these phenomena. The third feature is that its ability to handle the cases of large Re and high swirl rate where a lot of new and interesting phenomena happen, e.g.. multiple bubble breakdowns, periodic and chaotic breakdown solutions. All the previous axisymmetric simulations gave us at most two bubble breakdowns. The three bubble breakdown cases are reported in our simulation. Although Lopez (1990) reported the existence of periodic axisymmetric solutions in his simulation of in ternal flows in a closed cylindrical container, no similar reports have been found in other simulations. Therefore, the existence of periodic solutions in swirling jets, more specifically in the slender vortex flows, is a quite interesting phenomenon. The last one is that this numerical code is a parallel code with several features of flow visualization tool subroutines embedded, e.g.. the streaklines tracing program, the passive scalar solver. It can be scaled well to 64 CPUs and it's suitable for large simulations on parallel machines with more than 20 CPUs. 1.7 Outline of the Study The primary goals of this thesis are trying to improve the fundamental under standings of the vortex breakdown phenomenon and its mechanisms for the bubble breakdown, spiral breakdown and the transition between them. 1 1 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Some specific questions we’re trying to seek the answers are: (1) W hat’s the mechanism to cause the axisymmetric bubble breakdown? Any correlation with some recent vortex breakdown theories, like wave theory or sta bility analysis? (2) W hat is the transition mechanism of bubble to spiral breakdown observed also in our simulations? Chapter 2 describes the m athem atical model of the problem. Chapter 3 presents a numerical method (Verzicco & : Orlandi 1996) with no boundary conditions for pressure. Chapter 4 presents all the necessary testing results for validating the code, and also compares our time-dependent solutions with the same time-independent solutions of other numerical N-S simulations. Chapter 5 presents solutions for Re = 200 and discusses in detail all the features of variable changes in this vortex flow to demonstrate some im portant aspects of the axisymmetric bubble vortex breakdown. Also, by tracing the axial velocity changes inside the vortex core, the propagations of axisymmetric waves are observed inside the vortex core. As the swirl rate increases, the upstream moving part of the waves becomes more obvious. That kind of upstream wave will move against the incoming flow from nozzle and bring up the bubble type vortex breakdown eventually. Chapter 6 presents some comparisons of results with experimental observations. Chapter 7 compares the results with some theoretical works. Chapter 8 verifies some vortex breakdown criteria based on the present numerical results. Chapter 9 lists some periodic solutions, which have been observed only recently in a rotating closed cylindri cal container and discusses their behaviors from the view of dynamical systems. Chapter 10 focus on the full 3D simulations and discussions for the transition from bubble to spiral breakdown. 12 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 2 M athematical Model 2.1 Governing Equations The axisymmetric form of the Navier-Stokes equations in cylindrical coordinates (r, 0, r) can be rewritten with the quantities qg = v$,qr = rvr and q. = vz (Verzicco & : Orlandi 1996). The continuity equation is 1 dqr 1 dqe , d q s ; * + ; a e + a 7 = 0' {la) In terms of q, the momentum equations in conservative form become Dqe _ 1 < 9 p §SL\ _ S i i L d2 (l < > , & 3l , - Dt rdO Re r dr r d r ’ r 2 r2 d02 dz2 r3 8 0 J > Dqr _ dp 1 . 2d(l< > i /]L\ Dt r dr Re dr r dr r2 d92 d z2 r d o ’' ( } _ dp 1 A d dq: I & 2q- , d 2g ., ' D _ l _ > 1 ~ ' » 2 >1212 > 1 _ 2 J- with Dt dz Re r dr dr r2 dO2 d z2 Dqe _dqe_. 1 drqeqr cfy| dqgq.- Dt dt r2 dr r d9 dz Dqr dqr 8 q; ^ , d , q6qr ^ dqrqz ■ > = ~TT + — + ) + —j d n Dt dt dr r c)Q r dz 13 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Lengths have been non-dimensionalized by vortex core radius S. velocities by the free-stream axial velocity {qz)0c. time by S/(qz)^,. and pressure by p{q:),^o after subtraction of pco, the uniform static pressure far from the vortex (Grabowski Berger 1976). The core Reynolds number is defined as Re = 2.2 Boundary Conditions The dimensionless boundary and initial conditions will be prescribed to mimic some experimental conditions as closely as possible. In dimensionless coordinates, the computational domain extends from z = 0 to ; = L in the axial direction, and from r = 0 to r = R in the radial direction. Therefore, if we choose the core radius at r = 0 as J, conditions are specified on the boundaries of 0 < r < /?. 0 < z < L where R and L are much larger than one unit. ra R qr *0 Figure 2.1: Boundary conditions Notice the boundary condition specified on r = 0.0 < c < L is for axisymmet rical cases only. For full 3D simulations, no boundary condition is specified on this line. 14 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. For the test cases we specify the same boundary condition as in the paper by Grabowski k Berger (1976). At the upstream boundary z = 0, the velocities are specified as function of r (see Figure 2.1). qr{r) = 0; 0 < r < R. qe(r) = Vr( 2 — r2), q:{r) = a + (I — a )?’2(6 — 8r + 3r2); 0 < r < 1, qg{r) = V/r qz{f) = 1 ? 1 < r < R- At the axis, r = 0,0 < z < L, For symmetric cases: qr = 0. qg = 0, (2.1a) (2.16) At the radial boundary, r = Ft. 0 < z < L, qr = 0, qg = V/ fl, q: = 1 (2.1c) At the downstream boundary, z = L,Q < r < R, a simple radiation boundary condition (Orlanski 1976: Verzicco k Orlandi 1996) is applied to each variable < 7 , vection speed of large-scale structures. As mentioned by Salvetti (Salvetti 1996), automatically satisfying the global conservation of mass. It has been used success fully in some simulations of spatially developing free shear layers, coaxial jets, etc. Also, it is used well in large eddy simulations of turbulent confined co-annular jets and turbulent flow over a backward facing step (Akselvoll 1995). For the determi nation of the constant C, it was pointed out in (Salvetti 1996) that the value of C is not critical to the solutions. Later in our solution checking procedures we'll confirm this. As indicated by Grabowski Berger (1976). the upstream condition at c = 0 was chosen to approximate the experimentally-measured velocity in vortex cores (2.2 ) where C is assumed to be a constant velocity, and this velocity represents the ad- this assumption leads to simple boundary conditions that have the advantage of 15 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 2.2: Axial velocity at inlet. Dashed line: a = 1.0; solid line: a = 2.0; dash- dotted line: a = 0.5. Figure 2.3: Azimuthal ve locity at inlet. Dashed line: V ’ = 1.0: solid line: V ' = 0.S944; dash-dotted line: V = 0.63. Figure 2.4: Azimuthal vor ticity at inlet. Dashed line: ci = 1.0: solid line: a = 2.0: dash-dotted line: a = 0.5. Figure 2.5: Axial vorticity at inlet. Dashed line: V = 1.0: solid line: V ’ = 0.S944: dash- dotted line: V ' = 0.63. 16 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. such as those of trailing vortices. V is the specified circumferential velocity at the core edge, it is equal to the circulation around the core after non-dimensionalized by 27r(qz)006. The cubic form of qg simulates the solid body rotation near the core center, and a smooth transition to irrotational flow at the core edge. The parameter a is the ratio of the axial velocity at the core center to the velocity in the free-stream. So setting a greater or less than one yields jet-like or wake-like axial velocity profile. Uniform axial flow results when a = L.0( Figure 2.2 through Figure 2.5). In this model, swirl number is defined as Note that when a = 1.0, we'll have S'0 = V. 17 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 3 Numerical M ethods Choosing a goocl numerical method is extremely im portant in simulating this kind of highly nonlinear problem. From our experiences the vorticity streamfunction method is good for a lot of closed container. For pipe flows or unbounded flows, esp. the later ones, this method is quite complicated and. in general, tend to yield less accurate numerical solutions than the corresponding primitive variable method (Orzag. 1971). Also vorticity streamfunction method is more limited to axisymmetric cases and cannot handle generally 3D problems. So we choose the primitive variable method to solving the full N-S equations. Although we need to solve an equation with a pressure related term , by using FFT technologies in a recent scheme (Verzicco k. Orlandi 1996), we can avoid the pressure boundary conditions required by some traditional primitive methods and the performance of the code is also enhanced. 3.1 Spatial Differencing The spatial discretization of the system of Equs. (1) is performed by a centered finite differences, second order accurate scheme. The staggered grid system is employed here. The variables are located on a staggered grid with the pressure at the center of the cell and velocities < /, at cell surfaces. Using of staggered grids, together with the choice of the variables < /,. simplifies the discretization at the axis (r = 0). where the governing equations are singular mathematically. As mentioned IS R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. by Verzicco (Verzicco 1996), these singularities could be the main reason that few examples of DNS are found for cylindrical coordinates, compared with many examples of second-order finite difference schemes in cartesian coordinates in the literature. Indeed, the advantage of using staggered grids is that only the variable qr is evaluated at the point j= l( r = 0); since at that point qr = 0 by definition, it's not necessary to discretize the equation for qr at r=0. The spatial difference and the treatm ent of the singularity at the axis is presented in detail in Appendix A. 3.2 Temporal Discretization The time-advancement of Equ. ( 1) employs a fractional-step method extensively described in (Kim Moin 19S5; Rai & Moin 1991; Verzicco & : Orlandi 1996). Verzicco’s approach will be summarized briefly at here. At the first step an inter mediate nonsolenoidal velocity field q is provisionally calculated, then it’s projected into a solenoidal one by introducing a scalar quantity $ to which the pressure is related. The third-order Runge-Kutta scheme is used to advance the momentum equations in time, in which the nonlinear terms are computed explicitly and the linear terms implicitly. Using an approximate factorization technique (Beam 1976), the inversion of a large sparse m atrix is approximated by inversion of tridiagonal matrices. At each substep of Runge-Kutta scheme the scalar $ is calculated. The computation of $ implies the solution of an elliptic equation for which trigonomet ric expansions are applied in the axial direction. This approach is very efficient to calculate $ directly and gives the solenoidal velocity field within machine roundoff errors. Details are presented in appendix B. 3.3 Coordinate Transform in the Radial Direction Vortex breakdown is essentially located in the vortex core, so that good resolution is required in and near the core. A nonuniform grid in radial direction was used 19 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. to cluster the points in the vortex core region by adopting the following transfor mation suggested by Grabowski (Grabowski & Berger 1976): y(rj) = ^ n ( l + rj/b) (3.1) here rj is the radial physical coordinate of the point j; and yj is the computational coordinate corresponding to a uniform distribution of the iV 2 points in a [0,1/ 2] interval. The values of a and b are determined from the specified value of R and the desired number and distribution of grid points in the radial direction. A typical grid system in the following test cases can be specified by setting R = 10, and in region 0 < r < 1 the number of grid points to be 17 out of total 41, then the values a = 6.79970,6 = 0.34531. to 9 I,:,'.... !' :|-:|.I i ! : ........................ ............ :: i: i ;l!::h ' ' l.ii -y : .......... .!|: IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIM j i i ; i i m i n i i m i [ i i ] i 'i i i i i i i ! i i i [ ! ; i ! i i i ! ] ! [ i [ i i i ] n i i i : i i ; i i i i i i i i i ! i i i i i i i i i i i i i i i i i i i i i i i i i i i . i i i ! i i i i i i i : i ! i i i i i i i i i ; i i i i i i ; i ! i ! H i i i i i i i i i i i n i ! i i i |[ i i i i i i i i |[ i i i i i i i i i i [ i i i i i i i i i i i i i i l l l l l l i J t ! H I I I I I I ] [ 1 l i l l l ! l l l l l l i ! ] l ] ! l l l l l l | ] l l l [ l l I I I ] l l i : i l l l l l l l l l l l l l ! l ! l l l ! l l i : i ) l ! l l l l l l l l l l l i n H I I I I I ! ! l l l ! l l l l l l l l l l 1 l ! l l l l [ U ! l l l l l l [ ! U l l l l l l l l l l l ! l l l l l [ l l l l l l l l l l l H I I ! l l l ! l l i i iiim iiitiim iiiim im m iiiM iiiiiiiitiiiiiiH iiiiiiin ijiiiiitm iH iiiiiiiiiiiiim ijn iiiiiitiiiU iin iiiJ iiim iiiiiiu ijiiiitiim iiiiiiim iiiiiiiM iiiiiiiiiiiiiiiiiiiM iitiiiiim i ii mi m u iiiiim mu nun iiiiiiin iiiii iiiM tiiiiiu ii nun iin im iiiiiiiiiiitiiiiiiiiiH iiiiiiim iiiiiiiiM m n iiin iiiiiiiiiiiu iiiiiiiiiiiiiiiii .. 10 12 Z (193pts) Figure 3.1: A typical grid in r-z plane The grid must be uniform in the z direction because of the use of cosFFT in solving the pressure-related Poisson equation which will be discussed in detail later. This requires a relatively large number of grid points in z direction to get the resolution close to that in r direction (Figure 3.1). 20 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 01535302024848532323232348235323482323482353485348234823015301532348230200000000000200532353230201020002534848232301000002 Chapter 4 Validation of Computer Code In this section we present some results of the vortex breakdown with the same flow conditions in Grabowski & : Berger's paper. Comparisons with numerical and ex perimental results obtained by different methods are shown to prove the validation of the computer code. 4.1 Check Computational Domain Size Appropriate locations for the radial and downstream boundaries. R and L, were determined from a sequence of solutions at Re=200. These showed that R could be set at 10 and L at 20, since the use of the larger values changed the solutions insignificantly (see Figure 4.1 to Figure 4.4). This also agrees with Grabowski & Berger's (1976) parameters of domain size. However, for solutions of higher Re or swirling rate, these values may not be large enough. Choosing smaller L value will affect the region where the interesting phenomena are still going on and also will enforce the downstream boundary condition prematurely in that region. So for some higher Re or swirling rate, the values of L and R should always be checked to get the unchanged solutions. 21 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 4.1: Stream function contours of two domains for case Re = 200, a = 1.0, V = 1.0. Grid size for upper figure: (257, 51); for lower figure: (194,41) Figure 4.2: Stream function contours in the same region of two domains; case param eters Re = 200, q = 1.0, V ' = 1.0. upper figure: por tion of domain [25,12.5] and grid [257,51]; lower figure: portion of domain [20,10] and grid [194,41]. (/’ -levels in the figures are not identical. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. > z Figure 4.3: Axial velocity at r = 0.031 vs. z. Solid line: axial velocity in domain [20,10], dashed line: axial velocity in domain [25.12.5]. Case parameters Re = 200. q = 1.0, V = 1.0. Z Figure 4.4: Axial velocity at r = 0.031 vs. z in a small re gion. Solid line: axial veloc ity in domain [20.10], dashed line: axial velocity in do main [25,12.5]. Case param eters Re = 200. a = 1.0. V = 1.0. 4.2 Check Spatial Resolution The result of Re=200. a = 1.0, V = 1.0 was used to check the sensitivity of breakdown solution to spatial steps for fixed other parameters. Results from two different resolution grids shows that a 193*41 grid for computational domain [20.10] is a good one which has enough spatial resolution for vortex breakdown problems. From now on. all calculations for Re = 200 will stick to this grid system. 4.3 Check time step Similar to the check of spatial step, the case of Re = 200, a = 1.0, V ’ = 1.0 was used to check the sensitivity of solution to time step. Three different time steps t = 0.05, f = 0.02 and t = 0.01 were used. Results demonstrated that converged breakdown solution is not sensitive to the time step up to dt = 0.05. To demonstrate that this code is also accurate for unsteady solutions, the tran sient behaviors of this breakdown case was used again. Figure 4.11 shows the 23 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1 0 9 0 t 2 3 4 5 6 7 8 » 10 t l 12 t ] U IS 16 >7 16 I t Figure 4.5: Stream function contours of two spatial res olutions for case Re = 200, a = 1.0, V = 1.0. Upper fig ure: grid size [257,51]; lower figure: grid size [193.41]. Domain sizes are the same. Figure 4.6: Stream function contours in the same region of two spatial resolutions; case parameters Re = 200. a = 1.0, V ’ = 1.0. Up per figure: portion of grid size [257,51]; lower figure: portion of grid size [193,41]. Domain sizes are the same. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. z » z Figure 4.7: Axial velocity at r = 0.031 vs. z. Solid line: axial velocity for grid size 193*41. dashed line: axial velocity for grid size 257*51. Case parameters Re = 200, Q = 1.0. V = 1.0. Figure 4.8: Axial velocity at r = 0.031 vs. z in a small region. Solid line: axial ve locity for grid size 193*41, dashed line: axial velocity for grid size 257*51. Case parameters Re = 200, a = 1.0. V = 1.0. comparisons of the axial velocity time history at point (1.77.0.139) for 0.01 and dt = 0.05 respectively. Results also show that the time step can be set to dt = 0.05 for unsteady calculations. 4.4 Check Solution Convergence To make sure the vortex breakdown solutions are really the steady ones, the fol lowing procedures were used. The error tolerance, which is the maximum flow field velocity change checked at every 5 non-dimensional time units, was set to 10- 4. 10-6 and 10_lt} respectively, the three respective results were compared with each other (see Figure 4.12 to Figure 4.14). It shows that the solutions are really steady ones and the convergency criterion can be set to 10-4 . 25 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 4.9: Axial velocity at r = 0.031 vs. z. Symbol o: axial velocity for dt = 0.05. symbol +: axial velocity for dt = 0.02,0.01. Case pa rameters Re = 200. a = 1.0, V = 1.0. Figure 4.10: Axial velocity at r = 0.031 vs. z in a small region of [0,20]. Symbol o: axial velocity for dt = 0.05, symbol +: axial velocity for dt = 0.02,0.01. Case pa rameters Re = 200. q = 1.0, V = 1.0. Figure 4.11: The time history of axial velocity at point (1.77.0.139). Symbol o: for dt = 0.05. symbol + : for dt = 0.01. Case parameters Re = 200. a = 1.0. I = 1.0. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. N > 0.5 - 0 5 Z Figure 4.12: Axial velocity at r = 0.031 vs. z. Symbol o: axial velocity for error tolerance = 10-4 checked at every 5 nondimensional time units, symbol +: axial velocity for error tolerance = 10“ 1 0 checked at every 5 nondimen sional time units. Case parameters Re = 200, a = 1.0. V = 1.0. Figure 4.13: Axial velocity at r = 0.031 vs. z in a small region. Symbol o: axial velocity for error tolerance = 10-4 checked at every o nondimensional time units, symbol +: axial velocity for error tolerance = 10- l° checked at every 5 nondimensional time units. Case pa rameters Re = 200. o = 1.0. I = 1.0. 27 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. e rro rs le-4 error * 1 e-6 % > i ! -20 5.5 Z S 35 65 Figure 4.14: Axial velocity convergence history. Symbol +: error tolerance = 10-4 ; symbol o: error tolerance = 10-6; solid line : error tolerance = 10-11. All errors are checked every 5 nondimensional time units. Case parameters Re — 200, a = = 1.0, V = 1.0. 28 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4.5 Check the Flow Divergence Because all the solutions are obtained by simulating time-dependent, incompress ible flows, checking the flow divergence is an im portant issue of flow incompress ibility. Our simulations with FFT for the Poisson equation can make sure the flow divergence is around machine round-off error even at each intermediate step (see Figure 4.15). 3.5 Q 1.5 0.5 500 200 300 350 400 450 100 150 250 T Figure 4.15: Time history of the flow divergence. Case parameters Re = 200. a = 1.0. V ’ = 1.0. 4.6 Check Downstream Boundary Condition It is mentioned in the section of boundary condition that the constant C ' used in Equ. (2.2) is not critical to the solution. To confirm that the same simulation as above was carried out for four different C values of 0.5. 0.6, 0.7 and 0.8. Results are basically unchanged(see Figure 4.16 and 4.17). 4.7 Check the Total Head on the Symmetric Line Last we want to check the solution behaviors at high Reynolds numbers. Since most of theoretical works are based on the inviscid assumption, therefore, the total head H = Pc/p + V 1 j'l should be nearly constant on the symmetrical line for high 29 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. N > 0.S - 0 .5 z Figure 4.16: Axial velocity at r = 0.031 vs. z. Symbol o: axial velocity for C = 0.6. symbol +: axial velocity for C = 0.5.0.7 and O.S. Case parameters Re = 200, Q = 1.0. V = 1.0. z Figure 4.17: Axial velocity at r = 0.031 vs. z in a small region. Symbol o: axial velocity for C = 0.6, symbol +: axial velocity for C = 0.5,0.7 and O.S. Case parameters Re = 200. a — 1.0. V = 1.0. 30 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Re. Figure 4.18 shows the variations of H with three high Re on the axis. Clearly, the results demonstrate that H will be closer to be a constant value as Re becomes higher. V =0.75 0.8 0 .7 O.S 0.5 0 .4 0.3 0.2 Figure 4. IS: Variations of total head H = Pc/p + V ’2/2 on the symmetric line with Re for the fixed swirl rate V = 0.75. For comparison purpose the three lines are aligned into one point at entrance. 4.8 Comparisons with Previous Results Figure 5.1 to Figure 5.67 show the six solutions for Re=200, a = 1. Cases of V = 0.63, 0.85, 0.8944, 1.0 and 1.095 were calculated and compared with the results from Grabowski k Berger’s paper (1976). Results up to Y=0.S944 are almost identical even for the axial velocity profile versus z near the axis. For the cases with vortex breakdown at higher swirling rates like Y=l.O and 1.095. from the axial velocity profile, the breakdown position and size are nearly identical too. 31 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. except the small differences in the expansion region after the breakdown. These small differences will be amplified for the cases with higher Re and swirling rate. The reason for that difference could be that the first order approximations for the derivative boundary conditions at the axis was used by Grabowski Sc Berger (1976). Unfortunately the vortex breakdown happens exactly on the axis for the symmetrical cases. So it is quite obvious to demand a second order accurate central difference scheme which keep its accuracy on and near the axis. The method adopted here has this kind of advantage. In Grabowski Sc Berger’s paper the largest Reynolds number at which solutions could be obtained over a reasonably wide range of a and V was Re=200. Beyond that, converged solutions could not achieved based on central difference scheme. So several calculations at higher Re were performed using upwind differences to approximate the convective terms in the momentum equations, however, as indi cated by Grabowski Sc Berger, these solutions are only suggestive of the actual flows since the numerical viscosity introduced by upwind difference is usually hard to be measured. The scheme used here is based on the staggered grid system and has finite- volume like behaviors so that it can keep mass, momentum and energy conserva tions better than the regular central difference scheme which is usually based on unstaggered grid system (Verzicco Sc Orlandi 1996). Our numerical experiments also demonstrate the robustness of this scheme for high Reynolds number and swirl rate simulations. So two more converged solutions of Re=500 and Re=1000 with q = 1.0, V=0.S944 were obtained. Note that we still need to follow the above checking procedures to make sure we get the true solutions (see Figure 4.19 and Figure 4.20 for Re=500). It Is interesting to notice a second vortex breakdown region after the first breakdown (see Figure 4.21). For the solution of Re=1000, the third breakdown region is emerged out (Figure 4.22). This phenomenon was not observed in Grabowski Sc Berger’s high Re solutions. The explanations for this will be shown later. The further comparisons are the cases of axial velocity behaviors of several Reynolds number under a low and a moderate swirling rate. Figure 4.23 shows 32 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. the results at various Reynolds number for V=0.63 which agree quite well with Grabowski’s solutions even for Reynolds number high up to 1000. However, the near breakdown results for V=0.S944 are quite different from Grabowski’s at Re=500 and Re=1000 while almost identical for Re < 200 (Figure 4.24). » > Z > z Figure 4.19: Axial velocity at r = 0.031 vs. z for three different domains. Dashed line: domain size [20,10], grid size [193,41]; Solid line: domain size [35,10]. grid size [361.41]; dash-dotted line: domain size [40,10], grid size [385.41]. Case parameters Re = 500, a = 1.0, V ' = 0.S944. Figure 4.20: Axial velocity at r = 0.031 vs. z for two different spatial resolu tions. Dashed line: do main size [20,10], grid size [193,41]; Solid line: do main size [20,10], grid size [257,51]. Case parameters Re = 500, q = 1.0, I ’ = 0.8944. 33 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9 2 4 * I 10 12 14 19 I I 20 Z 0 2 4 I I 10 12 14 I I I I 20 Z Figure 4.21: Stream func tion contour with super imposed axial velocity at r = 0.031 vs. z. Do main size: [35,10]; grid size: [361,41]. Case parameters Re = 500, q = 1.0, V = 0.8944. Figure 4.22: Stream func tion contour with super imposed axial velocity at r = 0.031 vs. z. Do main size: [50,10]; grid size: [501,41]. Case parameters Re = 1000. ci = 1.0. V = 0.8944. > z > I Figure 4.23: Comparisons of axial velocities at r = 0.031 vs. z for various Re. Case parameters a = 1.0. V ’ = 0.63. Figure 4.24: Comparisons of axial velocity at r = 0.031 vs. z for various Re. Case parameters a = 1.0. V " = 0.S944. 34 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 5 Axisym metric Solutions of R e = 200 Solutions with various combinations of Re and V were obtained with the intent that each solution should be as informative as possible to dem onstrate some basic features of vortex breakdown and compare further with the previous numerical solutions. All the solutions with Re = 200 can provide a comparison with Grabowski L Bergers (1976) results and the existing theory of Mager (1972). The development of a vortex breakdown, in terms of a sequence of steady state solutions with in creasing amounts of swirl, will be described here. The solutions were all obtained with Re = 200. a = 1.0 so that the axial velocity distribution at z = 0 was uni form. The swirling parameter V was set at 0.63,0.80,0.85.0.8944.1.0 and 1.095. respectively. In most cases of Re = 200, little of interest occurs outside the region 0 < r < 2, 0 < z < 6, which is a small region of 0 < r < 20, 0 < z < 10, the computational domain for Re = 200. So most of the figures for Re = 200 will display various aspects of the solution in this small region only. 5.1 Basic Features of Solutions The solution for V = 0.63 is presented in Fig. 5.1 through Figure 5.10. Itrs a steady solution that will have similar converging behaviour as in case Re = 200. o = 1 .0 .1 = 1.0 Figure 4.14. 35 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.3 displays axial velocity profiles at a number of axial locations, and Figure 5.5 shows the variations of the axial velocity at several values of r with distance downstream. The slight retardation of the axial velocity within the vortex core reduces the axial velocity from 1 at r = 0 to about O.SS at r = 6. Figure 5.4 shows the swirling velocity profiles at the same axial locations of Figure 5.3, and Figure 5.6 displays the variations of swirl velocity at the same fixed radial positions as in Figure 5.5. The swirl velocity at r = 0 is always zero, so it is not labeled in the Figure. A reduction in the swirl velocities in the core due to the diffusion of vorticity away from the axis is quite evident in these figures. For example, the swirl velocity at r = 0.463 decreases from 0.55 at r = 0 to 0.434 at r = 6. Figure 5.7 presents the variations with z of three quantities: the core pressure near axis pc: the radius of the vortex core. rc, which has been defined by Grabowski & Berger (1976) as the radial distance where the circulation is decreased to 0.95V; and the maximum swirl velocity in the core. VmaI. Since the pressure in the core will be less than the ambient pressure, the nondimensional pressure {p — p ^ /p W ^ . will always be less than zero. In this case pc varies very slowly with z. increasing from about -0.63 at z = 0 to -0.50 at c = 6. Diffusion of vorticity away from the core center is responsible for both the increasing of the core radius from about 0.90 at c = 0 to 1.18 at ; = 6, and the reduction of the maximum swirl velocity from 0.695 to about 0.554. Figure 5.8 through Figure 5.10 present the three vorticity components of the flow field. Note even though no breakdown happens for this swirling rate, the azimuthal vorticity has a large negative region. Brown & Lopez (1990) argued that the production of a negative azimuthal component of vorticity is a necessary condition for vortex breakdown. If the negative azimuthal vorticity component near axis is strong enough, it will cause the retardation of the axial velocity on axis and eventually bring it to zero. In Grabowski Berger’s (1976) paper, an explanation for the decreasing of axial velocity within the core was made by approximating the radial momentum equation. Here is the summary of his approaches. 36 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. For slowly expanding cores in the limit of high Reynolds numbers, the following approximation can be used dp l’2 /- n JT = P— - 0-1) Or r The above equation expresses the balance of a particle’s centrifugal acceleration with the restraining pressure force. It is also a part of the '‘quasi-cylindrical" approximation of the N-S equations (Mager 1972). Integration with respect to r from the axis to radial infinity, followed by differentiation with respect to z, yields dpc f vd v . — - -dr Opc f 0 0 v Ov 1 7 = - p L 7 5 T where pc is the pressure at the center of the core. Since v is always positive, and viscous diffusion of swirl makes d v/d z < 0, there results a positive axial pressure gradient at the axis at which the retardation of the axial velocity happens. Figure 5.1 shows the axial velocity near the axis for the entire length of the com putational domain. Note that the deceleration of the axial flow does not continue forever. As v and —d v/d z decrease, the rate of change of the pressure will also decrease, so that the adverse pressure gradient becomes smaller with increasing z. Viscous forces are then able to overcome the pressure gradient and accelerate the fluid to the free stream axial velocity eventually. 2 - 1 . 5 - 1 ' M > 0 .5 - 0 - _05l -------- 1 -------1 -------1 -----1 ------ 1 --------1 -------1 -------1 ------ 1 ------ 1 0 2 4 6 8 10 12 14 16 18 2 0 Z Figure 5.1: Axial velocity at r = 0.031 Vs. Z. Case param eters Re = 200. n = 1.0. \ = 0.63. Domain and grid sizes are the same as in case I = 1.0. 37 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.2: Stream function contour. Case parameters Re = 200. a = 1.0. V = 0.63. Domain and grid sizes are the same as in case V = 1.0. c t - . . . . . . . . . OS > ! ; ; ; ; ; ; J \ ! j ! ! ! I ■ 01 1 1 _1. ...I ■ 1, ■ I ii,„. .. I ill L,1 i L l 1. 1 1 I I 1 1 .. 0 1 2 3 4 1 * Z Figure 5.3: Axial velocity profiles. Re = 200, a = 1.0, V = 0.63. Domain and grid sizes are the same as in case V = 1.0. ) ■ 1 1 1 -vr ■ ■ ■ r u m m °0 t 2 1 4 S 4 Figure 5.4: Swirl velocity profiles. Re = 200, a = 1.0, V = 0.63. Domain and grid sizes are the same as in case V = 1.0. R eproduced with perm ission of the copyright owner. Further reproduction prohibited without permission. Figure 5.5: Axial velocity variation with z at fixed ra dial positions. Case Re = 200, q = 1.0, V = 0.63. Do main and grid sizes are the same as in case V = 1.0. (a) r=0; (b) r=0.463; (c) r=0.697 (d) r= 1.0 Figure 5.6: Swirl velocity variation with z at fixed ra dial positions. Case Re = 200. q = 1.0, V = 0.63. Do main and grid sizes are the same as in case V ' = 1.0. (b) r=0.463: (c) r=0.697 (d) r = 1.0 Pc - — - rc ♦---- Vnu* c 2 Figure 5.7: Core pressure at r = 0.031, core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V = 1.0. Case Re — 200. a = 1.0. V = 0.63. Figure 5.8: Azimuthal vor ticity contour. Re = 200, a = 1.0, V = 0.63. Domain and grid sizes are the same as in case V = 1.0. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.9: Axial vorticity contour. Re = 200, a = 1.0, V ’ = 0.63. Domain and grid sizes are the same as in case V = 1.0. Figure 5.10: Radial vorticity contour. Re = 200, a = 1.0, V = 0.63. Domain and grid sizes are the same as in case V = 1.0. F ig u re 5.11: Z-Y side view of vortex line. Re = 200, a = 1.0, V = 0.63. Domain and grid sizes are the same as in case V ' = 1.0. I Figure 5.12: X-Y side view of vortex line. Re = 200, a = 1.0, V = 0.63. Domain and grid sizes are the same as in case V = 1.0. 40 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. > -0 2 -0 3 0 2 t Figure 5.13: Vortex line (nor malized view). Re = ‘ 200, a = 1.0, V = 0.63. D om ain and grid sizes are the sam e as in case V = 1.0 . F igu re 5.14: Vortex line (real size view). Re = ‘ 200, a = 1.0. V = 0.63. Domain and grid sizes are the same as in case V = 1.0. 41 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Increasing V from 0.63 to 0.85 results in the solution presented in Figure 5.15 through Figure 5.24. The stream surface expansion in Figure 5.16 now can be seen clearly. The radii of the surfaces near the axis reach a maximum around r = 5. and the axial velocity on the axis is quickly decelerated to a minimum of 0.23 also near z = 5. The swirl velocity in the core( Figure 5. IS and Figure 5.20 ) after decreasing initially in a manner very similar to the previous cases, reaches a minimum near r = 5, and then begins to increase. The acceleration of the swirl velocity near the center of the core continues for a short distance, and then gives way to final deceleration, which eventually will lead to total dissipation. Figure 5.21 shows that after a rapid initial increase, the pressure reaches a constant value and begins to decrease slightly from c = 5. The dynamic situation may be explained as follows. A fluid particle enters the domain at r = 0 which is close to the axis with axial velocity and angular momentum. As it proceeds downstream, the amount of retardation of the axial flow is larger than the previous cases since both v and d v/d z are initially larger. The continuity requires that a significant amount of fluid move radially outward. This movement is apparent in the stream contour plot. Since the angular momentum of the fluid particle is conserved as it moves radially except for the viscous effects, its angular velocity must decrease. From Figure 5.20 a slowly rotating region near the axis is apparent around z = 4 to c = 6. These effects bring a much larger increase in the pressure at the axis. For larger z, however, v and dv/d z are very small near the axis, the adverse pressure gradient also becomes small, and the shear forces are able to accelerate the axial flow. Then continuity requires an inward motion of fluid toward the axis, and, again assuming the angular momentum conservation for a particle, so the slight increase in the angular velocity at fixed radial positions will occur in Figure 5.20. The increasing angular velocity near the axis requires that the pressure on the axis decrease slightly, as can be seen in Figure 5.21 around r = 6. Beyond that point the effect of viscous diffusion once again becomes significant, reducing the swirl velocity near the axis and leading to an adverse pressure again which will decrease the rate of acceleration of the axial flow. This can be seen in Figure 5.15. from about : = 7 to ; = S. After that, as in previous cases. 42 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. the adverse gradient decreases with r as v and d v/d z decrease so that the axial velocity will eventually be accelerated by shear forces to the free stream velocity. N > 0 .5 - 0 .5 2 Figure 5.15: Axial velocity at r = 0.031 Vs. Z. Case param eters Re = 200, q = 1.0, V ’ = 0.85. Domain and grid sizes are the same as in case V = 1.0. Domain and grid sizes are the same as in case V = 1.0. The solution obtained with V = 0.8944 is displayed in Figure 5.29 through Figure 5.38. A small, well-developed bulge in the stream surface near _ = 3 is readily apparent in Figure 5.30. A similar stream surface expansion were obtained experimentally by Sarpkaya(1971). The axial velocity decrease (Figure 5.31 and Figure 5.33) in the core is consid erably more rapid than those in the previous cases, and nearly results in stagnation at 2 = 3 where the axial velocity reaches minimum of 0.01. The acceleration of axial flow outside the core as it passes over the bulge is also clearly seen in Figure 5.33. The swirl velocity profiles show, as in the previous cases, the decrease in the swirl velocity due to viscous diffusion and the reduction induced by the stream surface expansion. This reduction near the axis is particularly substantial. As Figure 5.34 demonstrates, the swirl velocity at r = 0.463 decreases from 0.76 at r = 0 to 0.22 at : = 3. As shown in Figure 5.35. the pressure increases steeply from -1.24 at : = 0 to -0.75 at z = 2. and then becomes essentially constant within the bulge region. 43 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2 ,---------------------------------- 1 ---------------------------------- 1 ---------------------------------- r I.S F < r 1 0.5 F pi i :■ l '------------------------i 0 1 2 3 4 5 6 Figure 5.16: Stream function contour. Case parameters Re = 200, a = 1.0. V = 0.85. Domain and grid sizes are the same as in case V = 1.0. t 2 a « i « Figure 5.17: Axial velocity profiles. Re = 200, a = 1.0, V = 0.85. Domain and grid sizes are the same as in case V = 1.0. Figure 5.18: Swirl velocity profiles. Re = 200, a = 1.0. V = 0.85. Domain and grid sizes are the same as in case I' = 1.0. 44 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5 z 5 2 Figure 5.19: Axial velocity variation with z at fixed ra dial positions. Case Re = 200, q = 1.0, V = 0.85. Do main and grid sizes are the same as in case V = 1.0. (a) r=0; (b) r=0.463: (c) r=0.697 (d) r=1.0 Figure 5.20: Swirl velocity variation with z at fixed ra dial positions. Case Re = 200. q = 1.0. V = 0.85. Do main and grid sizes are the same as in case V ’ = 1.0. (b) r=0.463: (c) r=0.697 (d) r=l.O i 0 : 2 c 2 Figure 5.21: Core pressure at r = 0.031, core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V = 1.0. Case Re = 200. o = 1.0. r = 0.85. Figure 5.22: Azimuthal vor ticity contour. Re = 200, q = 1.0, V ' = 0.85. Domain and grid sizes are the same as in case V ’ = 1.0. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.23: Axial vorticity contour. Re = 200, a = 1.0, V ’ = 0.85. Domain and grid sizes are the same as in case V = 1.0. Figure 5.24: Radial vorticity contour. Re = 200, a = 1.0, V = 0.85. Domain and grid sizes are the same as in case V = 1.0. V -O •0 ► Z F ig u re 5.25: Z-Y side view of vortex line. Re = *200, a = 1.0, V = 0.85. Domain and grid sizes are the same as in case V = 1.0. X Figure 5.26: X-Y side view of vortex line. Re = ‘ 200, a = 1.0, V = 0.85. Domain and grid sizes are the same as in case V = 1.0. 46 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R « -2 0 0 . V « O .S O 03 02 >• - 0.1 -0 4 04 0 2 -0 2 Figure 5.27: V ortex line (nor m alized view ). Re = 200, a = 1.0, V = 0.85. D om ain and grid sizes are the sam e as in case V = 1.0 . Figure 5 .2 8 : Vortex line (real size view ). Re = 200, a = 1.0, V = 0 .85. D om ain and grid sizes are the sam e as in case V — 1.0. 47 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2 < r < 4. with a value of about -0.64. The contraction of the stream surfaces then causes the pressure to decrease to -0.73 at : = 5.5 as the angular velocity of fluid near the axis increases also. After that, for larger c viscous effects become dominant and begin to decrease the swirl velocity and. once again, pressure near the axis is increased slightly. Figure 5.35 also shows the rapid increase of rc between r = L and z = 3. which is due to more rapid radial convection of core vorticity that also decreases Vmax quickly. Figure 5.29 presents the variation of the axial velocity on the axis. The initial decrease, as explained above, is due to the large adverse pressure gradient, and the first acceleration which begins at c = 3 is initially due to viscous forces and the slightly decreased pressure gradient. The second deceleration, beginning at about z — 5.5, is the result of the second adverse pressure gradient. Following is another acceleration of flow which is dominated by viscous forces again. 1.5 N > 0.5 - 0 .5 Z Figure 5.29: Axial velocity at r = 0.031 Vs. Z. Case param eters Re = 200, a = 1.0, V = 0.S944. Domain and grid sizes are the same as in case V * = 1.0. Figure 5.43 through Figure 5.52 display a vortex breakdown solution obtained with V ’ = 1.0. The small stream bulge of the previous solutions has now devel oped into a closed bubble of recirculating fluid, similar in shape to some of those discovered by Sarpkaya (1971a). The streamlines within the bubble have small '48 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2 .5 1 0.5 0 5 6 2 3 4 0 1 Figure 5.30: Stream function contour. Case parameters Re = 200, a = 1.0, V = 0.8944. Domain ancl grid sizes are the same as in case V = 1.0. Figure 5.31: Axial velocity profiles. Re = 200, a = 1.0, V = 0.8944. Domain and grid sizes are the same as in case V ' = 1.0. Figure 5.32: Swirl velocity profiles. Re = 200. a = 1.0, V = 0.8944. Domain and grid sizes are the same as in case I =1.0. 49 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.33: Axial velocity variation with z at fixed ra dial positions. Case Re = 200. q = 1.0, V = 0.8944. Domain and grid sizes are the same as in case V ’ = 1.0. (a) r=0: (b) r=0.463; (c) r=0.697 (d) r= 1.0 Figure 5.34: Swirl velocity variation with z at fixed ra dial positions. Case Re = 200, o = 1.0. V = 0.8944. Domain and grid sizes are the same as in case V = 1.0. (b) r=0.463: (c) r=0.697 (d) r = 1.0 o t 3i 1 01 0 t 2 S 3 4 6 Figure 5.35: Core pressure at r = 0.031, core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V = 1.0. Case Re = 200. a = 1.0. V = 0.8944. Figure 5.36: Azimuthal vor- ticity contour. Re = 200, a = 1.0, V = 0.8944. Do main and grid sizes are the same as in case V = 1.0. 50 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1. e 0. 2 3 9 « 0 t 4 Z Figure 5.37: Axial vorticity contour. Re = 200, a = 1.0. V ' = 0.8944. Domain and grid sizes are the same as in case V = 1.0. Figure 5.38: Radial vorticity contour. Re = 200, a = 1.0, V ’ = 0.8944. Domain and grid sizes are the same as in case V = 1.0. Z x Figure 5.39: Z-Y sid e view o f vortex line. Re = ‘ 200, a = 1.0, V = 0.8944. D om ain and grid sizes are the sam e as in case V = 1 .0. Figure 5.40: X -Y side view o f vortex line. Re = ‘ 200, a = 1.0, V = 0.8944. D om ain and grid sizes are th e sam e as in case V = 1.0 . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.41: Vortex line (nor m alized view ). Re = '200, a = 1.0, V = 0.8944. D om ain and grid sizes are the sam e as in case V = 1.0. Figure 5 .4 2 : Vortex line (real size view ). Re = '200, a = 1.0, V = 0.8944. D om ain and grid sizes are the sam e as in case V = 1.0. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. negative values and the outer streamlines, after expanding over the bubble, neck down behind it and expand slightly again. The axial velocity plots (Figure 5.45 and Figure 5.47 ) show the very rapid deceleration near the axis upstream of the bubble. Again like in Figure 5.43. the axial velocity decreases first, then increases and decreases again before its gradual return to free stream conditions (Figure 5.43). Figure 5.46 and Figure 5.48 show that due to the stream surface expansion the swirl velocities within and near the bubble are very small. For example, at r = 0.463, the swirl velocity decreases from 0.87 initially to 0.108 around c = 1.78. After this, it increases to 0.56 around z = A due to the streamlines contraction behind the bubble. Figure 5.49 shows a pretty steep increase in the pressure near the axis ahead of the bubble. The pressure near the axis within the bubble is nearly constant due to the very slow rotation inside the bubble. As in the previous cases, the streamline contraction results in a pressure decrease, beginning around z = 2, which assists the viscous forces to accelerate the axial velocity. Beyond z = 3.5. the pressure increases again with the diffusion of axial vorticity due to shear stresses. N > 0.5 - 0 .5 Z Figure 5.43: Axial velocity at r = 0.031 Vs. Z. Case param eters Re = 200, q = 1.0, V = 1.0. Domain size: [20,10]; grid size: [193,41]. The solution for a large, more vigorous breakdown, obtained with I = 1.095. is presented in Figure 5.5S through Figure 5.67. 53 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2 1.5 1 0.5 0 5 6 3 0 1 4 Figure 5.44: Stream function contour. Case parameters Re = 200. a = 1.0, V ' = 1.0. Domain size: [20.10]; grid size: [193,41]. Figure 5.45: Axial velocity profiles. Re = 200, a = 1.0, V = 1.0. Domain and grid sizes are the same as in case V = 1.0. Figure 5.46: Swirl velocity profiles. Re = 200. a = 1.0. V ' = 1.0. Domain and grid sizes are the same as in case T = 1.0. 54 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. > 2 Figure 5.47: Axial velocity variation with z at fixed ra dial positions. Case Re = 200. q = 1.0, V = 1.0. Do main and grid sizes are the same as in case V = 1.0. (a) r=0; (b) r=0.463; (c) r=0.697 (d) r= 1.0 z > I Figure 5.48: Swirl velocity variation with z at fixed ra dial positions. Case Re = 200, o = 1.0. V = 1.0. Do main and grid sizes are the same as in case V ’ = 1.0. (b) r=0.463; (c) r=0.697 (d) r= 1.0 2 X Z Figure 5.49: Core pressure at r = 0.031, core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V ’ = 1.0. Case Re = 200. a = 1.0. I' = 1.0. Figure 5.50: Azimuthal vor ticity contour. Re = 200, a = 1.0, V ' = 1.0. Domain and grid sizes are the same as in case V = 1.0. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.51: Axial vorticity contour. Re = 200, a = 1.0. V = 1.0. Domain and grid sizes are the same as in case V = 1.0. Figure 5.52: Radial vorticity contour. Re = 200. a = 1.0, V ' = 1.0. Domain and grid sizes are the same as in case V = 1.0. Z > X Figure 5.53: Z-Y side view of vortex line. Re = ‘ 200, a = 1.0. V = 1.0. Domain and grid sizes are the same as in case V = 1.0. Figure 5.54: X-Y side view of vortex line. Re = *200, a = 1.0, V = 1.0. Domain and grid sizes are the same as in case V = 1.0. 56 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5 .5 5 : Vortex line (nor m alized view ). Re = 200, a = 1.0, V = 1.0. D om ain and grid sizes are th e sam e as in case V = 1.0 . Figure 5 .5 6 : Vortex line (real size view ). Re = 200, a = 1.0, V = 1.0. Dom ain and grid sizes are the sam e as in case V = 1.0. 57 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.57: Upper figure: axisymmetric bubble followed by spiral breakdown (from Sarpkaya 1971). Lower fig ure: Comparison of stream function contours (from Grabowski < S j Berger 1976). The lower are uniformly scaled. 58 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The bubble is larger than the case V ' = 1.0 and has moved further upstream. The streamline expansion around the bubble is now very rapid, and the contraction behind the bubble is followed by a second pronounced expansion, beginning at about c = 3. Sarpkaya (1971a) reported that the fluid did not enter through the front of the experimentally obtained bubble, but instead passed over it and entered from the back. Then, after mixing turbulently inside, fluid exited into a second core-like region behind the bubble. While turbulence is not included in this study, the similarity between Figure 5.59 and his experimental photo (Figure 5.57 is obvious. Based on the geometric agreement. Grabowski & : Berger (1976) also mentioned that the location of spiral breakdown roughly coincides with calculated secondary stream surface divergence. N > 0 .5 - 0 .5 Z Figure 5.5S: Axial velocity at r = 0.031 Vs. Z. Case param eters Re = 200, a = 1.0, V = 1.095. Domain and grid sizes are the same as in case V = 1.0. It’s interesting to note that when we increase the swirl rate to value V = 1.20. the breakdown bubble starts to move away from the axis, and leave only one point attached to the axis( Figure 5.72 and Figure 5.73 ). W ith higher V as 1.25, the breakdown bubble now is fully away from the axis, there is a vortex ring structure. Notice in both cases, the expansion after the first bubble or ring is very strong and there is always a second bubble breakdown! Figure 5.7 1 and Figure 5.75 ). As mentioned by Grabowski Berger, the second expansion might be the axisymmetric counterpart to the spiral breakdown which is always observed behind 59 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2 .5 t 0.5 0 0 2 3 5 6 1 4 Figure 5.59: Stream function contour. Case parameters Re = 200. a = 1.0. V = 1.095. Domain and grid sizes are the same as in case V = 1.0. Figure 5.60: Axial velocity profiles. Re = 200, a = 1.0, V = 1.095. Domain and grid sizes are the same as in case V = 1.0. Figure 5.61: Swirl velocity profiles. Re = 200, a = 1.0, V = 1.095. Domain and grid sizes are the same as in case V = 1.0. 60 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.62: Axial velocity variation with z at fixed ra dial positions. Case Re = 200, a = 1.0, V = 1.095. Domain and grid sizes are the same as in case V = 1.0. (a) r=0: (b) r=0.463; (c) r=0.697 (d) r=1.0 i 2 Figure 5.63: Swirl velocity variation with z at fixed ra dial positions. Case Re = 200, o = 1.0, V = 1.095. Domain and grid sizes are the same as in case V * = 1.0. (b) r=0.463; (c) r=0.697 (d) r= 1.0 > 2 0 Z ' p O T O * ts / At!----- 00704 ------ ----------- rj 1 / PJW " _ - i - - - - - - I \ ~ Z - — 1 . ___________________ as R W K ' f A i> - - - ........................ ---------------------------------------- 1 2 3 Z * 3 3 Figure 5.64: Core pressure at r = 0.031. core radius and maximum swirl velocity vs. z. Domain and grid sizes are the same as in case V = 1.0. Case Re = 200. a = 1.0. V = 1.095. Figure 5.65: Azimuthal vor ticity contour. Re = 200, q =’ 1.0. V = 1.095. Do main and grid sizes are the same as in case V ' = 1.0. 61 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.66: Axial vorticity contour. Re = 200. a = 1.0, V ' = 1.095. Domain and grid sizes are the same as in case V = 1.0. > 2 • ( I 1 '2 '* ’• 'I t Figure 5 .6 8 : Z-Y side view o f vortex line. Re = 200, a = 1.0, V = 1.095. D om ain and grid sizes are the sam e as in case V = 1.0 . Figure 5.67: Radial vorticity contour. Re = 200, a = 1.0. V ' = 1.095. Domain and grid sizes are the same as in case V ’ = 1.0. Figure 5 .6 9 : X -Y side view o f vortex line. Re = 200, a = 1.0. V ' = 1.095. D om ain and grid sizes are th e sam e as in case V = 1. 0 . Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 5.70: Vortex line (nor malized view). Re = ‘ 200, a = 1.0, V — 1.095. Domain and grid sizes are the same as in case V = 1.0. Figure 5.71: V ortex line (real size view ). Re = 200. a = 1.0, V = 1.095. D om ain and grid sizes are the sa m e as in case V = 1.0 . an axisymmetric breakdown in experiments. If the flow in this second expansion were not stable to spiral disturbances, which of course will not exist in this 2D simulations because of the axisymmetric assumption, such spiraling breakdown might then develop. In chapter 10, the same profile will be carried out again for 3D simulation, and the spiral breakdown in the second expansion region is indeed observed. From the above qualitative discussion of the solutions for Re = 200, we can ar gue that the emergence of the adverse pressure gradient is critical for the formation of vortex breakdown. Later we will show a simple breakdown criterion proposed by Mahesh (1996) is also in this spirit. 5.2 The Effects of Swirl It is clear from the list of the solutions presented in the cases of Re = 200. o = 1.0 that the appearance of a large adverse pressure gradient in the vortex core, with a corresponding large decrease in the axial velocity, is critically dependent on the 63 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R e= 2 0 0 . V =1.20 1.5 I 0 .5 Z Figure 5.72: Streamlines contour. Case parameters Re = 200. q = 1.0. V = 1.20. Convergence history of c a s e R e*200, V o l.2 0 0.1 0 .08 > 0 .0 4 0.02 4 0 0 700 800 1000 200 300 500 600 900 100 T Figure 5.73: Convergence history of maximum velocity change. Case parameters Re = 200, a = 1.0, V ‘ = 1.20. R e= 2 0 0 , V = 1Z 5 1.5 r tr 0 .5 z Figure 5.7-1: Streamlines contour. Case [parameters Rt = 200. q = 1.0. V = 1.25. 64 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Convergence history ot case Re=20O. V-1.25 0.1 0.08 0.06 > 0 .0 4 0.02 600 700 800 900 1000 400 S00 100 200 300 T Figure 5.75: Convergence history of maximum velocity change. Case parameters Re = 200. a = 1.0, V = 1.25. magnitude of the swirl which is introduced at r = 0. Figure 5.76 shows the minimum axial velocity, (Vj)m, changes with V. As V increases from zero, where no axial retardation will occur, d{Vs)m/d V becomes increasingly negative so that for small V. small changes in V yield only small changes in but for larger V, particularly near V = O.SO, small changes in V result in very large changes in (V':)m. It was first noticed by Grabowski k Berger (1976). Their explanation can be summarized as follows. For small values of V, the effect of swirl diffusion is to increase the pressure on the axis and thus to cause a small amount of axial flow retardation. The small radial flux away from the axis which is required by continuity, has little effect on the swirl velocity distribution though. Therefore, the swirl velocity field is basically decoupled from the rest of the flow field and affects it only through the small adverse pressure gradient. For larger swirl, however, this decoupling is no longer possible. The axial velocity retardation arising from viscous diffusion of the swirl is large enough to require a significant amount of radial outflow and this will decrease significantly the swirl velocities near the axis and at the same time introduce large values of dVf'dz, thus supporting a much larger adverse pressure gradient and an increased reduction of the axial velocity. As V becomes larger than O .SO . this kind of coupling or nonlinear behavior is much more obvious. 65 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0 .7 0.6 0.5 0 .4 0.3 0.2 0.1 - 0.1 0.6 0 .7 1.2 Figure 5.76: The first minimum axial velocity on the axis vs. V 5.3 Uniqueness of Solution R e = 200 Recently, non-uniqueness has been observed when Reynolds number is large enough, and is manifested as the abrupt formation of the reversed flow as the vortex swirl strength is increased through a critical value. Beran & Culick (1992) presented numerical solutions of a swirling flow in a pipe using the axisymmetric and steady Navier-Stokes equations. They found that there exist two limit points of the incoming swirl level, u;,which connect three distinct branches of steady-state solutions. Along the solution branches, if u increases up the first limit point, this branch describes near-columnar flows along the pipe. Another solution branch is between the first limit point and the secondary limit point and the swirling rate is reduced. The last branch, with increasing swirling rate, describes a vortex breakdown flow. Although one of the main points of the present study is to try to understand more about the process of the bubble breakdown, it's still worth checking the solution uniqueness in this free vortex problem. By using continuation method 66 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. tin i R0 3 2 0 0 , V sO .8944 0.9 0 8 0.7 0.6 £ 0.5 T f > 0 4 0.2 0 .t Figure 5.77: Comparison of results using different initial conditions. Case parameters Re = 200. a = 1.0, V = 0.8944. that uses previous result for the next step calculation, we calculated the above cases again using the solution of V = 0.80, which is far away from the critical region, as the start point. Figure 5.77 shows the comparison of results using two different approaches. Results are nearly identical, which means the solution is unique. Also, Figure 5.78 shows a tim e history to obtain solution of V = 0.8944 using continuation method. Clearly, by utilizing the solution below V = 0.8944 as initial condition, very small disturbance is introduced initially near inlet region, comparing the relatively strong initial disturbance using direct simulation. To trace the process of vortex breakdown, we checked the velocity changes again as we did above. This time, only small disturbance was observed to develop and trapped inside the vortex core and form a finite amplitude standing wave near inlet region. As tim e goes on, the amplitude of that standing wave subsides. Eventually, the flow will converge to a new structure. Figure 5.79 again demonstrate that this process mainly happens inside the vortex core. The continuation method is carried on to a case with full bubble breakdown, V ’ = 1.095. and the result is compared with the one from direct simulation. Both results are identical as shown in Figure 5.80. The axial velocity changes in Figure 67 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. o V fiitH/‘(t’1) V=0.8944(Continuation Method) • 0.0 1 0 .0 2 -0.03 0.04 T=20 -0.03 hO.04 1 0 1 2 .5 15 1 7 .5 2 . 5 20 T=30 7 . 5 10 1 2 .5 15 1 7 .5 |0-6 l [____ Ta70 0 -------------------------- r 0 . 01 ! - 0 .0 2 -0.03 hO.04 . 2 . 5 5 7^5 2. 1-0 12 ♦ 5 15 1 7 .5 20 Figure 5.78: Time history of axial velocity change at r = 0. Case parameters Re = 200. o = 1.0. V = 0.S944. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0V 5i^V ^t1 * Snapshots at T=40 for V=0.8944(Continuation Method) - 0 .0 2 • rO.03 ■ | j-0.04 • | j 0 2 .5 5 7 .5 10 12. S 15 1 7 .5 20 j -o.oi 1 - 0.02 i -0.03 -0.04 10 1 2 .5 1 7 .5 15 r=1 -o.oi 1 - 0 .02 j-0.03 1-0.04 10 1 2 .5 1 7 .5 15 20 r=2 -o.oi 1 - 0 .0 2 0.03 1-0.04 Figure 5.79: Snapshots of axial velocity change at several different po sitions. Case parameters Re = 200. a = 1.0. V = 0.S944. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. o: Direct Simulation ♦: Contfuatlon Mathod 0.9 OS 0.7 0.4 0 2 Figure 5.80: Comparison of results using different initial conditions. Case parameters Re = 200, a = 1.0. V ’ = 1.095. 5.81 show again the accumulation of small disturbances inside the vortex core and the upstream movement of finite amplitude waves. Besides continuation method, we also perturbed the steady above solutions and carry the calculations to new steady solutions. Results show that the above solutions are unique at this relatively low Reynolds number 200. Recent studies for the pipe flow problems also show that the solution non-uniqueness will only happen at relatively high Reynolds number that is larger than 360 (Tromp & : Beran (1996), Beran Culick (1992), etc.). Although the solutions of Re = 200 are unique for this free vortex problem, we do observe the non-uniqueness for high Reynolds number in the present problem too. Figure 5.83 shows a case with Re = 3000. The new steady state from perturbed solution of same V = 0.80 shows some difference about the axial velocity at r = 0. This observation also verifies the existence of non-uniqueness for the high Reynolds number cases. 70 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. V=1.095(Contlnuatlon Method) o ------------------------------- j •o.oi f- : -0.02 [ | -0.03 [ j-0.04[__________________________ i _ i _ I I 0 2 .5 5 7 .5 10 1 2 .5 15 1 7 .5 20 j T=20 rO.oi 0 .0 2 --0.03 -0.04 10 1 2 .5 1 7 .5 15 20 T=40 -0.01 r0.02 -0.03 10 1 2 .5 1 7 .5 20 15 T=50 r O . O l -0.03 -0.04 2 .5 Figure 5.81: Time history of axial velocity change at r = 0. Case parameters Re = 200. a = 1.0. V = 1.095. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. j Q vj ^ v* (t* 1 * Snapshots at T=40 V=1.095(Continuation Method) -0.01 • 0.02 15 1 2 .5 10 1 7 .5 20 j-0.01 1 r O . 0 2 -0.03 1-0.04 1 2 .5 15 lO 1 7 .5 20 V .(t)-V I(t-1) o . 6 i i -o.oi 0 . 0 2 -0.03 0.04 2 .5 r=1 7.5 10 1 2 .5 15 1 7 .5 20 ! O.ilr i ° -0.01 - 0 .0 2 -0.03 | 1-0.04 2 .5 r=2 7 . 5 ^ 10 1 2 7 5 15 1 7 7 5 20 Figure 5.S2: Snapshots of axial velocity change at several different po sitions. Case parameters Re = 200. a = 1.0. V ’ = 1.095. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. RO3000, V«0.*0 (nonuniqua wkitlon) ♦: Simulation w/o dtaturbanca o: Simulation w dtatuibanca 05 Figure 5.83: Comparison of steady state results using dif ferent initial conditions. Case parameters Re = 3000. Q = 1.0, V = 0.80. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 6 Some Comparisons with Previous Studies It is worthwhile to compare the present solutions with some previous studies. 6.1 Comparisons with Experimental Observations 6.1.1 Breakdown Response to Increase in Swirl Level The previous breakdown solutions of Re = 200 with V = 0.S944 to V ’ = 1.095 are consistent with the following experimental observations (Harvey 1962; Sarpkaya 1971). When swirl is sufficiently large, only the bubble form is seen. Continued increase in swirl causes upstream movement of the bubble until it reaches the upstream boundary of the apparatus. 6.1.2 Swirl Angle As mentioned earlier, there is very little quantitative experimental data regarding vortex breakdown. The quantitative information which does exist has, for the most part, been obtained from photographs taken of flows, after the introduction of dye or smoke into the flows. In this way, Harvey (1962) found that the swirl angle. d > = tan~l(v/w), ahead of the breakdown was never larger than 50.5°, and in a later similar experiment Sarpkaya (1971a) found that it was never larger than 51°. Theoretical predictions of o. vary from 45* (Squire. I960) to 62.5* (Bossel 1967). depending on the velocity profiles used for analysis. 74 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. The maximum swirl angle ahead of breakdown for Re = 200 of V ’ = 1 .0 and V = 1.095 varies from 47.4° to 49.9°. Therefore, there is reasonable agreement in terms of the swirl angle between the experimental results, theoretical predictions and the present numerical solutions. Grabowski & : Berger (1976) had a similar conclusion. 6.2 Comparisons with Theoretical Approaches 6.2.1 Quasi-Cylindrical Approximation The theoretical approach to breakdown, that it is a boundary layer separation-like phenomenon, predicted by the failure of the quasi-cylindrical approximation, is presented here to compare it with numerical results. For initial velocity profiles of the same form as used here, Mager (1972) provided a very precise criterion for the appearance of the discontinuities, which he assumed were the results of such failures, in his integral solutions of the quasi-cylindrical approximated N'-S equations. He defined a parameter 91 which is related to the invariant flux of axial momentum deficiency in the vortex core, the swirl ratio. 9\ is determined from the initial upstream condition. Mager then showed that if the initial conditions, as determined from the values of a and V, are such that 9X is less than a certain value, say 9", then the solution is a continuous one; if, however, 9\ > 9\, a discontinuity will arise, marked by the appearance of the infinite gradients regardless of the sub- or supercriticality of the upstream flow. The value of 9\ was determined to be -0.163989 which is independent of Re because of the high Re assumption in the quasi-cylindrical equations. Values of 9X at breakdown for various Re are shown in Figure 6.1. It is ap parent that the 9X for high Re at breakdown is around -0.1822 which is close to theoretical predicted value 9\ = — 0.163989. This also verified the argument that Mager’s results should be applied to high Reynolds number flows. The differ ence between numerical and theoretical values can be explained by the assumption made by Mager that axial gradient is smaller than the radial gradient so that some components in the fully N-S equations can be eliminated. However, as we 75 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. - 0.1 - 0.12 4 I H " -0 .1 6 -0 .1 8 - 0.2 4000 1000 2000 Re 3000 0.9 0.86 > 0.84 0.82 - 0 ---- 0.8 3000 4000 1000 2000 Re Figure 6.1: The breakdown parameters vs. Re showed previously, large axial gradients do appear in the solutions that do not exhibit stagnation. Second, Mager's assumed profiles may not be general enough to describe the flow approaching breakdown as mentioned by Hall (1972). The good agreement between this theory and our numerical solutions show that this theory is good for predicting the onset of breakdown, however, as by the inherent shortcoming of this theory, it will not provide any information after breakdown. 6.2.2 Wave Theory As mentioned earlier, vortex breakdown is thought of as a large amplitude wave motion (Leibovich 1973, 1990), analogous in some sense to a gas dynamical shock wave (for its relation with small amplitude waves). The case of Re = 1000, V ’ = 0.82 is chosen for comparison with some features of vortex breakdown predicted by wave theory. The reason to choose relatively high Reynolds number and the swirl rate V = 0.82 corresponding to the onset of breakdown, is because of the inviscid and critical flow assumptions in the theory, so the case will be close to the critical flow in theory. Figure 6.2 through Figure 6.5 are the time evolving details of this case. Lei bovich (1973) predicts a model of vortex breakdown with the weakly nonlinear theory approach. In his steady vortex breakdown (figure 13 of Leibovich) . the 76 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. breakdown bubble is a well-defined cell of recirculating fluid which is symmetrical around the tube center-line. The final solution of our case is quite similar to this model. Figure 6.2: Stream- function contours for t = 500,700,720 and 740. Case parameters Re = 1000, V = 0.S2. Domain size: [50,10]: grid size: [501,61]. Figure 6.3: Stream- function contours for t = 750.760,7S0 and S00. Case parameters Re = 1000, V = 0.82. Domain size: [50,10]: grid size: [501,61]. In the latter paper (Leibovich 1990), Leibovich predicted the emergence of the recirculation region based on the relationships between a fully nonlinear standing periodic wavetrain and solitary waves and the underlying columnar flows. Com paring the present time evolving solutions with the result in his paper (Figure 6.6) , it is quite clear that both numerical and theoretical bubbles develop almost sym metrically with respect to fore and after. The numerical bubble's mean position moves upstream as its size becomes larger which is consistent with the experimen tal observation, while the theoretical bubble only expands along its mean position, even though the similarity between these two is still striking. The similarities be tween numerical and theoretical solutions support the theoretical approaches from 77 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. to t, 1 0 1 0 t to Figure 6.4: Stream- function contours for t = 820,S40,S50 and 900. Case parameters Re = 1000, V = 0.S2. Domain size: [50,10]; grid size: [501,61]. Figure 6.5: Streamfunction contours for t = 1000,1300,1600 and 2500. Case parameters Re = 1000, V = 0.82. Do main size: [50.10]; grid size: [501,61]. 78 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. the view of wave theory which may describe more detail structures about vortex breakdown. 6.3 Absolute/Convective Instability of Solutions R e = 200 Recently interpretation of the occurrence of vortex breakdown in vortices in the light of the instability studies is a quite active field (Lim & : Redekopp 1997; Lessen, Singh & Paillet 1974; Tsai & : Widnall 1980; Delbende, C'homaz Huerre 1996). As generalized by Delbende, Chomaz Huerre (1996), the now widely used con cepts of absolute/convective instabilities (A I/CI) can be connected to the notions of supercriticality and subcriticality (Benjamin 1962). In this framework, the na ture of the instability is determined by examining the linear impulse response of a given basic flow for large time. If the response wavepacket is convected away from the source location, the flow is said to be convectively unstable, i.e. supercriti cal. If it contaminates the entire medium both upstream and downstream from the source location, the flow is said to be absolutely unstable, i.e. subcritical. Based on our steady numerical solutions as basic flows, we can test the impulse responses of the flow fields to determine the absolute and convective nature of the instability. Loiseleux. Chomaz & Huerre (1998) also studied the linear instability of Rankine vortex with axial flow, their A I/CI transition provided a reasonable estim ate of vortex breakdown onset. However, they found that the A I/CI criterion predicts a transition to a predominant helical mode m = — 2, while fully m ature vortex breakdown is primarily associated with mode m = 0 (axisymmetric bubble breakdown) or m = ±1 (asymmetric spiral breakdown) (Sarpkaya 1971). They mentioned that to get results close to experimental observations, the AI/C'1 anal ysis on the profiles measured before and after vortex breakdown onset should be undertaken. The code we are using now can also simulate the transient procedures of breakdown, esp. around onset ( see Figure 6.2 through Figure 6.o ). So based on the velocity held from our numerical solutions, the AI/CI analysis might shed more light on the nature of vortex breakdown. 79 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0 a i t Figure 6.6: Streamlines of solitary-wave solutions and the emergence of the recirculation region (From Lei bovich Kribus 1990). (a) A/^0 o = 0.80. (b)A//r0o = 0.70. a small recirculation bubble appears.(c)A//(00 = 0.001. a large recirculation bubble, (dj Detail of the bubble in (c). S O Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Swirl will also affect the dynamics of a wide variety of free shear flows such as jets, wakes and mixing layers. Delbande, Chomaz and Huerre (1996) presented a novel procedure which determines the absolute/convective instability properties of the Batchelor vortex by direct numerical simulation of the linear impulse response. The analysis of the large time asymptotic response leads to the determination of the complex wavenumber and frequency along each spatio-temporal ray. They found that axisymmetric jets and wakes experience a transition from convective to absolute instability for moderate swirling levels. Below in the chapter about 3D simulations, we'll see that before the occurrence of the spiral breakdown, the 3D process will basically repeat the 2D case result first, and only after the axisymmetric breakdown zone is more or less established in the swirling flow, the three dimensional disturbances appear inside the breakdown zone and its wake( Rusak, etc. 1998). So it will be very interesting to examine the 2D Re = 200 solution instabilities to understand more about the evolution process from bubble to spiral breakdown. This would be a very crucial step for the bubble to spiral breakdown transition mechanism. Look as back again at the steady solutions of Re = 200, we found at each downstream position the axial velocity and swirl velocity profiles more or less look like the Batchelor’s vortex profiles described by Delbende, et al.(1996). So at here we made a rough assumption that our profiles are in the Batchelor’s vortex profile family to make the analysis easier at this point. By assuming the Batchelor’s vortex profile, we can calculate the two parameters a and q defined by Delbende, et al.(1996) as followings: (t’.-)co (tf*)(r = l) a = — .q = ----- --------- At'; A t’; Where (i>;)oo is the free-stream axial velocity and (u;)c is the centerline axial ve locity. AU; = (u;)oo — (^i)c measures the velocity difference. (ue)(r = 1) is the swirl velocity at the vortex core radius. Based on the definition of external flow parameter a and swirl parameter q (Delbende. et al. (1996)). we can get the AI/CI transition curve in a-q plane for each solution of Re = 200. 81 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 6.7 show the a-q curve for solution V = 0.85. The upper plot shows the distribution of parameter a along the downstream direction. Several distinct position are marked as .4, B ,C and D. The lower plot shows the a-q curve. Com pare this plot to the superimposed figure from Delbende. et al. (1996), Clearly, the whole a-q curve is in the convective instability region of Delbende's figure. 1.75 m = -1 (Delbende, etc. 1996) 1.25 0.75 -3 -1 -0 .5 -2 .5 -2 -1 .5 Figure 6.7: a-q parameter plane. Case parameters Re = 200, a = 1.0. V = 0.85. At V = 0.8944, notice the most dangerous point B is inside a little bit the AI region (Figure 6.8). Notice another distinct point D which corresponds to the second expansion region also moves to the AI region quite obviously. Last, when we increase swirl ratio to the V = 1.095, still the most dangerous point is B , which corresponds to a position near or inside the breakdown bubble. However, we also notice the point D , which corresponds to the expansion position after the bubble, moves inside the AI region very quickly (Figure 6.9). Later in the same parameter simulation in 3D studies, we found out that it is this D point residing region that brings up the three dimensional small disturbances that eventually will cause the occurrence of spiral breakdown. 82 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. tiT) i «i w w o i m I 1 .75 m = -1 (Delbende, etc. 1996) 1.2 5 0.7 5 -1 .5 - 0 .5 -3 -2 .5 -2 -1 Figure 6.S: a-q parameter plane. Case parameters Re = 200. a = 1.0. V = 0.8944. The studies of this kind of a-q transition curves reveal that moderate amount of swirl strongly promote absolute instability, esp. for the helical mode m = — I. which is the mode of spiral breakdown. Thus this study help us understand the process of spiral breakdown in the 3D simulations. The above study is based on the assumption that the profiles at each down stream position are in the family of the profiles specified by Delbende. etc. (1996). To get rid of this assumption, the same numerical procedure should be implemented to examine the solution flow field response of the linear impulse. S3 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1.75 m = -1 (Delbende, etc. 1996) 1 .5 1.25 0.75 -1 -0 .5 -1 .5 -3 -2 .5 -2 Figure 6.9: a-q parameter plane. Case parameters Re = 200. a = 1.0, V ' = 1.095. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 7 A M echanism for Axisymmetric Vortex Breakdown Based on the above basic observations, the mechanism about the process of the bubble breakdowns will be brought up at this point. Since the observation of vortex breakdown, a lot of efforts were made from both theoretical and numerical ways to explain it. However, most of previous analysis of swirling flows can't provide clear insight into it(Rusak Wang 1996). Wang & Rusak (1997) have recently presented a novel theoretical approach which describes the axisymmetric vortex breakdown process based on the unsteady Euler formulation in a straight pipe of finite length. They argue that this process is a result of the interaction between disturbances propagating upstream and the incoming flow when the swirl ratio is near or above the critical level. In an effort to try to explain axisymmetrical vortex breakdown from numerical way, the axial velocity changes inside and outside the vortex core were traced for the above simulation cases. Figure 7.1 shows the tracing of axial velocity changes every 1 time unit at the position r = 0 which is the center of the vortex core. The caption V(t) — V(t — 1) means the measurement of axial velocity changes at every time unit on the axis. At T = I, the initial disturbance has already developed near the inlet boundary, as time goes to 10. the amplitude of the initial disturbance becomes larger. This disturbance development meanwhile also changes the flow structure. The propagating waves thus break into two parts: the "fast moving” part is washed So Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. avvav quickly by the convection of the flow: the "slow moving"’ part is trapped inside the vortex core, however in this case, it eventually will be washed away by the flow convection too. As time goes on, the amplitude of the disturbance subsides. The flow will converge to a new state which is still the columnar structure at this case. Figure 7.2 tries to demonstrate that the disturbance will mainly propagate inside the vortex core. Four positions are traced, r = 0 is the vortex core center, r = 0.5 is the position in middle of vortex core, r = I is the vortex core edge and r = 2 is the position outside the vortex core one more length unit. At T = 40, a snapshot is made at all the above positions, clearly, the disturbance has the largest am plitude near vortex core center. Similar figures (7.3. 7.4 ) are shown at the case V = 0.8944 which is the critical state. Now in this near breakdown state, the slow moving part of the disturbance is trapped inside the vortex core and also move slowly upstream from T = 10 to T = 100. After that, it becomes a standing wave at a position near inlet, with the amplitude subsiding. Clearly, it’s the accumulation of small disturbance propaga tion inside the vortex core that change the base flow structure from supercritical flow to subcritical one where the upstream propagation wave can stand against the incoming flow from the inlet. For the case with a full bubble breakdown V = 1.095. again, in Figure 7.5 the upstream movement of part of the disturbance is quite clear from T = 10 to T = 100, eventually this relatively strong upstream part will interact with the incoming flow from the inlet and form the standing wave quite near inlet. Figure 7.6 shows again that the whole process that brings up the axisymmetric bubble breakdown mainly happens inside the vortex core. 86 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. V=0.85(before breakdown) T=1 - 0 .0 1 - 0 . 0 2 -0.03 -0.04 10 12.5 15 17.5 20 T=30 -o.oi - 0.02 -0.03 -0.04 10 12.5 15 17.5 20 T=100 -o.oi - 0 .0 2 -0.03 -0.04 12.5 15 17.5 Figure 7.1: Time history of axial velocity change at r = 0. Case pa rameters Re = 200. q = 1.0. I" = 0.S5. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. o.Oi 0 •0.01 -0.02 ■0.03 ■0.04 C V > (M ) Snapshots at T=40 for V=0.85 (direct simulation) | ---------- rsO------------------------------------------------ j 2.5 5 7.5 10 12.5 15 17.5 20 j V ,(» ) I 0.01 1 0 ■0.01 ■0.02 -0.03 -0.04 ■ V ,(t-1) r=0.5---------------------------------------------- (> 2.5 5 7.5 10 12.5 15 17.5 20 V .O ) o.Oi 0 -0.01 j-0.02 j-0.03 -0.04 -V ,(t-1) r=1 ! ( > 2.5 5 7.5 10 12.5 15 17.5 20 ' v.(t> 0.01 0 -0.01 -0.02 -0.03 -0.04 ■ V ,(t-1) r=2 • * . _ . . . . . . . f (> 2.5 5 7.5 2 M 12.5 15 17.5 20 | Figure 7.2: Snapshots of axial velocity change at several different posi tions. Case parameters Re = 200. a = 1.0. I = O.So. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. V=0.8944(near breakdown) T=1 -0.01 - 0 .0 2 -0.03 -0.04 15 14 12.5 17.5 20 T=10 -0.01 - 0 .0 2 -0.03 -0.04 14 12.5 15 17.5 T=30 -0.01 - 0 .0 2 -0.03 -0.04 12.5 15 17.5 20 V.dhV.tM) 1 .0 1 r 0 0 -0.01 - 0 .0 2 -0.03 -0.04 T=100 2.5 12.5 15 17.5 20 Figure 7.3: Time history of axial velocity change at r = 0. Case pa rameters Re = 200. a = 1.0. V = 0.S9-14. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. V 6 ( W - 1) Snapshots at T=40 for V=0.8944 (direct simulation) -----------------r=0----------------------------------------- o o l-o.oi - 0 .0 2 -0.03 -0.04 2 . 5 5 7 . 5 10 1 2 .5 15 1 7 .5 20 r O .o i -0.03 rO. 04 10 1 2 .5 15 1 7 .5 20 0.01 0 .0 2 0.03 0.04 2 . 5 7 . 5 10 1 2 .5 15 1 7 .5 20 r=2 o.oi -0.04 Figure 7.4: Snapshots of axial velocity change at several different posi tions. Case parameters Re = 200. a = 1.0. V ’ = 0.S944. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. v.(t) 0.01 0 -0.01 -0.02 -0.03 -0.04 !v^t' 1) V=1.095 (breakdown) T=1 V / " 2.5 5 7.5 10 12.5 15 17.5 20 V.(t) 0.01 0 -0.01 -0.02 -0.03 -0.04 ■v,(t*i) 1=10 — ) 2.5 5 7.5 lO 12.5 15 17.5 20 V.(t) 0.01 0 -0.01 -0.02 -0.03 -0.04 -V.tM) \ / " ~ \ T=30 2.5 5 7.5 10 12.5 15 17.5 20 j v,(t) 0.01 0 -0.01 -0.02 -0.03 -0.04 -V,(t-1) ! _ ___ T=100 0 2.5 5 7.5 Z lO 12.5 15 17.5 20 Figure 7.5: Time history of axial velocity change at r = 0. Case pa rameters Re = 200. o = 1.0. V " = 1.095. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ) Snapshots at T=40 for Vsl .095 (direct simulation) 0.01 - 0.0 2 -0.03 -0.04 1 2 .5 15 1 7 .5 20 i t=0.5 -0.01 |-0.03 1-0.04 1 0 15 1 7 .5 1 2 .5 V.(*)-V,(t-1) 0.01 l 0 i j-0.01 ■ ■0.02 |-0 .03 '-0.04 2 .5 r= 1 7 .5 10 1 2 .5 15 1 7 .5 20 V.(t)-VI(t-1) 0.6l i 0 -0.01 - 0 .02 -0.03 -0.04 2 .5 r=2 7 . 5 2 127 5 15 1 7 7 5 20 Figure 7.6: Snapshots of axial velocity change at several different posi tions. Case parameters Re = 200. o = 1.0. V = 1.095. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 8 Verifications of Other Breakdown Criteria In the previous sections of comparing with experimental and theoretical works, some breakdown criteria, like swirl angle, the predicted onset of breakdown by Mager (1972). are already checked with the finding of good agreement. However, those breakdown criteria are not general enough and. moreover, it is quite com plicated like the theoretical prediction. Recently, more and more researchers have focused on simple breakdown criteria with only the upstream inflow parameters (Brown k Lopez 1990, Mahesh 1996. etc.). Brown k Lopez’s (1990) criterion will not give any information about the onset of breakdown in the present simulated case. Because the initial azimuthal vorticitv in this case is equal to zero, the breakdown criterion relation > 30 will be satisfied automatically for all our cases. Billant. Chomaz k Huerre (1997) proposed a simple breakdown criterion in the spirit of Brown k Lopez's (1990) work which argues that a necessary condition for vortex breakdown is that H' > 0, where H is the total pressure head defined as P/p + + v'g + t>j). However, the criterion is only suitable for flows at rest at infinity, which cannot be applied directly to our cases. Mahesh (1996) proposed a very simple criterion that works for both incom pressible and compressible inviscid flows, which states that vortex breakdown only occurs if the axial pressure rise exceeds the upstream streamwise momentum flux, which is ;»x — p . ; i > p.jl 1 in our cases. Here py is the pressure at infinity. p0 and L'0 the pressure and the axial velocity at upstream station. The good thing about this universal criterion is that one can derive a criterion for a particular problem 93 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. in this spirit. So based on his criterion we can derive the breakdown criterion for our problem as That is r™ it J O ___ r U* > 1.0 or j r > 0.775 O where v$ is the swirl velocity and V is the measurement of swirl strength as we used in the previous chapters. Compared it with the numerical critical value of V/U0 = 0.S11 from §7.1. this breakdown criterion works well in our case. 94 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 9 Other Axisym m etric Solutions with High R e 9.1 A dynamical systems representation of the flow 0.5 0.4 0.3 0.2 3 0.1 ~ V sO .8944 V =1.0 - 0.1 - 0.2 1 000 1 2 0 0 1 4 0 0 1 600 1 600 2 0 0 0 T 8 0 0 200 4 0 0 6 0 0 Figure 9.1: The history of the axial velocity at point (2.5,0.2297). Re = 500, a = 1.0, V = 0.83. Domain size: [50,10]; grid size: [501,61]. V V e carried on the calculations by increasing the Re to 500 with various swirling rates V. All the solutions for V up to 1.5 are steady (see Figure 9.1) and no periodic or chaotic phenomenon happened below this Reynolds number. Here we will try to use a dynamical system to represent the flow and to un derstand some topological behaviors from another view (Lopez 1990). Following is the summary of his approach. In axisvmmetrical cases, from Lagrangian view. 95 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. an axially centered circular ring of fluid remains axially centered and circular as it evolves in time. The only effect of swirl is to rotate the ring about the axis of symmetry. Rather than following the motion the particles perform in three dimen sions, it’s sufficient to study the motion of ring intersection points in an arbitrary meridional plane (r, c), where (r.O .z) is the usual cylindrical coordinates. The motion of the intersection point then can be described as 1 d t( r , z ,t) . ldi/>(r.x,f) r = Vr = t , - = V. = --------r------ r dz r or where i/>(r, x. f , ) is the Stokes function. By the coordinate transformation p = r2/'2 (Benjamin 1962), the above system becomes d ll' . On' P ~ ~ l h ' = ~ V which constitutes a Hamiltonian system of one degree of freedom. The velocity field (/?, i) is divergence free and hence area preserving in the canonical (p.z) phase space. Based on the above dynamical system about the particles, Lopez defined a Poincare map to exploit the temporal periodicity (r(<0),x(f0)) !-> • (r(f0 + r), z(tQ + r)). So the Poincare map maps the position of ring intersection points in (r,z) plane at tim e t0 to their new locations one period later. The image of a point x under one application of the map is denoted by F (x). A point x. such that F {x) = x, is a fixed point of the Poincare map. All possible maps can be obtained by choosing 0 < t0 < t and the topology or phase portrait of the map is independent of the choice made. Lopez also redefined the streaklines in the axisymmetrical flows as the positions of intersection points of rings rather than the positions of particles. In this sense, a point on a streakline is always mapped to another point on the same streakline and a streakline pattern is a phase portrait of the Poincare map. 96 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9.2 Topologies of the steady flows for R e = 500 Let’s analyze the steady results from Re = 500 (see Figure 9.2 to Figure 9.13 ). For the steady swirling flows, the only fixed points of the map which correspond to particles which never move in (r, 6, z) space, are the stagnation points on the axis of symmetry and some boundary points. All other hyperbolic and elliptic fixed points of the map occur off the axis where the azimuthal velocity is non-zero and they hence are periodic points of the steady flow. The time taken to complete one orbit is determined by the local azimuthal velocity. The development with increasing V of the steady flow’s topologies under Re = 500 in the area not far from axis (i.e. r < 2) is listed here. For the steady cases the streamfunction contours give a phase portrait of the Poincare map, of which they are invariant curves. For V < 0.835. there is no fixed point of Poincare map. For 0.835 < V < 0.8944, the topology of the Poincare map consists of two hyperbolic fixed points (also called saddle points of the dynamical systems) on the axis and one elliptic fixed point(also called center point) off the axis. Figure 9.14 through Figure 9.17 is the schematics of this. There are four special invariant curves originating or terminating at each fixed point. In the terminology of dynamical systems, these are the stable and unstable manifolds. Points on the two stable manifolds. IF3, asymptotically approach towards the fixed points under repeated applications of the map. Points on the two unstable manifolds, IT'1 1 , asymptotically move away from the fixed points under repeated applications of the inverse map. For V > 0.8944, the topology of the Poincare map consists of four hyperbolic fixed points on the axis, each of which has two stable and two unstable manifolds, and two elliptic fixed points off the axis. A separatrix connecting the two hyperbolic fixed points on the axis, corresponding to the two stagnation points of the flow, is an unstable manifold for one and a stable manifold for the other. As shown in the above figures, there are two bifurcations of V in the topology of Re = 500 flows. Around V = 0.835, flow pattern changes from no bubble region to single bubble region inside: and after V " = 0.87. two bubble regions begin to form and more fixed points of Poincare map show up. At I ’ = 1.05. two bubble regions move closer and the two fixed points near the entrance move towards each other 97 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. and eventually will form a single critical point, while the two saddle points in the other bubble will separate farther. This new critical point is called nonhyperbolic region is called an elliptic sector. After V ’ = 1.1, flow pattern becomes another one which has one recirculation ring, one bubble and a new saddle point near the ring. Although the flow pattern is the same, there are slightly different details. Note the level ip = 0.00093, which pass around the ring before V = 1.1, now pass through the ring after V = 1.1. From the descriptions above, changes in the swirl rate cause the flow field pattern changes by generating more hyperbolic points, e.g. saddle or center, with the birth of a bubble or a ring. The basic flow structure for Re = 500 is one or two bubbles, or a bubble and a ring. critical point in the terminology of dynamical systems, and the whole first bubble t O 1 2 U 1 6 1 1 2 0 Z 1 S S 1 a s Figure 9.2: Streamfunction contour Re = 500, a = 1.0. V '" = 0.83. Domain size: [50,10]; grid size: [501,61]. Figure 9.3: Streamfunction contour Re = 500. q = 1.0, V = 0.835. Do main size: [50.10]; grid size: [501,61]. 98 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R ta5 0 0 .V « O .S S RtaSOO.V>O.S7 Figure 9.4: Streamfunction contour Re = 500, a = 1.0, V = 0.85. Domain size: [50,10]; grid size: [501.61]. Figure 9.5: Streamfunction contour Re = 500, a = 1.0, V '' = 0.87. Domain size: [50,10]; grid size: [501.61]. RtaSOO. VaQ.S9*4 RtaSOO, V»0.W Figure 9.6: Streamfunction contour Re = 500, a = 1.0, V = 0.8944. Do main size: [50,10]; grid size: [501,61]. Figure 9.7: Streamfunction contour Re = 500, a = 1.0, V = 0.95. Domain size: [50,10]; grid size: [501,61]. 99 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A *-500. V at.O feaflC O . V a t.0 5 Figure 9.S: Streamfunction contour Re = 500, a = 1.0, V = 1.0. Domain size: [50.10]; grid size: [501,61]. Figure 9.9: Streamfunction contour Re = 500, a = 1.0, V = 1.05. Domain size: [50,10]; grid size: [501.61]. n — 500. V a t 1 fe aflO O .V a t.lS Figure 9.10: Streamfunction contour Re = 500, a = 1.0. V = 1.1. Domain size: [50,10]; grid size: [501,61]. Figure 9.11: Streamfunction contour Re = 500, a = 1.0, V = 1.15. Domain size: [50,10]; grid size: [501,61]. 100 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. o 2 « < a io ts 14 te ic 20 o 2 4 a a 2 10 12 14 ie ta 20 2 Figure 9.12: Streamfunction contour Re = 500, a = 1.0, V — 1.2. Domain size: [50,10]: grid size: [501,61]. Figure 9.13: Streamfunction contour Re = 500, a = 1.0, V ' = 1.3. Domain size: [50,10]; grid size: [501.61]. w — joo. v-o.cs Ri.Mft v-ots ............................................ j ■ * O J ! i J „ s , ! / \ " i M c , s , s , A t-r-7, —. ffrrt . . . . 1 Figure 9.14: Schemat ics Re = 500, q = 1.0, V = 0.835. Domain size: [50,10]; grid size: [501,61]. Figure 9.15: Schematics Re = 500, Q = 1.0, V = 0.85. Do main size: [50,10]; grid size: [501.61]. 101 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. fW floo. v-aar « W M 0 .V - 0 tt* 4 Figure 9.16: Schematics Re = 500. a = 1.0, V = 0.S7. Do main size: [50,10]; grid size: [501,61]. Figure 9.17: Schemat ics Re = 500, a = 1.0. V = 0.8944. Domain size: [50.10]; grid size: [501,61]. 102 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 9.IS: Schematics Re = 500. q = 1.0, V = 0.95. Do main size: [50,10]; grid size: [501.61]. Figure 9.19: Schematics Re = 500. q = 1.0. V = 1.0. Domain size: [50,10]; grid size: [501.61], x 2 X Z Figure 9.20: Schematics Re = 500, q = 1.0, V = 1.05. Do main size: [50,10]; grid size: [501.61]. Figure 9.21: Schematics Re = 500, a = 1.0, V = 1.1. Domain size: [50,10]; grid size: [501.61]. 103 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. IU .6 0 0 V .1 IS (W -M 0 .V .I2 £ 0 » X 2 Figure 9.22: Schematics Re = 500, a = 1.0, V = 1.15. Do main size: [50,10]; grid size: [501.61]. o. 0 s 0 0 2 Figure 9.23: Schematics Re = 500, a = 1.0. V = 1.2. Domain size: [50,10]; grid size: [501,61]. R e= 5 0 0 , V » 1 .3 0 .0 3 5 1 .0 1 0.8 0.6 1.0193 0 .0 1 9 3 j 0 .0 1 0 .4 0.2 14 Figure 9.24: Schematics Re = 500. ft = 1.0, V = 1.3. Domain size: [50.10]: grid size: [501.61]. 104 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9.3 Topologies of the steady flows for R e = 1000 The above results for Re = 500 demonstrate the interesting and complicate phe nomenon of swirling flow under high Re number. To further exploit its behaviors. Re was increased to 1000 and more interesting results including periodic solu tions have been obtained. The steady solutions are investigated first here. All the Reynolds numbers in the following cases are set to be 1000 and the swirling rates are changed slowly from V = 0.81 to V = 1.05. All the solutions for V < 1.05 are steady ones (see Figure 9.25 through Figure 9.33). The flow patterns can be analyzed through their topologies. Around V = 0.81 the flow pattern experienced first bifurcation from no reverse flow region to a new one bubble structure which is constructed of two saddle points and one center point rotating in the clockwise direction. Around V = 0.835. the flow pattern bifurcates to another one which had four saddle points and two center points, and formed two bubble regions. Increasing V ’ to 0.85. the two bubbles become larger and both move towards the entrance while the distance between them becomes smaller. Around V = 0.875 the flow bifurcates again and now contains three bubbles, which are formed by six saddle points and three center points rotating in the same direction. This is a quite interesting result. To our knowledge, it’s the first time to observe the three bubble structure in numerical simulations. Increasing V moves all bubbles towards entrance, gets them closer each other. The two bubbles near entrance shrink on the axis while the third one expands. Notice that at V = 1.03 the first bubble's two saddle points already collapsed to a single point on the axis. Figure 9.35 through Figure 9.43 plots the schematics of the topologies for the steady ones under Re = 1000. Above all, the changes of flow pattern for high Re is associated with additional bubbles, which bring more critical points to the flow fields. We will see in next section that there also is ring and bubble structure for high swirl rate at Re = 1000, which is similar to the cases of Re = 500. 105 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 9.25: Streamfunction contour Re = 1000, a = 1.0, V = 0.81. Domain size: [50,10]; grid size: [501,61]. Figure 9.26: Streamfunction contour Re = 1000, a = 1.0. V ’ = 0.82. Domain size: [50,10]; grid size: [501,61]. Figure 9.27: Streamfunction contour Re = 1000, a = 1.0, V = 0.835. Do main size: [50,10]; grid size: [501,61]. Figure 9.28: Streamfunction contour Re = 1000, a — 1.0, V = 0.85. Domain size: [50,10]; grid size: [501,61]. 106 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. f W ta o a . V -O J7S p — m o . v - g a e Figure 9.29: Streamfunction contour Re = 1000, a = 1.0, V = 0.875. Do main size: [50.10]; grid size: [501,61]. Figure 9.30: Streamfunction contour Re = 1000, a = 1.0, V = 0.95. Domain size: [50,10]; grid size: [501,61]. z z Figure 9.31: Streamfunction contour Re = 1000, a = 1.0, V = 1.0. Domain size: [50,10]; grid size: [501,61]. Figure 9.32: Streamfunction contour Re = 1000, a = 1.0, V = 1.03. Domain size: [50,10]; grid size: [501,61]. 107 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. z Figure 9.33: Streamfunction contour Re = 1000, a = 1.0, V = 1.05. Domain size: [50,10]; grid size: [501,61]. Figure 9.34: Power spectra of the axial velocity at (2.5,0.2297) for V as indi cated. The spectra have been taken over 20 < t < 4000 using 21 7 samples. 108 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. te a ia o o . v- o j i 2 Figure 9.35: Schematics Re = 1000, a = 1.0, V ' = 0.81. Do main size: [50,10]; grid size: [501,61]. c Figure 9.36: Schematics Re = 1000, o = 1.0, V = 0.82. Do main size: [50,10]; grid size: [501,61]. 109 R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. Ito-tflOO. Figure 9.37: Schematics Re = 1000, o = 1.0, V = 0.S35. Do main size: [50,10]; grid size: [501,61]. fWiaoa i a* a t a.4 a * ° 0 2 S. 4 St • I St « s , 1 2 1 4 1 « I* » Figure 9.38: Schematics Re = 1000, a = 1.0, V = 0.85. Do main size: [50,10]; grid size: [501,61]. fe.iaoo. v«ojn c Figure 9.39: Schematics Re = 1000, a = 1.0, V = 0.875. Do main size: [50,10]; grid size: [501,61]. fttelOOO. V «4M 44 S Ts7s3 ioStiT Figure 9.40: Schemat ics Re = 1000, a = 1.0, V = 0.8944. Domain size: [50,10]; grid size: [501,61]. 110 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. <to-iaoa v*i o c Figure 9.41: Schematics Re = 1000, Q = 1.0, V = 1.0. Do main size: [50,10]; grid size: [501,61]. A a » 1 0 0 0 .V * 1 .0 3 s Z Figure 9.42: Schematics Re = 1000, a = 1.0, 1/ = 1.03. Do main size: [50,10]; grid size: [501,61]. R»tO OO , V>1.05 0.8 0.6 < E 0.4 0.2 1.00093 Z Figure 9.43: Schem atics Re = 1000, o = 1.0, V = 1.05. Do m ain size: [50.10]; grid size: [501,61]. I ll R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9.4 Topologies of temporally periodic flows for R e = 1000 In recent papers about confined vortex breakdowns, Lopez (1990) gave some cases for time-periodic solutions at high swirling rates which also set the Reynolds num bers to be high values in his definition. This is a really interesting phenomenon for vortex breakdown. Here in our unconfinecl problem, we have also observed the periodic behaviors for high swirl rate. Figure 9.44 to Figure 9.46 show the results for case of V = 1.05. The flow field is still steady. However strong oscillations are showing up. It took a long time to completely dampen these oscillations. This is the last steady case we can get for Re = 1000. It indicates the onset of periodic flows. Power spectra of the axial velocity at point (2.5.0.2297) for 0.8944 < V < 1.05 are shown in Figure 9.34. It verifies that the dominant frequency is zero for all steady flows. Finally, when the swirling rate reached V — 1.11, the solution is obviously periodic and the length of a period is about Tp = 42.5 (see Figure 9.47). Figure 9.49 to Figure 9.61 list the changes of instantaneous streamfunction contours in one period. Based on these, the changes of flow topology are plotted in Figure 9.62 through Figure 9.73. During a typical period, the flow structure experiences the following changes. The bubble and ring sizes shrink and expand periodically. From t = 3880 to t = 3892, the first ring (nearest to entrance) shrinks, then it expands until t = 3916. After that, it shrinks slowly again to the same size as t = 3880. The second ring also experiences the same size change with phase shift from the first one. The bubble size also changes slightly, but not as obvious as the sizes of those rings. With the size changes the positions and the critical points of rings and bubble move periodically. The movement of bubble is accompanied by the movements of critical points along the axis, while the movements of rings are with saddle points and center points away from the symmetrical line. The movements of the saddle points inside the rings cause the coalescing and separating of instantaneous streamlines. To demonstrate that, a typical streamline v = 0.00093 is traced in 112 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. one period. From t = 3884 to t = 3888 the curve of x b = 0.00093 inside the bubble breaks into two parts, and at t = 38S8 the saddle point forms in this curve. After that, the broken curve merges again and forms a smooth curve passing around the rings and bubble. Increasing V up to 1.3, solutions are still periodic (see Figure 9.74 and Figure 9.75). However, from the power spectrum density (PSD) plots (Figure 9.78) more high frequency components shows up to play roles. More high frequency compo nents will come out as the V value increases. This higher frequency eventually will bring the flows from periodic oscillating state to chaotic state. The solutions for two high rotation rate V = 1.5 and V = 2.0 basically are chaotic ones (see Figure 9.76 and Figure 9.77 ). The coalescing and separating of streamlines inside rings and bubble with fixed frequency characterize the periodic behaviors of high Re flows with large swirl rate. The increasing of swirl rate will bring more high frequency components into the flow field. 113 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Axial Vel. change at point (z,r)=<2.5.0.2297) -0 .0 5 N > - 0 .t 5 - 0.2 3300 1300 1700 2500 2900 3700 1 0 0 500 900 2100 t Radial Vel. change at point (z.rM 2.5.0.2297) -0 .0 5 > -0 .1 5 - 0.2 1300 1700 2500 2900 3300 3700 100 500 900 2100 Figure 9.44: The convergency history of the axial and radial velocity at point (z,r)=(2.5,0.2297). Case parameters Re = 1000, q = 1.0, V = 1.05. Domain size: [50,10]; grid size: [501,61]. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Maximum axial Vel. change 0.8 0.4 0.2 1500 2000 2500 3000 3500 4000 t x iq'*2 Maximum divergency nistory 2 2 500 1000 1500 2000 t 2500 3000 3500 4000 Figure 9.45: The convergencv history of maximum veloc ity change and maximum local divergency checked every 5 nondimensional time. Case parameters Re = 1000. a = 1.0. V = 1.05. Domain size: [50.10]; grid size: [501,61]. 2 , — i ------------------1 --------- r 0 2 4 6 8 10 12 14 16 18 2 0 Z Figure 9.46: Stream function contour at t = 4000. Case parameters Re = 1000. a = 1.0. V = 1.05. Domain size: [50.10]: grid size: [501.61]. R eproduced with permission of the copyright owner. Further reproduction prohibited without permission. ' 4 port IX/H 2 S 02397) > i at p o rt U./W 3 4 0 n * n t Figure 9.47: The convergency history of the axial and radial velocity at point (z,r)=(2.5,0.2297). Case parameters Re = 1000, a = 1.0, V = 1.11. Do main size: [50,10]; grid size: [501,61]. Figure 9.4S: The period of the axial and radial velocity at point (z.r)=(2.5.0.2297). Case parameters Re = 1000. q = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501.61]. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. ^ z 2 Figure 9.49: Stream func tion contour at t = 3880. Case parameters Re = 1000. Q = 1.0, V = 1.11. Do main size: [50,10]; grid size: [501,61]. Figure 9.50: Stream func tion contour at t = 3884. Case parameters Re = 1000, Q = 1.0, V ' = 1.11. Do main size: [50,10]; grid size: [501.61]. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0 2 i « • 10 12 <4 1* t« z Figure 9.51: Stream func tion contour at t = 3888. Case parameters Re = 1000. Q = l.o, V = 1.11. Do main size: [50,10]; grid size: [501,61]. M M Figure 9.53: Stream func tion contour at t = 3896. Case parameters Re = 1000, a = 1.0, V = 1.11. Do main size: [50,10]; grid size: [501,61]. Figure 9.52: Stream func tion contour at t = 3892. Case parameters Re = 1000, a = 1.0, V = 1.11. Do main size: [50,10]; grid size: [501,61]. M H O Figure 9.54: Stream func tion contour at t = 3900. Case parameters Re = 1000, Q = 1.0, V = 1.11. Do main size: [50,10]; grid size: [501,61]. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 9.55: Stream func tion contour at t = 3901.25. Case parameters Re = 1000. a = 1.0. V = 1.11. Do main size: [50,10]; grid size: [501,61]. Figure 9.56: Stream func tion contour at t = 3904. Case parameters Re = 1000. a = 1.0. V = 1.11. Do main size: [50,10]; grid size: [501,61]. Figure 9.57: Stream func tion contour at t = 3908. Case parameters Re = 1000, a = 1.0, V = 1.11. Do main size: [50,10]; grid size: [501,61]. Figure 9.58: Stream func tion contour at t = 3912. Case parameters Re = 1000, Q = 1.0, V = 1.11. Do main size: [50,10]; grid size: [501,61]. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 9.59: Stream func tion contour at t = 3916. Case parameters Re = 1000, a = 1.0, V = 1.11. Do main size: [50,10]; grid size: [501,61]. Figure 9.60: Stream func tion contour at t = 3920. Case parameters Re = 1000, a = 1.0, V = 1.11. Do main size: [50.10]; grid size: [501,61]. (=3922.5 0 2 « 6 8 10 12 1* 16 18 20 Z Figure 9.61: Stream function contour at t = 3922.5(same as t = 3880). Case parameters Re = 1000, a = 1.0, V' = 1.11. Domain size: [50,10]; grid size: [501,61]. 120 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. c I c z Figure 9.62: Schematic at t = 3880. Case parameters Re = 1000, a = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. Figure 9.63: schematic at I = 3884. Case parameters Re = 1000. q = 1.0, V = 1. 11. Domain size: [50,10]; grid size: [501,61]. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. s Z Figure 9.64: schematic at t = 3888. Case parameters Re = 1000, a = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. c z Figure 9.65: schematic at t = 3892. Case parameters Re = 1000, q = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. M M S Z Figure 9.66: schematic at t = 3896. Case parameters Re = 1000, q = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. M H O K Z Figure 9.67: schematic at t = 3900. Case parameters Re = 1000, o = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. 1 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. K 2 Figure 9.68: schematic at t = 3904. Case parameters Re = 1000, q = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. c 2 Figure 9.69: schematic at t = 3908. Case parameters Re = 1000, o = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. c 2 s Figure 9.70: schematic at t = 3912. Case parameters Re = 1000, o = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. Figure 9.71: schematic at t = 3916. Case parameters Re = 1000, q = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. 123 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. c z s z Figure 9.72: schematic at t = 3920. Case parameters Re = 1000, ct = 1.0, V = 1.11. Domain size: [50,10]; grid size: [501,61]. Figure 9.73: schematic at t = 3922.5. Case parame ters Re = 1000, ft = 1.0, V — 1.11. Domain size: [50,10]; grid size: [501,61]. D M V * . C M n p « p o r t U J X i 1 0 » * > ) i s I Figure 9.74: The period of the axial and radial velocity at point (z,r)=(2.5,0.2297). Case parameters Re = 1000, a = 1.0, V = 1.20. Domain size: [50,10]; grid size: [501,61]. A M M ftt cfanptfparairrvdtoavn 5 t m o m a o Figure 9.75: The period of the axial and radial velocity at point (z,r)=(2.5,0.2297). Case parameters Re = 1000, ft = 1.0, V = 1.30. Domain size: [50,10]; grid size: [501,61]. 124 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. A t t f i V tt t f w n » * « l p o r t ( l » H 2 5 0 2 2 9 7 ) 02f > I 3 I Figure 9.76: The period of the axial and radial velocity at point (z,r)=(2.5,0.2297). Case parameters Re = 1000, a = 1.0, V = 1.50. Domain size: [50,10]; grid size: [501,61]. - - * 4 t 5 Figure 9.77: The period of the axial and radial velocity at point (z.r)=(2.5.0.2297). Case parameters Re = 1000. a = 1.0, V = 2.0. Do main size: [50,10]; grid size: [501,61]. 125 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. R e * io o o .v » i.u 0 0.1 0 .2 0 .3 0.4 0.5 R e » l0 0 0 ,v » 1 .2 R e«1000.V -1.5 a . - i o 0.1 0 .2 0.3 0.4 0.5 Frequency -1 0 -1 5 -20 0.2 0.4 0.5 0.1 0.3 0 -10 -1 5 -20 0.2 0.4 0.1 0.3 R e«1000.V «2.0 Frequency Figure 9.78: Power spectra of the axial velocity at (2.5,0.2297) for V as indicated. The spectra have been taken over 20 < t < 4000 using 21 7 samples. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 10 Parallel 3D Simulations of Vortex Breakdown In this chapter, vortex breakdown is examined again for the same parameters we used in axisymmetric cases. Typically, we focus on the case Re = 200, swirl strength V > 1.095, which are cases of strongly-swirling flows. Tromp (1995) observed strong 3D reverse flows only for strongly-swirling channel flows. At Re = 200, the three-dimensional breakdown is still laminar, so the main features of the flow are easily discernable. The boundary and initial conditions are the same as we used in 2D cases, except that the axisymmetric boundary condition on r = 0 is removed. P a ra lle l s p e e d u p o n T 3E (S D S C ) . . . 1 1 1 1 N«1 I I 1 I C|: G rid ^ lz * 6S’fl* S 1 3 I I - I . I I I _ l. I L i i i _ i. i i i _ i. i I ’ V N - a 1 * • 1 \ 1 1 1 \ 1 1 1 1 1 1 1 i X i i ' v i i i ---------- 1 ----------- r V - i n j — i----------- r y j m t s -------- r ---------1 ------------ ---------- j----------- ' _ b ^ _ | ----------- j-------- r 1 • 1 — 1 1 1 1 — 1 1 1 1 1 1 1 1 1 1 X n m | ! ^ J c - 3 2 ! 1 1 1 1 \ | 1 1 1 1 'V 1 1 ^ i s N - « 4 1 1 1 ^ 1 1 1 1 1 1 1 _______1 _______ 1 _______ j_______ 1 _______ 1 ____ X * i ----------- 1 1 1 1 1 _______ 1 _______ 1 _______ 1 _______ 1 _______ 1 _____ 3 .5 03 1 3 2 2 3 lo g (N ) 3 3 4 3 Figure 10.1: Code Scalability. 127 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. It’s well known that 3D DNS is quite CPU time and memory consuming. For example, for a typical 2D case of Re = 200, V = 1.095 with grid size 193 x 41, the typical CPU tim e in our local parallel machine SGI 0rigin2000 with 64 CPUs is about 2 hours for one CPU. However, carrying the same case to 3D simulation with grid size now 193 x 41 x 61, the CPU time is about 19 days if just a single CPU is used! So it’s quite necessary that the numerical code should be parallelized to utilize more efficiently the resources of some parallel machines which are more popular and easier to be accessed than before. A message passing interface (MPI) technology is employed here. This MPI technology recently has become popular and m ature with the advantages of hiding machine dependent parts and leaving the Fortran programmer the high level, machine independent function-like calls. The code written with MPI is portable for different parallel machines. The present numerical code is tested in several parallel machines, including SGI 0rigin2000 parallel machine with 64 CPUs and 16GB shared memory. IBM SP2 with 30 CPUs and 12SMB per CPU distributed memory, HP Convex with 16 CPUs and S GB main shared memory. When porting the code to different platform machines, only a few lines about machine-dependent part in the header file need to be changed. Figure 10.1 shows a test result for two different grids on a parallel machine T3E in San Diego Supercomputer Center, demonstrating that present numerical code is well scalable to large CPU numbers. 10.1 Transition from bubble to spiral breakdown and its mechanism Since the first observation of vortex breakdown phenomenon, a lot of research efforts have been put on this from theoretical, experimental and numerical ways. However, most of numerical works (Kopecky & Torrance 1973, Grabowski Berger 1976, Spall, Gatski Ash 1990) are based on the axisymmetry tim e independent assumption. Only recently, time dependent 3D simulations (Tromp 1997, Lucas 1997. etc.) were reported and some (Lucas 1997) are quite successful in simulating the bubble, spiral breakdowns and the transition between them. Still, because of 128 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. the CPU tim e and memory consuming in 3D simulations, the 3D vortex breakdown mechanism is not studied well. By utilizing the parallel computing technique and the highly accurate numer ical scheme (Verzicco Sc Orlandi 1996), we have a chance to examine 3D vortex breakdown carefully and hope to shed some light on the mechanism about spiral breakdown and the transition from bubble type. Here, the same parameter Re = 200. a = 1.0. V = 1.095 used in 2D case is carried out to 3D simulation. Figure 10.2 shows a comparison between 2D and 3D about dv, which is the maximum velocity change in temporal domain. There is an interesting point to note in this result. During 0 < T < 600. it’s quite obvious from this figure that the 3D case basically repeats the 2D process and develops the bubble type breakdown first as happened in 2D case. After T > 600, the three dimensional disturbances come up and bring flow field to a new structure. Notice that, in Figure 10.2. the instability happens very quickly for the 3D case after T = 500. To understand more about the procedure of this instability, instead of tracing the maximum velocity change in the whole flow field, the temporal fluctuations of the three velocity components dvr, dvz, dvg at c = 1. r = 0.0153, 0 = 0° are investigated. The results in semi-logarithm form are plotted in Figure 10.3. From results, it’s clear that, in a small region near T = 500, the velocity changes grow almost exponentially as function of dimensionless time. T. Figure 10.4 to Figure 10.6 show the time sequence of streaklines in the flow field. Initially, a circle of particles are placed on the position near inlet center at every 36 degree. As time goes on, we first observe the development of bubble breakdown as in 2D case. A short tim e after the bubble type breakdown is fully developed, the three dimensional disturbances occur in the wake region of the bubble. At T = 600, small disturbances are observed in the bubble wake region. After T = 650, this kind of small disturbances develop into a large finite amplitude movement and bring up the spiral breakdown. Figure 10.7 to Figure 10.9 show the trace of axial velocity changes near r = 0 line. From this figure, before T — 600. the disturbance wave movements have 129 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. RO200, Vx1.095, oz1.0 o: 20, ♦: 30 0 . 1 s 0.05 •0.05 1500 1000 500 Figure 10.2: Temporal fluctuation of the maximum ve locity change. Case parameters Re = 200, a = 1.0, V = 1.095. similar behavior as in 2D case. Basically, the downstream moving part is washed away by the convection of the main flow and the upstream moving parts move towards the inlet while subsiding in amplitude. Between T = 200 and T = 600, the flow field is in a relatively steady state with the developed bubble breakdown near inlet region. At T = 600, we can see that some small disturbances occur first in the region around Z = 5. Recall the instability analysis about the 2D solution for the same parameters from Chapter 6. The secondary most dangerous point D is in the region around Z = -5 , the 2D instability analysis is consistent with 3D observation. After T > 600, small disturbances develop into a large finite am plitude wave movement in most of the streamwise region, revealing that a new flow structure occurs in the flow field. Recalling the local A I/CI analysis in the Chapter 6 again about the axisymmet ric solution of Re = 200, V = 1.095, in the finite local AI region around Z = 5, we observed the spiral breakdown system exhibiting self-excited oscillations at specific 130 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. I p n g I n ox mi i I cx ti i n p « i f i Figure 10.3: Temporal fluctuation of the three velocity component changes at 2 = 1. Case parameters Re = 200, a = 1.0, V = 1.095. 131 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. frequencies shown in the figures 10.10 to 10.13. It’s globally unstable and in the later tim e period like 1500 < T < 2000. the system even exhibits some irregular behaviors (see Figure 10.11 and Figure 10.13). From our simulation of this vortex breakdown case, we suspect that this kind of self-excited resonance is due to the presence of a finite region of local AI, which in turn is formed after the axisym metric bubble breakdown. In this sense, it demonstrates the im portant role that the axisymmetric breakdown plays in the emerging of the spiral breakdowns. 132 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 10.4: Time history of streaklines, T =5, 50 and 400. Case pa rameters Re = 200, a = 1.0, V = 1.095. The gaps in the lowest figure are due to the attractions of the rear stagnation point. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. m Figure 10.5: Time history of streaklines, T=500, 600 and 650. Case parameters Re = 200, a = 1.0, V = 1.095. 134 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. |[u A B J Z Z it ~ m , i Figure 10.6: Time history of streaklines, T=750, 850 and 1000. Case parameters Re = 200, a = 1.0, V = 1.095. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. V * (t> - 0.01 0 -O.Ol -O .02 -0.03 -O . 04 V * (t' 1) 3 D ,V = 1 .0 9 5 i i i 0 2.5 5 7.5 1« 12.5 15 17.5 20 V.(t) 0 .01 0 -0.01 -0 .02 L o .03 -0 . 04 *V,(M) j C 2.5 5 1 .5 1-0 12.5 15 17.5 20 { 0.01 0 -0 . 0 1 -0 .02 -0.03 -0 .04 ( 2.5 5 7.5 1 < ) 12.5 15 17.5 20 v.(t) 0.01 0 -0.01 K0.02 |-0 . 03 1-0.04 -V z (t-1) 1 < ) 2.5 5 7.5 7 1 < > 12.5 15 17.5 20 ...... ............... Figure 10.7: Sequence of axial vel. changes (a), Case parameters Re = 200, a = 1.0. V = 1.095. 136 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3D , V =1.095 T=100 -o .o i -0 .02 -0 .03 0 .04 1« 1 2 . 5 15 1 7 . 5 20 T=200 - 0.01 ;-0 . 02 -0 . 03 -0 . 04 1 7 . 5 1 4 1 2 . 5 15 20 j V,(t)-V,(M) -0 .01 1 2 . 5 15 1 7 . 5 20 V.(D-V,(M) o .o i 0 -0 .01 -0 .02 (-0 .03 -0.04 I T=620 2 . 5 1 2 . 5 1 5 1 7 . 5 20 Figure 10.S: Sequence of axial vel. changes (b), Case parameters Re = 200. a = 1.0, V = 1.095. 137 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3D, V=1.095 T=700 - 0 . 1 2.5 lO 15 17 . 5 12 . 5 20 T=9oa i-0 . 1 17 . 5 2.5 15 12 . 5 T=1000 - 0 .1 17 . 5 2.5 1« 12 . 5 15 20 T=1500 - 0 . 1 2.5 15 17 .5 12 . 5 20 Figure 10.9: Sequence of axial vel. changes (c). Case parameters Re = 200, a = 1.0, V - 1.095. 138 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. s s Figure 10.10: Power spectra density of the axial velocity at ; = 5. It’s sampled by 21 8 samples over 550 < T < 1000. Figure 10.11: Power spectra density of the axial velocity at : = 5. It’s sampled by 2IS samples over 1500 < T < 2000. fe* < K l« M t i W A M K M M Figure 10.12: Power spectra density of the axial velocity at s = 10. It’s sampled by 21 8 samples over 550 < T < 1000. Figure 10.13: Power spectra density of the axial velocity at z = 10. It’s sampled by 21 8 samples over 1500 < T < 2000. 139 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. To reveal more of the physical meanings of the streaklines shown above, several flow physical parameters are examined here. The lower part of Figure 10.14 shows vortex lines generated near inlet center at T = 1000. The similarity between the streaklines and vortex lines is quite streaking. Later we’ll see that particles are also concentrated in the larger vorticity magnitude regions(Figure 10.16). Also, we checked the pressure contour in the flow field for the plane of 9 = 0° and 9 = ISO0 and superimposed the streaklines into this field. It can be seen from Figure 10.15 that the particles basically concentrate in the region with lower pressure. By checking the vorticity m agnitude field as shown in Figure 10.16. the particles are found to be concentrated in the regions where the magnitude of vorticity is relatively large. To reveal more details of streaklines, we also solve the equation about the passive scalar $ and we choose the Peclet number as 1000. while Reynolds number is 200. Figure 10.17 to Figure 10.20 show the comparisons between the streaklines and the similar concentration contours at certain level. From the above observations, axisymmetric breakdown is an im portant step in the evolution of the three-dimensional finite amplitude disturbances and the unsteady 3D vortical flows. The same idea is expressed in the recent theoretical works by Rusak, Wang and W hiting (1998), though so far most of their works are based on the axisymmetry assumption. 140 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 10.14: Upper figure: Streaklines at T = 1000, lower part: Vortex lines at T = 1000. Case parameters Re = 200, a = 1.0, V = 1.095. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 00)I i iwii 1-o.ot ao> a t? o .m a s s a u a sa a <2 o.7t o j o a te o a t toy i t s 1 .2 s % y t 1 * 3 > 5 3 rat 1 7 0 T=1000, pressure contour of 6=0° and 0=180f l T=1000, streaklines Figure 10.15: Upper figure: pressure contour of plane 0 = 0° and 180° at T — 1000, lower part: streaklines at T = 1000. Case parameters Re = 200, q = 1.0, V = 1.095. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. flO l 1 I M lH W I laato >190 t.40o t eio tjao a.oao m » 2.449 z u a a w 1079 3 .2 0 3.* n i t o s ia ia T=1000, vorticity m agnitude contour (6=0° and 180°) o o n >M IW >1 T=1000, streaklines x Figure 10.16: Upper figure: vorticity magnitude contour of plane 0 = 0° and 180° at T = 1000, lower part: streaklines at T = 1000. Case parameters Re = 200, a = 1.0, V = 1.095. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. r. c=0.909 Figure 10.17: Passive scalar concentration at T = 50 T k . X Figure 10.18: Passive scalar concentration at T = 300 144 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 10.19: Passive scalar concentration at T = 700 5 Streakings Figure 10.20: Passive scalar concentration at T = 1500 145 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10.2 3D case of R e = ‘ 200, a = 1.0. V’ = 1.2 Some experiments (Sarpkaya (1974), Escudier (1988). etc.) have also observed some other forms of vortex breakdowns like double helix, etc., depending on the swirl strength of the inlet flows. At this 3D case, we increased the swirl strength V to a larger value of V '" = 1.2. Figure 10.21 to Figure 10.23 give a tim e sequence of streaklines of the flow field. Clearly, at T = 350. the second expansion region after the bubble breakdown is larger and longer than in case V ’ = 1.095. Until T = 475. the simple spiral breakdown is developed with a relatively small bubble near inlet. After T = 475, the double helix breakdown occurs and at T = 500. this kind of breakdown form is quite obvious. The mechanism for the occurrence of this kind of breakdown is not quite clear yet and still under investigation. Based on the observations of numerical simulation, the only comment we can make so far is that the wake-like region after the bubble breakdown is the crucial region where the complicated three- dimensional disturbances will be generated and develop into large finite amplitude movements. 146 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 10.21: Time history of streaklines, T=5, 100 and 350. Case parameters Re = 200, a = 1.0, V = 1.2. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3 Figure 10.22: Time history of streaklines, T=400, 450 and 475. Case parameters Re = 200, a = 1.0, V — 1.2. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. lE T iW .tm Z li Figure 10.23: Time history of streaklines, T=500, 550 and 600. Case parameters Re = 200, a = 1.0, V = 1.2. 149 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 11 Simulations of Jet-like Profiles Previously the simulations are based on the previous C'o-flow studies; the jet-like profiles have not been used yet. To mimic experiments of Liang k Maxworthy (1999, private communication), here we pick up a new set of initial profiles for azimuthal and axial velocities as followings: V*(r) = V > ( l - 0 . 5 ( e r / ( ^ ^ ) + l)), H(r) = (1 -0 .5 (er/(^ -y ^ ) + 1)). where function erf(r) is the error function, the two parameters H and S control the swirling strength and shear thickness respectively. A typical initial profile for both velocity and vorticity is shown in Figure 11.1 and Figure 11.2. As usual, the parameter V measures the swirling strength. In this chapter, vortex breakdown is examined again for both axisymmetric and 3D cases. 150 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Th« initial v*toc«y profiM at mtat 0 8 02 0.4 V e 0.7 0.8 0 9 OS 0.2 0 3 0 4 V 2 Figure 11.1: The initial velocity profiles at inlet The initial wioctty profiltf at inlat 4 3 1 O '— -2 .5 2 2.5 0.5 1.5 0 1 •2 -1 .5 -0 .5 - 1 O h 4 1 i i[ | j r 0 0 —-------e------1 5 1s :> z o>e Figure 11.2: The initial vorticity profiles at inlet Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 11.1 A xisym m etric cases o f Re = 200 Three cases of Re = 200 are investigated here. Figure 11.3 shows the streamfunc- tion contour in part of computational domain for medium swirling rate V = 1.2. Basically, it’s a steady solution. Increasing the swirling strength to V ' = 1.4, the solution reaches the steady state very slowly, even we did a long time calculation (about 20.000 steps compared with only a few hundred steps for the C'o-flow jets :n previous studies) (Figure 11.4). Figure 11.5 shows the velocity convergency time history. At V = 1.5, a quasi-steady vortex breakdown solution is obtained. Figure 11.6 and 11.7 show the results. R e=200. V=1.2 0.5 z Figure 11.3: The streamfunction contour for case Re = 200, V = 1.2 152 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Re=200.V=1.4 H 0.5 Z Figure 11.4: The streamfunction contour for case Re = 200. V = 1.4 C on vergen ce history ol velocity, Delta=0.2. V=1.4 * 0.6 ■ o 0.4 0.2 0.2 0.4 0.6 0.8 T 4 X 10 Figure 11.5: The tim e history of velocity convergence, case Re = 200, V =1.4 R e= 200.V = 1.5 0.5 Z Figure 11.6: The streamfunction contour for case Re = 200. V = 1.5 153 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Convergence history ot velocity. Delta=0.2. V=1.5 0.8 ^ 0.4 02 0.2 0.4 0.6 1.2 T 4 xIO Figure 11.7: The time history of velocity convergence, case Re = 200, V — 1.5 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 11.2 A xisym m etric Solutions o f H igh Re Recalling previous studies about C ’ o-flow jets profiles that high Re solutions exhibit either multiple bubble breakdowns or periodic behaviors, the solutions for current jet-like profiles are worth being examined again. For cases of Re = 500. Figure 11.9 and 11.8 show the solutions for a medium swirling strength V ’ = 1.2. Basically the solution doesn’t have the vortex break down, and it doesn't reach steady state in this case. Increasing to swirling strength V = 1.5,the quasi-steady status doesn't exist anymore (Figure 11.10). After long time calculation (0 < T < 20,000), the flow field still has the conical structure and extents to the upper boundary (Figure 11. 11). For higher Re number like Re = 1000. the same param eter of I ’ = 1.5 is carried on again for long time calculation. The similar conical structure is observed also (Figure 11.13 and Figure 11.12). Velocity con vergen ce history, D elta= 0.2.V =1,2,B e= 500 0.8 0.4 1.2 0.2 0.4 0.6 0.8 1.4 1.6 T x 10* Figure 11.8: The time history of velocity convergence. Case Re = 500, V = 1.2 155 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Re=500. V=1.2 1.5 0.5 Z Figure 11.9: The Streamfunction contour. Case Re = 500. V = 1.2 Velocity co n v erg en ce history. D elta= 0.2.V =1.5,R e=500 0 0 .2 0.4 0 .6 0 .8 1 1.2 1.4 1 .6 1.8 2 T x 10* Figure 11.10: The tim e history of velocity convergence, case Re = 500. V = 1.5 156 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Re=500, V=1.S O C 4 Z Figure 11.11: The Streamfunction contour. Case Re = 500. V = 1.5 Velocity c o n v erg en ce history. D eH as0.2,V 3l.5,R e=lO O O 1 - q! -------------- t L_ _ 1 _ 1 _ . L i t I 0 0 ^ 0 .4 0 .6 0 .8 1 1.2 1.4 1.6 1.8 2 Figure 11.12: The time history of velocity convergence, case Re = 1000, V = 1.5 157 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 11.13: The Streamfunction contour. Case Re = 1000, V = 1.5 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 11.3 3D simulation of case R e = 200, V = 1.5 We carry out the simulations to 3D to verify some of the observations of vortex breakdown about the Co-flow jets profiles. So far, the case of Re — 200. V = 1.5 is focused at this point. Figure 11.14 shows a comparison between 2D and 3D about d i\ which is the the maximum velocity change of flow field in temporal domain, it's clear from figure that 3D case basically repeats the 2D process and develops the bubble type breakdown first as happened in previous modeling profile case. After T > 1400. the three dimensional disturbances come up and bring flow field to a new structure. Unfortunately, the code cannot handle the strong instabilities for this kind of flow, the simulation couldn't be carried on beyond T = 2100. To explain this blowing up of the code, we did more tests from reducing time step to decreasing spatial step into the half or even quarter of the original ones, still the code blows up around T = 2100. There could be one possible reasons for this, which is that the effects from boundaries, esp. the upper one is not negligible any more. Recalling from the calculations before T = 200(Figure 11.15). the flow experiences strong sudden expansions towards the upper boundary where we assume the free flow condition there, once the 3D disturbances occur, the flow quite possibly might experiences the strong expansions again and bring down the strong boundary effects and blows up the code. The above result is basically consistent with our previous modeling profile stud ies in the sense that 3D instabilities comes out from 2D quasi-steady flow structure. However, instead of just forming a bubble breakdown, the current 3D flow expe riences a conical-like structure first, then shrinks to a much smaller bubble type breakdown. Figure 11.15 shows the time sequence of the streaklines for the 3D simulations and demonstrates the above observations. 159 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 0.1 0 .0 9 0 .0 8 0 .0 7 ’ - G 0 .0 6 0 .0 5 0 .0 4 0 .0 3 0.02 0.01 -0.01 2000 1000 1 5 0 0 5 0 0 T Figure 11.14: Temporal fluctuation of the maximum ve locity change. Case Re = 200, V = 1.5 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Figure 11.15: 3D simulation. Time sequence of streak- lines, T=30, 75, 200 and 800. Case Re = 200, V = 1.5 161 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 12 Conclusion • Based on Verzicco Orlandi’s (1996) scheme, developed a high resolution, time-dependent 3D N-S solver with parallel computing technology. The nu merical code is validated by comprehensively checking and comparing pro cedures. The results demonstrate that this code maintains good resolution on and near the symmetry line, where equations are singular mathematically and some traditional finite difference methods degrade. The present code has some obvious advantages over regular primitive variables, makes no need for pressure boundary conditions which are hard to be given correctly. • W ith current relatively simple model, some features of vortex breakdown, e.g., bubble and spiral breakdowns, have been recognized in the present sim ulations. The simulated features are consistent with some of the experimental observations and theoretical works. • Some breakdown criteria are investigated by present simulations. A simple vortex breakdown criterion proposed by Mahesh (Mahesh 1996) is examined and good agreement is found between theoretical prediction and the asymp totic numerical value at high Reynolds number. • A mechanism for axisymmetric vortex breakdown is confirmed to explain the process of bubble breakdown phenomenon. From our simulations of axisymmetric and full 3D cases, the bubble breakdown is found to be the result of disturbances propagating upstream inside the vortex core. Also. 162 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. the local instability analysis about the axisymmetric solutions of Re = 200 reveal that there exists a finite local AI section at the wake region of bubble breakdown, which will bring up the three-dimensional disturbances later and causes the occurrence of the spiral breakdown in full 3D simulations. • In 3D simulations, above a critical swirl number, the axisymmetric bubble- type solution becomes unstable, and the transition to a temporally periodic, helic vortex breakdown is observed. This, together with the finite absolute instability region observed behind the bubble breakdown, indicates the pos sible existence of a global mode (Hurre Sc Monkewitz (1990)). • From our simulations, axisymmetric breakdown is an important step in the evolution of the three-dimensional finite amplitude disturbances and the un steady 3D vortical flows. The same idea is expressed in their recent theoret ical works by Rusak, Wang and Whiting (1998). • The periodic axisymmetric solutions are observed for higher Reynolds num ber cases. A new view of the dynamic systems proposed by Brown Sc Lopez (1990) is used to explain the topological changes to try to explain the nonlin ear behaviors that happened at high Reynolds number. For the steady cases of high Re such as 500 or 1000, the effects of swirl rate can be explained as the changes in the flow field where more dynamical points, e.g. saddle or center point, are shown; while for the periodic solutions the change of instantaneous streamfunction pattern is characterized by the coalescing and separating of some streamlines inside the rings and also the oscillating of the mean position of the rings and bubble. • A set of jet-like profiles which are closer to experiments are studied. And the results are basically consistent with the studies from modeling profiles, esp., the 3D simulations will repeat the 2D results first, then come up the 3D instabilities. 163 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 12.1 Future Works • More simulations on the jet-like profiles more close to experimental studies. The jet-like profiles we used in Chapter 11 still need to be verified with experiments. • It would be more interesting and accurate that more efforts are put on the instability analysis about the axisvmmetrical solutions. So far, we do local analysis on the velocity profiles at every downstream station and assume they’re Batchelor's vortex profile-like. It should be better to studies the propagation of disturbance in the whole field, instead of at local station, then we can verify it with real 3D simulations. 164 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reference List [1] Althaus. V V . k Krause. E.. "Flow visualization of flow with concentrated vor- ticity,” EC Contract SCl-0212. Progress Report Dec. 1990. [2] Althaus. Y V . k YY'eimer. M., "Review of the Aachen work on vortex breakdown." IUTAM Symposium on Dynamics of Slender Vortices. 1998,331-344. [3] Akselvoll. K., and \Ioin. P.. "Large Eddy Simulation of Turbulent Confined Coannular Jets and Turbulent Flow over a Backward Facing Step." Thermo sciences Div.. Dept, of Mechanical Engineering, Report No. TF-63. Stanford Univ.. CA. Feb. 1995. [4] Beam. R. M.. and YVarming. R. F.. J. Comput. Phys. 22. 1976. p.87. [5] Benjamin T. B. Theory of the vortex breakdown phenomenon. J. Fluid Mech.. 1962, 14:593-629. [6] Benjamin T. B. Some developments in the theory of vortex breakdown. J. Fluid Mech., 1967, 28: 65-84. [7] Beran, P. S. k Culick. F. E. C.. “The role of non-uniqueness in the development of vortex breakdown in tubes.’ ’ J. Fluid Mech. 242. 1992. 491-527. [S] Billant, P., Chomaz, J. M. k Huerre, P., “Experimental study of vortex break down in swirling jets,” J. Fluid Mech. 11, 1997. [9] Bossel H. H.. "Vortex breakdown flowfield,” Phys. Fluids 12, 1969. 498-508. 165 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [10] Bruecker, CH. Sc Althaus. W. ” Study of vortex breakdown by particle tracking velocimetry(PTV), Part I: Bubble type vortex breakdown modes,” Exps. Fluids 13, 1992, 339-349. [11] Bruecker, CH. Sc Althaus, V V . ”Study of vortex breakdown by particle track ing velocimetry(PTV), Part 3: Time-dependent structure and development of breakdown modes,” Exps. Fluids 18. 1995, 174-186. [12] Brown, G. L. Sc Lopez, J. M., "Axisymmetric vortex breakdown. Part 2. Physical mechanisms,” J. Fluid Mech. 221, 1990, 553-576. [13] Delbende. I., Chomaz, J.-M ., Sc Huerre, P. "Absolute/convective instabili ties in the Batchelor vortex: a numerical study of the linear inpulse response,” Submitted to JFM , Dec. 1996. [14] Escudier. M. P., Bornstein, J. Sc Maxworthy. T., "The dynamics of confined vortices,” Proc. Roy. Soc. Lond. A 382, 335-360. [15] Escudier, M. P. Sc Zehnder, N., “Vortex-flow regimes.” J. Fluid Mech. 115, 1982, 105-121. [16] Escudier, M. P., “Vortex breakdown: observations and explanations,” Prog. Aerospace Sci. 25, 1988, 189-229. [17] Faler, J. H. Sc Leibovich. S., “Disrupted states of vortex flow and vortex breakdown,” Phys. Fluids 20, 1977, 1385-1400. [18] Gaxtshore I. S., “Recent work in swirling incompressible flow,” NRC Can. Aero. Rep. LR-343. 1962. [19] Gatski, T. B. Sc Spall, R. E. “Numerical studies of vortex breakdown: from helices to bubbles,” Fourth Int. Symposium on Comp. Fluid Dynamics, 1991, vol. I, 418-423. [20] Gelfgat, A. Yu., Bar-Yoseph, P. Z. Sc Solan, A., “Stability of confined swirling flow with and without vortex breakdown,” J. Fluid Mech. 311. 1996. 1-36. 166 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [21] Grabowski. V V . J.. and Berger. S. A.. “Solutions of the Navier-Stokes Equa tions for Vortex Breakdown,” J. Fluid Mech., Vol. 75, Part 3, 1976, p. 525. [22] Hashimoto, H. A, usoliton on a vortex filament.” J. Fluid Mech.. 1972, 51:477- 485. [23] Hall, M. G., Vortex breakdown. Ann. Rev. Fluid Mech. 4, 1972. 195-217. [24] Hall M. G., " A new approach to vortex breakdown,” Proc., Heat Transfer and Fluid Mech. Institute, 319-340, 1967. Stanford Univ., Press. [25] Harvey J. I\.. "Some observations of the vortex breakdown phenomenon,” J. Fluid Mech. 14, 1962. 195-218. [26] Huerre P., Monkewitz P. A.. "Local and global instabilities in spatially devel oping flows.” Annu. Rev. Fluid Mech. 22. 1990. 473-537. [27] Hopfinger, E. J., Browand, F. K. and Gagne. Y., "Turbulence and waves in a rotating Tank,” J. Fluid Mech. 125. 1982. 505-534. [28] Kida, S. A. uvortex filament moving without change of form,” J. Fluid Mech.. 1981. 112: 397-409. [29] Kim, J., and Moin, P., ^Application of a Fractional-Step Method to Incom pressible Navier-Stokes Equations.” J. C ’ om ptut. Phys. 59, 1985, p. 308. [30] Kopecky, R. M., and Torrance, I\. E., "Initiation and Structure of Axisym metric Eddies in a Rotating Stream.” Computer & : Fluid, Vol. 1, 1973, p.289. [31] Kribus, A. & c Leibovich, S., ^Instability of strongly nonlinear waves in vortex flows,” J. Fluid Mech. 269, 247-264. [32] Krause, E. Shi, X. & Hartwich, P. M., “Computation of leading edge vortices,” AIAA Paper no. 83-1907, 1983 [33] Lambourne. N. C. & : Bryer, D. VV.. “The bursting of leading-edge vortices- some observations and discussion of the phenomenon." Aeronautical Research Council R M 3282. 1-36. 167 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [34] Lambourne. N. C. . uThe breakdown of certain type of vortex.” Aeronautical Research Council(G. B.), Rep. No. ARC-CP-915: NPL-AERO-116. [35] Leibovich. S., uVVeakly non-linear waves in rotating fluids," J. Fluid Mech.. 1970. 42:803-822. [36] Leibovich. S. Randall. J. D. A amplification and decay of long nonlinear waves. J. Fluid Mech., 1973, 5S:481. [37] Leibovich, S., "The Structure of Vortex Breakdown." Ann. rev. Fluid Mech.. Vol. 10. 1978, p. 221. [38] Leibovich, S. Ma. H. Y.. "Soliton propagation on vortex cores and Hasimoto soliton,” Phys. Fluids. 1983, 26: 3173-3179. [39] Leibovich. S. & Kribus. A.. "Large amplitude wavetrains and solitary waves in vortices.” J. Fluid Mech. 216. 1990. 459-504. [40] Lessen. H.. Singh. P. J. & Paillet. F. "The stability of a trailing line vortex. Part 1: Inviscid theory.” J. Fluid Mech. 63. 1974. 753-763. [41] Lim. D. V V . Redekopp, L. G.. "Absolute Instability Conditions for Variable Density. Swirling Jet Flows.”. 1997 (??) [42] Loiseleux T.. Chomaz, J. M. Huerre P.. "The effect of swirl on jets and wakes: Linear instability of the Rankine vortex with axial flow," Phys. Fluids 10, 1998, 1120-1134. [43] Lopez, J. M.. "Axisymmetric Vortex Breakdown Part 1. Confined Swirling Flow.” J. Fluid Mech, vol 221. 1990, p.533. [44] Lopez, J. M.. uOn the bifurcation structure of axisymmetric vortex breakdown in a constricted pipe,” Phys. Fluid 6, 1994, 3683-3693. [45] Lucas, K., “Numerical investigation of three-dimensional vortex breakdown." Ph.D thesis. Stanford university. [46] Ludwieg, H.. "vortex breakdown." Deutsche Luft-und Raumfahrt. Rep. 70-40. 168 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [47] Lundgren, T. S. & Ashurst. V V . T.. “Area vary waves on curved vortex tubes with application to vortex breakdown,'’ J. Fluid Mech. 200, 1989, 283-307. [48] Mahesh, K., "A model for the onset of breakdown in an axisymmetric com pressible vortex,'1 Phys. Fluids 8. 1996, 3338-3345. [49] Mager, A., “Dissipation and Breakdown of a Wing-Tip Vortex," J. Fluid Mech., Vol. 55, Part 4, 1972, p. 609. [50] Marshall. J. S., “A general theory of curved vortices with circular cross-section and variable core area,’ " J. Fluid Mech. 229, 1991, 311-338. [51] M artin, J.E., and Meiburg, E., "On the Stability of the Swirling Jet Shear Layer,'1 Physics of Fluid, vol. 6, No. 1. 1994. [52] Maxworthy T., “on the structure of concentrated, columnar vortices.'1 1971, Astronautica Acta., Vol. 17. 363-374. [53] Maxworthy T.. Mory M. Hopfmger E. J.. "Waves on vortex cores and their relation to vortex breakdown.” 1983, Proc. AGARD Conf. on Aerodynamic of Vertical Type Flows in Three Dimensions, AGARD C'PP-342, paper 29. [54] Maxworthy, T.. Hopfmger, E. J. and Redekopp. L. G. Wave motions on vortex cores. J. Fluid Mech.. 1985, 151: 141. [55] Maxworthy T., "Waves on vortex cores,” 1988, FI. Dyn. Res. 3, 52 [56] Pritchard, W. G.. "Solitary waves in rotating fluids.” J. Fluid Mech., 42, 1970,61-83. [57] Orlanski, I., “A Simple Boundary Condition for Unbounded Hyperbolic Flows,” Journal of computational Physics 21, 1976, p.251. [58] Panda, J.. and McLaughlin, D. K., “Experiments on the Instabilities of a Swirling Jet,” Physics of Fluid, vol. 6. No. 1, 1994. [59] Rai. M. M.. and Moin. P.. ‘Direct simulations of Turbulent Flow Using Finite- Difference Schemes." J. C'omptut. Phys. 96. No. I. 1991. p. 15. 169 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [60] Randall J. D. and Leibovich S.. “The critical state: a trapped wave model of vortex breakdown,” J. Fluid Mech. 58, 1973, 495-515. [61] Rusak, Z. “Axisymmetric swirling flow around a vortex breakdown point,” J. Fluid Mech. 323,1996, 79-105. [62] Salvetti, M. V.. Orlandi P.. and Verzicco, R.. “ Numerical Simulations of Transitional Axisymmetric Coaxial Jets,” AIAA Journal. Vol. 34. No. 4, April 1996. [63] Sarpkaya, T., “On Stationary and Traveling Vortex Breakdowns,” J. Fluid Mech., Vol. 45. Part 3. 1971, p. 545. [64] Sarpkaya, T.. “Effect of the adverse pressure gradient on vortex breakdown," AIAA J., 12. 1974, 602-607. [65] Shtern, V.. Borissov, A. & Hussain F.. "Vortex sinks with axial flow: solution and applications,” Phys. Fluids 9. 1997, 2941-2959. [66] Sliyy, V V .. Thkur, S.. and Wright. J.. “Second-Order Upwind and Central Difference Scheme for Recirculating Flow computation,” AIAA J., Vol. 30. No. 4. 1992, p. 923. [67] Spall, R. E., Gatski, T. B. & Grosch, C. E. “A criterion for vortex breakdown.” Phys. Fluids 30. 1987, 3434-3440. [68] Spall, R. E., Gatski. T. B. & : Ash, R. L. “The structure and dynamics of bubble-type vortex breakdown.” Proc. R. Soc. Lond., 1990, A429. 613-637. [69] Squire, H.B., Imperial College, London. Dep. of Aero. Rep. no. 102, 1960. [70] Tromp J. C., “The dependence of the time-asymptotic structure of 3D vortex breakdown on boundary and initial conditions,” Ph.D thesis, Air Institute of Technology, Department of the Air Force, VVright-Patterson Air Force Base, Ohio. July 1995 170 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. [71] Tromp J. C.. & Beran P. S.. uThe role of nonunique axisymmetric solutions in 3D vortex breakdown,” Phys. Fluids Vol. 9, No. 4, April 1997, 992-1002. [72] Tsai, C.-Y. < V VVidnall, S. E. uExamination of a group-velocity criterion for breakdown of vortex flow in a divergent duct,” Phys. Fluids 23, 19S0. 864-870. [73] Verzicco. R., and Orlandi, P.. "A Finite-Difference Scheme for Three- Dimensional Incompressible Flows in Cylindrical Coordinates,” J. C'omptut. Phys. 123, 1996, p. 402. [74] Wang, S. & : Rusak, Z., “The effect of slight viscosity on a near-critical swirling flow in a pipe,” Phys. fluids. Vol. 9. No. 7. 1997, 1914-1927. [75] Wang, S. & Rusak, Z.. "The dynamics of a swirling flow in a pipe and tran sition to axisymmetric vortex breakdown,” J. Fluid Mech. 340. 1997. 177-223. [76] Wang, S., Rusak, Z. W hiting C. H.. "The evolution of a perturbed vortex in a pipe to axisymmetric vortex breakdown,” J. Fluid Mech. 366. 1998, 211-237. 171 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Chapter 13 Appendices 13.1 Appendix A: Spatial Finite Differencing The centeral difference scheme based on staggered grid applied to the momentum Equs. (lb) and (lc). The computational region 0 < z < L. 0 < y < 0.5. 0 < 0 < 2- in which the solution is to be obtained is overlaid with a grid system of N3 points in the z direction. N2 point in the y direction and Nl points in the 9 direction. Notice that only r direction is mapped from physical range 0 < r < F t to computational range 0 < y < 0.5. The location (O.y.z) of the grid point in the center of computational cell is given by the indices (ic+1/2, jc + 1/ 2. kc+ 1/ 2) (see Figure 13.1). qg(y,9,z), qr(y,9,z) and q:(y,9.z) are represented by the function qg(ic.jc.kc), q?(ic,jc,kc) and q”(ic,jc,kc) defined at each grid point (ic,jc + 1/2,kc + 1/2), (ic+ 1/2, jc, kc + 1/2) and (zc+ l / 2 ,j c + 1/2. A rc) which is located on the surface of the computational cell, here n means tim e level. p(y,9,z) is represented by the function pn(ic.jc, kc) defined at the center of the cell (ic.jc, kc). So the equation for qg is discretized around (ic.jc + 1/2, kc+ 1/2), qr equation around (ic -t-1/2,jc,kc+ 1/2) and q. around (ic + 1/2,jc + 1 /2, A rc). Centered differences are used for all spatial derivatives. For example, for a typical convective term like (1 /r2)j^(rqgqr) at work point (ic,jc+ 1 /2 ,kc + 1/2) in the qg momentum equ., it can be discretized as following: 1 d(rqeqr) , ~ I ] fr I.C .JC + 1 /2 .A .-C + L /2 - 172 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Iz A x i a l R a d ia l A z im u th a l Figure 13.1: A typical computational cell [( ^'^lO^lr ) |i'c,jc+l,/:c+l/'2 ( r<7 fl?r) |ic.jc,fcc+l/2 ’ IC.JC+1 / 2 ,kc+ 1 /2 ( .4. 1a ) where (?*) l.c.jc+u-c+i/2= 0.5 * [ < 7 e( ic .jc + 1, A rc) + qa(ic. jc, A rc)] (<7r) |,c j c + i .J tc + i / 2 = 0.5 * [qP(ic - 1, jc + 1. A rc ) + qr{ic.jc + I. A rc)] ^ ^ M |.c.Jc./:c + i/2 = 0.5 * [ < 7 fl(i'c, j c - 1, A rc) + r/fl(ic. jc. A rc)] ( 9 r ) |icjc.*c+L/2= 0.5 * [qr{ ic - 1, jc.Arc) + gr(/c, jc. A rc)] and F \ic,jc+i/2M+i/2= (ic'Jc+ 1/2, A rc + 1/2). (.4.1c) here, (dy/dr)(ic,jc + 1/2, Arc+ 1/2) is the transform function in r direction at point (ic.jc + 1/2, A rc + 1/2), 8y is the grid spacing in y direction. Similarly other spatial derivative terms can be obtained in this way. 13.1.1 Treatment o f the axis The discretization of the region around r=0 is the most im portant feature that is quite different from other schemes. As indicated by Verzicco (Verzicco & Orlandi 1996). the difficulties in its accurate representation are the main reason for the few 173 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. simulations in the literature. This has been extensively described in Verzicco k Orlandi’s paper and here we just summarize his treatm ents to show how to avoid the singularities in the Equs. (lb) to (lc). Because of qr = 0 at r = Oby definition the equ. for qr can be discretized for all j c > 2 without any difficulties. The equations for q: and qg at j c = | require the evaluations of the radial derivatives around r=0. For the convective terms the advantages of qr = 0 and r = 0 at jc = l avoid the evaluations of qg and q. at jc = i( r = 0), thus these convective terms can be discretized without any approximations. For the viscous radial derivative of Equ. qz (.4.2) the advantage of r t = 0 avoids the approximation of ^ | j C =i with a reduced accu racy at r= 0. The viscous derivative for Equ. qg is (.4.3) Again, because of r t = 0 no approximations have been done and the derivative maintains its second-order accuracv at the axis. 174 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 13.2 Appendix B: Temporal Finite Difference A fractional-step method and a low storage third-order Runge-Kutta scheme were employed for the time-advancement of Equs. (lb) and (lc). This was extensively described in Verzicco & Orlandi’ paper (Verzicco & Orlandi 1996), here just sum marize his approaches. For the details of the followings, please look at their paper. The nonsolenoidal velocity field q can be evaluated as: = [l t H‘ + - a,G iPl + a,(A> g + A ir + A ,-.-)-y^ ]. (0.1) Where Hi contains the convective terms and those viscous terms with a single velocity derivative or without derivative. 6 ’ ,, .4,0, .4r and .4,, denote discrete differential relations for the pressure and viscous terms respectively, whose explicit expressions are the same as in the spatial finite diff. in Appendix A. For the explanation of i, / and q, please look at their paper. a;, 7/ and pt are the coefficients of the tim e advancement scheme and the same values as in the paper (Rai,1991) which are 7i =8/15,72 = 5/12,73=3/4. Pi = 0,p2 = —17/60, p3 = -5 /1 2 . (0.2) c*i = 7i + pi = S/15,ct2 = 72 4-/92 = 2/15,03 = 73 4- P 3 = 1/3. The implicit treatm ent of the viscous terms eliminates the numerical viscous stability restriction which is particularly severe for low Reynolds number flows and near boundaries where stretched meshes are used. However, this also requires the inversion of large sparse matrices to obtain solutions. These can be reduced to three tridiagonal matrices by employing a factorization procedure with error 0(8t3) (Beam & Warming 1976). By using the increment Sq, = <7, — q\ Equ. (B .l) can be written in delta form, giving (1 - /M * )(l - ^-4,,)(1 - &(Air)Sqi = <^[7/H[ + piH[ 1 — o/G’ ,y 4- o/(.4,0 4- -4,r 4- .4 ,-;)< ■ /{ ]. {BA) 175 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Where 3i — aiSt/2. After obtaining the interm ediate velocity field q, the continuity equation is then enforced on each cell by introducing a scalar variable $ to map it to the solenoidal velocity field q. qj+l = -a ,S tG i* l+l, (BAa) where r 1 d . d _ 3 0 , = 7 W ’ - G' = r d~r‘r' $ is calculated from (see Appendix C ' in detail) L ^ +l = - ^ ;D i- q , (BAb ) Q/Ol where D, and L = D, • ( S ', are the following discrete differential relations „ i a , 1 5 . 5 _ Di = r ov r or oz I d 2 I d d d2 7 * !W + ~ r ’ d r J" d^ + 5?' ( ^ Then the pressure can be calculated as pl+l = p‘ + $ '+ l - (ctiSt/ (2Re))L4>l+l. (BAc ) 176 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 13.3 Appendix C: Solving Poisson Equ. about $ using FFT To solve the Poisson Equ. (B.4b) about pressure-related term $ L$l+l = -KrDi ■q = Q(ic,jc. kc) (C .l) Q;0£ where 1 d2 i a d d2 r2 d&2 + r d r d r + dz2' A solution to the above Equ. and the corresponding boundary conditions can be easily obtained using Fourier transform methods (Kim & Moini985). Let * V i — 1 iV 3 — l -r.m-ikc—-) - 2 ti.c ${icjc,kc) = 11, ${lJc,m)cos[- - H e - ' ,V l (C*.2) /=0 m=0 * ^ for ic = 1,2. • • •. ;Vt, jc = 1.2, • • •, N 2, kc = 1,2, • • •, A’ 3. Substituting into Equ. (C .l) and using the second order central difference on staggered grid for all the derivative terms, we obtain 1 d 9 $ * » ~rFr{ r f r ] ~ = QiLjC-m )’ (C,3) where k'{ = 2[1 — cos(2tt • I/ . \ \ )] and k'm = 2[1 — cos(- • m/A'’ 3)] are the modified wave numbers in Fourier frequency domain. For each set of the wave numbers, the above tridiagonal system of equations can be easily inverted, and $ is obtained from Equ. (C.2). Note that $ is determined to within an additive constant. So after getting the solution, the constant can be subtracted from the whole field. To determine the Fourier transform of $ (/,jc , m) using the standard periodic FFT packages instead of constructing a special routine for cosFFT transform, the following auxiliary function $ ’(ic,jc,kc) can be assumed as . 11 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. { complex(^(ic, jc, A rc), 0), 1 < kc < N$ — 1 (C.3) com pie j:(0. 0), A rc > ;V 3 The inverse Fourier transform is /V, 2-<V 3 : 2re-kc m ; 2ir tc I £ V (icjc,k c)e' 2s* e vi (C-4) ic = l fcc=l for I = 0,1, • • •, iV | — 1, m = 0,1, • • -,2 • 'V 3 — 1. According to the definition of $'(ic.jc. kc), Equ. (C.4) becomes .v , .v , $ ’ (/. jc.m) = < & ;ea/ + = j r 5^ $(ic,ic.A-c)e,^ _ e,_^ - (C.5) ic=L fcc=I where .V, ;V 3 , i M v— ' . . . . T T ' k c - m , i 2~jc 1 $ re a ( = 2 - 2 - $ { l C. JC, k c ) c O s [ — ]e ^'ma3 = E E ${ic.jc, kc)sin[‘' fe m ]e,l77r- \' ic= l kc=L * 3 so easily we can have * ■ m ■ 7;, - , .it • m • 4, ^ rea/co5[— = j r $[ic,jc. A . ,c)cos['r m — —]e* vi (C.6) ic = l fcc=l i 3 The right part of Equ. (C.6) is exactly the inverse Fourier transform of Equ. (C.2), So $ (/,jc , m) can be calculated as ${l,jc, m) = $ ;eaic° 4 " - ^ " ■ ■ ] + ^'ma3sin[ - - ^ — -] {C.7) - ’ 3 ^ 3 Which < & ’ can be obtain easily through Equ. (C.4) using standard periodic FFT. ITS Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. -4.2 -4.6 -5.4 -5.6 - 6 1 - - 2.2 - 1.7 - 1.6 -1.5 -1.9 Figure 13.2: Accuracy checking for Poisson Eqn. about To test the above approach, a simple axisymmetric test case was constructed. Assuming boundary conditions for Equ. (C .l) are 0$ n d < t > — — = 0. r = 0 and r = 1. or d$ — = 0, r = 0 and z = 1, (C.8.1) and Q(r, z) = z2( I - z)2 • 2(2 - 9r + Sr2) + r2(l - r)2 • 2(1 - 6c + 6c2) (C.S.2) Under these assumptions, Equ. (C.l) has a exact solution $ e(r. c) as 179 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. <M r,c) = r2( l - r ) * . r 2( l - r ) 2 (C.S.3) Solving the above Equs. using FFT and comparing the numerical solution with exact solution to get the maximum error, vvedl get the accuracy checking plotting(Figure 13.2). The slope of the line in Figure 13.2 is about 2 which is the accuracy for the second order central difference for the above equations. ISO Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 13.4 Appendix D: The Transport Equation for Passive Scalar dc i .1 d d c s i d2c d2c. ,oiN Dt ~ R e x Sc r dr dr } + r 2 dO2 + dzl ^ ( ] Where C is the passive scalar to be solved, and D C _ _ d £ IdqrC 1 dqeC_ dqzC Dt dt + r dr ^ r dO dz 181 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.

Asset Metadata

Creator Chen, Peiji (author)

Core Title Numerical simulations of swirling flows

Contributor Digitized by ProQuest (provenance)

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag engineering, aerospace,OAI-PMH Harvest

Language English

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c16-59730

Unique identifier UC11338335

Identifier 3017990.pdf (filename),usctheses-c16-59730 (legacy record id)

Legacy Identifier 3017990.pdf

Dmrecord 59730

Document Type Dissertation

Rights Chen, Peiji

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA