"15868 tlllllllllll llllllllllllitllllllfill 3 1293 01686 3270 This is to certify that the dissertation entitled Relaxation Model for Homogeneous Turbulent Flows presented by S teven M. Parks has been accepted towards fulfillment of the requirements for 1 Ph.D. Chemical Engineefing degree in duke? Major professor Date M 2' ) ‘qfi 7 MS U is an Affirmative Action/Equal Opportunity Institution 0— 12771 LIBRARY Michigan State University —v— '— rv PLACE IN REI‘URN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. I DATE DUE DATE DUE DATE DUE ISEPi2zaeomo “13291 12007 [int ighpfifi :,.I~J we WWW.“ RELAXATION MODEL FOR HOMOGENEOUS TURBULENT FLOWS by Steven M. Parks A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Chemical Engineering 1997 ABSTRACT RELAXATION MODEL FOR HOMOGENEOUS TURBULENT FLOWS by Steven M. Parks Turbulent closure for the Reynolds stress is investigated for homogeneous turbulent flows. The basic structure of the transport equations for the turbulent kinetic energy and turbulent dissipation is adopted from the standard k-e model of turbulence. An integral analysis of the Kman-Howanh relation for isotrOpic turbulence yields a prediction for the destruction-of-dissipation coefficient in the e-equation. This coefficient is applied to the problem of decaying isotropic turbulence and is shown to be able to quantitatively reproduce several data sets for the flow available in the literature. A Green’s function analysis of the equation for the fluctuating velocity and a subsequent smoothing approximation yields a turbulent pre-closure relating the Reynolds stress to mean field quantities and an unclosed quantity: the turbulent pre-stress, which is related to both pressure fluctuations and the fluctuating Reynolds stress. An isotropic pre- stress closure assumption is introduced and applied to the problem of homogeneous shear. This closure is found to guarantee realizability a priori and yields a non-zero primary normal stress difference. Subsequent extension of the closure to an objective, anisotropic pie-stress accounts for both stress relaxation effects as well as a non-zero second normal stress difference for homogeneous shear flows. The case of homogeneous shear flow in a rotating frame of reference is examined due to its qualitative similarity to a flow in an inertial frame which also includes swirl and/or streamline curvature. The predominant result is that the combination of shear and rotation qualitatively changes the nature of the flow. For intermediate relative rotation rates, both I: and e are growing without bound. For large absolute relative rotation rates, turbulent production is cut-ofl‘ and both I: and a decay to zero. to Maya iii ACKNOWLEDGEMENTS First and foremost, I would like to express my gratitude to Professor Charles Petty for his guidance and continuing support throughout this entire research. Additional thanks are due to Dr. Syed Ali and Professor John Foss for their helpful suggestions during the early period of the research. Finally, heartfelt appreciation is expressed to my graduate student colleagues, Maya Murshak, Klaus Weispfennig, and G. Dale Wesson, for their valuable discussions and constant encouragement. I would also like to acknowledge my families for their unending support and source of inspiration: Mom, Dad, Violetta, Isaiah, Eric, Mike, Isaak, and May, thanks to you all. Special thanks are given to F.E.C. Laboratories and Michael Goergen for the unique opportunities ofi‘ered to me, allowing me to pursue both my professional and academic goals. I would like to recognize my friends and colleagues at the Lehrsruhl fiir Stramungmechanik in Erlangen, Germany, for making my personal growth there as important as the academic. Great thanks go to Professor Cameron Tropea, Suad Jakirlic, Marco Marengo, and Jochen Volkert. The author would like to thank the Department of Chemical Engineering and the Hydrocyclone Development Consortium at Michigan State University for financial support during this research. Further thanks are extended to the DuPont corporation, the DeVlieg foundation, and the DAAD (Deutsche Akademische Austauchdienst; German Academic Exchange Service) for additional funding. iv TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES LIST OF NOTATION CHAPTER 1. INTRODUCTION 1.1 Background 1.2 Objectives 1.3 Methodology 2. THE DECAY OF ISOTROPIC TURBULENCE IN AN INERTIAL FRAME 2.1 Introduction 2.2 Local Analysis of the Kannan-Howarth Equation 2.3 2.4 2.5 2.6 2.7 Global Analysis of the Karman-Howarth Equation Approximate Analysis of Dissipation and Production Integrals Parameter Estimates Transient Isotropic Decay Conclusions 3. ALGEBRAIC PRE-CLOSURE THEORY FOR THE REYNOLDS STRESS 3.1 3.2 3.3 3.4 3.5 3.6 Introduction Pre-Closure Theory Isotropic Pre-Stress Theory Parameter Estimates Results and Discussion Conclusions 4. ANISOTROPIC PRE-STRESS THEORY FOR HOMOGENEOUSLY SHEARED FLOWS 4. l 4.2 4.3 4.4 Introduction AnisotrOpic Pre-Stress Theory Anisotropic Homogeneous Decay Homogeneously Sheared Turbulence Page V111 14 15 21‘ 21 29 33 35 42 50 61 70 76 81 83 92 94 94 103 4.5 Asymptotic Homogeneous Shear 4.6 Transition States for Homogeneously Sheared Turbulence 4.7 Conclusions HOMOGENEOUSLY SHEARED TURBULENCE IN A ROTATING FRAME OF REFERENCE 5.1 Introduction 5.2 Pre-Closure Theory 5.3 Coriolis Redistribution of Turbulent Energy for Homogeneous Flows 5.4 Isotropic Pie-Stress Theory for Rotating Homogeneous Shear 5.5 Anisotropic Pre-Stress Theory for Rotating Homogeneous Shear 5.6 Conclusions CONCLUSIONS RECOMMENDATIONS 7.1 Further Study 7.2 Further Application APPENDICES an encore n—r O Navier—Stokes Equation in a Rotating Frame of Reference Reynolds Stress Equation in a Rotating Frame of Reference Turbulent Dissipation Equation in a Rotating Frame of Reference Objectivity of Convected Trme Derivatives Tables of Referenced Data Relation of Homogeneous, Isotropic Turbulent Production to the Velocity Derivative Skewness Asymptotic State for Homogeneous Shear Computer Program Listings H.l Program TRANS . BAS H.2 Program IPS_ROT. BAS Properties of the Eigenvalues of the Reynolds Stress LIST OF REFERENCES vi 106 109 120 124 124 128 130 132 152 161 163 168 168 170 175 177 179 181 183 197 199 201 201 206 209 211 TABLE 2.1 3.1 3.2 4.1 El E2 E3 E4 E5 E6 E7 E.8 E.9 LIST OF TABLES Summary of Parameters Used in the Analysis of the “(H-Equation Parameter Estimates for the IPS-Theory Predictions of Asymptotic Statistical Properties for Homogeneously Sheared Turbulence Parameter Estimates for the APS-‘Iheory Summary of Isotropic Decay Data; Batchelor and Townsend [1948]; M=0.635cm (grid mesh) Summary of IsotrOpic Decay Data; Comte-Bellot and Corrsin [1971]; U0=10m/s Summary of Isotropic Decay Data; Sirivat and Warhaft [1983]; M=2.5cm (grid mesh) Summary of Isotropic Decay Data of Speziale er a1. [1987]; Rea = 35.1 Summary of Isotropic Decay Data Bardina er al. [1985]; Rea = 45.4 Homogeneous Shear Flow Data (Tavoularis and Kamik “A”, 1989) Homogeneous Shear Flow Data (Tavoularis and Kamik “D”, 1989) Homogeneous Shear Flow Data (Tavoularis and Kamik “G”, 1989) Homogeneous Shear Flow Data (Tavoulan’s and Kamik “TC”, 1989) 15.10 Homogeneous Shear Flow Data (Tavoularis and Karnik “L”, 1989) 13.11 Return to Isotropy Data: Choi and Lumley [1984] (Plane Distortion) E.12 Return to Isotropy Data: Choi and Lumley [1984] (Axisymmetric Expansion) E.13 Return to Isotropy Data: LePenven et al. [1985] (Positive Third Invariant) 84 85 110 184 185 186 187 188 189 190 191 192 193 194 195 196 LIST OF FIGURES FIGURE 1.1 L-Diagram for Realizable Anisotropic States 1.2 Anisotropic States for Homogeneous Shear Flows 1.3 Closure Methodology 2.1 Relaxation States for Decaying Isotropic Turbulence 2.2 Literature Data for the Skewness and Correlations for the Destruction-of-Dissipation Coefficient 2.3 Qualitative Behavior of the Integrands I 1, [2, and I 3; Estimation of the Cut-ofi‘ Integration Distance EC 2.4a. Prediction of the IKH-Equation for the Skewness Using a Prescribed Function for the Destruction of Dissipation 2.4b Parametric Study for the Effect of Re‘ on the Resulting Skewness Prediction 2.4c Parametric Study for the Effect of a on the Resulting Skewness Prediction 2.5 Predictions of the IKH-Equation for the Destruction of Dissipation Coefficient 2.6 Predictions of the IKH-Equation for the Production of Dissipation Using a Prescribed Function for the Destruction of Dissipation 2.7a Model Computations for the Isotropic Decay Data of Batchelor and Townsend [1948] (U = 150 cm/s, Re a = 7.58) 2.7b Model Computations for the Isotropic Decay Data of Batchelor and Townsend [1948] (U = 643 cm/s, Re a = 42.6) viii ll 20 24 28 37 43 47 48 49 51 52 2.7c Model Computations for the Isotropic Decay Data of Batchelor and Townsend [1948] (U = 1286 cm/s, Re a = 83.7) 53 2.7d Model Computations for the Isotropic Decay Data of Comte-Bellot and Corrsin [1971] (U = 10 m/s, Rea = 354) 54 2.7c Model Computations for the Isotropic Decay Data of Comte-Bellot and Corrsin [1971] (U = 10 m/s, Rea = 769) 55 2.7f Model Computations for the Isotropic Decay Data of Sirivat and Warhaft [1983] (U = 340 cm/s, Rea = 139) 56 2.7g Model Computations for the Isotropic Decay Data of Sirivat and Warhaft [1983] (U = 630 cm/s, Re a = 262) 57 2.7b Model Computations for the Isotropic Decay Simulations of Speziale, et al. [1987] (Re a = 35.1) 59 2.7i Model Computations for the Isotropic Decay Simulations of Bardina, et al. [1985] (Re 0 = 45.4) 60 2.8 Predicted Transient Decay of the Turbulent Kinetic Energy as a Function of the Initial Turbulent Reynolds Number 62' 3.1 Asymptotic State for Homogeneously Sheared Turbulent Flows 65 3.2 Energy Simplex with Transition States for Homogeneously Sheared Turbulence 67 3.3 Anisotropy Invariant Diagram with Transition States for Homogeneously Sheared Turbulence 69 3.4 Anisotr0py Invariants Predicted by the [PS-Theory for Homogeneously Sheared Turbulence 80 3.5 3.6 3.7 3.8 Estimation of the Relaxation Group for Asymptotic Homogeneously Sheared Turbulence 82 The Effect of the Development Time on the Relaxation Group for Homogeneously Sheared Turbulence 86 The Effect of the Development Time on the Distribution of Kinetic Energy for Homogeneously Sheared Turbulence 87 The Effect of the Development Time on the Shear Component of the Normalized Reynolds Stress for Homogeneously Sheared Turbulence 91 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 5.1 5.2 5.3 5.4 5.5 5.6 5.7 Relaxation to Isotropic, Homogeneous Decay 102 The Effect of the Development Time on the Relaxation Group for Homogeneously Sheared Turbulence 112 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7): Components of the Normalized Reynolds Stress 114 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7): Components of the Anisotropic Pre-Stress 1 15 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7): dHyz/dg 116 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (”It = 0.7): dez/dE 117 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7 ): Anisotropy Invariants 121 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7): Energy Simplex 122 Schematic for Homogeneously Sheared Turbulence in a Rotating Frame 125 Coriolis Redistribution of Turbulent Kinetic Energy for Homogeneous Flows 133 The Effect of Rotation on the Turbulent Relaxation Time for Asymptotic Homogeneous Shear (IPS-Theory) 137 The Effect of Rotation on the Components of the Asymptotic Reynolds Stress for Homogeneous Shear (IPS-Theory) 138 Distribution of the Energy Components for Homogeneously Sheared Turbulence 139 Transient Response of the Relaxation Group for Homogeneously Sheared Turbulence in a Rotating Frame (IPS-Theory) 141 Relative Growth Rates of Turbulent Statistics in the Asymptotic Regime (IFS-Theory; N; = 1.0, m5 = —o.5) 143 5.8 Relaxation of the Reynolds Stress Components in the Asymptotic Regime (IPS-Theory; N; = 10 9/5: -0. 5; R" = R (No 52/5), cfi Eqs. (5.10)-(5.14)) 144 5.9 Relative Growth Rates of Turbulent Statistics in the Decay Regime (IPS-Theory; N; = 1.0; 9/5 = + 1.0 and +2.0) 145 5.10 Relaxation 0fth Reynolds Stress Components in the Decay Regime (IPS-Theory; = 1.0; m5 = +1.0 and +2.0; 13" = 13m" n/S) , cf Eqs. (5.10)-(5.14)) 146 5.11 Anisotropy Invariants for the Asymptotic States of Homogeneously Sheared Turbulence in a Rotating Frame (IPS-Theory) 149 5.12 Anisotropy Invariants for the Decaying States of Homogeneously Sheared Turbulence in a Rotating Frame (IPS-Theory) 150 5.13 The Effect of Rotation on the Turbulent Relaxation Time (NR = ‘R‘C S) for Asymptotic Homogeneous Shear (APS-Theory) 154 5.14 The Efl'ect of Rotation on the Components of the Asymptotic Reynolds Stress for Homogeneous Shear (APS-Theory) 156 5.15 The Efl‘ect of Rotation on the Asymptotic Pre-Stress Anisotropy for Homogeneous Shear (APS-Theory) 157 5.16 Distribution of the Energy Components for Homogeneously Sheared Turbulence in a Rotating Frame (APS-Theory) 158 5.17 Asymptotic Anisotropy Invariants for Homogeneously Sheared Turbulence in a Rotating Frame (APSoTheory) 159 5.18 Transient Response of the Relaxation Group for Homogeneously Sheared Turbulence in a Rotating Frame (APS-Theory; N0- = 1.0; H”: 0) 160 7.1 Qualitative Sketch of a Hydrocyclone Centrifugal Separator 171 LIST OF NOTATION English Symbols II III w 0" "6" II)» I" 1“ 2w 2 a.“ we we (3. >9 an (I!) "N H: N y—n Second invariant of 9 Third invariant of _b Parameter in objective time derivative; parameter in theoretical correlation for CD (Re); empirical parameter for the estimation of d[ln (11)] /d[ln (Re) ] Matrix relating the Reynolds stress and the pre-stress in the pre-closure Anisotropy dyadic of the normalized Reynolds stress Empirical parameter for the estimation of 13/11 Double longitudinal velocity conelation Double normal velocity correlation Parameter in theoretical correlation for C D (Re) Coefficient for turbulent destruction of dissipation Coefficient for turbulent production coefficient of dissipation Coefficient for relaxation time scale of Green’s function Coefficient for linear strain term in phenomenological model for pre-stress Coefficient for relaxation term in phenomenological model for pre-stress Eddy viscosity coefiicient Viscous destruction of dissipation Turbulent Deborah number One-dimensional energy spectrum Fluctuating force Turbulent pie-stress Anisotropic portion of pre-stress Unit dyadic ” Integral pr0perty of the turbulence; I 1 = IBLLJE 0 ”N are- in, m Q re mg 11:5 Integral property of the turbulence; I2 = Iii-3%;- (B L L) d§ 0 Integral property of the turbulence; I3 = I-é—deé Fluctuating Reynolds stress Turbulent kinetic energy Energy spectrum wavenumber Linear convective-diffusive operator Memory function for the Green’s function Empirical parameter for the prescription of turbulent destruction Turbulent relaxation group; ratio of relaxation time of the Green’s function to characteristic time scale of the mean field Instantaneous pressure; parameter in theoretical correlation for CD (Re) Fluctuating pressure Mean pressure Turbulent production of dissipation due to vortex stretching Ratio of the relaxation time scale of the pre-stress to the growth rate of the turbulent kinetic energy or homogeneous turbulence Normalized Reynolds stress Proper orthogonal, time-dependent frame transfonn ation dyadic Redistribution of Reynolds stress components due to the isotropic portion of the pro-stress Turbulent Reynolds number Empirical parameter for the prescription of turbulent destruction Parameter in theoretical conelation for CD (Re) Magnitude of velocity gradient d (uz)/ dy Symmetric part of mean velocity gradient Velocity derivative skewness Time xiii Triple longitudinal velocity correlation LLL u Instantaneous fluid velocity r_¢' Fluctuating fluid velocity (3) Mean fluid velocity (g'g' Reynolds stress (___ Antisymmetric part of mean velocity gradient Greek Symbols (1 One-half the trace of the pre-stress [3 Linear strain parameter in phenomenological model for the pie-stress A Minimization parameter c Scalar turbulent dissipation 5 Permutation triadic 11 Kolmogorov length scale A Taylor microscale; relaxation time scale for the pre-stress A, Eigenvalues of the characteristic equation for the transient development of homogeneously sheared turbulence Molecular viscosity v Kinematic molecular viscosity Vt Eddy viscosity § Dimensionless separation distance for turbulent correlations; dimensionless evolution time for homogeneous shear flow EC Cut-off separation distance for approximation of turbulent integrals p Density ID Characteristic timescale for destruction of dissipation 1P Characteristic timescale for production of dissipation TR Relaxation time for the Green’s function xiv (3) Mean molecular stress S_2 Frame rotation rate 2 Frame rotation dyadic Subscfiptr/Superscriptr a Asymptotic value e Experimentally observed value 0 Initial value or limiting value at zero Reynolds number no Limiting value at infinite Reynolds number ~ Dimensional quantity XV CHAPTER 1 INTRODUCTION 1.1 Background Turbulence Closure Problem The ability to predict the low order statistical properties of turbulent flows is an important area of research. Although the exact equations governing the instantaneous pressure and velocity fields of constant density, constant viscosity, Newtonian fluids are known, the field equation for the mean velocity (1'. e. the Reynolds equation) is statistically unclosed. With E = (u)+r_a' andp = (p)+p'. 3(5) 7+(§)°V(y>+291\(§) = Q - (2 _V[ 152 _ (_ a5) 2 (_ "-"]+vv2-v. (2.!» (1.1) The statistical correlation p (r_¢'r_¢') is the Reynolds stress and accounts for the transport of momentum associated with the fluctuating field. In the above equation, the mean fields and the fluctuating fields are relative to a frame of reference rotating at a constant angular velocity £2. 9 is related to the temporal connection between the inertial and non-inertial frames by the following expression Qa—Q.QT = . 9. (1.2) It“ In Eq. (1.2), £2 is the rotation dyadic, g. is the permutation triadic, Q is an orthogonal dyadic operator (to. g . QT = _I), and g is the time derivative of g. The orthogonal operator 9 defines an arbitrary, time-dependent coordinate transformation: 9 = vx’, (1.3) where V is the gradient operator in the rotating frame of reference and f. is the position vector in the inertial frame. For constant density fluids, the mean velocity field also satisfies the continuity equation: V- (y) = 0. (1.4) Eddy viscosity type models for the Reynolds stress presume that the mean velocity. field and the second order correlation (u'g’) are related by the following model 21: (143') = - 3:1 - 2V,<§>r (1.5) where k denotes the kinetic energy per unit mass of the fluctuating field, 21:505-330. (1.6) and v‘ is a turbulent eddy viscosity coefficient which depends on the local statistical state of the turbulence. (.38) is the mean strain rate, which is traceless for constant density fluids. This results in a diffusion-type transport model which permits estimates of the mean field behavior, provided a model for the scalar eddy viscosity can be specified. Turbulent flows can be computed directly by solving the continuity and the Navier- 3 Stokes equations. Because turbulent fluctuations are three-dimensional, the spatial grid required for direct numerical simulations is also three-dimensional. The computational grid must be fine enough to resolve all the physically significant scales of motion. Even the “simplest” turbulent flows require a four-dimensional discretization (one ill time and three in space) whose resolution is directly dependent upon the Reynolds number. As the Reynolds number increases, the grid must become finer because the Size of the small scales decreases. This reduction in the small scales can be seen by comparing the Kolmogorov and integral length scales for isotropic turbulence (Tennekes and Lumley, 1972; Hinze, 1975): ( 3/ )1/4 I". = 1+. “Rafi/4, (1-7) In the above expression, the dissipation e is defined by £‘=’V(Vt_¢': (Vg'fl), (1.8) and the integral length scale is given by the integral of the two-point longitudinal velocity correlation r = j (u'L(_J_r, t) u'L(§+f, r) )dr. (1.9) 0 Although direct numerical simulations (DNS) of turbulence are possible, they have only been performed for simple geometries at relatively low Reynolds numbers. For instance, Kim et al. (1987) have reported DN S results for fully-developed channel flow at Reynolds numbers less than 10,000 based on the channel half-width. Such simulations provide important statistical information which is difficult to observe experimentally. However, the use of direct numerical simulation for practical engineering flows with 4 complex geometries and/or high Reynolds numbers would require computational power beyond what is presently available. Thus, improved models capable of predicting the low order statistical properties for a wide class of turbulent flows remains an intense area of research. Turbulence models can yield physical insight into the behavior of the mean field provided the important qualitative properties of the Reynolds stress are retained. Thus, turbulence modeling has several distinct advantages over direct simulation. For instance, because of the intrinsic three dimensional and unsteady nature of turbulent fluctuations, direct simulation cannot take advantage of simplifying statistical properties. For example, a flow which is statistically stationary must still be simulated as a transient. Furthermore, if a specific flow geometry possesses symmetry planes and/or axisymmetric features, the flow must still be simulated as fully three dimensional. Thus, a Reynolds stress model which is not overly complex will almost certainly be computationally faster (and thereby lessiexpensive) than a direct simulation of the same flow. Realizability The term realizable as applied to Reynolds stress models relates to whether or not the eigenvalues of the statistical correlation (u'u’) are non-negative. An important aspect of this concept relates to invariant mapping, as discussed by Lumley [197 8]. Specifically, the invariants of the anisotropic stress, defined by (1414') b= ' '= 2k :1, (1.10) le are computed and cross plotted to form a ”-11! phase plane of permissible turbulent states. The invariants of b are defined as 11: tr (9 - 5) (1.11) and III; tr (_b - -b) . (1.12) By definition, the first invariant of the isotropic stress is zero (Iztr (:13) = 0). Therefore, at least one of the anisotropic eigenvalues associated with b is negative. Because II is the sum of the squares of the anisotropic eigenvalues, the second invariant is always positive; however, III can be either positive or negative, depending on the local dynamic state of the turbulent field. There are two key mathem atical properties of the Reynolds stress which restrict all turbulent States to a subset of the ”-111 plane. These properties are: (l) the trace of the Reynolds stress is twice the total kinetic energy of the turbulence (see Eq. (1.6)); and, (2) the Reynolds stress dyadic is positive semi-definite, i.e., z - (u'u') - 3 2 0, (1.13) for any non-zero _z in 153. These two conditions imply that the anisotropic stress is traceless and its eigenvalues are bounded by 2/3 and -l/3 (see Appendix I). Lumley [197 8] showed that all realizable turbulent states map into a quasi-triangular domain in the [HI] phase plane, illustrated by Figure 1.1. The portions labeled A, C and E are the vertices of the domain, while B, D and F are the boundaries. Each of the portions are listed in an accompanying table and described mathematically with either the coordinates of the vertices or the equation for the boundaries. If a coordinate pair (1!, III) falls on or within the boundaries of the L-diagram (where “L” stands for “Lumley”), then this anisotropic state may be associated with a realizable Reynolds stress. On the other hand, no realizable anisotropic state can fall outside the L-diagram. Non-turbulent, albeit realizable, states A II A B C Ila-"(99) F ”litre-2.2) D—> E III P Description of Points on Anisotropy Simplex I Invariant" I Eigemnluu of the Reynolds Stress Name A 1 Component B 2 Component Anisotropic II = 2/ 9 + 211! 1-x x 0 s x s 112 c 2 Component Isotropic 1/6 L -1I36 1/2 112 u Oblate Axisymmetric m = -6 (II/6) ”I x x 1-2: 1/3515 112 E 3 Component Isotropic 0 l 0 1/3 1/3 113 F Prolate Axisyrnmetric III = '-6 (II/ 6) 3’4 x x 1-2: 0 s x 5 1/3 -4 Notes:(1) “n-Component” means that the Rdnolds stress has “n” non-zero eigenvalues. (2) “Axisymmetric” two of the eigenvalues of the Reynolds stress are equal. Invariants 0f the anisotropic stress. Figure 1.1 L-Diagram for Realizable Anisotropic States 7 may fall within the L-diagram. Further a priori, model—independent restrictions of the anisotropic turbulent states within the L-diagram remains an open question for research. The L-diagram will be employed in this work to represent the local state of turbulence associated with different turbulent models with different experimental data sets. k-e Turbulence Model The most commonly applied turbulence model is the Boussinesq approximation (defined by Eq. (1.5)). The k—e theory uses an eddy viscosity related to the kinetic energy and dissipation associated with the fluctuating field: 2 v, = Cuk /e. (1.14) The universal coefficient Cu is estimawd from plane jet data to be about 0.09 (Launder and Spalding, 1972). The first models of this type date back to Prandtl (see Speziale, 1991) and use empirical, algebraic models for the specification of the characteristic turbulent scales. Later models have moved away from this and solve differential equations to specify the turbulent scales. These are generally termed two-equation models, as they require two differential equations for the scale-determining parameters. The most prevalent example of this approach to Reynolds stress modeling is the be model of Hanjalic and Launder (1972), hereinafter denoted as the standard k-e model. Transport equations for the turbulent kinetic energy (1:) and the turbulent dissipation rate (a) are used to compute turbulent time and length scales at each point in the spatial domain: at . .. __ . 2 b7 _ -(r_u_r).V(u) e+v [[v+6k]Vk]’ (1.15) De 8 , , . _ 82 v, ‘5‘- - -CP-k-(l-‘l-‘ ).V(y) CDI+V [[V+O’c]ve], (1°16) 8 While there are many variants on the two-equation model, the standard k-e and its low Reynolds number derivatives (see Patel, et al., 1984, for an overview) have become the most widely tested and applied engineering turbulence models. The continuing use of the standard k-e model is rooted primarily in its relative simplicity combined with its ability to provide some acceptable predictions. Some of the computational advantages of the k-e model include the fact that it is simple to implement into existing flow predictor codes and that there is not much additional computational demand inasmuch as only two additional scalar transport equations need be solved simultaneously with the Reynolds equation. Additionally, the k- 6 model has been found to be computationally stable and robust, being largely independent of the initial guess for field variables. Because of these advantages and despite its limitations, the k-e model has become the standard against which other models are judged. The k-e model is, therefore, well documented in a variety of flows and summarizes a wealth of information. In this respect. any new model ought to be at‘least comparable to the he model in some basic test flow cases. One of the many criticisms of the Boussinesq approximation lies in the fact that the eddy viscosity coefi'rcient is a scalar-valued function that only depends on the mean field implicitly through I: and 8.. Therefore, at a fixed position in the flow field, Eqs. (1.5) and (1.10) imply that = ___B , (1.17) for all combinations of base vectors 5“ and 513' Eq. (1.17) may fit one of the cross correlations of the anisotropic stress for simple shear flows, but cannot explain the anisotropic distribution of kinetic energy among the components of the fluctuating velocity. Indeed, Kitoh [1991] has shown experimentally that Eq. (1.17) is inconsistent with the observed behavior of the anisotropic stress for the decay of swirling flow in a 9 pipe and concluded that the Boussinesq model was an inadequate closure for complex flows with strong curvature effects. Simple Shear Flows While eddy viscosity models typically make good predictions for two dimensional mean shear flows (Hanjalic, 1994), they do have many limitations. For example, the Boussinesq approximation does not predict a normal stress anisotropy for fully developed channel flows or for homogeneous shear flows. A simple shear flow is defined by d(uz) V9!) = dy 5,52, (1.13) where S t d(uz)/dy. For this flow, the Boussinesq model predicts an equipartition of turbulent kinetic energy among the fluctuating components of the velocity, which is clearly unphysical. Tavoularis and Corrsin [1981], Tavoularis and Kamik [1989], as well as Gibson and Kanellopoulos [1987] have all observed normal stress anisotropies for homogeneously sheared turbulent flows. ' Additionally, the Boussinesq model is unable to explain the transient relaxation effects observed in return-to-isotropy experiments or the reorganization of turbulent energy associated with the transient development of homogeneous turbulent shear flows. For simple mean shear flows, the anisotropic stress associated with the Boussinesq model can bewrittenas s=-(%)e1ea+e» The invariants of the anisotropic stress for this situation follow directly from Eq. (1.19), viz... 10 I = "(5) = 0, (1.20) 2 11 = "(5.5) = (Ci/2) (1:3) , (1.21) and III = tr (5- - g) = 0. (1.22) Experimental studies (see, esp., Choi and Lumley, 1984; LePenven, et al., 1985) clearly Show that turbulent flows which are initially anisotropic gradually relax towards an isotropic state when the mean strain rate is removed. The Boussinesq model. however, predicts that the Reynolds stress responds instantaneously to a sudden change in (3'). Therefore, the intrinsic memory effects observed in low-order statistical properties of turbulent flows cannot be represented by Eqs. (1.5) and (1.19). Experimental data for homogeneous shear flows (see, esp., Tavoularis and Kamik, 1989) also exhibit a transient relaxation to a strongly anisotropic asymptotic state. figure 1.2 illustrates the trajectory of the relaxation process on the L-diagram. For simple shear flows, the accessible realizable states permitted by the Boussinesq theory all lie on the line III = 0 and OSIISZ/9. Clearly the Boussinesq model cannot form the basis of a universal closure theory for the Reynolds equation inasmuch as the Boussinesq states are too restrictive. It is noteworthy that for simple shear flow the Boussinesq model is realizable for III = 0 and 0 $11 = (cf/2) (ii—")2 5 2/9. (1.23) With Cu = 0.09, the above inequality for the Boussinesq model implies that the ratio of turbulent to mean field time scales must be restricted to IIEtr(=b-=b) IIIEtr(_b--b-§) (1/6, —l/36) Boussinesq States Isotropic State \ Figure 1.2 Anisotropic States for Homogeneous Shear Flows (2/9, 0) A11 V (II, III) = (2/3,2/9) Esp. Data - Homogeneous Shear O Asymptote Transient Approach to Asymptote III ‘ r 12 (:83) < 7.4. (1'24) 8 Thus, application of the Boussinesq model to turbulent flows for which k/e » 7/ S would require Cu to be a function of Sk/e in order to maintain realizability. Hanjalic [1994] has also noted the poor performance of eddy viscosity type models in flows with streamline curvature, especially in swirling flows. This is probably related to the misrepresentation of the underlying mechanism which distributes the kinetic energy of the turbulence among the three components of the fluctuating velocity. A primary source of velocity fluctuations in turbulent flows is the convective coupling between the mean field gradient and the fluctuating field, via, 11' - V (a). If (g) has the following local structure for a fully developed swirling flow (g) = (u9)(r)§9. (1.25) then the local mean field gradient can be written as the sum of the following symmetric and antisymmetric strain rates V (11) = 1nd. r(‘(_u0>)[‘,39 +e e +--— (r(u9)) [59 e e05} (1.26) 9' ' 2 r dr Thus, the convective coupling between V (a) and r_r' has two distinct contributions: u' (u ) ° ° (1.27) . . d e -V <14) - u r(;;)s,- s,~ The r-component, or cross-stream component, of Eq. ( 1.27) provides a means to shift the energy into the pressure field, and, thereby, increase the highly anisotropic distribution of turbulent kinetic energy among the three components of the fluctuating velocity. The coupling between the fluctuating field and the strain associated with streamline l3 curvature does not arise in simple shear flows in inertial frames inasmuch as S +S V(u)= — y_e_ez +ezey+ zygz—ggzey], (1.28) which yields 'oV<>- '(i< >) (129) r_4 u — uy dy uz ez. . Only streamwise fluctuations induced by the mean field occur for simple shear flows. However, as implied by Eq. (1.1), velocity fluctuations in non-inertial frames are produced by convective coupling between the mean field gradient and the rotational strain rate associated with the frame (see Appendix B). Thus, for simple shear flows in a rotating frame of reference, velocity fluctuations occur by strearnwise convective coupling and, most significantly, by cross-stream convective coupling with the rotational strain: u'-[V(u)+£2] = u' (1(u)+(2)e -u’ Ga (130) - - = y dy z -z z -y‘ ' Eq. (1.30) assumes that the rotation dyadic is constant and given by 2 = “(5,521.3 ). (1.31) where (2 represents the angular velocity of the frame rotation about the x-axis. Thus, as previously noted by Speziale [1991] and others (Speziale, er al., 1990), an analogous streamline curvature coupling between the mean field gradient and the fluctuating field can be achieved by a simple shear flow in a non-inertial frame. This observation provides a means to study the efficacy of phenomenological turbulent models under relatively simple kinematic conditions with an anticipation that a turbulence model which performs l4 adequately for simple shear in a non-inertial frame may provide a reasonable approach to flows with strong streamline curvature in inertial frames. The ad hoc basis for this conjecture stems from a comparison between the velocity gradient in an inertial fiarne (see Eq. (1.28)) and the velocity gradient in a rotating frame, viz... V <10 = fiestas] + (i + “Isa-5.9.]- “32> Thus, the mean field gradient associated with streamline curvature in an inertial frame is mathematically analogous to the effective mean field associated with a simple shear flow in a rotating frame provided 59"}(903) and (29%. r r I' 1.2 Objectives This research examines a new class of phenomenological models for the Reynolds stress. The study supports a long-standing goal of turbulence research to achieve a practical, albeit physical, statistical closure of the Reynolds equation for the mean field (see Eq. (1.1)). The proposed approach addresses some of the limitations of eddy viscosity type models, while maintaining some of their advantages. The following five issues are central elements of this investigation: (1) realizability; (ii) primary and secondary normal stress differences; (iii) spatial and temporal relaxation, or memory, phenomena; (iv) frame dependence; and, 15 (v) universality. The Boussinesq model for the Reynolds stress (Eq. (1.5)) which provides closure for the ubiquitous k—e theory (Hanjalic and Launder, 1972; Speziale, 1991) as well as for the sub- grid scale model for large scale eddy simulations (Bardina, er al., 1985), misrepresents some aspects of the above five physical issues. Items (i)-(iv) are fundamental and necessary for any closure hypothesis which aims for some limited form of universality. However, the possibility of a low-order, universal closure theory for turbulent flows remains a Speculative, albeit desirable, goal. This research is limited to the development of a closure for the normalized Reynolds stress t I: < . '). (1.33) Previously published experimental data for homogeneously sheared turbulent flows provides an extensive resource to guide the development (see, esp., Tavoularis and Kamik, 1989; Tavoularis and Corrsin, 1981; Gibson and Kanellopoulos, 1987). In order to develop some understanding of the potential applicability of the proposed approach to flows with streamline curvature, the new theory is used to predict the distribution of asymptotic states for homogeneously sheared turbulence in a rotating frame of reference. This study complements a parallel development of this theory for fully developed channel flows by Weispfennig [1997]. 1.3 Methodology The theoretical approach to closure stems directly from the continuity equation and the 16 equation of motion for a constant density, constant viscosity, Newtonian fluid. An explicit dependence of the Reynolds stress on the reference frame (i.e. Issue (iv) in Section 1.2) is introduced by a direct analysis of the equation for the fluctuating velocity field relative to a rotating frame of reference (Chapter 3 and Appendix B). The development introduces one level of memory into the closure by using a smoothing approximation related to the relative relaxation of the space-time structure of the turbulence and a convective-diffusive Green’s function associated with the non-local transport of momentum fluctuations by the mean field. This strategy, which has been previously employed by Petty [1975] and others (see, esp., Hill and Petty, 1996) for turbulent mass transfer, leads to an algebraic pre- closure for the Reynolds stress which relates (r_r'_u') to the gradient of the mean field, to a scalar-field'relaxation time 1R associated with the temporal structure of the turbulence, and to a statistical correlation hereinafter referred to as the turbulent pre-stress (1R ([I'». The pre-closure theory in a rotating frame of reference has the following structure 33’; (2'13) -5__1 = cm). (1.34) where the operator .3 is defined by £5!+IR(V (u)+£2). (1.35) The fluctuating vector _f‘ is related to instantaneous fluctuations in the instantaneous Reynolds stress and to pressure fluctuations: f'a v . Egan - (5:50]. (1.36) An inertial frame pre-closure for (u’u’) follows directly from Eqs. (1.34) and (1.35) by setting!) = 0. 17 A phenomenological relaxation model is introduced for the pie-stress and this approach provides an additional means to further address the physical issues associated with Items (i) and (iii) in Section 1.2. The pie-stress is formally decomposed into an isotropic and anisotropic contribution: R0,):20t 1: _ _ -3— (1.37) "at _I+ where the anisotropic pre-stress :1 is both symmetric (1:?! = ET) and traceless (tr (£1) = 0). This strategy provides an explicit, self-consistent model to relate the isotropic coefficient a to the mean field properties inasmuch as Eqs. (1.34) and (1.37) imply that 2t: = trig]: (u'u')-z=i]. (1.38) Thus, no additional closure hypothesis is nwded to evaluate the isotropic contribution to the pre-stress. In Chapter 3, the theory is applied to homogeneously sheared turbulence under the assumption that the anisotropic pre-stress is unimportant. This isotropic pre- stress theory (IPS-theory) provides a theoretical limit case to compare with the anisotropic pre-stress theory (APS-theory) developed in Chapter 4 for nontrivial if. However, the IFS-theory provides a relatively simple extension of the Boussinesq model which incorporates major issues associated with Items (1), (ii), Crv), and (v). The inability of the IFS-theory to account for the well-documented phenomena associated with the return-to- isotropy partly motivates the extension to the APS-theory in Chapter 4. In Chapter 4, the following property is attached to the anisotropic pie-stress in order to address Item (v) in Section 1.2: The anisotropic pre-stress is an objective statistical property of a turbulent flow 18 and any phenomenological closure model for {:1 should be frame invariant. The above assumption has been applied incorrectly to phenomenological models for (r_r'r_r'), of which the Boussinesq model is a prime example as well as some algebraic theories proposed by Speziale [1987]. The assumption here is that the operator :1 introduced by the pro-closure provides the necessary frame-dependent effects manifested by the Reynolds stress. The isotropic pro-stress contains a self-consistent frame dependence inasmuch as the anisotropic pro-stress is traceless: tr (5?) = 0. (1.39) Thus, the above closure hypothesis related to the behavior of the anisotropic pre-stress implicitly supports the goal of developing a universal closure. The IPS- and APS-theory are used in Chapter 5 to predict the effect of frame rotation on the asymptotic states of homogeneous shear. The phenomenological coefficients introduced by the model are sealed with the turbulent kinetic energy It and the turbulent dissipation e. The large Reynolds number (ie. kz/vc » l) scalar transport equations are used to estimate the behavior of k and 6. However, the coefficients in the k-e equations are recalibrated using the IPS- and APS- closures for the Reynolds stress. The benchmark flows used to calibrate certain aspects of the theory include: - homogeneous isotropic decay, ° retum-to-isotropy, and - asymptotic homogeneous shear. Experimental data for these flows relative to an inertial frame are used to estimate model coefficients. 19 In Chapter 2, the classical problem of homogeneous isotropic decay is examined to determine the effect of the turbulent Reynolds number (kz/ve) on the decay process. An approximate integral analysis of the Karmén-Howarth equation is employed to relate velocity derivative skewness data to the dissipation coefficient which appears explicitly in the e-equation (see the coefficient CD in Eq. (1.16)). An extension of the APS-theory developed hereinafter to flows for which k2/ve at 1 will require a complementary low Reynolds number k-e theory to calculate the turbulent time scales needed for momentum transport Figure 1.3 illustrates the su'ucture of the theory and highlights the specific area of focus. Presently, the development is incomplete and requires further theoretical work. Applications to specific test flows (see the recommendations in Chapter 7) are essential to detect flaws in the proposed closure. However, the methodology developed as part of this research provides a clear and unambiguous framework to introduce further improvements in the proposed low-order closure of the Reynolds stress. 20 F luctuating Field Equations 3 _ 2 .__._ _. (gun-V vV)2 — r W!) I Figure 1.3 Closure Methodology V - u' = 0 Chapter 3 Preclosure k-e Eq , ns T r t _ Q '<‘_‘.‘f)’§ "qu Chapter-2 A A ChapterS Isotropic Pre—stress Reynolds 2a Stress ‘ tRQ’f') = -3—:I ’ (g'g') Chapter3 Anisotropic Pre-stress Mean Field Equations . V mm = 391+H (-a-+(u)-V-VV2)(u) = -v--2v2(IV(Vr_r)1E[V(vr_r')7]). (2.3) Although Eqs. (2.2) and (2.3) are formally exact, two unknown statistical correlations 21 22 appear in Eq. (2.3). These two terms correspond to the production of dissipation due to vortex stretching and the destruction of dissipation due to viscous effects, respectively. The modeling approach for the dissipation equation introduced by Hanjalic and Launder [1972] assumes that the destruction of turbulent dissipation scales with the local dissipation. Thus, the following statistical ansatz is used to represent the underlying physical effect associated with the terms on the right-hand-side of Eq. (2.3): §=_i co dt 71) where to represents a characteristic relaxation time for the dissipation. The conventional k - e theory further assumes that TD ~ k/e for large turbulent Reynolds numbers. Extension of this closure to low Reynolds numbers involves the introduction of an empirical, albeit universal, destruction of dissipation coefficient CD (Re) : de 22 — = -C R -— . 2. dt D( e) k ( 5) The utility of decaying isotropic turbulence is that Eqs. (2.2) and (2.5) provide a means of determining the closure coefficient CD (Re) . Much experimental (Batchelor and Townsend, 1948; Comte-Bellot and Corrsin, 1971; and Sirivat and Warhaft, 1983) and numerical (Bardina et al., 1985; Speziale et al., 1987; and Mansour and Wray, 1994) data are available for the analysis of this model parameter. Eqs. (2.2) and (2.5) earl be combined to yield the following equation for the turbulent time scale k/c: ditUc/e) = CD (Re) - 1. (2.6) Experimental data available for isotropic decay indicate that the time scale k/e increases 23 throughout the decay process, which implies that CD (Re) > 1. Analogously, Eqs. (2.2) and (2.5) can be combined to yield an equation for the turbulent Reynolds number Re ( s kz/ve), i = .. E * dt (Re) [C D (Re) 2] Rek. (2.7) Eq. (2.7) implies that the decay process can occur at a constant Reynolds number Ret defined by CD (R?) = 2. This critical Reynolds number Re‘ is stable if CD (Re) > 2 for Re < Re. and CD (Re) < 2 for Re > Re. . However, experimental measurements for isotropic decay indicate that Re decreases throughout the decay process. Therefore, the dissipation coeflicient CD (Re) , introduced by the closure model given by Eq. (2.5), must satisfy the following inequality: 1 < CD (Re) < 2. (2.8) With 1 500, see Comte-Bellot and Corrsin, 1971), Hanjalic and Launder [1972] estimated the decay coefl‘icient to be about two. This value was later refined to 1-92 by Launder and Spalding [1974]. Subsequent investigators introduced empirical functions to account for the effect of the Reynolds number on the coefl‘icient CD. For instance, Hanjalic and Launder [1976] proposed the following expression CD = 1.8[1-gexp (—[%]2)], ' (2.10) whereas Lumley [197 8] suggested 18 CD = 1.4+0.49exp[-J;;]. (2.11) Both of the above empirical expressions for CD employ the final decay coefficient estimated by Batchelor and Townsend [1948] (Le. CD -) 1.4 as Re 40). The high Reynolds number asymptote of Comte-Bellot and Corrsin [1971] is used in Eq. (2.10); Eq. (2.11) incorporates the limit CD —> 1.89 as Re -> oo. Both Eqs. (2.10) and (2.11) assume that the destruction of dissipation coefficient increases monotonically as the Reynolds number increases. 27 Mansour and Wray [1994] studied decaying isotropic turbulence by direct numerical Simulation. They computed the skewness of the velocity derivative as well as the destruction-of-dissipation coefficient, C D. The velocity derivative skewness is defined by (see p. 58 Monin and Yaglom, 1965): . =-<@: M (3:11”. and, as discussed by Tennekes and Lumley [1972], is an important statistical property related to vortex stretching (also see Monin and Yaglom, 1965): -2v([vr_r': (Vr_r' . V3371); R— - ail—55“,; (2.13) Eq. (2.13) will be derived from the classical Karmén-Howarth equation for isotropic turbulence in Section 2.2. For Re > 100, Mansour and Wray [1994] observed that S k 4- 0.3 to 0.4. For Re —) 0, S k = 0. Figure 2.2 shows the effect of Re on the velocity derivative skewness given by Mansour and Wray [1994] and Tavoularis et al. [1978]. The empirical representations of C D given by Eqs. (2.10) and (2.11) are also shown for comparison. The goal of this chapter is to develop a relationship between the destruction-of- dissipation coefficient C D and the turbulent Reynolds number. In the next section, the rate of change of dissipation will be formally identified as the difference between two statistieal properties, viz, the production of dissipation due to vortex stretching and the destruction of dissipation due to viscosity. For decaying isotropic turbulence, the dissipation rate exceeds the production rate, with the net result that the dissipation turnover time (615 Eqs. (2.4) and (2.5)) is positive: _1. e (n -r) E >0. (2.14) 1‘1) 28 EoBEooU consummma mo 534580 2,: .8 223.230 can Rogue—m one .8 Sea 8882: mm m: mm: 3.: u-q-«u q 4 —--q.—« u d —q-qudd4 d‘ $8: $53 \ 83: S53 pea 6:252 Tm— 4..— n.— o.— h.— w.— o.— ~.N ouswE mm m: m: m: 8: 3: -----EH d —--qq- d 1 ddqdqqq - 1 l—uqd-qqd 1 0° See: >83 93 58:32 I 33: .3 3 $3302; D D . .o a oofibmom _er>.m . o D I Lao DD a . l I I I two I a a D P II I . ca 0 J a D— a m. I a a a a men I D I D .3 1 n1 0 I Ice 0 I D . .. . o a PD 0 I h . 3 on a B one a a a a D U DO 0,0 29 The destruction coefficient .0 and the production coefficient P can be related to specific statistical pr0perties of the flow and possibly provide a means to understand the underlying behavior of the k - 8 model coefficient C D. It follows directly from Eqs. (2.4), (2.5), and (2.14) that CD (Re) = s(Re) —P(Re). (2.15) A better understanding of CD (Re) is essential for further improvements of the widely used It — 8 transport equations for turbulent flows. In Section 2.3, a semi-theoretical approach is used to relate CD and Re. Subsequently, an integrated form of the Karmén-Howarth equation (hereinafter referred to as the IKH- equation), which links the double and triple longitudinal velocity correlations, is employed to develop an approximate relationship between .8 (Re) and 1? (Re) . An empirical representation for a (Re) , which agrees with experimental results for high and low Reynolds number asymptotes, permits the velocity derivative skewness factor to be predicted from the IKH-equation. The adjustable parameters introduced by the approximation are determined in part from the experimental data summarized by Figure 2.2. In Section 2.6, the non-linear dynamic equations given by Eqs. (2.2) and (2.5) are solved in the time domain for an arbitrary, albeit realizable, initial state (kg, 80) . 2.2 Local Analysis of the Kurman-Howarth Equation The Karman-Howarth equation stems directly from the Navier-Stokes equation and provides the following fundamental relationship between the double and triple longitudinal velocity correlations for isotropic turbulence (see Karman and Howarth, 1938; Hinze, 1975; and, esp., p. 122 Monin and Yaglom, 1965), $3”) = 2v (%+;)%(Ru) + (r%+§)im' (2.16) In Eq. (2.16), ELL (r, t) represents the double longitudinal velocity correlation, Run, t) = (u'LQr, t) u'L(x+f, r) ), (2.17) and 7m (r, t) is the triple longitudinal velocity correlation, 71,10, r) = (trig, t) u'L(§+f, r) ). (2.18) Ru (r, t) is an even function of r (z||_r||) and I'LLL (r, t) is an odd function of r. Therefore, a Taylor series representation of the two correlations can be written as follows, - .. 321} 3B Bu = Bu(0,t)+-21;[ :1] r2+-41-'-[ 1"] r‘+... , (2.19) ' 3’ r20 ' 3" r=0 and - a i rm = i'[ ‘3”) r3+... . (2.20) 3. a, '80 An equation for the kinetic energy follows directly from Eq. (2.16) by setting r = 0: dB 0 10v 32%” 221 It LL( 9‘) ' 3’2 - (- ) r=0 Because ELL (0, t) = 2k/ 3 for isotropic turbulence, Eqs. (2.21) and (2.2) imply that (2.22) 31 With the Taylor microscale 1 defined by the following expression (Tennekes and Lumley, 1972): 3253 R (0, r) [ :L] = "—L-L—3—- , (2.23) 37‘ 7 = 0 2- Eq. (2.22) can be rewritten as I} 0, r e = 15v—L—L—(—) = 3315 . (2.24) v 2.2 2.2 Eq. (2.21) governs the rate of viscous decrease of the turbulent kinetic energy for an isotr0pic decay process. The Taylor microscale, or the energy dissipation length scale, is defined by the double longitudinal velocity correlation; It can also be related to the second moment of the velocity derivative (see p. 143, Monin and Yaglom, 1965). Eqs. (2.23) and (2.24) suggest that an equation for the energy dissipation also follows from the KArmén-Howarth equation (i.e. Eq. (2.16)) by differentiating twice with respect to r and then setting r = 0. This analysis yields (see p. 144, Monin and Yaglom, 1965) 33 as at 1[ :1”) = Ev[ 1"”) +1 [ 2“] . (2.25) d‘ 3r 0 3 3r , = 0 3 0r rs r=0 Eqs. (2.23) and (2.24) can be used to rewrite Eq. (2.25) as 2 ds 2 2? .. -(s-e)-k-, (2.26) where 4 a 1‘} pal._ 7‘ [ :1) (2.27) 15Bum) ar ,30 32 r=0 and P __ -7 Auk—e [firm] (2 28) - 3/2 3 ' 355- (Buton 3' The first term on the right-hand side of Eq. (2.26) causes a decrease in energy dissipation due to viscous damping as .8 is always positive, i.e., air [ “] >0. (2.29) 3r4 r=0 The second term on the right-hand side of Eq. (2.26) accounts for the “production” of dissipation by vortex stretching (see p.144 Monin and Yaglom, 1965; and, esp., p. 83 Tennekes and Lumley, 197 2). For isotropic turbulence, P is generally positive inasmuch a3i~ [ Lu] <0. (2.30) 3 at r-O The production term can also be written in terms of the velocity derivative using the following result (see Appendix F) )3 2 (again, = ‘G—Z')3>/[<@i)2>]mE-sr (2.31) .. 3/ (Ba, (0)) Eq. (2.26) shows that the net destruction of dissipation coefficient CD (Re) has two contributions (see Eq. (2.15)). In the next section, an integral analysis of the Kérman- Howarth equation is used to estimate the effect of the Reynolds number on the coefficient CD by using velocity derivative skewness data (see Figure 2.2). 33 2.3 Global Analysis of the Karman-Howarth Equation Integrating Eq. (2.16) from r = 0 to r = on yields an integrated form of the Karman- Howarth equation (IKH-equation): d I ~ "1 ~ "1 - E gaudr = 8v (1) :51r-(Bu) dr+ 4 g :TLudr . (2.32) The following conditions on the double and triple longitudinal velocity correlations have been used to obtain Eq. (2.32) al‘iu .. [8r ] = 0 and (Tm)r=0,~ = 0' r20,» With dimensionless velocity correlations defined as B , t Bu. 2 —_. LL (r ) (2.33) an (0, t) and I’ , r T iii—(Ll (2.34) Tm (0. 0 Eq. (2.32) can be rewritten as - . - 3/2 52”“- (0, t) 2.11] = £2.23” (0, t) 12 + 4 (Bu (0, t) ) 13. (2.35) The dimensionless integrals in the above equation are defined as follows: I1 = jaunt, (2.36) 0 I2 = 15%er d§,and (2.37) 0 I3 = jédeg, (2.38) 0 where g: r/7t. Because I: = 3131.110, r) /2, Eq. (2.35) may also be written as d _ SW: 32 3/2 5053.11) - T12+f-3:k 13. (2'39) With 3.2 = lka/ c, it follows directly from Eq. (2.6) that temporal changes in the Taylor microscale are governed by d2. lOve _ = c -1 _. 2.40 dt ( D ) 4k ( ) Therefore, because C D > 1 for isotropic decay in an inertial frame of reference, Eq. (2.40) implies that the integral microscale increases during the decay process. Note, however, that the rate of change in A decreases inasmuch as the ratio k/s —) co as k and s decay (see Section 2.1). Eqs. (2.40), (2.2) and (2.5) can be used to eliminate the Taylor microscale from Eq.(2.39) with the result that k d 8 I2 ’3 5Re C -3 +2- —I I =- -—+— —— . 2.41 ( D ) e[drn(1)] 5L, [Ii 3 l ( ) An important assumption regarding the decay process is that the integral I 1 depends only 35 on the instantaneous Reynolds number (kz/ve): 11 = 11 (Re). (2.42) The above universality hypothesis uncouples Eq. (2.41) from the initial conditions. Therefore, it follows from Eqs. (2.42) and (2.7) that 25F»: (11)] = 2 2.43 8 dt ( ) k dRedl"(Ir)] dln(ll) E[ d: dRe '2(CD-2)dzn(ke)' Eq. (2.41) can now be rewritten as dlnul) 8 ’2 13 ’Skel (CD-3)+2(CD—2)dln(Re) -— 3|:I—l+;; —3— . (2.44) The above equation is an integral property of the KArmén-Howarth equation and will be subsequently used to develop some understanding of the possible dependence of the destruction-of-dissipation coefficient CD as well as the velocity derivative skewness on the turbulent Reynolds number for 0 < Re < on. 2.4 Approximate Analysis of the Dissipation and Production Integrals The previously developed series expansions for the double and triple longitudinal correlations (see Eqs. (2.19) and (2.20)) can be rewritten in terms of §, .8, and P by using the definitions given by Eqs. (2.27) and (2.28): 253 1 4 B = l-- +— +... 2.4 36 and I 15 P 3 Lu Re l4§ ( ) Figure 2.3 illustrates the behavior of the three functions which determine the integrals I 1, 12, and I 3. As an approximation, the integrals in Eq. (2.44) will be estimated using a finite cut-ofl' separation distance related to the behavior of the dissipation integrand 12 (see Figure 2.3): §. 11 {nude .. gc 6§c+56§c+..., (2.47) g"1 20.9 .. -51 = _ 3 12 (j) gagwumg —§‘+168§c+"" (2.48) and g‘1 15 a = — —— 3 13 (j) Erma: ’Re ”get... (2.49) The cut-off separation distance §c, defined by Figure 2.3, scales with the Taylor microscale and the dimensionless viscous destruction coefficient, 9 14 — = —. 2. 0 gc ' 2. 53 ( 5 ) The experimental results of Tavoularis and Corrsin [1981] and the direct simulation results of Mansour and Wray [1994] show that the velocity derivative skewness S k < no as Re —) co; therefore, Eq. (2.13) implies that B . . . T LL / 181 three Taylor Series Terms ~\. , .1 E §\ .. \ 1l e JBLLdfi ’ o \ \ Osculating Parabola ; +§ Err ”/2 A 1&1“; ) l gag LL ‘ o l :5. is 4.1 a _ 1,4 [Egg—(Brad: o -l O >§ \ gr \ ~ 1 13 ~ I ETLLLdé o Figure 2.3 Qualitative Behavior of the Integrands of I 1, 12, and I3; Estimation of the Cut-off Integration Distance fie. 38 lim w~s r. (2.51) Re-Hn It follows from Eqs. (2.15) and (2.51) that limCD d—s >1 (2.52) Re—)0 Experimental data of Batchelor and Townsend [1948] at low Reynolds numbers imply that c; = 1.4; therefore, it follows from the above limits that so = 1.4. At large Reynolds number, the data of Tavoularis and Corrsin [1981] and the simulations of Mansour and Wray [1994] (see Figure 2.2) indicate that lim S =S°° ~O..4 (2.53) Re -) o- Therefore, the production of dissipation coefficient increases with the Reynolds number as follows (see the definition of 0’ given by Eq. (2.13)) 7 lim 0’ = —S°°JRe . (2.54) Re -+- 3J15 " Because 1< CD < 2, it follows that the asymptotic behavior of .B at large Reynolds numbers must be as follows lim 5-5” = s“./R_. (2.55) Ike-too CBS—7.]? S: where C; ~ 1.8 to 1.9 (see Section 2.1; Comte-Bellot and Corrsin, 1971; and Sirivat and Warhaft, 1983). It follows directly from Eq. (2.50) and the asymptotic behavior of the viscous destruction term that the integral cut-off parameter §c has the following limits at high and 39 low Reynolds numbers lim gc = J5 (2.56) Re —) O and lim §c = 0 (2.57) Re -) a- Each of the integral ratios appearing in Eq. (2.44) can be estimated in terms of Q. For instance, it follows directly from Eqs. (2.47) and (2.48) that 20.8 3 ’2 -§c 168§C+ T- l 3 , (2.58) 3 5 ‘ ic'zichic As Re —) co, Eq. (2.50) for EC implies that the ratio 12/1l asymptotically approaches —2/ 3 inasmuch as 2 -§ . (2.59) , ewe) m“ [3] = 451 1 a Red- 1___+0 _) 1503 2 11 For Re -) 0, it follows from Eq. (2.38) in the previous section that lim 1—2 -§(C"-3) (260) Re—)0 ll -8 D ’ . provided ’3/11 < on and dIl/dRe < .. to: Re -> 0. Therefore, with c; = 1.4, Eq.(2.60) implies that ’2 Km _ =_1. - (2.61) Re-)0 I1 40 Because 0 S QC 5 J2, a linear interpolation for 12/11 as a function of EC will be used between the two limits given by Eqs. (2.59) and (2.61). Therefore, I § _2 = _Z__l._£, (2.62) 11 3 3,]; whichcanalsobewrittenas I 2 = __2___1. _1‘_‘_. (2.63) 1—1 33105 It follows from Eqs. (2.47) and (2.50) that Il depends on the turbulent Reynolds number implicitly through EC and, consequently, .8. Thus, £311. _ [£32] £55 (264, I1 dRe s dRe 11 ds ' ' It follows from Eqs. (2.47) and (2.50) that 1 44(2). a...) Me) For Re —) co, Eq. (2.55) shows that B —) on; therefore, Q (13 fl 11 J dI1 1 lim —— = -—. 2.66 For Re-—)O, .8 ”1 1 lim —— = -—a. 2.67 R¢-)0[Il d9] 2 ( ) The parameter a remains as an adjustable universal parameter. A linear interpolation 41 similar to Eq. (2.62) is used for intermediate values of the Reynolds number. Therefore, .0 d1 1 l: gc] —— = -- 1+ [a-1]— , (2-63) I1 ds 2 J5 or, in terms of the viscous destruction coeflicient, 5 £111 1[ 14 _]O __ = —— 1+ 1 2.69 ,1 am 2 [a- 1 m < > The estimate of 13/11 follows from the approximations expressed by Eqs. (2.47) and (2.49) along with the integration cut-off distance given by Eq. (2.50). — — + ... I ('_"3 c 7?. .. 15R; , (2.70) 1 5c " 6.0:: + The limiting behavior of Eq. (2.70) is not apparent, so the lead term in the numerator rs factored and the remaining ratio is modeled as a universal constant: '1 Pb (2.71) J15Re3b The parameter b is a constant to be evaluated at the high Reynolds number asymptote. At ”NI N II infinite Reynolds numbers, Eq. (2.55) implies that (2.72) C" D 16 8 — - 2 = . . 2 15 15 (2 73) 42 This gives b = 0.0343 for C; = 1.83. 2.5 Parameter Estimates The asymptotic states (i.e. Re —> 0 and Re -) on) have been used to approximate the integral expressions found in Eq. (2.44). Eq. (2.44) may now be used to predict the intermediate behavior of CD (Re) . A parametric study of CD (Re) with respect to the empirical parameter in the representation for .8 (Re) will conclude this section. The selection of .8 (Re) may be judged as acceptable based on the resulting prediction of the velocity derivative skewness compared to experimental and siruulation data for S k (Re) . The form of 8 (Re) is specified by the following function: .8= (C;+7S;° J1?) (C; - Cg) exp [- (Re/Re. ) n] , (2.74) 35 where S: represents the high Reynolds number asymptote for the skewness. Re. and n are empirical parameters which control the transition for the destruction term between its high and low Reynolds number asymptotes. The parameters a, Re. and n are selected to reproduce the trend in St (Re) consistent with Eq. (2.44) and the foregoing approximations to I 1’ 12, and I3. The constants, which have been determined as a result of a parametric study, are a = -3.5, Re’ = 2.2 and n = 1.25. Figure 2.43 shows the behavior of S): (Re) predicted by Eq. (2.44). By design, the skewness assumes a value of 0.4 at large Reynolds numbers. Interesting is the prediction that a local maximum of St, "m = 0.52 occurs at Re ~ 5, which compares favorably with the data for the skewness: S k, ”3 ~ 0.5 to 0.6 at Re ~ 5 to 10 (Tavoularis er al., 1978; and Mansour and Wray, 1994). All asymptotic and empirical pararueters used to support this analysis are summarized in Table 2.1. Of the ten 43 0.7 0.6 - U L” D 3 D B D so ("3 r a L '— ' D U :1 0.5 - a 13.3 D I i3 0.4 - [:1 L] I S [i D k r CHI 93 0.3 1- I I 0.2 - I Reference D E] Tavoularis et al. [1978] 0.1 I Mansour and Wray [1992] 0.0 1 11111111 1 1111111l 1 11144111 1 1111111 lE-l 1E0 lEl lE2 1E3 Re Figure 2.4a Predictions of the IKH-Equation for the Skewness Using a Prescribed Function for the Destruction of Dissipation ( a = —3.5, Re 2 2.2 and n = 1.25). 44 Table 2.1: Summary of Parameters Used in the Analysis of the IKH-Equation Parameter Value Basis c; High Re data of Comte-Bellot and Corrsin [1972]; direct simulations of Mansour and Wray [1994] C; 1.4 Final decay period of Batchelor and Townsend [1948] s: 0.4 Data of Tavoularis er al. [1978]; direct simulations of Mansour and Wray [1994] (12/11) -1 Consistency with 6‘}, 0 (12/1 1) -2l3 Asymptotic behavior of series representations for '° integral quantities 3 (111 -ll2 Asymptotic behavior of series representations for [— —] integral quantifies II (1.8 a -3.5 Empirical parameter in (.8/11) (dIl/d.8) selected for consistency with data for skewness b 0.0343 Consistency with C3 Re. 2.2 Erupirical parameter in 3 (Re) selected for consistency with data for skewness n 1.25 Empirical parameter in .8 (Re) selected for consistency with data for skewness 45 parameters listed in Table 2.1, seven are derived directly from experimental data and/or from consistency conditions with the asymptotic behavior of various quantities. Only three quantities are left as free constants to reproduce the data for the velocity derivative skewness. The effects of the parameters a and Re. on the resulting skewness prediction are shown in Figures 2.4b and 2.4c. In each case, the parameters used in the computations are the same as those listed in Table 2.1, with the exception of the variable shown in the figures. Figure 2.4b indicates that the parameter Re. influences both the location (in terms of Re) and the magnitude of the local maximum in the skewness. Figure 2.4c shows that the parameter a also influences the magnitude of the local maximum in the skewness as well as the low Reynolds number behavior of the skewness. Figure 2.5 illustrates the semi-theoretical prediction of CD (Re) and Figure 2.6 shows the individual contributions of .8 (Re) and 0’ (Re) to the destruction-of-dissipation coefficient. The curve for .8 (Re) is the erupirical prescription described by Eq. (2.74). The curve for P (Re) is the resulting prediction of Eq. (2.44). Within the context of the modeled form of the IKH-equation, the local maximum in CD (Re) is a direct consequence of the transition of .8 (Re) between high and low Reynolds number behavior. For instance, if CZ — C; = 0 (with all other parameters being the same as given above), then C D (Re) is constant. As C; - C; increases, however, the non-monotonic behavior becomes more pronounced. At large Reynolds numbers, the growth of both production and destruction are pmportional to JR; and CD = .8 — P -9 C3 = 1.83. Given that the non-monotonic behavior of C D (Re) is supported by: (1) the analysis of the classical Karman-Howarth Equation; and, (2) the numerical siruulations of Mansour and Wray [1994], the following empirical representation of CD (Re) is proposed as an improvement on the previous expressions developed in the literature (see Figure 2.2), C"+C“‘(R /R )’ , CD: D D e eA [1+C [exp(-p{ln(Re/Rep)2})]:|. (2.75) 1+(Re/ReA)‘ 46 0.7 00 4 44111111 1 11111111 1 1 1111111 1 1 lllLll lE-l 1E0 lEl IEZ IE3 Re Figure 2.4b The Effect of Re“ on the Velocity Derivative Skewness (a = —3.5 and n = 1.25). 47 0.6 0.5 *- 0.4 - 0.3 - -2.0 0.2 - 0.11- 00 111111111 1 11111111 1 11111111 11111111 113-1 150 ”31 182 1E3 Re Figure 2.4c The Efi‘ect of a on the Velocity Derivative Skewness (Re. = 2.2 and n = 1.25). 48 1.9 Eq. (2.75) Eq. (2.44), IKH-equation 14 111111111 lllLLLlll 111111141 11 111111 lE-l 1E0 1E1 1E2 1E3 Re Figure 2.5 Predictions of the IKH-Equation for the Destruction of Dissipation Coefficient 49 / Asymptotic grow h oszandfl (~ Re) 01 1111111l 11111111l 11111111] 11111111 [El ”30 113] 152 1E3 Re Figure 2.6 Predictions of the IKH-Equation for the Production of Dissipation Using a Prescribed Function for the Destruction of Dissipation. 50 In Eq. (2.75), Re A determines the Reynolds number at which the decay coefi'icient approaches its high Reynolds number asymptote and s is a sharpness parameter for this approach. Equivalently, Re P determines where CD (Re) achieves its local peak and p relates to the sharpness of this peak. The quantity C. determines the value of CD at the local maximum. The following parameters have been selected in order to reproduce the local maximum and minimum in CD (Re) illustrated by Figure 2.5: Re A = 5, s = 0.4, Re}, = 1.5, p = 0.65,and C' = 0.054. 2.6 Transient Isotropic Decay With CD (Re) expressed by Eq. (2.75), Eqs. (2.2) and (2.5) may be employed to compute an isotropic decay process for any arbitrary initial conditions given by 1:0, 80, and Re 0 (kg/veg). These equations have been integrated using a fourth order Runge- Kutta integration algorithm (Carnahan et al., 1969; see also Appendix H). For an appropriate Specification of the initial conditions, this calculation can be used to reproduce the isotropic states which correspond to the experiments of Batchelor and Townsend [1948], Comte-Ballot and Cousin [1971], and Sirivat and Warhaft [1983] (see Figure 2.1). Figures 2.7a-2.7g present the results for the seven experimental data sets used in this chapter. The results are given as a crossplot of the turbulent kinetic energy vs. the dissipation. Eqs. (2.2) and (2.5) imply that . (2.58) 84% min- _1- CD The local slope of the computed line compared to the slope of the data points gives an indication as to the quality of the specification of C D (Re) . In general, the reproduction of 51 0.1 ' 0.01 Figure 2.7a Model Computations for the Isotropic Decay Data of Batchelor and Townsend [1948] (U = 150 cm/s, Reo = 7.58) k/ko 01 1 11111111 1 11111111 0.01 0.1 1 Figure 2.7b Model Computations for the Isotropic Decay Data of Batchelor and Townsend [1948] (U = 643 cm/s, Rea = 42.6) 53 01 1 11111111 1 11111111 0.01 0.1 1 Figure 2.7c Model Computations for the Isotropic Decay Data of Batchelor and Townsend [1948] (U = 1286 cm/s, Reo = 83.7) 00] 111114111 1 11111111 1 11111111 0.001 0.01 0.1 1 who Figure 2.7d Model Computations for the Isotropic Decay Data of Comte-Bellot and Corrsin [1971] (U= 10 m/s, Rea = 354) 55 1 1111111 1 11111111 0.1 L 0.01 0.1 1 Figure 2.7c Model Computations for the Isotropic Decay Data of Comte—Bellot and Corrsin [1971] (U= 10 m/s, Rea = 769) 56 01 l 1 1111111 I 1 1111111 0.01 0.1 l Figure 2.7f Model Computations for the Isotropic Decay Data of Sirivat and Warhaft [1983] (U= 340 cm/s, Rea = 139) 57 01 1 1 1111111 1 11111111 0.01 0.1 1 Figure 2.7g Model Computations for the Isotropic Decay Data of Sirivat and Warhafi [1983] (U = 630 cm/s, Rea = 262) 58 the data sets is quite good. Specifically, the initial decay stages of each of the seven data sets is well represented, although some sets show deviation from the computations at points farther along in the decay process. This is not surprising, however, as the degree of homogeneity and isotropy of the flow worsens at larger downstream positions in the wind tunnel (Comte-Bellot and Corrsin, 1971). Of particular note are the two cases which represent the extrema for the experimental Reynolds numbers. Figure 2.7c represents the highest Reynolds number data of Comte-Bellot and Cousin [1971] (Re a = 769). The slope of the calculation is slightly higher than that of the data, indicating that CD is relatively low, and perhaps approaches its asymptote too slowly. Figure 2.7a represents the low Reynolds number data of Batchelor and Townsend [1948] (Re a = 7.58). Here, the representation of the data is excellent throughout the entire decay process. The good agreement of the data sets which span two decades of Re a is an indication that Eq. (2.75) is a quantitatively good expression for CD. Figures 2.7h and 2.7i present the results of the same calculations, but compared with the direct simulations of Bardina et al. [1985] and Speziale et al. [1987]. The solid line shows the results of the computations using the data at tea/k0 = 0 as the initial conditions, whereas the dashed line shows the results using the simulation data at tea/k0 = 1 as the initial state for the solution of Eqs. (2.2) and (2.5). What is most noteworthy about these two simulations is the fact that the initial decay period is very poorly characterized, while the final decay is well reproduced. This, however, does not necessarily mean that the model for CD (Re) is poor. Rather, the limitations of direct simulations are brought into question. Huang and Leonard [1994] address this problem. They indicate that, during the initial stages of a direct simulation, the maximum resolvable wavenumber 75m“ of the turbulent energy spectrum (E (12)) is relatively small. Thus, although the larger, energy-containing eddies (i.e. low wavenumber I?) may be well characterized, the smaller, dissipative eddies (i.e. high wavenumber I?) are not. Without an accurate characterization of both the energy-containing and dissipative eddies, the decay 107cc 0.1 0.01 Figure 2.7h Model Computations for the Isotropic Decay Simulations of Speziale, et al. [1987] (Reg = 35.1) 60 0.1 0.01 Figure 2.7i Model Computations for the Isotropic Decay Simulations of Bardina, et al. [1985] (Reg 2 45.4) 61 transient is clearly suspect. Huang and Leonard explain, however, that in,“ increases during the course of the simulation. Therefore, at longer times, both types of eddies are accurately represented and the quality of the decay transient is much more reliable. Figure 2.8 presents the decay of the turbulent kinetic energy in the time domain for three different initial Reynolds numbers: Re a = 10, 100, 1000. The kinetic energy normalized by its initial value is plotted versus the dimensionless decay time (‘t = re 0/ k a). It is noted that the small differences among the three cases are manifested only in the long time behavior. The curves follow the trends expected from the nature of Eqs. (2.2) and (2.5): At higher Reynolds numbers, C D is higher and causes the turbulent dissipation to decay more rapidly, causing the turbulent kinetic energy to persist longer. 2.7 Conclusions The integral form of the KArmdn-Howarth equation for isotropic turbulence provides a semi-theoretical prediction for the destruction coefficient CD in the equation for the turbulent scalar dissipation. Integral expressions containing the double and triple velocity correlations are estimated using Taylor series representations of the velocity correlations and a characteristic cut-ofl' integration length. A modeling hypothesis is made which assumes that the ratio of the integral quantities is correctly represented, although the absolute magnitudes of the integrals may not be. Empirical parameters which are chosen to reproduce experimental data for the velocity derivative skewness in decaying isotropic turbulence predict that CD spans its upper and lower limits in a non-monotonic fashion. This is unlike previous approaches to this problem, which have simply bridged the two limiting cases with some monotonic, empirical function. The non-monotonic behavior CD (Re) is seen in some direct simulations of this flow. The resulting prediction for CD (Re) is able to reconstruct the 1.0 k/ko Re 2 1()()() 0 100 10 O] l 1 1 A l 1 1 1 I 1 0.0 1.0 2.0 3.0 4.0 5.0 Figure 2.8 Predicted Transient Decay of the Turbulent Kinetic Energy as a Function of Initial Turbulent Reynolds Number 63 various experimental and direct simulation data available in the literature. Within the context of the modeled form of the IKH equation, the non-m onotonicity of CD (Re) is a direct consequence of the fact that CD varies between two different upper and lower Reynolds number limits. For C3- (3; = 0, CD (Re) would be a constant function, while the local maximum in CD (Re) becomes more pronounced as C; - C; increases. The disparity between the calculations and the initial stages of decay in direct simulations illustrates that care must be taken when treating direct simulation data. Only the longer-time simulation data resolve both the large and small turbulent scales and provide meaningful data for comparison. CHAP'I'ER3 ISOTROPIC PRES-STRESS CLOSURE THEORY FOR HOMOGENEOUSLY SHEARED TURBULENCE 3.1 Introduction Homogeneously sheared turbulence has the interesting feature that the turbulent time scale k/e approaches a constant finite value as the flow develops (Tavoularis and Cousin, 1981; Gibson and Kanellopoulos, 1987; Rohr er al., 1988; and, Tavoularis and Kamik, 1989), although the kinetic energy I: and the dissipation 2 associated with the velocity fluctuations continuously increase as 2 —) co. Figure 3.1 slunmarizes the asymptotic state attained by this flow for which the shear rate S is approximately constant. The existence of this state requires that the transport equations governing the turbulent kinetic energy and the scalar dissipation rate satisfy the following model-independent condition ldk 1dr: lim -— = lim -— . 3.1 z-r~|:kdz] z-r-oLtdz] ( ) For homogeneous shear flows at high turbulence Reynolds numbers (It2 » vs), the k-e equations (Hanjalic and Launder, 1972) are (emu-35 = —:V—e. (3.2) 2 do (2.2.)‘VQQ a (uz) (0)2; - -CP 1P -CD-‘E;. (3.3) (uz) (0) is the axial component of the mean velocity along the centerline of the flow field. IP and 1'0 represent, respectively, turbulent time scales associated with the production of 64 :2 zS/ (0) D Q 22 ............................................................... (“2) (O) J d=V _ Cn-l akin-i e l'cp-l’ Eqs. (3.2) and (3.3) are also consistent with experimental data for homogeneous decay provided CD = 1.83 (see, esp., Comte-Bellot and Corrsin, 1971; Mansour and Wray, 1994; and Chapter 2 of this dissertation). The Reynolds stress can be written as the sum of an isotropic and an anisotr0pic stress (Speziale, 1991) '( 't_r') = 3359211), (3.5) where b = 97 and tr (2) = 0. For a strictly homogeneous flow with S = 0, the anisotropic stress 2 is zero. Because (r_t'r_r') is a non-negative operator, the eigenvalues associated with the Reynolds stress are non-negative (Schumann, 1977). Moreover, because tr(r_e'_u') = 2k>0, (3.6) at least one eigenvalue of (r_¢'_r_r') must be non-zero and positive. This means that the normaliwd components of the turbulent kinetic energy are confined to the energy simplex illustrated by Figure 3.2. The homogeneous shear data (Tavoularis and Kamik, 1989) summarized by Frgure 3.2 show that this flow produces a positive primary normal stress difference, RR 67 00:23.5... eoeeocm bmzoocowofiom .8 $35 5228:. .23 538% Queen NM 23w:— : .o .8 8 v ”V N ”museum wine—05D mo =o_woM|/A T... 3% oiobofl 8 ._ .8 u meanness 68 lim [(u'zu'z) - (u'yu'y) a 0.377, 2k 1 as well as a negative second nounal stress difference, lim [(uw’) — “in"? 2-0040. z -r o. 2]: The anisotropic stress b has two non-trivial invariants: II = tr (_b - _b) and 111 = tr (_b - b - b) . Lumley [1978] has shown that all realizable turbulent Sta-tes—(II, 11]) must fall .on_ or within a two dimensional domain illustrated by Figure 3.3. The experimental data for homogeneous shear approach the asymptotic state given by 111m” (”,1”) a (0.138, 0.0174) . For a constant density fluid, the Boussinesq model for the anisotropic stress is (Hinze, 1959 and Brodkey, 1967) 21.5 = -ve[V(g)+V(t_e)r] = -2v¢(__s). (3.7) Q) represents the mean strain rate dyadic and ve is a scalar valued eddy viscosity. For Ve > 0, Eq. (3.7) implies that the kinetic energy is irreversibly transferred from the mean field to the fluctuating field inasmuch as —(_u'_u'):V(u) = 2ve(§): (.3) >0 (3.8) for all turbulent flows. This feature partly justifies the use of Eq. (3.7) as an approximate model for the anisotropic stress. However, for homogeneous shear flows, Eq. (3.7) also predicts an equipartition of kinetic energy among the components of the fluctuating velocity as well as a zero third invariant for b (i.e., ”150). These results clearly contradict the experimental measurements summarized by Figures 3.2 and 3.3. Thus, as previously noted by Speziale [1991], the Boussinesq model qualitatively misrepresents the underlying mechanism associated with the flux of momentum due to velocity fluctuations. 111 = —6(11/6)3/2 69 111 = 6(11/6)3/2 Asymptotic State: (0.0174, 0.138) Region of Developing States: 2 < g < 26 III Isotropic State: (0, 0) Figure 3.3 Anisotropy Invariant Diagram with Transition States for Homogeneously Sheared Turbulence 70 The proposed pre-closure representation of the Reynolds stress supports a long- standing goal of turbulence research to achieve a practical statistical closure of the mean field equations and, thereby, compleruents other algebraic turbulent closure models for the anisotr0pic stress _b (see, esp., Speziale, 1991; Taulbee, 1989; Reynolds 1989; and, Hanjalic, 1994). The development of homogeneously sheared turbulent flows towards an asymptotic state provides a critical experimental test flow to partially guide the evaluation of the proposed theory. The objective of this chapter is to develop this new pre-closure for the Reynolds stress. In Section 3.2, an analysis of the goveming equations for the fluctuating velocity yields a relationship between the Reynolds stress, mean field quantities, and the pre-stress (see Eq. (1.37)). A simplifying assumption is made regarding the isotropic nature of the pre-stress in Section 3.3, which is found to directly impact Issues (i) and (ii) from Section 1.2. The asymptotic state of homogeneously sheared turbulent flows is used to calibrate universal model parameters and the transient approach to this asymptote is examined. 3.2 Pre-closure Theory The fluctuating velocity associated with a constant density, Newtonian fluid satisfies the following equation (Monin and Yaglom, 1965) <00») = 12' . (3.9) where 11"“! ° V+f . (3.10) and 71 f'a V . [%!+r_e'r_4' - (1_e't_¢'):l . (3.11) The convective-diffusive operator (2) depends on the mean velocity field (3) and the kinematic viscosity of the fluid, v: (.1) gait). (10' v -vV2. (3.12) Eq. (3.9) is exact and emphasizes that fluctuations in momentum (i.e. pg') are produced within the flow domain by (l) a convective coupling between the mean velocity gradient and the fluctuating velocity; (2) pressure fluctuations; and, (3) fluctuations in the instantaneous Reynolds stress, g'g' - (u'u'). Momentum fluctuations are transported by viscous fluctuations in the molecular stress and by mean field convection. An exact, albeit fouual, representation of the flucurating velocity can be written in terms of a Green’s function (Morse and Feshbach, 1953) associated with the linear differential operator (2). For statistically stationary flows in an unbounded domain, ‘ t t_t' (a. t) = - j diIdWG) (f, t [13,?) l." (3.1). (3.13) .... v For 0 S t -i at "it -§E"2/v, the Greens’ function is spatially peaked in a frame of reference moving with the local mean velocity; however, as t -3 -) no, viscous momentum transport causes the Green’s function to relax to zero over the entire spatial domain. For an unbounded spatial domain, the Green’s function satisfies the following integml property jdi/(G) (a. t [3,?) = 1. ' (3.14) v The analog of Eq. (3.13) for a passive scalar field has been previously used by Hill and Petty [1996] and many others. A fouual representation for the Reynolds stress follows by either pre- or post- multiplying Eq. (3.13) by the fluctuating velocity u' (it, t) and then founing the enseruble 72 average to obtain 1 (314') = - I dildWG) (a. t pg?) (g' (gt. 0 12' (52. 1)) = _. V t - I 4214mm: pg?) (11' (33):; (r. t) >. (3.15) ... v Eq. (3.15) relates the Reynolds stress to the Green’s function and to the symmetric space- time conelation def‘ med by (513:. 0511523)): (2' (5. t) 5' (5.?) )- WEN- (£15. t)['(§.?)) = .. _ r - _ _ - VQI) ~(t_¢'(§. 014' (3. OH (Hg. 0 55' (3. t) )- (3.16) All the components of (51' (5.1-)1! (it, t)) relax to zero for either It-ilrtH or '15—?" » I”, where I” (~k/e) and 1H (~k3/2/e) represent finite turbulent time and length scales, respectively. The physical hypothesis that I” <6. and tH< co partially motivates the use of a spatial smoothing approximation to simplify the non-local representation of the Reynolds stress given by Eq. (3.15). A spatial Taylor series expansion of 11' (S, 1) about (it, 3) gives it. (is?) = é.(§’;) + (5-5) ° vb. (59;) +551?) (5—3) :vvy (5,?) + (3.17) Inserting Eq. (3.17) and Eq. (3.14) into Eq. (3.15) gives the following representation for the Reynolds stress in terms of the spatial moments of the Green’s function t fl 0314') = - I at (Lt'(§.'i)l_¢'(§.t))+ 2; 53‘” (am-3) (3.18) -¢° i=1 73 where 53‘" =- [ Idiot—5t) (G) (a. t lice] - (V 5' (you (gt. t)) (3.19) V and 5 ‘2’ a g [Idiot—3) (5 -§) (G) (at. t I30]: (We (55);; (s. t) >. (3.20) v — The temporal correlation 1.1 (n) (it, t-?) involves the nth spatial moment of the Green’s function contracted with the nth order gradient of the fluctuating field 5'. For t = ‘i, all of the spatial moments of (G) (r_r, t [13,?) are zero; however, these moments become non- zero on a time scale associated with the viscous transport of momentum. The foregoing expansion of Eq. (3.15) exploits the ides that the Green’s function (G) (3:, t '3, 3) acts like a spatial delta distribution on a time scale for which turbulent correlations become uncorrelated. For turbulent flows at large Reynolds numbers (i.e. 1H « 111/ v), the spatial moments of the Green’s function are assumed to remain small over the finite time scale for which turbulent fluctuations are temporally couelated. This hypothesis motivates the use of a spatial smoothing approximation which reduces Eq. (3.18) to the following statistically stationary approximation (g'g' s-l (1v (5.?)3' (5. t) >41 = -J (515.0? (grinds. (3.21) 0 0 where ‘t E t-i. Eq. (3.21) gives a representation of the Reynolds stress in terms of the teruporal history of the turbulence. As previously noted, the autocouelation (_l_r' (3,?)14' (it, t)) must be a symmetric dyadic-valued '0perator. This follows directly from the fundamental representation of r_t' in terms of [1' , the source of momentum fluctuations. 74 Eq. (3.21) can be simplified further by assuming that there exists a scalar-valued memory function m (5,1) with a finite cut-ofl' time such that (if (13.?) g' (5. t) ) = (11' (at. t) g' (at. t) )m (gr. I) . (3.22) Thus, Eq. (3.21) can be represented by (3'27 = 4p (1114') = 43051:). (3.23) where the relaxation time In is a phenomenological coefficient dependent on the temporal structure of the turbulence: “CR 5 1m (x, 1:) d1 . (3.24) 0 For statistically stationary flows, the integral time scale TR depends on the local statistical properties of the turbulence k and e, the mean velocity gradient (S a I] V (_u)||), and the viscosity of the fluid. Dimensional reasoning suggests that 't = C S , (3.25) where CR is a dimensionless function of kS/e and the turbulent Reynolds number kz/ve. It follows from Eq. (3.10) that Eq. (3.23) can be written as (2'!) = -1nQ!'>-tnvr- (2'27 = -tn<2'1'?-tn<2'e')-V- (326) Algebraic equations for the turbulent correlations ([14’) and (u'f') follow by pre- and post-multiplying Eq. (3.13) by ['(x, t) and taking the ensemble average. The previously anployed spatial smoothing approximation for large turbulent Reynolds numbers and the use of the memory ansatz expressed by Eq. (3.22) yields the following result for the 75 symmetric correlation (If): we = «122-«Ree» V = -<2'[HRVT- (err. (327) In the above expression, IR is assumed to be the same relaxation time as defined by Eqs. (3.24) and (3.25). The idea that a single scalar relaxation function can be used to characterize the temporal structure of all low-order, space-time correlations is a significant unifying step towards a practical closure theory. Eqs. (3.26) and (3.27) can be combined to yield the following pre-closure representation of the Reynolds stress BHRVQ‘UT' (Leg). [.IHRVQOJ = rim"). (3.2s) Eq. (3.28) is called the pre-closure hereinafter because it relates the Reynolds stress to three aspects of the flow field: (1) the spatial gradient of the mean field; (2) the relaxation time TR; and, (3) the unclosed turbulent pre-stress, fig?) . For small dimensionless relaxation times: NR 31R" VQr)" « 1, (3.29) the pre—stress approaches the Reynolds stress. Like the Reynolds stress, the pre-stress is a non-negative operator, and, thereby, has only non-negative eigenvalues. Eq. (3.11) could be used to founally relate the turbulent pre-stress to other unknown statistical correlations involving spatial gradients of the fluctuating pressure and the divergence of fluctuations in the instantaneous Reynolds stress. Alternatively, a phenomenological theory for the pre-stress could be deveIOped analogous to Eq. (3.5) by first writing 112, (ff) as the sum of an isotropic pre-stress and an anisotropic pre-stress 1,2, ([f) a 2&1” 2kg , (3.30) 76 Like the anisotropic stress a, the anisotropic pre-stress g is both symmetric (g = LIT) and traceless (tr (H) = 0). The isotropic pre-stress coefficient a must satisfy the following normalization condition ring) = 2d = 2k+ 21,, (Leg): V(y) 4r?R [tr (V(g)T- (_u'u')- V(_r_r))]. (3.31) As the relaxation group N R increases, the isotrOpic and anisotropic parts of the pre-stress should make significant contributions to the distribution of kinetic energy among the components of the fluctuating velocity as well as to the shear components of the Reynolds stress. Eqs. (3.28), (3.30), and (3.31) provide an alternative closure strategy to previously developed approaches based on the anisotropic part of the Reynolds stress. In this paper, the efficacy of using Eq. (3.28) with (332) ll m 1| 11 Q will be explored. This assumption provides a closure model for the normalized Reynolds stress < '5') 2 (3.33) Hz ae- A non-trivial model for L! is presented in Chapter 4 of this work. Weispfennig [1997] has also considered non-trivial models for the anisotropic pre-stress. 3.3 Isotropic Pre-stress Theory For homogeneous shear and £1 = 9, Eqs. (3.28)-(3.30) reduce to the following set of equations for the components of the normalized Reynolds stress l Rxx = 3 + NR2 ’ (3.34) 1 W = 2 ’ (3'35) 3+NR 1+NR2 Rzz = 2 , (3.36) 3 + NR and -N R yz= 2. ea) 3 + NR The relaxation group N R is defined by N R I IRS. The isotropic pre-stress coefl‘icient a is given by 3 2. 3+NR 0t — = 3.38 k ( ) The above set of equations, hereinafter teuned the isotropic pre-stress (IPS-) theory, approaches the Boussinesq theory (see Eq. (3.7)) for NR « 1 inasmuch as r r . . . . Rn = Ryy = Rzz = 3, Ryz = -5NR, and a = k. The eddy vrscosrty coefficient forthrs limiting case is the same as the standard k-e theory, viz... 2 2 1: lim v = -C —. 3.39 NR—) 0 ‘ 3 R e ( ) For N R » l, the IFS-theory shifts the kinetic energy to the axial component of the fluctuating velocity with the result that Rn = a” = l/NR2 and Ru -+ 1. Eq. (3.37) predicts that the shear component of the Reynolds stress becomes inversely proportional to the relaxation number, Ryz = -1/NR, and Eq. (3.38) gives an: = 3m),2 for NR » 1. 78 Eq. (3.37) shows that a maximum in —1ryz (i.e. 1/./l_2) occurs at NR = J3. Furthermore, Eqs. (3.25) and (3.37) together with NR = IRS imply that . . 9404,) 40 -(uyuz) €11: dy (1) where 2C c = R . (3.41) 11 3+le2 Thus, the IFS-theory predicts that the eddy viscosity coefficient cu decreases as the relaxation group increases, provided CR is modeled as a universal constant. This non- linear dependence on the mean shear field is fundamentally different than the behavior presumed by the k-e theory, which employs a constant eddy viscosity coefficient. The IFS-theory and the Boussinesq theory also predict qualitatively different results for the invariants of the anisotropic stress. For homogeneous shear, 2 2 2 2 II = bu + by, + bzz + 2byz (3.42) and _r' 3 3 3 2 III — b1”r + by, + bzz-3bnbyz. (3.43) For the Boussinesq theory, Eqs. (3.7) and (3.39) imply that but = by), = bzz = 0, and byz = -NR/3. (3.44) Therefore, for this special case, Eqs. (3.42) and (3.43) yield II 2 2 B = 5117,, and 1118 = 0. (3.45) 79 Because the third invariant is zero, it follows directly from Figure 3.3 and Eq. (3.45) that realizable Boussinesq turbulent states are restricted to N R e [0, l] . It follows directly from Eqs. (3.5) and (3.34)-(3.37) that the IFS-theory yields the following results for the components of the anisotropic stress 2 NR bu = — 2 , (3.46) 3 (3 + NR ) 2 NR b” = 2 , (3.47) 3 (3 + N R ) bZZ = —2byy , (3.48) and -N fl = R 2. (3.49) 3 + NR The invariants II and III for this theory can be written as 2 II = — (3.50) [PS 2 3 3 + NR and 9 + 21th2 NR2 3 III [PS = 2 2 . (3.51) 9NR 3 +NR As the relaxation number goes to zero, Eqs. (3.50) and (3.51) give the isotropic pair (”,1”) = (0,0); as NR-Mo, (”,1”) = (2/3,2/9); and, for NR = 1, (II, III) = (1/6, 11/576) . Figure 3.4 shows the locus of realizable homogeneous shear states predicted by the IPS-theory for 0 S N R S co. The asymptotic state measured by 8O Locus Qf Invariant Domain [PS-States Boundary Asymptotic State: NR = 1 (0.0174. 0.138) Boussinesq Regime III Figure 3.4 Anisotropy Invariants Predicted by the [PS-Theory for Homogeneously Sheared Turbulence 81 Tavoularis and Kamik [1989] is also shown on Figure 3.4. Clearly, the IFS-theory represents a significant improvement over the classical algebraic Boussinesq theory for which III B = 0 for all homogeneous shear states (ch Figure 1.2). The spatial development in the relaxation group is determined by the transport equations for k and e. Eqs. (3.2) and (3.3) can be combined to yield the following non- linear, ordinary differential equation for N R dNR 71—? = 2NRRyz(CP- 1) + CR(CD— 1). (352) where the dimensionless development time g is defined by § = 25/ (:41) (0) . (3.53) and Ryz is given by Eq. (37). As é approaches infinity, Eq. (3.52) reduces to the asymptotic condition given by Eq.(3.4). ' 3.4 Parameter Estimates Although Figure 3.4 shows that the IFS-theory gives an improved prediction of the anisotropic invariants compared to the Boussinesq theory, the experimentally observed asymptotic state is nevertheless unattainable. Figure 3.5 illustrates a selection strategy for the asymptotic value of N R which minimizes the relative error between the asymptotic experimental state (112.1112) and the locus of attainable states consistent with the IPS- theory. Therefore, with ”-11: 2 III-III: 2 m A={ .H—A . In 11“ N; is selected to minimize A (see figure 3.5). This procedure yields N; = 0.945, ooeoifia. BEBE EmsooeoonquonEfi/x e8 955 nausea—om of .8 :ozmsfimm Wm 83mm a. WEE: I 3 ed mo . . a . q d % . 0.0 :2 m._ a; m5 _ . . . a i . . 3. 28m EBESAQEEmEthRM t o.— eenmz m5 Q a .c:\: Q— = wfld n . .. 2 88m. 922%??? min: a N n «2 mmaeamfiafik named L m.— 83 R; = 12‘;y = 0.257, R; = 0.486, and -R;z = 0.243. With N; = CR(Sk/8)a, the relaxation coefficient CR = 0.227 for (Sk/e) a = 4.16, as reported by Tavoularis and Kamik [1989] (see Figure 3.1). The existence condition for an asymptotic state (see Eqs.(3.4) and (3.52)) provides an additional equation for the model coefficient C P: CR(CD- 1) C=1 P (3.55) a 211/fenyz Eq. (3.55)implies that c, = 1.41 for CD = 1.83 and Nine:z = -o.229. Table 3.1 gives a summary of the model parameters for the IPS-theory. Table 3.2 gives the predicted properties of the asymptotic state for homogeneously sheared flows. 3.5 Results and Discussion Eq. (3.52) with R yz given by Eq. (3.37) was solved numerically using a fourth order Runge-Kutta algorithm (Carnahan, er al., 1969; see also Appendix H) with a variable time step. The three model parameters C1,, CD, and CR are assumed to be universal constants independent of Sk/e. Figure 3.6 shows that the relaxation group N R ( 5 1R5) monotonically approaches its asymptotic limit for a wide range of initial states (0 S N; S 10). The development time for N R a N; depends on the initial conditions. For an isotropic initial condition (i.e. N; = 0), the asymptotic state is approached for developmental times on the order of ten. However, for highly anisotropic initial conditions (N; = 10), the approach to the asymptote requires § ~ 30. Figure 3.7 shows the response of the normal components of the normalized Reynolds stress subjected to an isotropic initial state (5 = {/3 and N; = 0). Like the relaxation group, an asymptotic state is attained for §-10 to 30, depending on the initial turbulent state. This result is comparable to the development times observed experimentally by Tavoularis and Cousin [1981] and is sensitive to the initial state of turbulence. The 84 Thble 3.1: Parameter Estimates for the IPS-Theory Parameter Estimate Basis TT Isotropic Decay (Comte-Bellot and Corrsin, 1971; Mansour and Wray, 1994; see Chapter 2) C}: 1.41 Existence condition for k-e equations CR 0.227 Optimization of anisotropy invariants Table 3.2: Predictions of Asymptotic Statistical Properties for 85 Homogeneously Sheared Turbulence Property [PS-Theory Experimental NR I 0.945 WAT kS/e NIAl 4.16 Rx: 0.257 0.236 R», 0.257 0.196 Ra 0.486 0.568 -R),z 0.243 0.165 II 0.153 0.138 III 0.0162 0.0174 * Data point not applicable to this entry. 86 10 ; . 5.0 10.0 t 3.5 1.5 1 r ’ ' Asymptotic State N; = 0.945 0.50 NR 0.15 0.1 f b N; = 00 001 I 11111111 1 1 11111Ll 1 1 11411” 0.1 1 10 100 Figure 3.6 The Effect of the Development Time on the Relaxation Group for Homogeneously Sheared Turbulence 87 0.6 r 00 1 1 1111111 1 1 IILUL] 1 10 100 Figure 3.7 The Efi'ect of the Development Time on the Distribution of Kinetic Energy for Homogeneously Sheared Turbulence ( — [PS-theory; experimental data (Tavoularis and Kamik, 1989): D R”; A Ryy; 0 R22) 88 asymptotic state for the IPS-theory is as close to the data as allowed by the optimization of the model parameter CR (see Figure 3.5). For an anisotropic pre-stress (g at (:1), the invariants of the Reynolds stress anisotropy can be reproduced exactly (see Chapter 4). The IFS-theory clearly does not predict a second normal stress difl‘erence, and the predicted primary normal stress difference, R; — R3), = 0.229, is significantly smaller than the experimental estimate of 0.372 (see Figure 3.2). These predictions, however, can easily be improved by using a phenomenological theory with a non-trivial anisotropic pre- stress. For instance, it follows directly from Eqs. (3.28), (3.29), and (3.30) that 10: R11 = g-k-+Hxx , (3.56) R -1“ H 357 yy-EZ+ yy’ (' ) and R —1°‘ H 2NR 1112 358 7.2-33+ zz- Ryz_ RR)?" (°) Eqs. (3.56) and (3.57) stern directly from the pre-closure and imply that an anisotropic pro-stress is a necessary condition for a non-zero second normal stress difference. Eqs. (3.57) and (3.58) also show that the primary normal stress difi'erence arises from two distinct physical efl’ects: (l) the primary normal difference sz—Hyy; and, (2) the convective coupling between the transverse velocity fluctuations and the mean field gradient. d(uz) L‘" VQ‘) = u'y dy fz' (3.59) This fluctuating field is directly responsible for the two statistical correlations (u'zu'y)S and (u'yu'y)S, which, as indicated above, make the Reynolds stress anisotropic even if the pre-stress is isotropic. 89 The normalization of the pre-stress (Eq. (3.31)) gives the following expression for (it/k (see Eq. (3.31)) a _ 2 .1; - 1+2NRRyz+NRRyy. (3.60) For an isotropic pre-stless, Ryz = —NRR”, and Eq. (3.60) reduces to a — 1 N 361 7" " + RRyz’ ( ' ) With this result, Eqs. (3.56)-(3.58) can be re-written in canonical form: Rn = %(1-@) , (3.62) 1 R” = 3(1-0) , (3.63) and 1 R22 = §(1+20), (3.64) where was N - (u'zu'y)S _. RRyz _ 4R 2k , (3.65) It follows from Eq. (3.25) that the above expression for 12 can be rewritten as - C R (u'zu'y)S 28 ° (3.66) Thus, the ratio of production to dissipation of turbulent kinetic energy (see Eq. (3.2)) determines the redistribution of energy among the velocity fluctuations. It follows from 90 Eq. (3.4) (see also Eq. (3.55) and Table 3.1) that CR CD—l lim 0 = — = 0.229. (3.67) Eqs. (3.62)-(3.64) predict a non-zero primary normal stress difference and a zero second normal stress difference: R -R =0 , R —R =0. (3.68) For the IFS-theory, Figure 3.7 and Eqs. (3.62)—(3.64) show how a redisuibutes the energy produced by the coupling between the shear component of the Reynolds stress and the mean gradient. The experimental data of Tavoularis and Kamik [1989] show the same qualitative trend as the theory but, clearly, R yy - R xx :9 0. Figure 3.8 illustrates the possibility that the shear component of the Reynolds stress may not approach its asymptotic state monotonically. The figure also shows that -Ry'z has a maximum value at NR = J3. Because N; = 0.945 < J3, the monotonic behavior of the relaxation group towards its asymptote causes the shear stress to relax monotonically for initial states characterized by N; < J3. On the other hand, for N; > J3, -Ryz increases to a local maximum and then decreases to its asymptotic value for development times greater than ten. If N; > 2.7, the transient shear component of the Reynolds stress overshoots its asymptotic state at some finite time (see Figure 3.8). This complex transient response may have the appearance of a quasi-asymptotic condition, but Figure 3.8 shows that R yz requires § > 30 to attain its ultimate asymptotic state for highly anisotropic initial conditions. The maximum in —Ryz occurs because of two competing physical processes. For small values of N R, the IFS-theory approaches the Boussinesq (or gradient) transport regime for which (u'yu'z) ~ S; for large values of the relaxation group, the shear component of the Reynolds stress approaches the so—called equilibrium transport regime: 91 0.31’ b d Asymptotic State 0.2 '- C -Ryz 0'31 b d 0.1 - as» c ”Ryz L a ,,‘. P a 00 2 NR 4L A 6 O 1 11111111 1 lllllLLl 1 11111111 01 1 10 100 Figure 3.8 The Effect of the Development Time on the Shear Component of the Normalized Reynolds tress for Homogeneously Sheared Turbulence (a: N; = 0.2;b: N; = 3;c: N; = 5.0;d: N; = 0.945, —R‘;Z = 0.243) 92 2k 28 lim . . =__=__. 3.69 Sun yu Z) ”R CRS ( ) Eq. (3.69) shows that, in the equilibrium regime, the production of turbulent energy, -(u'yu'z)S, is proportional to the dissipation of turbulent energy. This non-linear transport regime arises because the energy in the transverse fluctuations (u'yu'y) decreases significantly as N R increases due to the redistribution of energy by a: (see Eqs. (3.63) and (3.68)). Thus, the [PS-theory for homogeneously sheared turbulent flows yields a shear- thinning eddy viscosity (see Eqs. (3.40) and (3.41)) which bridges the gradient transport regime (N R -) 0) with the equilibrium transport regime (N R —-) no): -(u'yu'z) 2CR k2 ”W = W“? = Mi? (3.70) For some simple shear flows, the relative time scale N R may span a wide range of values; Therefore, both transport regimes may occur in the same flow field, albeit at difl'erent spatial locations. 3.6 Conclusions As a direct consequence of the pre-closure theory given by Eq. (28), a non-zero primary nounal stress difference obtains regardless of the closure hypothesis for the pre- stress; however, a non-zero second normal stress difference requires the existence of a non-trivial anisotropic pre-stress. For positive I: and e, the IFS-theory predicts realizable turbulent states for 0 S N R S co, whereas realizable states for the Boussinesq theory occur only for OSNRS 1. This theoretical result obtains for the IPS-model because the pre-closure relates the Reynolds stress to a pre-stress having only non-negative eigenvalues. For the special case of an isotropic pre-stress, all the eigenvalues of the pre-stress are all equal to a/3k (> 0). The IFS-theory with CR interpreted as a universal constant predicts two distinct 93 transport regimes bridged by an effective eddy viscosity which depends on tRS. For “IRS « l, a gradient transport regime occurs (i.e. (u'yu'z) ~S); whereas, for 1R5 » 1, an equilibrium transport regime occurs (i.e. (u'yu'z) ~ US). The time required for the turbulence to achieve an asymptotic state is strongly dependent on initial conditions. The developmental times needed to reach an asymptotic state based on the IPS-theory agree qualitatively with experimental observations. CHAPI'ER4 ANISOTROPIC PRE-STRESS THEORY FOR HOMOGENEOUSLY SHEARED TURBULENT FLOWS 4.1 Introduction In Chapter 3, an algebraic pre-closure theory for the Reynolds stress from an analysis of the equation of motion for statistically stationary turbulent flows was developed. A spatial smoothing approximation and the use of a memory ansatz for turbulent temporal correlations were key elements in the development of the following pre-closure theory for the Reynolds stress 542's»: = with . (4.1) where g -.-=. [3+tRV(g)] . (4.2) Eq. (4.1) relates the Reynolds stress to the gradient of the mean velocity, the relaxation time IR, and the turbulent pre-stress. For large turbulent Reynolds numbers (i.e. k2 » vs), the relaxation time is assumed to scale with the characteristic eddy tumover time k/s, 1: TR = CR; ’ (4.3) where It represents the kinetic energy of the turbulent fluctuations (2k 501' . 3'» and 8 represents the dissipation of turbulent kinetic energy (a §v(V l_¢': (V l_a') T)). The model 94 95 coefficient CR is assumed to be independent of turbulent and mean field statistical properties. The turbulent pre-stress I: ([f) depends on statistical correlations related to pressure fluctuations and fluctuations in the instantaneous Reynolds stress, I EEK-(WE) . (4.4) With f defined by . p' E O —I , 4.5 f V [p33] ( ) it follows that the correlation (ff) can be expressed as ' .. 11' E.’ 11' . (m - + +((V-{)V(%))+((V-{)(‘73)). (4.6) Although a closure for the pro-stress could be developed by analyu'ng the statistical correlations appearing in Eq. (4.6), the approach employed in Chapter 3 was based on an alternative strategy which incorporated a direct decomposition of the pre-stress into isotropic and anisotropic components: 2 2“! 2w 4 TR I") 3 ?: + = . ( .7) The anisotropic pre-stress 2k_l‘_I is symmetric and traceless. Because tr(u'r_¢') = 2k, Eq.(4.1) requires that the isotropic coefficient a be determined by the following equation 20: = ring) = 21r+2tR(r_t'_a');V(g)+t§rr(V(u)T. (tag). v0.» . (4.3) 96 Eq. (4.8) provides a means to relate the isotropic coefficient or to the mean field and the specific closure hypothesis for the anisotropic pre-stress fl . In Chapter 3, the implications of an isotropic pre-stress (IPS-) theory, for which H = 0, were examined. Eqs. (4.1) and (4.7) were applied to a homogeneously sheared turbulent flow (i.e. S a dIV—8 , (4.9) (Z de IV 8 049(0)); 2 —CP 1P _CDT—D , (4.10) where I P : I D = k/e. ('p and CD are constants independent of the local state of turbulence. The [PS-theory predicts a positive first normal stress difference and a shear thinning eddy viscosity coefficient which bridges the more traditional gradient transport regime for which (u'yu'z) ~ S with an equilibrium transport regime for which (u'yu'z)~ 1/ S . However, the [PS-theory erroneously predicts that the second normal stress difference is zero for homogeneously sheared turbulence. Moreover, the algebraic pre-closure theory with an isotropic pre-stress cannot explain retum-to-isotropy experiments for homogeneous turbulent flows (Choi, 1983; Choi and Lumley, 1984; and, LePenven, etal., 1985). Therefore, the purpose of this chapter is to further demonstrate the utility of Eqs. (4.1) and (4.7) by using an anisotropic pre-stress (APS-) model. In Section 4.2, an APS-model for the pre-stress which incorporates a phenomenological relaxation process consistent with the return-to-isotropy phenomena is introduced in accordance with Issue (iii) of Section 1.2. Retum-to-isotropy data are used to determine the phenomenological relaxation parameter. The extension of the closure model to an anisotropic pre-stress also serves to address the limitation of the IFS-theory, which had a zero second normal stress 97 difference (see (ii) in Section 1.2). The asymptotic state of homogeneous shear is used as in Chapter 3 to determine the APS-model constants. The non-algebraic nature of the APS- theory does not permit an a priori evaluation of realizable turbulent states; however, the realizability of the transient computations are verified a posteriori. 4.2 Anisotropic Pre-stress Theory Eq. (4.1) can be generalized to a class of non-inertial frames rotating at a constant angular velocity relative to an inertial frame by replacing the frame-dependent V(u) operator with V. (g')+ £2, where g. is the instantaneous velocity in the non-inertial frame and 52 represents the anti-symmetric temporal connection between the inertial and non-inertial—frames (Bird et al., 1977): a— . g, (4.11) "D where g is an orthogonal dyadic-valued operator (i.e. g - QT = I) and g is the time derivative. Thus, the pro-closure theory for the Reynolds stress contains an explicit dependence on the frame of reference through the operator :1 defined by Eq. (4.2). Unlike the Reynolds stress, the theory deve10ped here assumes that the anisotropic pre-stress £1 (= 2kg) is an objective property of the motion associated with the mean field and should not depend on the reference frame. Whence, the anisotropic pre-stress in a non-inertial frame can be related to the anisotropic pre-stress in an inertial frame by the following expression (Mase and Mase, 1992) T. (4.12) II at it ll lug II it 1110 With this hypothesis, a frame indifferent closure model is proposed for Iii. Thus, an APS- closure theory combines the pre-closure equation for the Reynolds stress (see Eqs. (4.1) 98 and (4.7)) and the following linear relaxation model for the anisotropic pre-stress {1+ hut—g (21:501] = ms). (4.13) For large turbulence Reynolds numbers, the phenomenological parameters A and B are assumed to scale with k and e: k e and p=uC-. an) CB and CA are universal model coefficients. In the above equation, the mean strain rate dyadic (_S) is defined by 29)=V@H«V@»K (4M) and the mean vorticity dyadic (2') is defined by 2(3) = V(y)- (V0017. (4.17) (if) can also be written as T-£1-§-(3)-a[(§)-§+£1-(§)]- (4.19) With a = 0, Eq. (4.19) reduces to the corotational Jaurnann derivative. For a = 1 and a = -1, Eq. (4.19) yields the upper and lower convected derivatives of Oldroyd, respectively. Appendix D demonstrates that the operator 51(3) is objective for —oo < a < +oo. Physically, the Jaumann derivative represents the—temporal changes in the pre-stress relative to a frame of reference moving with the local mean velocity and rotating with an angular velocity equal to the mean vorticity. Note that Eqs. (4.13) and (4.19) maintain the symmetry and the contraction properties of the anisotropic pre—stress, i.e., and (4.20) “:1 11 11:1 "37‘ ”It 11 O 4.3 Anisotropic Homogeneous Decay For homogeneous turbulence with no mean shear, the anisotropic pre-stress £1 equals the anisotropic stress 2kb and Eqs. (4.1), (4.8), and (4.13) imply that (g'g')s -2§k:1+§ (4.21) and For this flow, the kinetic energy is governed by the following equation (see Eq. (2.2)) 195 = -5 . (4.23) 100 For anisotropic homogeneous decay, the relaxation coefficient Ct. associated with the anisotropic pre-stress can be estimated from retum-to-isotropy data (Choi, 1983; Choi and Lumley, 1934; and, LePenven, et al., 1985).W1th 5'} = 2kg, it follows from Eq. (4.22) that -— (H) = _——g. (4.24) Eq. (4.24) shows that the pre-stress retums to an isotropic state provided C). < 1. For homogeneous, shear-free flows, the anisotropic pre-stress _I_I is equivalent to the anisotropic stress b: (5'!) 1 - —I . 4.25 2]: 3'— ( ) 11:1; 113. With the second invariant of the anisotropy tensor defined as II = tr (_b - b) , (4.26) it follows directly from Eqs. (4.24) and (4.25) that 5:1(11) = 2(l-CA) II . . 1»: dt C3. (4 27) Combining Eqs. (4.23) and (4.27) yields d(II) _ 2(l-Ct.)I_I dk .. ———C1 k (4.28) If C71 is taken to be a constant, it then follows from Eq. (4.28) that 101 II = II (1‘— )", (4.29) where 2 (l — C1) n 5 _C—_ . (4.30) 1 Eq. (4.29) provides a convenient means for comparing the above theory with retum-to- isotropy data. In Eq. (4.29), kc and Ila are reference values taken as the first data point for- which the gradient of the mean field has been effectively removed. Figure 4.1 shows the retum-to-isotropy data of Choi and Lumley [1984] and of LePenven et al. [1985], which are summarized in Tables E.11-E13 of Appendix E. In Figure 4.1, only data with a positive third invariant (III > 0) are chosen, as they more closely represent the homogeneous shear flow of interest. In order to model C). for a wider class of flows (Issue (v) in Section 1.2), C). could be viewed as a universal function of the invariants II and III (cf Sarkar and Speziale, 1990 and Lumley, 1978). However, in this work, C). is assumed to be a constant. The solid line in Figure 4.1 indicates the trajectory of the linear APS-theory in the [Mt phase plane for C1 = 2/ 3 (i.e. n = 1). It is apparent from Figure 4.1 that most of the long time decay data are consistent with a decay exponent of unity. The long time data may deviate from the solid line due to the fact that the initial reference state is not truly homogeneous and/or shear free. The data in Figure 4.1 can also be correlated by using a phenomenological closure model which assumes that the pressure-strain rate conelation in the second-order moment equation for the Reynolds stress is proportional to the anisotropic stress (Launder et al., 1975), .S. —— 0)) 103 4.4 Homogeneously Sheared Turbulence For homogeneously sheared turbulence, v (10“) 432 (E) - dy fygza ( ° ) and (1.4"!) 5= 2]: = Rxxfxf +R e ye +Rzzfzfz +Ry e y_e +Rzy_ez_ey. (4.33) The pre-closure theory (see Eqs. (4.1), (4.2), (4.7), and (4.8)) applied to homogeneous shear yields the following relationships between the components of the Reynolds stress and the components of the pre-stress: 101 Rn = -3- I + anr , (4.34) R - ”+11 435 yr " ‘3'; yr ’ (' ) 201 R = 1--—+Hzz 4. 6 a 3k ( 3) and Ryz = -NRRyy+Hyz. (4.37) In Chapter 3, the above set of equations were examined for the special case H= 0 (IPS-theory). Eqs. (4.35) and (4.36) imply that the first normal stress difference rs given by 104 R -R =1—%+H —H . (4.38) The second normal stress difference follows by subtracting Eq. (4.34) from Eq. (4.35): R” — Rn = H” — H“ . (4.39) Eq. (4.39) shows that, if R” —Rn :6 0, then the pre-stress must have an anisotropic component which has a second normal pre-stress difference. Eq. (4.38), however, shows that the isotropic portion of the pre-stress causes a primary normal stress difference in the Reynolds stress even if the anisotropic part of the pre-stress is zero. Thus, as previously noted by Parks et al. [1997] (also see Chapter 3), the primary role of the isotropic pre- stress is to redistribute the kinetic energy of turbulent fluctuations among the velocity components subject to the normalization requirement that tr (5) = 1. For homogeneous shear, Eq. (4.8) reduces to or Z = 1+2NRRyz+N2RRyy , . (4.40) where Sk NR-a-tRS = CR‘E . (4.41) The phenomenological role of the isotropic pie-stress is clearly indicated by Eqs. (4.34)- (4.36). On the other hand, the underlying statistical aspects of the flow which cause a are expressed in Eq. (4.40) and the contraction of Eq. (4.6). The APS-theory presented here assumes that the anisotropic pre-stress is governed by the phenomenological model given by Eq. (4.13). With {:1 = 2kg, Eq. (4.13) can be written in dimensionless form as 105 “meets-idea] = we («42> where the temporal operator _ll! is defined by Eq. (4.19) with :1 replaced by __I_l. In the aboveequation, le sl-— . 4.43 4 km ( ) and 1: _ka — — Vk. 4.44 D: '81 +(u-) ( ) For homogeneous shear, the normal components of the anisotropic pre-stress satisfy the following set of ordinary differential equations (1an e‘ d§ xx'EaD‘rHyz = 0, (4.45) dHyy e, dg’ Hyy+(§ +1)DeH c, yz= (4.46) and D dH“+(1+)H + “—1 D H -—0 447 er}? 4 :2 3 er yz " - ( ' ) Note that the sum of Eqs. (4.45)-(4.47) preserves the anisotropy property, tr (:1) = 0. In the above equations, the dimensionless development time is defined as E; = zS/(0) . (4.48) The turbulent Deborah number De, introduced in the dimensionless formulation 106 compares the phenomenological relaxation time for the anisotropic pre-stress with the characteristic time of the mean field: De 525: C E. (4.49) e The parameter q (see Eq. (4.43)) compares the characteristic relaxation time for the anisotropic pre-stress with the characteristic tumover time of the turbulent h'netic energy. The off-diagonal components of the anisotropic stress satisfy the following ordinary difi‘erential equations dey e, da’ 0 , (4.50) dez e! dgz =,0 (4.51) and dHyz Det CDC, 1 [CS 1 d: 701,, -16Iy )+ (sz+Hyy) = 438. (4.52) 4.5 Asymptotic Homogeneous Shear As the dimensionless development time increases, the above set of equations predict the existence of an asymptotic state provided 168’ Its lim — = (—) (on. (4.53) Thus, with q -> q“ < co, De, -) De,“ < co, and (kS/e) a < co, the above set of equations imply that (see Appendix G) H; = 0, H; = 0 107 DeaHa .. = 35'- ‘ ". (4.54) 3 1+qa DeaHa W -(1+3) ‘ ", (4.55) 3 l+q DeaHa H“ = (1—3)——‘—’—‘. (456) a 3 1+q“ and C (1+q“)Dc° H“ = B ‘ . (4.57) yz 2C a 2 a 2 2 1(1+q) +(Det) (l-a /3) It follows from Eq. (4.43) and (4.9) that le a —— = -1 , 4. 8 4 1“” CM ) (5) where 4’ represents the ratio of production to dissipation of turbulent kinetic energy: -(u' u' )S r s —’—?—. (4.59) 8 Eqs. (4.9) and (4.10) imply that CD-l glimP=lPa=C 1' (4.60) 409 P- Thus, the asymptotic value of q can be related to the model parameters in the e-equation and to the relaxation coefficient C1: 108 Ct (Co " Cpl lim = a: . (4.61) gel] 4 CP-l It follows directly from Eqs. (4.54), (4.55), and (4.57) that a negative second normal stress difference develops for the APS-theory provided CBC). > 0 (assuming -J3 < a < 1) inasrnuchas a a Ha Ha DefHa Rid-Rn: yy- n=-(1-a)l a yz: +9 _ C C 2 (12“) B ‘ (“is"). (4.62) a 02 a2 2 (l-l-q) + (Det) (l—a /3) The APS-theory at large turbulent Reynolds numbers contains six phenomenological coefficients: CD, CD CR, C13’ C13 and a. Isotropic, homogeneous decay requires CD = 1.83 (see Chapter 2). Anisotropic, homogeneous decay (i.e. retum-to—isotropy) requires C1 = 2/ 3 (see Figure l). ‘ The statistical properties of asymptotic homogeneous shear measured by Tavoularis and Kamik [1989] can be used to estimate Cp CR, CB’ and a. For instance, as § —> on, (I?) = 4.16, (4.63) a and ng—R; = —0.040. (4.64) The invariants II and III associated with the anisotropic stress _b have also been measured. These parameters are given by (also see Eq.(4.30)): II = tr (:b ~ £7) (4.65) 109 and III = tr (£2 - _b - b) . (4.66) For g—m, the data of Tavoularis and Kamik [1989] imply that 11: = 0.138 and 111: = 0.0174. The existence condition for an asymptotic solution (see Eq. (4.60)) can be rewritten as kS CD ‘ 1 —2R" (—) = . (4.67) ’1 e a CP—l A consistent set of model parameters (CB CR, CB’ a) can be identified for which Eqs. (4.38), (4.63), (4.64), and (4.67) are satisfied exactly and which also exactly reproduces the experimental values for the asymptotic normalized Reynolds stress components. It is found that C? = 1.60, CR = 0.271, CD = 0.174,and a = —2/3.Tab1e4.l summarizes these parameters 4.6 Transition States for Homogeneously Sheared Turbulence The component equations for the anisotropic pre-stress depend on the development of the time scale k/c. Because the mean strain rate is constant (S is a constant) the k- and 8- equations can be combined into a single equation for the dimensionless relaxation group (see Eq. (3.52)): (Wk 7&- = ZNRRyZ(CP-l) +CR(CD-1)’ (4568) where N R = CRSk/e. The development of the Reynolds stress towards an asymptotic state can be calculated by solving Eq. (4.68) along with Eqs. (4.45)-(4.47) and (4.50)— (4.52). The transient calculations assume that an initially isotropic, homogeneous 110 Table 4.1: Parameter Estimates for the APS-Theory Parameter Estimate Basis CD 1.83 Isotropic Decay (Comte-Bellot and Corrsin, 1971; Mansour and Wray, 1994; see Chapter 2) Cp 1.60 Existence condition for k-e equations CR 0.271 , , . Reproduction of asymptotrc state for Cg 0.174 homogeneous shear (Tavoularis and Kamik, a 413 1989) C; m Retum to isou'opy data (III > 0; Choi and Lumley, 1983; LePenven et aL, 1985) 111 turbulence is subjected to an instantaneous increase in the mean shear. Therefore, 0 for g < 0 N = (4.69) N; for§ = 0. Because N R = 0 for § < 0, the Reynolds stress and the pre-stress are isotropic; however, once N R = N2, the Reynolds stress and the isotropic part of the pre-stress respond instantaneously and attain a state consistent with the IPS-theory: £1 for §< 0 (4.70) 11 :5 II R0 for§ = 0, where 13" is the normalized lPS-Reynolds stress for N§>0 and _H" = 0. Thus 5° satisfies 4;" ~g= ———I. (4.71) The operator :10 is defined by Eq. (4.2) with NR = 117;. Thus, Eq. (4.70) defines the initial state of turbulence for which the pre-stress is isotropic and the Reynolds stress has a degree of anisotropy commensurate with the initial relaxation parameter N2. Because 5" = 9, it follows from Eqs. (4.50) and (451) that ny = 0, and H” = 0 for §>0. The transient behavior of the components of the anisotropic pre-stress together with the relaxation group N R were determined by numerical integration using a fourth-order Runge-Kutta integration algorithm (Camahan, et al., 1969; see also Appendix H). Figure 4.2 shows the transient response of the relaxation parameter N R as a function of its initial value (with 0 < N; < 10). Unlike the IPS-theory, the approach to the asymptote is not necessarily monotonic; as the trrrbulence nears the asymptotic state, oscillation in ll2 10: 10.0 5.0 3.5 \ 1 1.5 ........... g N; = 1.127 0.50 NR 0.15 0.1 r N; = 0.0 001 1 11111111 1 1111141L 1 1 1.1.11] 01 1 g 10 100 Figure 4.2‘ The Effect of the Development Time on the Relaxation Group for Homogeneously Sheared Turbulence 113 the Reynolds stress components (cfi Figures 4.3 and 4.4) can cause the first term on the right-hand-side of Eq. (4.68) to generate oscillations in dNR/dE. Approximately ten developmental time units are required to reach the asymptote, with more being required for larger values of NOR. Figures 4.3 and 4.4 present the transient development of the turbulence from an initial state characterized by N; = 0.7 . This initial value of N R is selected to closely represent initial states (in terms of the invariant pair (II, 111)) from the experimental dataThese figures present the individual components of the normalized Reynolds stress and the anisotropic pre-stress, respectively. The development times are on the order of ten to twenty units. Compared to the IFS-theory, the salient difference is an exact representation of the asymptotic values of the Reynolds stress through a non-zero secondary normal stress difl'erence. From Eqs. (4.34)-(4.37), it is apparent that oscillations in transient development of the Reynolds stress components are directly related to similar oscillations in the components of the pre-stress anisotropy H. As such, it is instructive to investigate the equations governing the transient behavior of the anisotropic pre-stress as well as the individual terms which contribute its transient evolution. Eqs. (4.52) and (4.47) may be recast as the following by dividing each term by Det: dH yz 1 p (1+q) 1 a —_ =-—---———-—-—-H -- H -H -- H +H 4.2 d§ 2c2t Det fl 2( 22 yr) 2( n it) ( 7) and (4.73) dez _ (1+q)H _ (g_ l . d§ De, zz 3 ’7 Figures 4.5 and 4.6 show the contributions to the development of dHyz/d§ and dez/d§, respectively, for the simulation shown in Figure 4.4. In Figures 4.5 and 4.6, the solid line represents the net time derivative and the thin lines represent the individual contributions to the net time derivative. From Eq. (4.72), it is evident that an initially 114 0.7 r l 0.6 - * R“ = 0.568 1' ZZ 0.5 r 0.4 r Ry. 0.3 r ij = 0.236 0.2 - Ryy = 0.196 b —R"_ = 0.165 _lz 0.1 - l 00 1 11111111 1 11111411 1 11111‘111 0.1 1 E 10 100 Figure 4.3 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7): Components of the Normalized Reynolds Stress 0.2 ~ H‘z'Z = 0.152s 0.1 r W = 0.0561 Hit l 0.0 H; = —().0564 -0 1 . H‘y’v = —().0964 _02 #1111111] 1 11114111 1 1 111111] 0.1 1 10 100 c‘. Figure 4.4 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7): Components of the Anisotropic Pre-Stress. 116 0.2 - 1% 2 C). 0 1 ’ dHyz (Ia dH yz dé a A —5 (”a +Hyy) 0.0 - VAV l _ ( 0:!) H” l -0.1 - 1 H H —E( zz— w) _02 J 11111111 1 11111111 1 llllllJJ 0.1 1 E, 10 100 Figure 4.5 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7): dHyz/dé. 117 0.15 l l 0.10 - a (l-§)Hyz 0.05 - (1H 22 dé /\ dHZZ 0.00 v=~ (m -0.05 * _(1+q)H Del 12 l _010 1 11111111 1 11411111 1 1(111111 0-1 1 10 100 Figure 4.6 Transient Response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7): dHZZ/dé. 118 isotropic turbulence (i.e. _H = 0) would remain isotropic in the absence of the linear strain-coupling coefl‘icient—( CB).and the linear relaxation coefficient (CA). Since both of these coefficients are modeled as constants, the term gCfl/ C1 provides a constant source of anisotropy in the pre-stress, causing the shear component of If to develop more rapidly than the normal components (cfi Figure 4.4). Once a signifith component Hyz has developed, it serves as a means to develop normal pre-stress anisotropies (cfi Eqs. (4.45)— (4.47)). From Figure 4.5, it is seen that, once appreciable normal pre-stress anisotropies have been developed, they serve to balance the anisotropy source term iCB/ C1 for the component Hyz. For the normal pre-stress anisotropies, the ultimate asymptotic value obtained represents a balance between a growth term (proportional to H yz) and an exponential decay term (proportional to the component itself). Figure 4.6 shows this representative behavior for the component H zz' The source of the oscillations may be determined by investigating the governing equations for the anisotropic pre-stress components in the following matrix form: 11' (4.74) 11 nu, ne'- 4. 1e- where the tensors and vectors in Eq. (4.74) are defined as follows: ldHyy/dg = dez/dE , (4.75) _dHyz/d§_ 3: 1 (4.76) I a- 111 m b a 0 (4.77) 119 and a 0: 0 “-3-1 9;. 0 o -;+t. (4.78) -a+1-a-1 a l. 2 2 .. In Eq. (4.78), a is the parameter associated with the convected time derivative (see Eq. (4.19)) and (1 denotes the term -(1+q)/Det. Note that the component H” does not appear in the above equations as it is not independent (due to the anisotropic pre-stress being traceless) and has been eliminated from the system of equations. Mth the system of equations expressed in the form given by Eq. (4.74), the solution to the characteristic equation for the mauix 9 yields the eigen values of :13. Positive eigenvalues indicate that the asymptotic solution to Eq. (4.74) is not stable, while the reverse is true for negative eigenvalues. In the event of imaginary eigenvalues, the solution exhibits oscillatory behavior (Camahan et al., 1969). The characteristic equation for the matrix 3 is: det (g - 1:!) = 0. (4.79) Eqs. (4.78) and (4.79) yield the following cubic equation: (a-X)[(a-?~)2-¢(a)l = o. (4.80) where the parameter 4) (a) denotes the following collection of terms «a=(—:-)(2:—1—l+(;+1)(;) By inspection, the roots to Eq. (4.80) are 120 3.1: or 7.2 = ot-JS . (4.82) 3.3 = ore/6 Since a < 0 (i.e. q > —1) for all points during the transient simulation, this represents a contribution to a stable asymptotic state. For -./3 < a < 4/3, 4) (a) < 0, meaning that 2.2 and 3.3 are a complex conjugate pair, as the value for a selected in these simulations is -2/ 3. Thus, since the real parts of the eigenvalues are all negative, the asymptotic state is stable. However, since two of the eigenvalues have non-zero imaginary parts, the approach to the asymptote will exhibit oscillations. The isotr'Opic and anisotropic parts of the pre-stress represent two different responses of the turbulence to an external force. The isotropic portion of the pre-stress causes an instantaneous response to a mean shear, immediately reorganizing the Reynolds stress to an anisotropic state. Subsequently, the components of the anisotropic pre-stress relax towards their asymptotic values as illustrated on the invariant phase plane by Figure 4.7 and on the A-hyperplane by Figure 4.8. On a very short timescale, the turbulence relaxes towards the isotropic state, but then reverses direction and approaches the asymptote. As it nears the asymptote, the oscillations in the components of the pre-stress (see Figure 4.4) result in the two trajectories in Figures 4.7 and 4.8 approaching the asymptote in a contracting orbit. 4.7 Conclusions The presented anisotropic closure for the pre-stress extends the predictive capabilities of the IPS-theory for the case of homogeneous shear due to the interaction of a linear dependence on the mean strain rate and a frame invariant relaxation effect. Both effects are required to generate a non-zero second normal stress difference as well as an improved primary normal stress difference. The fact that experimentally observed values of the 121 0,20 - Locus of IFS-Turbulent States—\ 0.15 L . / I] APS-Asymptote ’ APS—Initial State 0.10 - APS-Transient 0.05 - l Invariant Boundary 000 _ 1 1 1 J 1 l 0.00 0.01 0.02 0.03 [11 Figure 4.7 Transient response of Isotropic Turbulence to a Sudden Increase in the Mean Strain Rate (N; = 0.7): Anisotropy Invariants 122 538% $55 A We 0 wzv 3mm Embm :82 2: E 388:— :ocesm m 9 ooze—spear 3358— we 3:88. 6388... we came a? c 5.8 . an: ex _ ceasefimm? EflmemrTmaQIIUp V/l 8.3m we: “cougar—Mme? 1 . %.a \ v/ moemem ma: team» meafizm ,.. .. . 80583.44. made. v/ 8mm— .axm Am a .J u Ase be .445 ___ m a .3 123 second normal stress difference are negative requires that the phenomenological parameter C13 > 0. The choice of a convected time derivative in the phenomenological model for the anisotropic pre-stress is required in order to have a non-trivial second normal stress difference in the absence of non-linear terms in the mean strain rate. Moreover, an interpolated convective derivative (with a _=_ —2/ 3) is specifically chosen to represent the finite memory associated with the anisotropic pre-stress (Issue (iii), Section 1.2), inasmuch as it is able to exactly reproduce the experimentally observed asymptotic Reynolds stress components. The turbulent pre-stress exhibits two qualitatively different means of reacting to an external driving force, such as mean shear. The isotropicportion of the pre-stress causes an instantaneous reorganization of the turbulence in response to the imposition of the mean shear. This state is equivalent to the state predicted by the [PS-theory at a given value of the relaxation parameter. Note that this rapid effect results in a non-zero primary normal stress difference, although the second normal stress difference remains zero. Conversely, the anisotropic portion of the pre-stress reacts slowly to an imposed mean shear. The time scale for this reaction is related by the turbulent Deborah number, De,. The anisotropic pro-stress causes the second normal stress difference and partly influences the primary normal stress difference. As the IPS-theory is realizable for all NR given a positive turbulent kinetic energy and dissipation, the initial state of these transient calculations is always realizable. The turbulence predicted by the APS-theory remains realizable during the entire transient approach to the asymptotic state. Similar to the [PS-theory, the developmental times for the turbulent statistics to reach their asymptotic values depend on the initial conditions and are comparable to those seen experimentally CHAPTERS HOMOGENEOUSLY SHEARED TURBULENCE IN A ROTATING FRAME OF REFERENCE 5.1 Introduction ‘ Homogeneously sheared turbulence relative to a rotating frame of reference provides a critical test for Reynolds stress closure theories. Figure 5.1 shows a schematic of this flow. The mean velocity gradient is constant within the non-inertial frame: (5.1) V (y) = 55,32. where S is a constant and both 5’ and _ez represent mutually perpendicular unit vectors in the non-inertial frame. For this problem, turbulent fluctuations couple with the mean velocity gradient and the rotational tensor, def’med by (5.2) 1.3 = £19 = meg-5,5,). where g) (.=. 951) is the rotation vector, 5: is the permutation triadic, and Q is the scalar rotation rate. The impetus for investigating flows within a rotating frame of reference stems from the following observations drawn in Section 1.1: the superposition of a frame rotation upon a simple mean shear generates a cross- stream turbulent production term analogous to the cross-stream turbulent production which arises in inertial frame flows with streamline curvature. 124 125 y A d(uz) = S Jr 24 19 F low Properties: V (y) = $5ng 19 = 951 Q = Q(e e -e e ) = -y-z -z-y Figure 5.1 Schematic for Homogeneously Sheared Turbulence in a Rotating Frame. 126 Previous authors (Spen'ale and Mhuiris, 1989; Speziale, et al., 1990) have investigated rotating homogeneous shear flows using two-equation models and second order closures. They report that asymptotic homogeneous shear states exist for a finite range of the relative rotation rate (Q/ S ), -—oo < (fl/S) m." < Q/S < (SI/S) m < co. (5.3) The asymptotic states are attained for large development times (t E z/ (“2) (0)) and have the following limiting behavior lim (k, e, g) —> (on, on, constant) . (5.4) t-too For (SI/S) < (fl/S)”. and (fl/S) > (SI/SM“, the flow changes qualitatively inasmuch as the turbulent kinetic energy It and the dissipation e decay rather than grow, i.e., ‘lim. (1:, e, g) -> (0, 0, oo) . (5.5) The standard k-e model of turbulence is frame indifferent due to the fact that no explicit frame-dependent terms appear in the k- and e-transport equations (see Appendix B and C) and the fact that the basic Boussinesq approximation for the Reynolds stress is frame-indifferent. Therefore, the k-e model incorrectly predicts no influence of Q on the low-order statistical properties of homogeneous shear. However, the He theory does predict a mean field dependence on the rotation because of the Coriolis terms in the Reynolds equation (see Eq. (1.1)). Moreover, second-order closure models for the Reynolds stress (see Appendix B) also include explicit Coriolis effects and, thereby, account for the influence of Q on low-order flow statistics. However, closure models for 127 the statistical correlations which appear in the second-moment equation for (5'1!) are often assumed to be frame indifferent. For example, the [RR-model of Launder et al. [1975] relative to a rotating frame of reference reduces to the following equation for statistically stationary rotating homogeneous shear flows: d(u'u') <0) 2; + (yap- [1700+ 231+ [1700+ 23]" <_'t_r> = 2 e --3-g=I-Cli=b +Ca (<1!!!) ° [V <9) +2] + [‘7 +§]T- (5'1”) - § (2'22” <10!) (56) The convective coupling of the Reynolds stress with the mean velocity gradient and the rotation dyadic on the left-hand-side of Eq. (5.6) arises naturally with a change of frame. The terms on the right-hand-side account for two distinct physical processes: (1) an isotropic destruction of the Reynolds stress due to turbulent dissipation of energy; and, (2) a redistribution of energy among the components of the fluctuating velocity. Because (u'u'):__f=l = 0 and tr (b) E 0, a contraction of the two redistribution terms is identically zero. Thus, these terms do not explicitly influence the energy balance (cf Eq. (1.15)) which follows directly from the trace of Eq. (5.6): (uz)(0) int - u') = -2(u'u'):V (u)—2e. (5.7) dz - - - - - Eq. (5.7) has the same form and physical interpretations as its inertial frame counterpart. Clearly, changes in the kinetic energy of the fluctuating field for homogeneous shear occur because of turbulent production and turbulent dissipation. The last two terms on the right-hand-side of Eq. (5.6) account for the slow and the fast redistribution of energy due to the pressure-strain rate correlation. 128 Application of the LRR-model to retum-to-isotropy experiments in an inertial frame (see Section 4.3) reduces to the same dynamic model as the APS-theory. With a Rotta coefficient C1 = 3.0, the LRR-predictions and the APS-predictions (for CR = 2/ 3) are indistinguishable when applied to the relaxation of homogeneous turbulence to an isotropic state. However, these two closure strategies predict significantly different responses in non-inertial frames. In this chapter, the APS-theory is applied to rotating homogeneous shear flows and the results provide a basis to compare the LRR-model and the APS-model (see Section 5.5). The standard k-e transport equations are employed to determine the turbulent time scale (i.e. k/e) for both the LRR- and APS-models. However, recent research (see Speziale et al., 1987) suggests that the scalar dissipation equation should also have an explicit Coriolis effect. Earlier DNS results for isou'opic decay by Bardina et al. [1985] and by Speziale et al. [1987] wem to support this view. However, this fundamental question is not addressed in this dissertation. Instead, the explicit Coriolis effects in the operator 3 (see Eq. (4.2)) are evaluated and compared with the LRR-model. 5.2 Pre-closure Theory The turbulent pre-closure for rotating homogeneous shear is the same as that expressed by Eqs. (4.1), (4.2) and (4.7), with the exception that the velocity gradient in the inertial frame (cf Eq. (4.2)) is replaced with the sum of the mean velocity gradient in the non-inertial frame and the rotation dyadic (see Section 4.1). Thus, the turbulent pre- closure for rotating homogeneous shear is QT. (QM-1:1 = 201 _ 71+ 2151, (5.8) 129 where §=!+1R(V (9+2) = I+NR[(1+%)-€y—ez-S—;Ez§y] . (5.9) The trace of Eq. (5.8) utilizing Eq. (5.9) yields the following expression for the isotropic portion of the pre-stress: o‘-121v 921v: 1 92111211 510 Z ' + RRyz+ 3" RRzz+ +37 R Y)’ ° (' ) Eqs. (5.8) along with (5.9) and (5.10) gives the following component equations for the pre-closure: R — 1“ 1 xx - 32+Hxx’ (5.1 ) lot a o 2 2 R” — a; +Hyy+2§NRRyz- (E) NRRZZ’ (5.12) lot a o 2 2 1122 = 3? +Hu-2—S—Nkkyz—2Nkkn-(l +3) NRR”, (5.13) and n o o a [-1+§(1+E)N;](—Ryz) = Hyz+§NRRu- (1+?)NRR”. (5.14) Note that Eqs. (5.10)-(5.l4) reduce to Eqs. (4.40) and (4.34)-(4.37) for Q = 0. As was observed in Chapter 3, the isotropic portion of the pre-stress (Eq. (5.10)) serves to redistribute energy among the normal components of the Reynolds stress. Each of the last three terms on the right-hand-side of Eq. (5.10) is found in the normal component equations with equal magnitude, but opposite sign. In Chapter 3, it was noted 130 that the term ZNRRyz shifts energy from R” and R n to the streamwise component Ru. Conversely, the final two terms in Eq. (5.10) cause a redistribution due to coupling with the frame rotation to remove energy from R” and R zz’ which subsequently transfers it to the axis of rotation, i.e. to er The term 2 (Q/S) NRR),z also arises in the component equations for R” and R 22’ but with opposite sign. These relate directly to the source/sink terms for velocity fluctuations due to coupling with the rotation dyadic (see Eq. (1.30)). Because they are of opposite sign, energy transfers from R” to Ru (or vice versa, depending on the direction of rotation). The mathematical nature of the operator fl is important in that, if its determinant were zero, then there would be no algebraic connection between the Reynolds stress and the statistical correlations responsible for the pre-stress (see Eq. (4.7)). For homogeneous shear flows, the determinant of g is det (g) = 1 +1112?- (1 + 552). (5.15) From Eq. (5.15), it is apparent that det (i) could only possibly be negative for -l < Q/S < O. In this range, it follows directly that NR < 2 is a sufficient condition for det (g) > 0. For the asymptotic state associated with a particular value of Q/ S , it would be possible to verify a posteriori that det (1:1) > 0. However, for a transient calculation in the range -1 < Q/S < 0, there would exist some upper bound above which an initial condition for N R could not be specified: 9 Q -l/2 NR’W=[(-§)(1+—)] 22, for -1 Ill . N +1122 =;I+1:RQ e ez -e ey]. (5.17) - .1— From Eq. (5.10), it follows directly that the isotropic portion of the pre-stress is given by a-IQ = 1+ (239) 2 (Ryymu). (5.18) For isotropic turbulence in an inertial frame, 3 = 9 and a/k = 1, since 2 = :1. As 0 changes from zero (either positive or negative), the anisotropic pre-stress remains zero in the absence of any mean field deformation process. However, this isotropic part of the-pre- stress changes according to Eq. (5.18) and the operator 5 depends on the rotational rate 0. Therefore, with S = O, Eqs. (5.8) and (5.17) lead to the following component equations for the pre-closure: 1 R1): = adv/k, (5.19) 2 _ I Ryy-r- (1R0) Ru - ga/k, (5.20) 1912+ (2129) 2Ry y = got/k, (5.21) and Ryz = 0. (5.22) Eqs. (5.20) and (5.21) combine with Eq. (5.18) to give 132 1 R = R = -—-—- . (5.23) ’y 12 3+ (1R9)2 This implies that 1+ (11,0)2 )3 = (5.24) 3+ (cko)2 since Rxx+Ryy+Rzz = 1. Thus, the noninertial pre-closure theory for the Reynolds stress predicts a redistribution of turbulent kinetic energy among the components of the fluctuating velocity in the absence of any mean field deformation. As the characteristic time scale for the turbulence becomes large compared to the time scale for rotation (i.e. Itkfll » l), the turbulence approaches a one-component state (cfi Figure 1.1), with the turbulent energy aligned along the axis of rotation. This redistribution process is symmetric about 1R0 = 0 and is summarized in Figure 5.2. It is also noted that Eqs. (5.23) and (5.24) express a realizable turbulent state for all IRQ 5 [-oo, oo] . The goal of this chapter is to apply the noninertial pre-closure theory to homogeneously sheared turbulence (see Issue (iv), Section 1.2). The asymptotic states for both the IPS- and APS- theories are examined along with the effect of rotation on the qualitative nature of the turbulence. The symmetry about (2 = 0 predicted by Eqs. (5.23) and (5.24) will be broken by the presence of S > 0. 5.4 Isotropic Pre-Stress Theory for Rotating Homogeneous Shear Similar to the approach taken in Chapters 3 and 4, the consequences of the IFS-theory are considered prior to evaluating the APS-theory. Thus, for the following analysis in this 133 I 0.8 0.6 1 XX j I 0.4 T Y)” RZZ 0.2 00 llllLLlll 11111111L 1111111 '0.01 0.1 1 10 100 TRQ Figure 5.2 Coriolis Redistribution of Turbulent Kinetic Energy for Homogeneous Flows 134 section, H = 0. For the problem of rotating homogeneous shear, no new phenomenological parameters are introduced, so the results may be computed directly from the parameters outlined in Table 3.1. The asymptotic states for the IFS-theory are determined by solving the initial value problem for the relaxation group N R at a fixed value of Q/S. Since the evolution equation for NR does not explicitly introduce frame-dependent effects, the governing equation is the same as Eq. (4.68): dNR 1E- : 2NRRyz(CP—1)+CR(CD—l) . (5.25) If asymptotic states exist (i.e. dNR/d§ -) 0), then it follows directly from Eq. (5.25) that N RR” does not depend on Q/ S inasmuch as: lim (-NRR ) = CR(CD_1) = 0230 (526) g..- fl 2(CP—l) ' ' ' The above numerical result assumes that the IPS parameters listed in Table 3.1 are universal. With RM = l-Ryy-Rzz and g = Q, Eqs. (5.ll)-(5.l4) may be reduwd to three coupled algebraic equations for R”, R”, and R”. During the transient calculations, the coefficient matrix for equations relating these three components is inverted using a Gauss-Jordan elimination technique (Camahan et al., 1969). Although there is no explicit dependence on Q/ S in the equation for NR, the asymptotic state N; does depend on Q/ S due to the Reynolds stress components. The initial state used to compute the approach to the asymptotic states is the inertial frame asymptote: N; = 0.945. (5.27) The initial value problem was solved using a fourth order Runge-Kutta integration scheme (Camahan et al., 1969; see also Appendix H). 135 With H = 0, Eqs. (5.12) and (5.13) may be expressed as follows: C C [R ] 1 Q/S 11 12. ’2 = ( )-2NR ( ) (528) R yz ’ ' (C21C22) Rzz l 1+Q/S where (211 = €22 = 2, (5.29a) Q 2 C12 = l+N§(—§) , (5.2915) and Q 2 c21 = 1+N§(1+—S-) . (5.29c) For g —) no, N RR yz approaches 0.230 as indicated above. Therefore, the determinant of the coefficient matrix in Eq. (5.28) must be non-zero: det(£) = 4-[1+N§(%)2][1+Ni(l+%)2]#0 . (5.30) For a fixed value of Q/ S , there exists an N; such that the condition det (9 = 0 is met; this identifies states for which asymptotic states do not exist. Thus, the limits on Q/ S for which an asymptotic state cannot exist is determined by the intersection of N; and N2. Both N; and N; depend on ms. The solution to Eq. (5.30) for det (g) = 0 is (111;)2: (”4” 2[-11J1+ 124’ 2], (5.31) 2(1+Q/S) (1+¢) where 136 = (1+Q/S) 2 ¢‘[ (II/S) ] Note that only the positive root from Eq. (5.31) is considered, as N; (SI/S) is a real- valued quantity. Figure 5.3. presents the results of the IPS-asymptotic states for N2. These results are qualitatively similar to other second-order closure models. The IPS-theory predicts unbounded growth of k and 8 over the range -l.26 < 0/ S < 0.26. A similar behavior is predicted by the second-order closure of Launder et al. [1975] for —0.11 < Q/S < 0.39. Outside of this range, the flow changes character; both It and 8 become decaying functions and N R —-) no (see below). The minimum value of the relaxation group is N; m." = 0.834, which occurs at Q/S = —0.5; the maximum value of the relaxation group is N; m a 1.27, occurring at 9/5 5 0.26 and ms 5 -l.26. Figures 5.4 and 5.5 complement Figure 5.3 in that they show the components of the asymptotic Reynolds stress as a function of Q/ S and the corresponding behavior on the turbulent energy simplex. Figure 5.4 shows the redistribution effects due to the isotropic pre-stress and the rotation coupling outlined in Section 5.2. As IQ/SI increases, the terms (Q/S) zNiku and (1 + Q/ S) 2N2Ryy grow, causing the isotropic pre-stress to transfer turbulent energy to the rotation axis. Similarly, the redistribution term arising due to rotation (~ (SI/S) N RR”) causes a net transfer of energy to Rzz from R” when Q/S > 0 and vice versa when (2/ S < 0. Both figures indicate that the solution is symmetric about the point 9/ S = -0.5. Specifically, the following properties are noted: Raul/S) = R3 (-Q/S— 1) , (5.32) Ryy(Q/S) = Rzz (-Q/S- 1), (5.33) 137 N‘RW = 0.834 1.2 ~ 2 .9 M 1.0 _ N; 0 45 0'8 P \M N; "m :— 1.27 l ‘ 0.6 ~ NR z 0.26 0.4 F max I l l ' Decay Regime | Asymptotic Regime : Decay Regime I k —-> 00 l k —> O 0.2 '- I I I a —> 00 I 8 —> 0 I k/e ——> constant ' k/e —> 00 I I 0 O 1 l l 1 l l 1 l -2 0 -1.0 0 0 1 0 Q/S Figure 5.3 The Effect of Rotation on the Turbulent Relaxation Time (N R 5 IRS) for Asymptotic Homogeneous Shear (IFS-Theory) I38 0.6 ' Ry)" W R I I 22 \f\‘ I I I I I l 0.4 - I | I I I I Rn l \. I Rij ’ I I I I /0\ I M I 0.2 - | my, : l I I I I ' I I Symmetry I . Line _ (9) V1.26 : (€53) z0.26 min I max I I I I 00 1 I l l l I J J -20 -r0 00 10 Q/S Figure 5.4 The Effect of Rotation on the Components of the Asymptotic Reynolds Stress for Homogeneous Shear ([PS-Theory ) I39 X" -1.25 Q Isotropic Point \ \ \ \ \ \ Symmetry Line\ \ Figure 5.5 Distribution of the Energy Components for Homogeneously Sheared Turbulence in a Rotating Frame (IPS-Theory). 140 and Ryz (Q/S) = Ryz (-Q/S— 1) . (5.34) It follows directly from Eqs. (5.32)-(5.34) and the property R x x + R y y + Rzz = 1 that the anisotropy invariants share the same symmetry property: ”(Q/S) = II (—Q/S - 1) , (5.35) and III (Q/S) = III (-Q/S — l) . (5.36) These results stem from the isotropic pre-stress assumption and the fact that g (9/5) = {Hz/s — 1) . (5.37) Eq. (5.37) notwithstanding, subsequent introduction of an anisotropic pre-stress will be seen to destroy the symmeuy properties possessed by the IPS theory. Figure 5.6 shows the transient behavior at infinite Reynolds numbers of N R for three values of ms. The initial conditions correspond to N; = 1. The results illustrate that two different flow regimes occur: an asymptotic regime (described by Eq. (5.4)) and a decay regime (described by Eq. (5.5)). For 0/ S = -0.5, N R approaches a constant value in about ten dimensionless time units. For 0/ S = 1.0 and -2.0, however, N R grows monotonically without bound. The relative growth rates of the turbulent quantities k, a, and Re follow directly from Eqs. (3.2) and (3.3): I41 5.0 r I 4.0 I Q/S = +1.0 and -2.0 3.0 11111 NR —> oo §—)oo NR 2.0 1.0 52/5 2 -0.5 alim NR = 0.834 0.0 1 l 1 l 1 l 1 I 0 5 10 15 20 Figure 5.6 Transient Response of the Relaxation Group for Homogeneously Sheared Turbulence in a Rotating Frame (IFS-Theory). 142 __ = .. -1 , 5.38 1: dr CR ( ) 1 2C R .295 _._ __£N_£_YE_C , (5.39) t: dt CR D and 321‘: = — 2(2—CP)NRRy Re (It CR ‘ - (2 — CD) . (5.40) In the above equations, ‘0 a k/e. In the absence of turbulent production (i. e. N RR yz = 0) Eqs. (5.38)-(5.40) describe the decay problem presented in Chapter 2. If an asymptotic state obtains (see Eq. (5.26)), then all three of the above parameters approach the same limit given by r t r C -C (.1292?) = (.2515) = (.2213) = D P, (5.41) 1: dt 8 dr a Re dt 0 CP-l For the IT’S-parameters listed in Table 3.1, the dimensionless asymptotic growth rate is 1.024. The relative growth rates in the asymptotic regime are shown in Figure 5.7. Figure 5.7 shows that, after approximately ten dimensionless time units, k, e, and Re achieve their limiting growth rates. Similarly, Figure 5.8 indicates that roughly the same amount of time is required for the Reynolds stress components to reach their asymptotic values. As expected, (N RR”) 0 is reproduced exactly (see Eq. (5.26)). The initial values for the Reynolds stress components are determined by the isotropic portion of the pre- stress; a specification of values for N; and Q/ S causes an instantaneous reorganization of an isotropic turbulent state to an anisotropic turbulent state defined by Eqs. (5.10)- (5.14) with :1 = 9. Figures 5.9 and 5.10 show the transient behavior of the turbulent time scales and the I43 2.2 - 2.0 1.8 IDds ?E 1.6 ~ 1.4 - Erik rDEk/s k dt 1.2 - Lee: CD-CP = 1.024 _ - CP— 1.0 1 I 1 I 1 I 1 I 0 5 10 15 20 Figure 5.7 Relative Growth Rates of Turbulent Statistics in the Asymptotic Regime (IFS-Theory; N; = 1.0, m5 = —0.5). 1.0 0.8 0.6 0.4 0.2 0.0 Figure 5.8 144 I I R R R“,R“ = 0.387 Y)’ 12 (—NRR ) = 0.230 yza r . R = 0.226 /R..r “ Relaxation of the Reynolds Stress Components in the Asymptotic Regime (IFS-Theory; N; = 1.0; Q/S = —0.5; R" = R (N2, Q/S) , cf. Eqs. (5.10)-(5.14)). — - I45 05 ' ID: k/e In (1R8 Re 71—!— 0.0 \\ Figure 5.9 Relative Growth Rates of Turbulent Statistics in the Decay Regime (IPS- Theory; N; = 1.0; Q/S = +1.0 and +2.0). 146 1.0 r RH for both (2/ S 0.4 ' RZZ (for (US = +1.0) R”. (for 52/5 = —2.0) 0.2 ’ R” (for Q/S = +1.0) RzZ (for Q/ S -— —2.0) 0 0 -NRRYE 1 1 1 1 1 r *1 0 5 10 15 20 E Figure 5.10 Relaxation of the Reynolds Stress Components in the Decay Regime (IPS- Theory; N; = 1-0; Q/S = +1.0 and +2.0; R” = 5(N2, Q/S) , cf Eqs. (5.10)-(5.14)). 147 Reynolds stress in the decay regime (Kl/S = +1.0, -2.0). k, e, and Re ultimately decay with time (i.e. negative growth rates at large development times). This occurs because N RR N -> 0 in the decay regime (as seen in Figure 5.10). The asymptotic decay rates for k, e, and Re follow directly from Eqs. (5.38)—(5.40) by setting N RR” = 0 These asymptotes are indicated in Figure 5.9 by the dashed lines. Although NR -) co in the decay regime, Ryz remains realizable and decays more rapidly than NR grows, such that the product is also approaching zero. Additionally, it is seen that the effect of rotation in the decay regime is to asymptotically shift the turbulent energy into the energy component along the axis of frame rotation: R x -> 1, (5.41) R” —> O, (5.42) and Rzz —-) 0. (5.43) For an isotropic pre-stress, it follows directly from Eqs. (5.11) and (5.41) that 9 —> 3 (5.44) k in the decay regime. Thus, Eq. (5.41) implies that (u'xu‘x) —> 21: -) 0. Eq. (5.12) with H n = 0 and Eq. (5.44) together with the observation that N RR” -—) 0 yield the following limiting behavior for R”, 2 R” ->1- (o/S) NfiRzz. (5.45) 148 Since R” -) 0, the product NRRzz must be approaching a constant in the decay regime in order to balance the first term on the right-hand-side of Eq. (5.45). Thus, in the decay regime for large values of N R, Eq. (5.45) reduces to 1 R —> ————. (5-46) “ (1)/$2111; Similarly, Eqs. (5.13) and (5.44) combine to give 2 2 Ru -) 1 — (1+Q/S) NRRyy (5.47) in the decay regime. Thus, for large values of N R, Eq. (5.47) implies that 1 . _, (5.48) R . ” (1 +Q/S)2Ni It follows directly from Eq. (5.14) that Ryz a 0 for the limiting behavior expressed by Eqs. (5.46) and (5.48). Thus, Schwarz’s inequality is satisfied. For decaying homogeneously sheared turbulence, Eqs. (5.46) and (5.48) provide a prediction for the distribution of turbulent energy in the plane of rotation: R 2 i = (1+Q/S) . (5.49) Ryy (fl/S)2 For the two relative rotation rates presented in the decay regime (i. e. +1.0, -2.0), the transient ratio of R zz/ Ryy approaches 4.0 and 0.25, respectively, which is in accord with Eq. (5.49). Figures 5.11 and 5.12 address the issue of realizability of the IPS-theory for rotating homogeneous shear. Figure 5.11 shows all of the states within the asymptotic regime (i.e. I49 0.4 r Trajectory Q/S a 0.25, -125 i b 0.15, -l.]5 : c 0, -1 0 3 _ -0.5, -0.5 d h : i" II 0.2 1 locus of IPS asymptotic states; parameterized by Q/ S 0.1 r , Shaded area: states inaccessible to Boundary ' asymptotic homogeneous shear. O 0 1 1 1 1 I 1 1 I 1 1 I -0.03 0.00 0.03 0.06 0.09 Figure 5.11 Anisotropy Invariants for the Asymptotic States of Homogeneously Sheared Turbulence in a Rotating Frame (IPS-Theory). 150 0.4 r [I Trajectory Q/S _: a 0.26, -l.26 if. b 0.5, -15 j C I, 2 4' d 3. 4 f 0.3 . i [I 0.2 . 0.1 - Invariant Shaded area: states inaccessible to Boundary asymptotic homogeneous shear. O 0 1 1 1 I 1 1 I 1 1 I -0.03 0.00 0.03 0.06 0.09 [I] Figure 5.12 Anisotropy Invariants for the Decaying States of Homogeneously Sheared Turbulence in a Rotating Frame (IPS-Theory). 151 —1.26 < Q/S < 0.26). Each point within the L-diagram is parameterized by at least one pair of (Q/ S, N R) , with the exception of the two shaded areas. These regions contain no IPS-trajectories which relax to an asymptotic state. Trajectories on the L-diagram represent lines of constant 52/ S, along which the various points are parameterized by N R. Due to the symmetry of the IPS-theory about Q/S = —0.5 (see, esp., Eqs. (5.35) and (5.36)), each trajectory corresponds to two values of the relative rotation rate. The trajectories for Q/ S and —Q/S — 1 are identical. All trajectories are realizable and are attracted to the highlighted line in the figure, which represents the locus of asymptotic states for -1.26 < 0/ S < 0.26. The arrows indicate the path that the anisotropy invariants take as they approach the asymptote. Figure 5.12 is the analog of Figure 5.11 for the decaying homogeneous shear states (Q/S>0.26 and Q/S<-1.26). For the decaying states, only one region of the L- diagram is inaccessible. As indicated in Figure 5.10, the decaying states all relax towards a one-component turbulence asymptote, with all of the turbulent energy aligned along the rotation axis. All states between the isotropic state and the onercomponent limit remain realizable. From the preceding results, it is apparent that a rotation of frame coupled with a mean shear has two significant effects. First, for small absolute relative rotation rates in the asymptotic regime, rotation causes two redistribution effects among the normal components of the Reynolds stress via the isotropic pre-stress: Rzz <—> R” and (R yy, Ru) —) Rn. In this case, however, the flow retains the same qualitative character it had in the inertial frame: both It and a grow without bound, but at the same relative rates. However, at larger absolute relative rotation rates, the effect of frame rotation drastically affects the turbulent structure. large rotation rates have the effect of completely cutting off turbulent production due to coupling of the mean shear and Reynolds stress (i.e. N R -) 0 rather than NRR),z -) constant). This causes the turbulence to decay and, R yz thereby, to transform into a one-component turbulence with R n -) 1. It is noteworthy that 152 either asymptotic or decay trajectories by the IPS-model remain within the L-diagram, although some of the L-states are only attainable in one of the two regimes. All states predicted by the IPS-theory are realizable. 5 .5 Anisotropic Pre-Stress Theory for Rotating Homogeneous Shear The extension of the model predictions for the APS-theory follows directly from the procedure outlined for the IPS-theory. The phenomenological model for the anisotropic pro-stress was developed in Chapter 4. Again, no additional unknown closure terms and/or phenomenological parameters arise in the rotating frame (relative to the inertial frame), so the parameters in Table 4.1 allow the complete description of the rotating homogeneous shear flow. Eq. (5.25) for the evolution of the relaxation parameter applies as it did previOusly. However, the asymptotic invariant for the product N RR” is now given by CR(CD- 1) (-N ) = RR” a 2 (CP— 1) = 0.186. (5.50) The previous approach for determining the asymptotic states by solving the initial value problem remains valid; however, instead of having one differential equation for N R coupled with a set of algebraic equations, additional differential equations are required to determine the anisotropic pie-stress Li. Since the model for g is frame invariant, the governing equations have the same form as before (see Eqs. (4.42), (4.45)—(4.47), and (4.50)-(4.52)). Specifically, the four components of the anisotropic pro-stress are governed by the following set of ordinary differential equations: (11! H 4’; 2 Der + (l +q) Hn-gaDe‘Hyz = 0, (5.51) NOTE TO USERS Page(s) were not included in the original manuscript and are unavailable from the author or university. The manuscript was microfilmed as received. 153 I | (9) s —1.56 (9) e 0.19 s . . s 0.4 - I I Decay I Asymptotic Regime I Decay Regime - Regimel I I k —> 00 I k —) O 02' I saw I 890 I k/e —> constant : k/e —> 00 I 0.0 I 1 I 1 I 1 1 I -2.0 -l .0 0.0 1.0 Q/S Figure 5.13 The Effect of Rotation on the Turbulent Relaxation Time (N R E IRS) for Asymptotic Homogeneous Shear (APS-Theory) 155 group. The boundaries of the asymptotic regime are shifted (-1.56 < Q/S < 0.19 in the asymptotic regime), but, more significantly, the symmetry about 9/ S = —0.5 no longer occurs. Figures 5.14 and 5.15 show how the components of the Reynolds stress and the anisotropic pre-stress are distributed in the asymptotic regime. Figure 5.15 indicates that, even though the phenomenological equation governing g is frame indifl‘erent, the implicit dependence of £1 upon the frame due to the coupling with the relaxation group causes the pre-stress anisotropy to vary with Q/S. From Eqs. (5.11)-(5.l4), it is seen that the anisotropic pre-stress serves to systematically increase the levels of Rzz and RN, while Ru and R y y are systematically decreased. These trends directly result in the destruction of the symmetry property that existed for the IPS-theory about the point 0/ S = -0.5 . Figures 5.16 and 5.17 show the asymptotic states for the APS-theory on the energy simplex and the L-diagram. The locus of asymptotic states on the energy simplex further illustrates the asymmetry of the energy distribution illustrated by Figure 5.14. Qualitatively, the trend of redistributing the turbulent kinetic energy into the component of the Reynolds stress along the rotation axis is seen again at the boundaries of the asymptotic regime. The broken symmetry is less pronounced on the L-diagram, although it is seen slightly due to the fact that points along the locus of asymptotic states are no longer parameterized by two values of Q/S; for the APS-theory, each point is characterized by a unique value of Q/ S . The trajectories for the APS-theory are not presented in the same detail as for the IPS- theory. The purpose of these calculations was to demonstrate the qualitative differences between the asymptotic and decay regimes. The transition between these two regimes is controlled by the behavior of the product N RR,z (see Eqs. (5.38)-(5.40)). Although the asymptotic growth rates of k, e, and Re have different numerical values for the two theories (a value of 0.374 vs. 1.024, sec Eq. (5.41)), the qualitative features of the flow remain unchanged. Figure 5.18 shows the effect of rotation on the relaxation group N R for 156 0.6 - : .M ' I I I I I I 0.4 ~ I I I I I I I I Rd | I I I I I 0.2 ~ I | I I ~ I I (9) z—l.56 (9) zO.l9 I | 00 I 1 I 1 1 1 1 I -2 0 -l.0 0 0 l 0 Q/S Figure 5.14 The Effect of Rotation on the Components of the Asymptotic Reynolds Stress for Homogeneous Shear (APS-Theory ) 0.2 0.1 -O.2 Figure 5.15 157 1 I Q/S The Effect of Rotation on the Asymptotic Pre-Stress Anisotropy for Homogeneous Shear (APS-Theory ) 158 V)’ Figure 5.16 Distribution of the Energy Components for Homogeneously Sheared Turbulence in a Rotating Frame (APS-Theory). 159 0.16 . -0.15 Experimental Data (Tavoularis/Kamik, 1989) 0.14 I. W locus of APS asymptotic states; parameterized by Q/ S RRMoch I] 0.12 ' Asymptotcs Q/S=0, 0.25 Q/ S = —l.50 0.10 . g -1.56 A k-e Model I Asymptote; I all as l l I I Invariant : Boundary I J 1 0.08 020 ' ' . 1 -0.005 ' 0.005 0.015 0.025 III Figure 5.17 Asymptotic Anisotropy Invariants for Homogeneously Sheared Turbulence in a Rotating Frame (APS-Theory). 160 51)” I QAS=440 411‘ Q/S = —2.0 311* NR 211' L0 I Q/S = —0.5 Flim NR = 0.828 0.0 1 I 1 I g L 1 l 0 5 10 15 20 Figure 5.18 Transient Response of the Relaxation Group for Homogeneously Sheared Turbulence in a Rotating Frame (APS-Theory, N; = 1.0; fl” = 0). 161 Q/ S = —2.0, —0.5, 1.0. As before, the stable trajectory approaches its asymptote in roughly ten time units, while the two cases in the decay regime exhibit unbounded growth. The lack of symmetry is also evidenced in this plot by the different growth rates of N R for the two cases Q/S = —2.0 and +1.0. 5.6 Conclusions Both the IPS- and APS-theories predict two distinct regimes for homogeneously sheared turbulence in a rotating frame of reference. In an intermediate range of relative rotation rates, i.e. (Q/S) m." < Q/S < (Q/S) MI, an asymptotic regime is observed in which both the turbulent kinetic energy and dissipation grow without bound, but at the same relative rates. In this regime, the isotropic portion of the pre-stress serves to redistribute the components of the Reynolds stress as a function of Q/S and II/RRyz approaches a constant independent of Q/ S . For Q/S< (Q/S) m." and Q/S> (Q/S) I'm, the two clusure theories examined predict a qualitative change in the nature of the flow. The presence of large absolute relative rotation rates completely eliminates turbulent production (i.e. N RR yz —> 0 for large development times, rather than N RR” -> constant). Thus, the turbulence becomes uniformly decaying, with all of the turbulent energy being aligned along the axis of rotation. The nature of the pie-closure operator A combined with the isotropic pre-stress assumption yields a symmetry property for the asymptotic state about (US = —0.5. Subsequent introduction of an anisotropic pre-stress with non-zero normal stress differences breaks this symmetry property. The IPS-theory is shown to predict only realizable states for -oo < Q/S < co and is seen to span the entire L-diagram of anisotropic stress invariant pairs. Certain regions of 162 the L—diaglam are only accessible to the asymptotic or decaying turbulent regimes. The asymptotic states of the APS-theory are verified a posteriori to be realizable. CHAPTER 6 CONCLUSIONS The concept of realizability was addressed for each of the homogeneous flows studied herein. Specifically, it was desired to construct a turbulence model which would only predict realizable turbulent states. Where possible, this quality was verified a priori. When the nature of the governing equations did not permit such an analysis, the realizability of the turbulent states was verified a posteriori. For the problem of isotropically decaying turbulence, the decay coefficient CD in the dissipation equation was selected to generate a realimble decay process for 0 < Re < co. Specifically, if 1 < CD (Re) < 2, both It and e are decaying and non-negative. In the final form of the model, CD (Re) ranges between 1.4 and 1.83, guaranteeing a realizable decay. Similarly, the algebraic nature of the pre-closure when applied with the isotropic [ire-stress assumption yielded a Reynolds stress tensor which was realizable for all values of the relaxation group 0 0.26 and Q/S < -l.26 for the IPS-theory), however, the turbulence changes to a uniformly decaying process. This exchange is caused by a concentration of turbulent energy along the axis of rotation, which effectively eliminates turbulent production, leaving only turbulent destruction effects. Additionally, the IPS-theory is found to yield realizable predictions for all —oo < Q/S < co, again owing to the character of the pre-closure which preserves the trait 1 ()6 of positive eigenvalues in the pre-stress. The qualitative features mentioned above apply equally to both the isotropic and anisotropic formulations for the pre-stress. Small quantitative differences exist, such as the exact values of Q/ S which delineate the boundaries between the asymptotic and decay regimes. For instance, the asymptotic regime for the APS-theory exists over the range —I.56 < 0/ S < 0.19, which is similar to the range for the [PS-theory given above. The salient difference between the two theories lies in the symmetry property exhibited by IPS-theory about Q/ S = -0.5 which is not present in the APS-theory. This is due exclusively to the fact that the pre-stress is no longer isotropic (i.e. :1 $9) and its components provide contributions to the Reynolds stress components which are not all equal. A major portion of this work has been dedicated to the development and evaluation of two closure hypotheses for the pre-closure theory developed in Chapter 3. The isotropic pre-stress (IPS-) theory was initially investigated because it introduced no additional closure hypotheses. This was then extended to an anisotropic pre-stress (APS-) theory in order to capture relaxation effects as well as a second normal stress difference for simple shear flows. In retrospect, however, each of the two closure theories has its own distinctive advantages and disadvantages. For instance, the IPS-theory represents a significant improvement over the traditional Boussinesq approach; in a simple shear flow, the IPS-theory is uniformly realizable and predicts a non-zero primary normal stress difference. Relative to the APS-theory, the algebraic character of the IPS-theory may make it more attractive for use in practical applications. Although a 3x3 matrix inversion must be performed at each discrete point in a given domain, no additional differential equations beyond the four required for the continuity equation and the equation of motion are needed. This is a computational advantage compared to both the APS-theory as well as other traditional second-order Reynolds stress modeling approaches, which require the solution of an additional six 167 differential equations in the most general case. Additionally, the algebraic nature of the IPS-theory lends itself more readily to theoretical analysis. For instance, it is possible to demonstrate a priori that the IPS-theory predicts only realizable states for homogeneous shear flows and that the entire Lumley diagram is accessible to the theory, with each point therein being parameterized by one or more pairs of values for (N R’ Q/ S) . In contrast the APS-theory expands upon the predictive capabilities of the IPS-theory by providing additional degrees of freedom in terms of model parameters which allow it to capture more physical effects. For example, the phenomenological model for the anisotropic pre-stress is able to reproduce both the relaxation effect seen in return-to- isotropy as well as the second normal stress difference exhibited in homogeneous shear flows. As mentioned above, however, this improved predictive capabilities come at the cost of additional differential equations for each of the six anisotropic pre-stress components. Not only would this greatly increase the computational burden for more complex flow simulations, but, from a theoretical standpoint, also makes analysis and interpretation of results more complicated. CHAPTER 7 RECOMMENDATIONS 7.1 Further Study The analysis of the KArman-Howarth equation involved the semi-empirical expression of several inteng properties of the double and triple longitudinal velocity correlations. While there is not the expectation that the absolute values of the integrals themselves are necessarily correct, the working hypothesis is that the integral ratios, as they appear in the analysis of the Kannan-Howarth equation, are reasonably well characterized. Either direct simulation or experimental measurements of isotropically decaying turbulence could be applied in order to estimate the extent to which these integral properties are adequately represented, provided that both the double (B u) and triple (Tm) longitudinal velocity correlations are accurately determined (cfi Eqs. (2.36)-(2.38)). With this simulation and/or experimental data. the validity of the modeled ratios 12/ I 1 (see Eq. (2.63)) and 13/11 (see Eq. (2.71)) which appear in the modeled form of the IKH- equation (Eq. (2.44)) may be verified directly, as the parameters Re, P, and .8 can all be derived from the velocity correlation data. The variation of the integral parameter I 1 with the turbulent Reynolds number (or, equivalently, within the context of this theory, the destruction coefficient .8) can be verified by integrating Eq. (2.69): J trial/11.0) = g I [1+ (a—l) IfilIdMfl). (7.1) J 168 169 Ultimately, the degree to which the integral approximations applied in the theory do or do not agree with experimental values will either lend support to this approach or suggest further avenues for research, in terms of the manner in which the integral ratios are modeled. The process of retum-to-isotropy could be further addressed, as it has been analyzed only for the case of a positive third invariant (the state observed for homogeneous shear in an inertial frame) and has been modeled as a universal constant. There exists a larger body of experimental data which may suggest that either: (1) the time constant for decay is a function of the local turbulent state such that C1. = C71 (II, III, Re) ; and/or, (2) the decay process should include both linear and non-linear relaxation effects (Sarkar and Speziale, 1990; Speziale, 1991). ‘ The dissipative time scale in the e-equation has been modeled as 10 cc k/e. This formulation is not capable of explaining the effect of a slower rate of decay of an isotropic turbulence in a rotating frame of reference. Specifically, direct simulations of this flow (Bardina et al., 1985; Speziale er al., 1987) indicate that an increasing rate of frame rotation has a pronounced effect in that it reduces the turbulent dissipation e, causing the turbulent kinetic energy I: to persist over a longer period of time. This flow is governed by a system of two ordinary differential equations, one each for k and e. The equation governing the turbulent energy (see Eq. (2.2)), however, is exact and has no explicitly frame-dependent terms. Therefore, this increased persistence of turbulent energy must be localized to the modeled equation for the turbulent dissipation (see Eq. (2.4)). It is thus suggested that the e-cquation be revisited and the characteristic time scale 1: D be modified to include frame-dependent effects. 170 7.2 Further Application The work contained herein was primarily aimed towards achieving Objectives (1)-(iv) in Section 1.2. That is, theoretical problems were selected such that each of these four issues could be isolated or built upon based on previous results within the dissertation. For example, the flows of isotropic decay and return-to-isotropy allow the determination of the model parameters CD and C , respectively, independent of all other model parameters. With this basis, the application of the APS-theory to homogeneous shear flow then allows the final calibration of model parameters to achieve the goal of accurately representing both the primary and secondary normal stress differences. Objective (v), i.e. universality, however, is not addressed directly in this work. This question can be answered through the further application of the APS-theory to other test flows with a subsequent comparison of the predictions with known results. Of course, the extent to which the mean field and turbulent quantifies agree with experimental values is important. However, other questions also arise. For example, for a given flow, what restrictions, if any, are required to ensure realizability or to ensure that the pre-closure operator 11 (see Eq. (4.1)) is invertible. Of course, an additional test of universality is the application of the theory to a practical application of engineering significance. Since this study was, in part, motivated by the poor performance of the Boussinesq theory in flows with streamline curvature, a cyclone would be an ideal case for testing the developed model. The cyclone as a centrifugal separator has a wide variety of separation and classification functions. A generic cyclone with its relevant dimensions is shown in Figure 7.1. A classical application for a cyclone is the separation of particulate solids from a gas or liquid medium. Another example is the use of a hydrocyclone for liquid/liquid separations. For this application, there generally exists a relatively small centrifugal 171 Feed: contains the continuous phase and the dispersed phase to be separated. DI—I Ir— |¢-——1.1 ’I‘ Ira—F Underflow: contains the continuous phase and a majority of any heavy dispersed phase. Overflow: contains the continuous phase and a majority of any light dispersed phase. Hydrocyclone Top View: Tirngentid entry creaes a strong swirling flow, generating a strong centrifugal force. {—Feed Stream Figure 7.1 Qualitative Sketch of a Hydrocyclone Centrifugal Separator 172 driving force due to the small density difference between the two liquid phases (typically ~ 0.2 g/ml). While it is generally true that the ability to accurately predict the interior flow patterns may aid in cyclone design, this capability is more important in the liquid/liquid separator. The large density difference in a solid/fluid cyclone can yield acceptable separation efficiencies for a wide range of operating parameters. Because this density difi‘erence is so small for the liquid/liquid application, however, it is important to be able to characterize the interior flow within the hydrocyclone in order to potentially maximize the separation driving force provided by the interior flow. There are several characteristics common to cyclone flows. In general, they are all high Reynolds number flows with a significant swirl component. This swirl flow is a Rankine-type vortex with a forced vortex structure along the cyclone core and a free-like vortex structure in the outer region. Although the flow may be adequately characterized as axisymmetric, an undulating axial velocity profile combined with toroidal recirculation zones makes for a fully three-component flow field. Despite the challenges that swirling flows present towards modeling, there exists a need for a practical engineering model which can accurately predict not only mean field quantities, but also turbulence quantities. An accurate turbulence model is instrumental in providing much of the information required to predict hydrocyclone separation efficiencies. Specifically, a converged solution to a turbulent flow problem provides both information concerning the mean field velocities and pressure as well as the six components of the Reynolds stress. The role of these terms becomes apparent when one examines the equation governing the acceleration of a dispersed phase particle (m a Lagrangian frame; Brodkey, 1967): d? "’r _ _ m4; = 3nud(y-g)-B:Vp+? ————vtVu . (7.2) In the above equation, it is the dispersed phase velocity and it the continuous phase 173 velocity. Similarly, md = (7rd3pd) /6 and me = (ndspc) /6 represent, respectively, the mass of a spherical volume of diameter d for the dispersed and continuous phase. Eq. (7.2) relates the acceleration of the dispersed phase particle acceleration to the Stokes’ drag force, the fluid pressure gradient, and the added mass effects. Eq. (7.2) does not include the Basset term, nor does it account for particle-particle interactions or drop coalescence/break-up/deformation. The application of Eq. (7.2) allows one to compute dispersed phase trajectories as a function of particle size which further permits the prediction of separation efficiencies given an inlet feed size distribution. Eq. (7.2) is written for the instantaneous particle and continuous phase velocities, but a turbulence model yields mean field quantities. There are two ways to proceed from this starting point. Integration of the particle trajectory equations several times for a given initial position and particle size would allow the formation of an ensemble mean collection probability. In this case, the instantaneous fluid velocity at a given location for a given component is computed in an ad hoc manner as the mean plus a random fluctuation proportional to the rms value of the fluctuations. This is the approach of Boysan, et al. (1982) and Hargreaves and Silvester (1990). An alternate method entails the Reynolds averaging of Eq. (7.2) in order to form the equation for the mean particle acceleration. In this case, there arises the following term which is not a result of the turbulence model predictions and must therefore be further modeled: (y- V y) = (y)- V (50+ (3' - V11)- (73) From Eq. (7.2), the role of the individual mean velocity components upon the separation is evident. Among the three mean velocity components, a good quantitative prediction of the tangential velocity (”9) is key, as the separation force is proportional to the swirl velocity, squared. The axial velocity (uz) is also influential in that it largely determines the residence time of a dispersed phase particle in the cyclone which correlates 174 with the time available for separation. The radial velocity (u r) is important in that it is the reference velocity against which the drift velocity is computed and can also have an effect on recirculation flows that can either work for or against the separation. With accurate computations of the complete mean velocity field, one can estimate dispersed phase trajectories as a function of particle size in order to evaluate a cyclone’s performance given a dispersed phase of a specified density and size distribution. Determination of the mean pressure distribution in the cyclone facilitates the computation of the power requirement for pumping the fluid through the unit operation. APPENDICES 175 Appendix AzNavier-Stokes Equations in a Rotating Frame of Reference The governing equations which are applied in turbulent flow situations stern directly from the instantaneous continuity equation and the Navier-Stokes equations. These are based on a mass and a momentum balance on an infinitesimal fluid element, respectively (Bird, et al., 1960). In a rotating frame of reference, the terms representing the centrifugal and Coriolis forces which arise may be referenced from Greenspan (1968): v.5 = 0 (A.1) and Du 2 -—'+20xa = -VP +vV u. (A.2) Dt - - ' ’ ‘ In Eq. (A.2), the substantial derivative is the instantaneous substantial derivative, i.e. D/Dt E a/3t+ it - V. In the Navier-Stokes equation (Eq. (A.2)), it is important to note that body forces (such as gravity) are not shown. Secondly, a homogeneous density is assumed, such that P, is actually the kinematic reduced pressure, defined by P = (’3’): (Ia-$6229) . (9x5). (A.3) Eqs. (A1) and (A.2) can be written in terms of mean field quantities by introducing the Reynolds decomposition. Explicitly, for the pressure and velocity, this decomposition is y = < )+u' (A.4) 176 and ’p =

+p'. (45> These expressions may be subsequently inserted into the instantaneous equations and then averaged as per Hinze (1975). This results in the final governing equations: the mean continuity equation V - (g) = 0, (A.6) and the Reynolds-averaged equation of motion DQ‘) 2 . . t +29x(g) = -V(Pr)+vV (y)—V-(1_¢t_4). (A.7) In Eq. (A.7), the term D/Dt now represents the substantial time derivative associated with the mean velocity field: D/Dt53/3t+ (g) - V. The last term in Eq. (A.7) is the Reynolds stress. These equations represent a system of four scalar equations for ten unknowns: mean velocity (3), mean pressure (1), and the symmetric Reynolds stress tensor (6). It is this unclosed nature of the governing equations which necessitates turbulence modeling. That is, given a phenomenological or theoretical expression for the Reynolds stress in terms of mean field quantifies, Eqs. (A.6) and (A.7) would be closed. In this example, a transport equation has been deve10ped for the mean velocity, but contains unspecified second-order moments. In general, were an equation for an arbitrary rim-order moment to be developed, it would contain unknown moments of order (n+1). The derivation for the transport equation for second-order moments is given in Appendix B. 177 Appendix BzReynolds Stress Equations in a Rotating Frame of Reference The derivation for the equation governing the Reynolds stress ((1_¢'1_¢')) begins with the equation for the fluctuating velocity (1_i'). This equation is formed by subtracting the Reynolds-averaged equation from the instantaneous (i. e. N avier-Stokes) equation: Du. pr 2 3.7-+29)“; = _V;+vv 51+V. (Etgt)_gtovgr_gr.v(y)o (B,1) The fluctuating velocity equation may then be rearranged with all of its terms on one side and written in general operator form: .£ (1_r') = 0. (13.2) The equation in this form is the source of both the equations for the Reynolds stress as well as the turbulent kinetic energy and the turbulent dissipation (see Appendix C for the dissipation). The second order transport equation is formed by taking the following moment: (gut (g') +4 (105') = 9. (13.3) The resulting equation is D32”! )+ 29 x (E'L‘Ii'l' 2 (9 X (5'50) T = -[(1_4'1_t')- V (19+ ((E'g'). V(E»1:I + —_c- V - (<_'g'g'>+ + 02143)") +VV2(5'I_4')- (B.4) The second term on the LHS of Eq. (B.4) represents a redistribution effect due to the 178 Coriolis terms in the equation for the velocity. On the RHS of Eq. (B.4), the terms are, respectively: (1) turbulent production due to coupling with the mean shear; (2) pressure- strain redistribution; (3) viscous dissipation; (4) turbulent transport; and, (5) molecular diffusion. While this research does not apply the Reynolds stress transport equation directly, it does serve as the basis for the k-equation. The turbulent kinetic energy is defined as 21: E (1_4' - 14'). (B.5) Thus, the contraction of the Eq. (B.4) yields the k-equation. 21: = _:v (y)—e—V- <§<2'>+ <12?» +vV2k- (23-6) Each of the terms on the RHS of Eq. (B6) are the scalar analogs of the terms in Eq. (B.4). An advantage of the k-equation is that the pressure-strain term, which is the focus of most second order modeling approaches, is identically equal to zero due to incompressibility. The only term which is unclosed is the diffusion term, being generally modeled with some form of gradient hypothesis (Hanjalic and Launder, 1972). It is also important to note that, although the Reynolds stress equation is frame dependent, there is no explicit appearance of the frame rotation rate in the kinetic energy equation. This is important, as it implies that rotation only serves to redistribute turbulent energy rather than to create or destroy it. 179 Appendix CzTurbulent Dissipation Equation in a Rotating Frame of Reference The second variable required as a turbulent scale-determining parameter is the scalar dissipation rate. Similar to the Reynolds stress equation, the dissipation equation is based on the equation for the fluctuating velocity. In this case, however, an alternate moment is formed (Speziale, 1991): 2V(V1_4'T:V.£(1_i')) = 0. ((2.1) This yields the following equation for the scalar dissipation rate: De 2 E = Pe-¢8+De+vV a. (C2) In Eq. (C.2), the fourth term on the RHS is the molecular transport term, while Pa, the, and De represent production, destruction, and turbulent diffusion of dissipation. Their exact expressions are P = -2v(Vr_¢'T- vg' +Vl_t' - vr_r°2):v (y)—2Vtr[(1_r'Vr_t'): (VVg')T] -2v([V1_¢I: (v_r_r' - Vr_¢') 7]), (c3) (D = 2v2tr(V Vg': (V Vg') T), (C.4) D = —vv. (re (Vr_¢'-V1_1'T))—2VV- (vp'-Vr_r'). (c.5) In the dissipation equation, there is no explicit frame dependance, providing that the dissipation is isotropic. In contrast to the Reynolds stress and k-equations, virtually all of the terms in the dissipation equation must be modeled. Typically, what is done is to I80 assume that each of the given processes in the dissipation equation are analogous to their k-equation counterpart, scaled by an empirical coefficient and a characteristic turbulent time scale (Hanjalic and Launder, I972). 181 Appendix DzObjectivity of Convected Time Derivatives The generalized convected derivative is (Denn, 1990; Joseph, 1990): 60 8A T EI93)aft+(y>-vr=t- -§-§- 11 1th u a. ' - - =QT. (D2) (S) is the symmetric, objective portion of the velocity gradient and (_lj’) is the antisymmetric, frame-dependent portion of the velocity gradient (Bird, et al., 1977): (>‘ =__Q-(=)-QTand (13.3) < >‘ =g-< ghgog. . 03.4) where Q is a time dependent, proper orthogonal (i.e. det (2) = +1) coordinate transformation tensor and Q its time derivative: ,9' T :1 (135) "(Q Eq. (DJ) in a rotating frame of reference is written as 5; . a . . . . E"; )=a—I(g )+ Wig. - Q Q A- Q +_Q- g Q’-Q-<§>-Q’]= aLQ 1<§>A+Ar<§» -Q’] (v.10) Combining Eqs. (D.6) - (13.10) yields 8; . 5 T s—tlfi ) = g I—(A)] Q (D11) 8:- and the convected time derivative Sa/St is thus objective for any choice of the parameter -°°231 II 1.76 0.151 -0.139 1.51 0.131 -0.119 1.33 0.126 41.108 1.19 0.116 -0.095 1.07 0.107 41.086 0.966 0.100 -0.077 0.884 0.097 -0.068 0.784 0.090 -0.067 0.704 0.094 -0.069 0.636 0.090 -0.062 0.580 0.086 -0.054 I Data not available; only normal components were measured. 195 Table E12: Retrrm to Isotropy Data: Choi and Lumley [1984] (Axisymmetric Expansion) 2HU2x103 822 1,33 b231 11 0.665 0.196 -0.117 0.613 0.189 41.112 0575 0.185 -0.110 0.541 0.175 .0103 0.510 0.177 41.101 0.483 0.171 -0.099 0.457 0.163 -0.091 0.427 0.163 -0.093 0.400 0.163 -0.094 0.376 0.160 -0.093 0.356 0.154 -0.091 0.339 0.150 -0.091 0.324 0.146 -0.083 I Data not available; only normal components were measured. _ 196 Table E.l3: Retum to Isotropy Data: LePenven et al. [1985] (Positive Third Invariant) 21/U2x 103 1,22 1233 0.238 0.183 41.045 0.229 0.181 -0.046 0.216 0.179 -0.049 0.207 0.174 -0.050 0.193 0.167 -0.053 0.182 0.164 -0.049 0.167 0.159 -0.052 0.159 0.158 -0.054 0.148 0.154 -0.052 0.139 0.154 -0.055 0.134 0.150 -0.053 I Data not available; only normal components were measured. 0.0532 0.0516 0.0480 0.0439 0.0426 0.0396 0.0387 0.0369 0.0364 0.0349 197 Appendix F: Relation of Homogeneous, Isotropic Turbulent Production to the Velocity Derivative Skewness Eq. (2.31) presents a relationship between the velocity derivative skewness and the Taylor microscale, along with properties of the double and triple velocity correlations. The velocity derivative skewness is The mean of the velocity derivative squared can be related to the second derivative of the double velocity correlation (see p.145, Monin and Yaglom, 1965): (GI—fl) = {325?} . (E2) r r=0 Eqs. (2.23) and (F.2) can be combined further to yield the following result, namely, ‘ (zen-=1 Similarly, the mean of the velocity derivative cubed is related to the third derivative of the triple longitrrdinal velocity conelation (see p.145 Monin and Yaglom, 1965): (8%.)3) = [33:2LL) 0. (E4) r r: It follows directly from Eqs. (El), (E3), and (E4) that 7.3 3,2 (33:?) =(@_: .3<))/[ (g: I2))]3/2§_SP (1:5) (Enron r=o 198 which is the result expressed by Eq. (2.31). Furthermore, the combination of Eqs. (ES) and (2.28) yield the relationship between turbulent production in homogeneous, isotropic turbulence and the velocity derivative skewness (cf Eq. (2.13)): w = ———S,JR_e. (E6) H9 Appendix GzAsymptotic State for Homogeneous Shear This appendix details the algebra required to express the asymptotic values of the anisotropic pre-stress components for homogeneously sheared turbulence (as given in Eqs. (4.54)-(4.57)). The asymptotic values for the normal components of the anisotropic pre-stress follow directly from Eqs. (4.45)-(4.47) by setting the derivatives with respect to § equal to zero: (1+ q“) H; —§ar)e,“1t1‘y’z = o, ((3.1) a a a (1+q )H;y+ (5+1)Det H; = O, (6.2) and (1+ q“) ng+ (g —1)Det"H;‘Z = 0. (G3) Rearranging Eqs. (G. 1)-(G.3) yields the relations given in Eqs. (4.54)-(4.56): De “H“ ”1‘. = a ‘ ". (6.4) 3 l + q“ 0 Def”; [fir-(1+3) ”q“, ((3.5) and De “H“ a : yz = 1 - — —. ((3.6) u ( 3) l + q“ Similarly, Eq. (4.72) combines with Eqs. (6.5) and (6.6) to give the expression for the asymptotic value of the shear component of the anisotropic pre—stress: 2(1) C a De a 2 De 0 O = 1J_£1_:_Z_).H‘z- ' H:z+-a— ‘ [1:2, (0.7) 2“; De,“ ’ (1+q“) 3 (1+q“) whence, upon rearrangement, becomes 11" C5 (1+ qa)De‘a yz = 2 (0.8) C1(1+qa)2+ (Det“)2(l-a2/3) ’ which is identical to the asymptotic result given in Eq. (4.57). 201 Appendix H:Computer Program Listings H.1 Program TRANS . BAS The program TRANS . BAS is written Microsoft QuickBasicT“. Its function is to integrate the initial value problem posed for the case of homogeneous shear flow. For the IPS-theory, the differential equation for N R is solved, as are the algebraic pre-closure equations for 5. With the extension of the model to include an anisotropic pre-stress, four additional differential equations are solved for the components of 51. Initial conditions as well as program control parameters are read as input through the file HOMG_SHE . INI. Transient output data is directed to a user-specified file and the final time step is output to the file SUMMARY . TXT. Variable Listing a convected derivative parameter a ak or/k B holding matrix for Gauss-Jordan inversion Cd , Cp C D’ C 1, cr, cb, c1 CR, CB’ C21 De relaxation group N R Deb, Del (CB/ZCRMI ; (CA/CR) NR Det dNR/d§ dt , fdt; Afi, variable step size parameter 9 Q/ S H , Hold anisotropic pre-stress _I1_I, holding variable for old values Ht: (lg/d: ‘- q, qq parameter q (see Eq. (4.43)); (l + q) /Det 202 R normalized Reynolds stress 5 ske S k/ a t , told dimensionless time é, holding variable for old time Program Listing REM full relaxation model with H + lambda'lIl/Dt - beta' (0/0t: Mixed Objective Derivative) REM full inversion of coefficient matrix to solve R given it and NJ REM transith calmlati on ' Variables DIM DIII DI)! DIM DIM DII DIM DIN DIM DIM 1H3, 3) AS DOUB 3(3. 4) AS DOUB H(3. 3) AS DOUB L8 LE LE Hold(3, 3) AS DWBLE dt, t, told, tdt AS 00181.8 'tine variables 0e. Deb. Del AS DWBLE Deold AS 00131.8 at A8 03081.8 eke as 00131.3 9 AS 00181.8 'reynolds stress 'transtornation tensor to be inverted 'anisotropic pre-stress 'anisotropic pre-stress ' Deborah nubers ' alpha lk 'Slt/eps '(hega/s (relative rotation rate) ' Tine derivatives or variables DI)! DI)! DI! DIM DIM DIM DIM DIM DIX DIN lit-(3. 3. 4) AS NOBLE Dot(d) A8 WUBL cr. cb. c1 AS 00181.8 Cd, Op AS nouar. q, qq AS 00181.8 a AS norms i. j. k. 1 A8 INTRIER a. n AS tour: but! as s'rnmo mdel AS swarm 'anisotropic pre-stress tine derivative 'Deborah number tine derivative 'rr nodal constants 'eps-eqn constants ' conv. derivat ive par-eter 'identitier tor IPS/APS node] ' Read initial values for H, De: OPE-8 'd:\data\rnsu\hang_she\hong-she.ini' FOR INPUT AS 01 'read model type (IPS/APS) LINE INPUT Ii, butts: aodels - ceasesuarruturrs, 3)) 'read initial H LINE INPUT .1. FOR i s 1 '10 3 nuts INPUT .1. "(1. 1): INPUT .1. “(1, 2): INPUT I1. H(i, 3) DEXTI 'read other para-stars LIB INPUT .1. LIFE INPUT .1. L138 INPUT I1. L118 INPUT .1. LIB INPUT .1. CLCBE .1 ' Specify nodal paraeters tufts: N115: 1:1!th NTTSr halts: SELECT CASE ”601$ C538 'IPS' er I .227. Ch I 0. C1 I 0. Cd 3 1.83. Op - 1.d1. INPUT .1. De INPUT .1. 9 INPUr Ii, dt. tdt, taax INPUT l1, mrint INPUi‘ 01. outtiles a I 0. CASE 'APS' CI I .2706. ch I .17‘2. cl I . 66666666666666660 Cd I 1.83. q; - 1.604. a I —. 6788999999999999. CASE ELSE PRINT 'No valid nodel epecitied; program aborted': BID BID SELECT ' Open data tile for output OPEN 'd:\data\lum\hang_she\' + (outtileS) FOR 0J'I‘PUT AS 01 PRINT .1. [SING 'MJM. ' Begin tine integration CLS t-O. DIO 'deternine initial Reynolds stress IF (”601$ I 'IPS') THEN “(1. 1) I 0I1H(1. 2) I 0O:l~l(1. ”(2, 1) I 09: "(2, 2) I 0.: 8(2, N(3. 1) I 0018(3, 2) I 0|: 1H3. 3) BID IF (DNS 1000 '4th order RK loop 00 ' output I? ((n / nprint) I INT(n / nprint)) 2): N2. 3) '1 13: PRINT 1:: De: 11(2) PRINT .1. 08116 “9.0.0... .: g PRINT ll. usrm 'DJIMI. mo IF 'aave old values Deold I DO Holdu, 1) I 11(1, 1) bid”. 2) I N”: 2) 8016(3, 3) - H(3, 3) Hold(2, 3) I 8(2. 3) told-t '1st integration step 3) 3) ': DI: .i R(10 'I "(10 I 0. I 0. I 0. max 1): 812. 2): M3. 3): R12. 3): 1): H12. 2): H13. 3): H(2. 3) 0812(1) I (2. ' DI ‘ R(2. 3) ' (Cp - 1.)) 4» (cr ' (Cd - 1“) De I Doold + .5. ' dt ' Det(1) IF (”691$ I 'APS') THEN Del - (cl lcr) ' 0e q.0el'(-20'R(2.3)-cr/0e) qq-(lI+q)/0el mu. 1. 1) I -Q ' H(1. 1) 4 (2| ' l / m(20 20 1) - 'm . "(2! 2) ":13: 3o 1) ' 'Qq . "(3r 3) ' “(20 3) ' + (Cl) / cl I Ht(2. 3. 1) I -qq ' 1H2, 3) H(1. 1) I Noldtl. 1) H(2. 2) I Nold(2. 2) 11(3, 3) I H01d(3. 3) 8(2, 3) I Nold(2. "(3, 2) I "(2' 3) 310 IF 3) 4 a + '0- .SI .5. .SI .5. 0 t O t dt (It (It dc - “(2. O O O O Ht(1. Nt(2. HI:(3. Ht(2. ‘deternine Reynolds stress for step i @918 1000 ' 2nd integration step 3)‘ 3|) ‘ “(2, 3) (I / 30 + 1!) (a / 3| - 10) 2.) + .50 ' H(3. 3) ' (-1. - a) + .50 " H(2, 2) ' (10 - a) 1. 1) 2, 1) 3, 1) 3, 1) Dat(2) I (2. ' D. ' R(2. 3) ' (Cp - 1.)) + (or ' (Cd - 19)) D. I Deold + .5. ' dt ' 0011(2) 1? (models I 'APS') THEN Del I (cl I or) ' De q I Del ' (-2| ' R(2. 3) - cr / De) qq - (1. + Q) I 001 Ht(1, 1. 2) I -qq ' Hti, 1) I (2| ' a / Ht(2. 2, 2) I -qq ' H(2, 2) - H(2, 3) ' Ht(3, 3, 2) I -qg ' H(3, 3) - H(2, 3) ' Htt2. 3. 2) I -qq ' H(2, 3) I (Ch / c1 / H(1, 1) I Hold(i, 1) I .50 ' dt ' Ht(1, H(2, 2) I Hold(2, 2) I .50 ' dt * Ht(2, “(3. 3) I Hold(3. 3) I .S| ‘ dt ' Ht(3. “(2, 3) I Hold(Z, 3) I .5| ' dt ' Ht(2. “(3. 2) I H(2. 3) END I? 'determine Reynolds stress for step 2 (13803 1000 '3rd integration step 3|) ' “(2. 3) (a / 3| I 1|) (a / 3| - 1|) 2|) 1. 2. 3. 3. + 2) 2) 2) 2) Det(3) I (2| ' De ' R(2. 3) ' (Cp - 1|)) I (cr ° (Cd - 1|)) De I Deold I dt ' 0et(3) IF (Ind-IS I 'APS') THEN Del I (cl / cr) ' De q I Del ' (-2| * 8(2. 3) - cr / De) aq-(1|+q) “t(1. 1. “t(2. 2. “t(3. 3. “(1. “(2. “(3. “(2, “(3. END IF 3. 1) 2) 3) 3) 2) 3) 3) 3) I “016(1. I “016(2. I Hold(3. I “016(2, I "(2, 3) / Del 1) 2) 3) 3) ++++ I -qq * H(1. 1) I (2| ' a 3) I -qq ' H(2. 2) - H12. 3) I -qq ' H(3, 3) - H(2, 3) I -qq 0 th. 3) + (ch I c1 ' “t(1. 1. 3) ' “t(2. dt dt dt * “C(3. 3. 3) dt 2. 3) ' “t(2, 3. 3) Idetermine Reynolds stress for step 3 00308 1000 '4th integration step / 3|) ‘ H12. 3) ‘ (a / 3| I 1|) ' (a / 3| - 1|) / 2|) I .5| ' “(3. 3) ' (-1| - a) I .S| ' “(2. 2) ° (1| - a) Det(|) I (2| ’ DI ' 3(2. 3) ‘ (Cp - 1|)) I (ct ' (Cd - 1|)) 1? (models I 'APS') THEN Del I (cl I or) ' De q I Del ' (-2| ' R(2, 3) - cr / De) qq I (ll I q) / Del “t(1. 1. “t(2. 2. Ht(3p 3' d) I -qq ' “(1. 1) I (2| ’ I / 3|) ' “(2, 3) 4) I -qq ' “(2. 2) - “(2. 3) ' (a / 3| I 1|) 4) I —qq ' “(3. 3) - “(2. 3) ‘ (a I 3| - 1|) “t(2. 3. d) I -qq ' “(2. 3) I (Cb / cl / 2|) I .5| ' “(3. 3) ’ (-1| - a) I .S| END IF 'dth-order RK-step n I n I 1 t I t I dt De I Doold I (dt I 6|) ' (Data) I 2| ' 00!:(2) I 2| ' Det(3) I Det:(d)) I? (”6015 I 'APS') THEN H11. 1) I Hold(l. 1) “(2, 2) I Hold(2, 2) H(3, 3) I Hold(3. 3) H(2, 3) I Hold(2. 3) “(3, END IF 2) - “(2. 3) O + + + (dt I 6|) (dt I 6!) (dt I 6|) (dt / 6|) 0 O O 0 (“C(1. ("C(2. (“C(3. ("C(2. 'deternine Reynolds stress for step 4 00808 1000 at I dt ' fdt traceH I “(1. 1) I “(2, 2) I “(3. 3) 1. 2. 3. 3. 1) 1) 1) 1) I 2| I 2| I 2| I 2| 0 O O 0 “C(1. “C(2, “t(3. “t(2. 1. 2) 2. 2) 3, 2) 3, 2) 0+++ 2| 2| 2| 2| ' “t(1. ' “t(2. ' “t(3, ‘ “t(2. .S| ‘ “(3. 3) ' (-1| - a) I .S| ' “(2. 2) ' (1| - a) ' “(2. 2) ' (1| - a) 1. 2. 3. 3. 3) 3) 3) I “t(1. I “t(2. I “C(3. I Ht(2. 205 I? (ABS(traceH) > .000001) 'iHI-Jl PRINT “)1 not traceless: '3 traoeH EXIT DO 340 IF LOOP UNTIL. (t > tmx) ' final output PRINT 1:: De: R(2. 2): R(2. 3) PRINT |1. USING '|||.|||. ': t: PRINT |1. [SING '||.|||||. ': De: PRINT |1. USING '|.|||||. '1 RM, 1): R(2. 2): M3. 3): R(2. 3): PRINT |1, (SING '|.|||||. '3H(1, 1): H(2, 2); H(3. 3); H(2, 3) CLCBE |1 ' Write asymptotic results to a eta-nary tile OPai 'd:\data\msu\hmg_she\sunary.txt' FOR APPEND AS |1 PRINT |1, 081K: '||.|||. '1 g: PRINT |1, [BIN '|||.|||. ': t; PRINT |1. usrno '||.|||||. ': De: PRINT II. 0811!: '|.|||||. ': R(1. 1): R(2. 2): R(3. 3)1R(2. 3): PRINT |1, usruc '|.|||||, '.H(1, 1); H(2, 2): H(3. 3): “(2. 3) CLCSE |1 END 1000 ' canpute R Iran H ° Comute initial 8 (coefficimt matrix) 8(1. 1) I 3| - ((1| I g) * D) “ 2| 3(1. 2) I 2| ' (g ' De) “ 2| 8(1. 3) I -(6| ‘ g I 2|) ' De R(2. 1) I 2| ' ((1| I g) ' De) A 2| R(2. 2) I 3| - (g ' 0e) “ 2| R(2. 3) I (6| ' g I 4|) ' De 8(3. 1) I (1. + 9) ' 0. 3‘3, 2) I -g . D. B(3. 3) I 1| - a ' (1| + 9) ' De A 2| ' Corpute initial 8 (known vector) 8(1, 4) I 1| I 3| ' H(2, 2) R(2. 4) I 1| I 3| ' H(3. 3) 3(3. |) 8 “(2. 3) ' Gauss-Jordan Elimination FOR k I 1 'IO 3 roe j - 4 '10 x crap -1 NR. 1) - 30!. 1) l 30:. k) NEXT) FOR 1 I 1 '10 3 IF (i o k) THEN FOR) I d m 1: STEP -1 3(1. 1) I 8(1. 1) - (8(1. k) ' 30!. 1)) EXT) END IF NEXT 1 NEXT k ' Note: post-inversion. 1st three colums are identity matrix, ' and the 4th colt-n is the solution vector. R(2. 2) M3. 3) NI. 1) R(2. 3) 0(1. 4) R(2. 4) 1| - R(2. 2) - R(3. 3) 3(3. 4) ' Verify realizability of R: I? (R(2. 3) " 2 > R(2. 2) ' R(3, 3)) MN PRINT 'Schumrs inequaltiy not satisfied.': END 310 IF IF(R(1. 1) < 0) THEN PRINT 'Rxx negative': END mo IP IP (R(2. 2) < 0) THEN PRINT 'Ryy negative“: END END IP 1? (NJ. 3) < 0) mm PRINT 'Rsz negative': END 206 END IF R(3. 2) I R(2. 3) RETURN H.2 Program IPS_ROT. BAS The program IPS_ROT. BAS is written in Microsoft QuickBasicm. It solves the algebraic pre-closure equations for _R given values for NR and Q/S. Program control parameters are read as input through the file IPS_STAT . INI. Output data is directed to a user-specified frle. Variable Listing The variables used are the same as listed in Section H.l. Program Listing REH IPS-model calwlations tor turbulent states as a tunction ct mega/S REH tull inversion of coefficient matrix to solve R givm H and N_R ' Variables DIH R(3. 3) AS DOUBLE 'reynolds stress DIN 8(3, 4) AS DOUBLE 'transtormation tensor to be inverted DIH H(3. 3) AS DOUBLE 'anisotropic pre-stress 0111 De AS mUBLB 'Relaxation Group 01)! Dei. De2 AS 00181.8 'Relaxation Group, initial/final DLH DeDel AS 00181.8 'delta Relaxation Group 01)! g AS 00181.8 'anega/S (relative rotation rate) 0111 cr, cb, cl AS 00181.8 'rr model constants 0114 Cd. Cp AS mUBLE 'eps—eqn constants 0111 i, j, 1:, 1 AS INTEGER 01)! m. n AS tom 01)! but! AS STRIN: DIN nstep AS INTESER 01)! model AS s'nunc 'identitier for IPS/APS model ' Read initial values for gamma: OPEN 'd:\data\msu\horng_she\ips_stat.ini' FOR INPUT AS |1 'read other parameters LINE INPUT |1. butts: INPUT N, g LINE INPUT |1, butts: INPUT |1. Del. De2, DeDel. {De LINE INPUT l1. buffs: INPUT |1, outfileS CLOSE |1 ' Specify model parameters ' CASE 'IPS‘ er I .227| Ch I 0| Cl I 0| Cd I 1.83| 0p I 1.41| a I 0| “(1, 1) I 0|: “(1, 2) I 0|: “(1, 3) I 0| “(2, 1) I 0|: “(2, 2) I 0|: “(2, 3) I 0| “(3, 1) I 0|: “(3, 2) I 0|: “(3, 3) I 0| ' men data file for output OPEN 'd:\data\msu\hang_she\' I (outfileS) FDR OUTPUT AS |1 PRINT |1, USING °||.|||. ': a DeIDei DO GOSUB 1000 PRINT De: R(2. 2): R(2. 3) PRINT |1. USING '||.||||. ': De: PRINT |1, USING '|.|||||. ') R(l, 1): R(2. 2): R(3, 3): R(2. 3) DeIDeIDeDel 03001 I DeDel ' fDe LOOP UNTIL (De > D02) CLOSE |1 END 1000 ' canpute R fran H ' Compute initial 8 (coefficient matrix) 8(1. 1) I 3| - ((1| I g) ‘ De) “ 2| 8(1, 2) I 2| ' (g ‘ De) “ 2| “(1, 3) I -(6| ' g I 2|) ' De 8(2. 1) I 2| ' ((1| I g) ' De) ‘ 2| 8(2. 2) I 3| - (D ‘ DO) “ 2| 8(2. 3) I (6| ‘ g I 4|) ' 0e 8(3, 1) I (1| I q) ' De B(3,2)I-g'0e 8(3.3)I1|-g'(1|Ig)'De‘2| ' Coupute initial 8 (known vector) 8(1. 4) I 1| I 3| * H(2, 2) 8(2. 4) I 1| I 3| ' H(3. 3) 8(3, 4) I “(2. 3) ' Gauss-Jordan Elimination FOR 1: I 1 '10 3 FOR 1 I 4 '10 1: STEP -1 30!. 1) I 30!. 1) / 30!. k) NEXTj FOR 1 I 1 '10 3 I? (i o 1:) THEN FOR 1 I 4 To 1: STEP -1 3(1. 1) . 3(1. 1) - (8(1. k) ' 8(k. 1)) NEXTj END 11? NEXT i NEXT 1: ' Note: post-inversion, 1st three columns are idmtity matrix. ' and the 4th colunn is the solution vector. 208 R(2. 2) I 3(1, () R130 3) . 3(2: 0) R(1. 1) I 1| - R(2. 2) - R(3. 3) R(2, 3) I 3‘3: ‘) Verify realizability of R: IF (R(2. 3) ‘ 2 > R(2. 2) ' R(3. 3)) THEN PRINT 'Schuarz inequaltiy not satisfied.': END END IF IF (R(l, 1) < 0) THEN PRINT 'Rxx negative': END END IF IF (R(2. 2) < 0) THEN pam'r 'Ryy negative': END END IF IF (R(3. 3) < 0) THEN PRINT 'Rzz negative': END END IF R(3. 2) - R(2. 3) RETURN 209 Appendix I: Properties of the Eigenvalues of the Reynolds Stress As mentioned in Chapter 1, the Reynolds stress (g'g') has the following properties: (5' 14') = 2k. (L1) and :0 (2'2) ' .z 2 0- V (12) More specifically, the trace of the Reynolds stress is twice the turbulent kinetic energy and the Reynolds stress is positive semi-definite. Since the normalized Reynolds stress is defined as (u'u') R E ’ ’ , 3 = 2k (1 ) Eq. (1.1) implies that the trace of the normalized Reynolds stress is unity: tr (5) = 1. (14) However, the trace of a tensor operator is also an invariant property of that operator which equals the sum of the eigenvalues: "(5) =1 =AR1+1R2+XRT (1.5) where AR i is the ith eigenvalue of the operator 5. The positive, semi-definite nature of the Reynolds stress expressed by Eq. (1.2) implies that each of the eigenvalues of the Reynolds stress (g'g') is non-negative. This character yields the following condition upon the values 1R i: 210 031mm. (1.6) since each of the eigenvalues is non-negative and the sum of the eigenvalues is one. The condition given by Eq. (1.6) yields information concerning the eigenvalues of the anisotropy tensor _h, defined as 9= MI. (L?) It follows directly from Eq. (1.7) that l 7m..- = ’m's’ (1.8) where Ab, i is the ith eigenvalue of the operator {2. Combination of Eqs. (1.6) and (1.8) yields the result stated in Chapter 1: ' l 2 ‘39:»,553- (1.9) LIST OF REFERENCES I. Bardina, J. Ferziger, and R. Rogallo, 1985, “Effect of rotation on isotropic turbulence: computation and modeling,” J. Fluid Mech., 154, 321. G. Batchelor and A. Townsend, 1948, “Decay of turbulence in the initial period,” Pmc. Roy. Soc., A194, 538. G. Batchelor and A. Townsend, 1948, “Decay of turbulence in the final period,” Pmc. Roy. Soc., A196, 527. J. Bennett, and S. Corrsin, 1978, “Small Reynolds number nearly isotropic turbulence in a straight duct and a contraction,” Phys. Fluids, 21(12), 2129. R. Bird, A. Armstrong, 0. Hassager, 1977, The Dynamics of Polymeric Fluids, volume 11, McGraw Hill. R. Bird, A. Armstrong, 0. Hassager, 1987, The Dynamics of Polymeric Fluids, volume I, McGraw Hill. F. Boysan, W. Ayers, J. Swithenbank, 1982, “A fundamental mathematical modeling approach to cyclone design,” Trans IChemE, 60, 222. R. Brodkey, 1967, The Phenomena of Fluid Motion, Addison Wesley. B. Carnahan, H. Luther and J. Wilkes, 1969, Applied Numerical Methods, John Wiley and Sons. K. Choi, 1983, “A study of the return to isotropy of homogeneous turbulence”, Ph.D. Dissertation, Cornell University. K. Choi and J. Lmnley, 1984, “Return to isotropy of homogeneous turbulence revisited,” in Turbulence and Chaotic Phenomena in Fluids, ed. T. Tatsurni, New York: North Holland. G. Comte-Bellot and S. Corrsin, 1971, “Simple Eulerian time conelation of full- and narrow-band velocity signals in grid generated, isotropic turbulence,” J. Fluid Mech, 48, 273. M. Denn, 1990, “Issues in viscoelastic fluid mechanics,” Ann. Rev. Fluid Mech, 22, 13. 211 212 M. Gibson and V. Kanellopoulos, 1987, “Turbulence measurements in a nearly homogeneous shear flow,” Conference on Transport Phenomena in Turbulent Flows: Theory, Experiment, and Numerical Simulation. H. Greenspan, 1968, The Theory of Rotating Fluids, Breukelen Press. K. Hanjalic and B. Launder, 1972, “A Reynolds stress model of turbulence and its application to thin shear flows,” J. Fluid Mech., 52(4), 609. K. Hanjalic and B. Launder, 1976, “Contribution towards a Reynolds stress closure for low Reynolds number turbulence,” J. Fluid Mech., 74(4), 593. K. Hanjalic, 1994, “Advanced turbulence closure models,” Int. J. Heat and Fluid Flow, 15(3). 178. J. Hargreaves and R. Silvester, 1990, “Computational fluid dynamics applied to the analysis of de-oiling hydrocyclone performance,” Trans IChemE, 68(A), 365. J. Hill and C. Petty, 1996, "Turbulent transport of a passive scalar field,” Chem. Eng. Commun, 152-153, 413. J. Hinze, 1959, Turbulence, McGraw-Hill. M. Huang and A. Leonard, 1994, “Power-law decay of homogeneous turbulence at low Reynolds numbers,” Phys. Fluids, 6(11), 3765. D. Joseph, 1990, Fluid Dynamics of Viscoelastic Fluids, Springer Verlag. T. von mm and L. Howarth, 1938, “On the statistical theory of isotropic turbulence,” Proc. Roy. Soc., A164, 192. O. Kitoh, 1991, “Experimental study of turbulent swirling flow in a straight pipe,” J. Fluid Mech., 205, 445. J. Kim, P. Moin, and R. Moser, 1987, “Turbulent statistics in fully developed channel flow at low Reynolds numbers,” J. Fluid Mech., 177, 133. B. Launder, D. Spalding, 1974, “The numerical computation of turbulent flows,” Comp. Meth. Appl. Mech. and Eng, 3, 269. B. Launder, G. Reece, and W. Rodi, 1975, “Progress in the development of a Reynolds stress closure,” J. Fluid Mech., 68(3), 537. L. LePenven, J. Gence, and G Comte-Bellot, 1985, “On the approach to isotropy of homogeneous turbulence: effect of the partition of kinetic energy among velocity components,” in Frontiers in Fluid Mechanics, ed. S. Davis and J. Lumley, Berlin: Springer Verlag. 213 J. Lumley, 1978, “Computational modeling of turbulent flows,” Adv. App. Mech., 18, 123. N. Mansour and A. Wray, 1994, “Decay of isotropic turbulence at low Reynolds numbers,” Phys. Fluids, 6(2), 808. G. Mase and G. Mase, 1992, Continuum Mechanics for Engineers, CRC Press. A. Monin and A. Yaglom, 1965, Statistical Fluid Mechanics: Mechanics of Turbulence, volumes 1 and 2, MIT Press. P. Morse and H. Feshbach, 1953, Mefltods of Theoretical Physics, part I, McGraw-Hill. S. Parks, K. Weispfennig, and C. Petty, 1997, “An Algebraic Pre-Closure Theory for the Reynolds Stress”, submitted to Phys. Fluids. V. Patel, W. Rodi, and G Scheuerer, 1997, “Turbulence models for near-wall and low Reynolds number turbulence”, AMA J., 23(9), 1308. C. Petty, 1975, “A statistical theory for mass transfer near interfaces,” Chem. Eng. Sci, 30, 413. W. Reynolds, 1989, “Towards a structure-based turbulence model”, Studies in Turbulence; editors T. Gatsld, S. Sarkar, and C. Speziale; p.76, Springer Verlag. J. Rohr, E. Itsweire, K. Helland and C. Van Atta, 1988, “An investigation of the growth of turbulence in a uniform mean shear flow,” J. Fluid Mech., 187, l. S. Sarkar and C. Speziale, 1990, “A simple nonlinear model for the return to isotropy in turbulence,” Phys. Fluids A, 2(1), 84. U. Schumann, 1977, “Realizability of Reynolds stress turbulence models,” Phys. Fluids, 20(5), 721. A. Sirivat and Z. Warhaft, 1983, “The effect of a passive cross-stream temperature gradient on the evolution of temperature variance and heat flux in grid turbulence,” J. Fluid Mech., 128, 323. C. Speziale, 1987, “On nonlinear H and k-e models of turbulence,” J. Fluid Mech., 178, 459. C. Speziale, N. Mansour, and R. Rogallo, 1987, “The decay of isotropic turbulence in a rapidly rotating frame,” Center for Turbulence Research, Proceedings of the 1987 Summer Program, (Stanford U.P., Stanford, CA, 1987) 205. C. Speziale and N. Mac Giolla Mhuiris, 1989, “On the prediction of equilibrium states in homogeneous turbulence,” J. Fluid Mech., 209, 591. C. Speziale, 1991, “Analytical methods for the development of Reynolds-stress closures 214 in turbulence,” Ann. Rev. Fluid Mech., 25, 107. C. Speziale, R. Raj, and T. Gatski, “Modeling the dissipation rate in rotating turbulent flows,” Studies in Turbulence; editors T. Gatski, S. Sarkar, and C. Spen'ale; p.129, Springer Verlag. D. Taulbee, 1989, “Engineering turbulence models”, Advances in Turbulence; editors W. George and R. Amdt; p.75, Hemisphere. S. Tavoularis, J. Bennett, and S. Cousin, 1978, “Velocity derivative skewness in small Reynolds number, nearly isotropic turbulence,” J. Fluid Mech., 88(1), 63. S. Tavoularis and S. Corrsin, 1981, “Experiments in a nearly homogeneous turbulent shear flow with a uniform mean temperature gradient,” J. Fluid Mech., 104, 311. S. Tavoularis and U. Kamik, 1989, “Further experiments on the evolution of turbulent stresses and scales in uniformly sheared turbulence,” J. Fluid Mech., 204, 457. H. Tennekes and J. Lumley, 1972, A First Course in Turbulence, MIT Press. K. Weispfennig, 1997, “Relaxation/retardation model for fully developed turbulent channel flow”, Ph.D. Dissertation, Michigan State University.