MISCHIGAN STAE I ”ll Nil/IIIHINHlllllllllllllllllllll 1293 01706 9927 I!!!” This is to certify that the dissertation entitled RELAXATION/ RETARDATION MODEL FOR FULLY DEVELOPED TURBULENT CHANNEL FLOW presented by Klaus Weispfennig has been accepted towards fulfillment of the requirements for Ph.D. degree in Chemical Engineering " / ‘ n d‘ Wkigh , M 3.1“ pro fessor Date/DQCQJV—A .2 1:) )l 0’37 MSU u an Allirmunw' Ar mm ’Iz'quul ()ppnrmnilv Institution 0-12771 LIBRARY MIchIgan State Unlverslty PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE ugvgale my“ jfiflBf ‘1 000 AUG 2 5 20“ 061812 1/98 c/CIRC/DateDue.p65—p.14 RELAXATION/RETARDATION MODEL FOR FULLY DEVELOPED TURBULENT CHANNEL FLOW By Klaus Weispfennig 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/RETARDATION MODEL FOR FULLY DEVELOPED TURBULENT CHANNEL FLOW By Klaus Weispfennig A representation for the Reynolds stress is developed in terms of a Green’s function associated with mean field convection and viscous transport of turbulent fluctuations. The finite memory of turbulent temporal correlations is used as an ansatz to develop a relationship between the Reynolds stress and a single-point self-correlation of an effective fluctuating force induced by pressure fluctuations and fluctuations in the instantaneous Reynolds stress. The theory gives an algebraic model as a pre-closure approximation relating the Reynolds stress to the velocity gradient and a statistical correlation termed prestress. A phenomenological relaxation/retardation closure model based on frame- invariant modeling relates the mean strain rate dyadic to the anisotropic part of this prestress. A new time scale is introduced to extend the theory to low Reynolds number regimes such as the near wall region in channel flow. The particular form of this time scale renders the theory realizable for simple shear flows (e.g channel/pipe flow, homogenous shear) independent of the solution of the transport equations for the kinetic energy It and the dissipation rate a. Two distinct regimes of momentum transport are identified, a gradient transport regime close to the channel center and an equilibrium type regime in the inertial sublayer. ~The application of the relaxation/retardation theory in the outer region of fully developed turbulent channel flow renders an overall good agreement with both numeriCal and experimental data. Developments for the near wall region are done and evaluated with respect to direct numerical simulation (DNS) data available. The associated transport equations for the turbulent kinetic energy k and the dissipation rate 8 are employed to determine the distribution of the turbulent time scale. These equations use modified terms for the turbulent transport derived using the same formalism as applied for the preclosure. The parameters introduced through the Reynolds stress model are determined using direct numerical simulation data. A parametric study of the associated model parameters manifest their choice. The parameters of the associated transport equations are recalibrated for the new theory using experimental data from various flow fields. A precursor to the relaxation/retardation theory has been investigated. This isotropic prestress representation yields a non-zero primary normal stress difference and a zero secondary normal stress difference. The wall region is characterized by a one-dimensional turbulence state in which all the energy is transferred into the streamwise normal component. The relaxation/retardation closure yields an energy partition which is in compliance with experimental and numerical data. The near wall region is characterized by an anisotropic two-dimensional state. Relaxation effects become important near the center of the channel for which both theories assume an isotropic state. Retardation retains energy in the spanwise component of the Reynolds stress allowing the existence of a nonzero secondary normal stress difference. In the low (turbulence) Reynolds number regime (i.e. near the wall) retardation prohibits the unbounded growth of this secondary normal stress difference. To my wife and best friend, Maria Victoria iv ACKNOWLEDGEMENTS I wish to express my gratitude to Dr. Charles A. Petty for his guidance and encouragement throughout this study and in preparation of this manuscript. Thanks are also given to Dr. John Foss, Dr K. Jayaraman and Dr Milan Miclavcic for their ideas and assistance during this research as members of my committee. I also appreciate the help of all my friends and all the students who supported me while accomplishing this work. I would like to thank my parents for their support and for making a study abroad possible. My wife, I would like to thank for all the love and patience during the preparation of this manuscript. I gratefully acknowledge partial financial support from the Michigan State University Foundation, the American Plastic Council and the Research Excellence Funds, State of Michigan. Partial support of this work by the Hydrocyclone Development Consortium at Michigan State University is appreciated. TABLE OF CONTENTS LIST OF TABLES viii LIST OF FIGURES ix NOTATION xi Chapter Page 1 INTRODUCTION 1 1.1 Motivation and Background 1 1.2 Objectives and Methodology 27 1.3 Summary Outline of This Research 29 2 REYNOLDS STRESS PRECLOSURE 1N INERTIAL FRAMES 30 2.1 Preclosure for the Reynolds Stress 30 2.2 Turbulent Relaxation Time 34 2.3 Properties of the Preclosure 36 3 PHENOMENOLOGICAL RELAXATION/RETARDATION THEORY 48 3.1 Motivation 48 3.2 Closure Hypothesis for the Prestress 51 3.3 Properties of the Preclosure 55 3.3.1 Approximation for Small Time Scale Ratios 57 3.3.2 Other Phenomenological Models 60 4 TRANSPORT EQUATIONS FOR TURBULENT KINETIC ENERGY AND DISSIPATION RATE 68 4.1 The Exact Equations 68 4.2 The Closure Hypothesis 70 4.2.1 Transport Terms 70 4.2.2 Production and Destruction Terms 78 4.3 The Final Model Equations 81 4.4 Discussion . 82 vi 5 MODEL CALIBRATION 88 5.1 Isotropic Decay 89 5.2 Return-to-Isotropy 93 5.3 Channel Flow 96 5.4 Asymptotic Homogeneous Shear 102 5.5 Equilibrium Region in Channel Flow 111 5.6 Near Wall Behavior 1 14 5.7 Summary and Discussion ' 121 6 MODEL PREDICTIONS FOR FULLY DEVELOPED CHANNEL FLOW 126 6.1 Parametric Study of the APS-Theory 127 6.2 Numerical Predictions for the Outer Region of Channel Flow 131 6.2.1 Model Formulation 131 6.2.2 Solution Strategy 135 6.2.3 Results 137 6.2.4 Summary of the Numerical Predictions 152 6.3 Mechanistic Influences of the New Closure Theory on the Reynolds Stress 155 7 CONCLUSIONS AND RECOMMENDATIONS 168 7.1 Conclusions 168 7.2 Recommendations 172 APPENDICES 177 A Direct Numerical Simulation Data 177 B Friction Factor Analysis 181 C Derivation of the Exact Turbulent Kinetic Energy Equation 184 D Derivation of the Exact Turbulent Dissipation Rate Equation 188 E First Order Approximation of the Equilibrium Region 193 F Near Wall Analysis 196 G Computer Program For Fully Developed Turbulent Channel Flow 200 LIST OF REFERENCES 210 vii 3.1 4.1 4.2 4.3 5.1 5.2 6.1 6.2 LIST OF TABLES Reynolds Stress Components in Terms of Isotropic and Anisotropic Prestress Contributions Modeling of the Turbulent Transport for the Turbulent Kinetic Energy Modeling of the Turbulent Transport for the Dissipation Rate Modeling of the Production and Destruction for the Dissipation Rate Equilibrium Homogeneous Shear/Summary of Equilibrium States Parameters for Anisotropic Relaxation/Retardation Model Individual Components of the Anisotropic Part of the Prestress in Terms of Contributions due to [3, A1 and X2 Reynolds Stress Components in Terms of individual Terms of the Anisotropic Part of the Prestress viii 58 83 84 85 109 122 158 162 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 1.10 2.1 2.2 2.3 2.4 2.5 5.1 5.2 5.3 5.4 LIST OF FIGURES Mean Velocity Profiles for Various Reynolds Numbers Mean Velocity Defect Profile Reynolds Shear Stress for Various Reynolds Numbers Spatial Distribution of the Eddy Viscosity Ratio Ratio of Bulk Average to Friction Velocity Energy Distribution Plane of the Turbulent Field in Channel Flow Realizability Diagram Turbulent Kinetic Energy in the Near Wall Region and in the Outer Region of Channel Flow Dissipation Rate from Numerical Simulations and from Experiments Shear Parameter kS/a Spatial Distribution of Re1 in the Near Wall Region Dependence of the Normal Stress Component Ryy and the Shear Component ~Ryz on cw° (based on DNS-data for 52395) a Normal Stress Component Ryy b Shear Stress Component -Ryz Energy Distribution Plane of the Turbulent Field in Channel Flow Spatial Distribution of Parameter 08 using DNS-Data for 5+=395 and both Estimates for ca Realizability Diagram for Isotropic Prestress Least Squares Analysis for cR and CB Normal Stress Difference Rxx-Ryy Energy Distribution Plane of the Turbulent Field in Channel Flow with the Anisotropic Closure for the Prestress Realizability Diagram for Anisotropic Prestress 11 12 15 16 18 21 23 25 41 42 42 43 45 47 98 100 101 102 5.5 5.6 5.7 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 Coefficient for Production Term in the 8-Equation Reynolds Stress in Asymptotic Homogeneous Shear Flow Production Terms P81+P82+P£3 Parametric Study within Energy Distribution Plane Parametric Study in the Invariant Plane Approach to Isotropy near Channel Center Turbulent Kinetic Energy for different Reynolds Numbers Dissipation Rate for various Reynolds Numbers Computed Mean Velocity Profiles for different Reynolds Numbers in Comparison to Experimental and Numerical Data Mean Velocity Defect Profile Shear Stress Distribution across the Outer Region of the Channel Spatial Distribution of the Turbulent/Mean Field Time Scale Ratio 6.10a Budget for the Turbulent Kinetic Energy 6.10b Budget for the Dissipation Rate 6.11 6.12 6.13 6.14 6.15 6.16 6.17 The Energy Partition among the Normal Reynolds Stress Components The Spatial Distribution of the Eddy Viscosity and its Coefficient Dependence of the Friction Factor on the Reynolds Number Decomposition of the Deviatoric Part of the Prestress into individual Components a Normal Component b Shear Component Influence of the Prestress on the Reynolds Stress Component Normal to the Wall Influence of the Prestress on the Reynolds Stress Component N ormal to the Wall ' Influence of the Prestress on the Shear Component -Ryz 108 110 118 128 130 132 138 140 141 143 144 146 148 149 150 152 155 160 160 163 164 165 E B "<7 u> llw 01,02 Cnn CM 0x2 CB cal, c£2 Ck»Ce 0?, CD Cs Del fd '5 ”'3 e” N 7r "H u: n:- NOTATION Second invariant of B Third invariant of B Prestress operator Anisotropy tensor Normalized anisotropy tensor Parameter in Rodi’s and Shih’s model Parameter in N isizima’s model Eddy viscosity coefficient Coefficient for relaxation term Coefficient for retardation term Strain rate coefficient Parameters in transport equation for e (i.e. standard model) Transport coefficients for turbulence model Parameters in transport equation for 8 (production, destruction) Parameter in Hanjalic’s model Turbulent Deborah number Friction factor Wall function employed by Nisizima and Hanjalic Wall function for eddy viscosity Fluctuating acceleration Wall function Prestress tensor Normalized prestress tensor Unit dyadic Turbulent kinetic energy Length scale in Speziale’s model Length scales for Green’s function decay, eddy size xi :1 Greek Letters Q m -< A419 7‘2 Convective diffusive operator Exponent for isotropic decay Pressure Production of kinetic energy, production tensor Production terms in e—equation Energy distribution parameter Normalized Reynolds stress Reynolds number Reynolds number based on Taylor microscale Reynolds number based on mesh size Turbulent Reynolds number Characteristic strain rate, strain rate dyadic (=O.5(Vy_+(Vu)T)) time Instantaneous velocity vector, fluctuating velocity vector Mean velocity vector Reynolds stress Friction velocity, bulk average velocity Vorticity tensor (=0.5(Vu—(Vu)T)) Position vector Kinetic energy of the prestress g Parameter in Rodi’s and Harlow’s model Dissipation rate Generic transport quantity Relaxation time, retardation time xii Gk, Ge Tm TP, TD 13:, TR 17w Karman number Dynamic viscosity Kinematic viscosity, turbulent kinematic viscosity (eddy viscosity) Parameter in Harlow’s model Density Parameter in Harlow’s model Turbulent Prandtl numbers in standard k-e model Parameter in Nisizima’s model Time scales for production and destruction Turbulent time, relaxation time for the Green’s function Wall shear stress Dimensionless variable indicating distances Symbols and Super/Subscripts xx, yy, 22, yz Tensor Components W <> Wall Ensemble average Normalized values (i.e with wall parameters) Fluctuating quantity Dummy variable, dimensional variable Reference value Friction velocity Parameter in Shih’s model Transpose of prestress operator A Gradient operator xiii CHAPTER 1 INTRODUCTION 1.1 Motivation and Background Motivation Turbulent motion contributes significantly to the transport of momentum, heat and mass. Swirling flows, which are typically high Reynolds number swirling flows, are used to separate materials with different densities. Gas cyclones have been used in cleaning flue gas (i.e. dust removal) and deoiling hydrocyclones have been used in removing oil from water. In order to develop a mathematical model capable of predicting the separation efficiencies of hydrocyclones based on the description of particle/droplet trajectories, it is necessary to devise a mathematical model which can quantitatively describe the motion of the continuous phase. Fluidized beds used for combustion or as catalytic reactors constitutes another example where turbulent motion contributes to transport processes. In combustion, turbulence enhances the mixing process and, thereby, allows a more efficient burnout of the fuels used. In catalytic reactors the turbulence contributes to better mass transfer rates of the constituents present in the bed. It is the interaction of the solid particles with the surrounding fluid which determines the level of turbulence and the efficiency of the overall process. Therefore, knowledge of the motion of the continuous phase is therefore essential for predicting efficiencies and to improve the design of the foregoing devices. 2 The presence of turbulence in almost all flows of practical interest has inspired many researchers to devise models which are capable of describing its motion. The direct computation of turbulent flow fields requires a three-dimensional spatial domain which is sufficiently large to capture the large scale motion of the flow but uses a computational grid which is sufficiently small to resolve the smallest scales of the turbulence. Since the small scales decrease in size with increasing Reynolds number (Tennekes and Lumley, 1972), the direct computation of even simple shear flows such as a two-dimensional channel flow is still subjected to low Reynolds number flows. Therefore, turbulence models need to be used for the computation of high Reynolds number flows. For most practical applications it is sufficient to know the influence of the Reynolds stress on the mean flow field. This Reynolds stress which arises from the nonlinearities in the Navier— Stokes equation appears explicitly in the Reynolds equation (Eq. (1.1)) and needs to be modeled: a+<9>v<fl>=71§V

+%V-(<;>-p<2'u'>). (1.1) The development of many turbulence models has been done using simple flow fields for which comprehensive experimental data for validation purposes are available. Statistically stationary, fully developed turbulent channel flow constitutes such a simple flow field. Many approaches to close the Reynolds equation by prescribing an algebraic relation for the Reynolds stress have made use of an eddy viscosity concept which models the Reynolds stress analogously to the molecular stress as being proportional to the mean strain rate dyadic (Eq. (1.2)). This approach is known as the Boussinesq approximation (1877, see Speziale, 1991). —=2v ——k1, (1.2) where v; is the turbulent or 'eddy' viscosity, <§_> is the mean strain rate dyadic and k is the turbulent kinetic energy. The eddy viscosity is commonly modeled in terms of local statistical turbulent parameters as V=C t u 2 k_, (1.3) e in which 8 is the dissipation rate. If the coefficient c,Ll is kept constant, it can be shown that in simple shear flows, such as channel flow, Eq. (1.2) in combination with Eq. (1.3) constitutes an unrealizable turbulence model for large values of kS/e where S is a characteristic strain rate defined as S=,/2<§>:<§>. (1.4) The realizability (i.e. the Reynolds stress is positive semi-definite) of this type of model relies on the solution of the associated transport equations for k and e to keep kS/e bounded. Thus, for turbulence models based on Boussinesq's approximation to remain realizable some empirical modeling approaches have been used. In wall-bounded flows the coefficient c” is commonly modeled as being a function of the distance from the wall and/or other parameters such as the turbulent Reynolds number Re, and has been used to provide that kS/e remains bounded, thus ensuring realizability for those type of flows. Eq. (1.2) points out another deficiency of the Boussinesq approximation. For vanishing diagonal components of the strain rate dyadic (which occurs for channel and pipe flow) the normal components of the Reynolds stress are proportional to k and isotropic. However, experimental investigations in those flows show that the equipartition of the turbulent kinetic energy does not exist (Laufer, 1951; Kreplin and Eckelmann, . 1979). The turbulent fluctuating velocity component normal to the wall is the smallest in magnitude and would therefore be overpredicted by Reynolds stress models employing Boussinesq's approximation. This leads also to the fact that turbulent mass transfer % —-—.4 4 normal to a wall due to velocity fluctuations of the normal components will be overpredicted as well. Background The mean flow field of statistically stationary, fully developed channel flow is characterized by the fact that the velocity field consists of only the streamwise component which depends solely on the normal coordinate as given by < 2 >=< u. > on. (1.5) Because turbulent statistical quantities do not vary in the x- and z-direction, the Reynolds equation as given by Eq. (1.1) can be simplified for this flow field. The corresponding equations for the streamwise and normal-to-the-wall components are given by a

a , , O=— 8y —$(p) (1.6) and em» d d d 0=— +—— Z -—— < ’ ’>. 1.7 82 dyfll dy ) dy(p 11,11. ) ( ) The gradient (with respect to z) of the mean pressure - since it does not depend on y - can be calculated from an overall force balance with the result that: 8

tw — =——. 1.8 dz 5 ( ) tw is the wall shear stress and 5 denotes the channel half width. Eq. (1.7) can be thus integrated with respect to y to yield the following form of the momentum equation: d dy I I _ y —p—tw(1—§). (1.9) #_—4_4 5 An expression for the variation of the pressure across the channel can be obtained by integration of Eq. (1.6) to

(y,z)=

(0,z)—p(y). (1.10) The Reynolds stress in fully developed channel flow consists of the normal components and the shear component arising from the correlation between streamwise and normal-to- the-wall velocity fluctuations and can be expressed as ’ ’ , ’ I I I I < u u >=< uxux > (y)§x§x+ < 11y“, > (Y)§y§y+ < uzuz > (”922. (1.11) + < u;u; > (y)§y§z+ < u;u; > (y)§:_zey All components depend solely on the normal coordinate y. This coordinate can be normalized with the viscous lengthscale v / u. to + yua- = . 1.12 y V ( ) where u. is the friction velocity defined as u‘ z "_w , (1.13) P The mean velocity gradient in the viscous near wall region can be expressed as d u? dy v. (1.14) Upon integration one obtain the well-known linear behavior of the mean velocity near the wall. u+ = y+. (1.15) The velocity is normalized with the friction velocity to u,:. (1.16) U. 6 The viscous sublayer for which Eq. (1.15) holds is observed experimentally for y+ < 5. The turbulent flow domain which covers most of the channel is made dimensionless by using the channel half width: gzy- (1.17) If the Reynolds number is high enough, there exists a region for which the influence of viscosity vanishes such that y+ » 1 yet E << 1 such that the flow field does not experience the influence of being in a confined domain (i.e. the channel half width does not yet enter the problem). The only length scale available is therefore the distance from the wall. The mean velocity gradient is thus expressed as d_u; (118) dy y’ which - upon integration - yields the well-known logarithmic velocity profile in this region (Tennekes and Lumley, 1972): u+=Aln(y+)+B. (1.19) Figure 1.1 shows some typical velocity profiles for different Reynolds numbers. It can be seen that for y+ S 30 all profiles seem to collapse onto a single curve indicating that the scaling with near wall values yields the prOper similarity whereas for y+ 2 30 the velocity profile can be well represented by Eq. (1.19) (center region excluded). The region between the viscous sublayer and the inertial sublayer is called the transition or ‘buffer’ region. In the inertial sublayer, the experimental data by Laufer (1951) for the mean velocity indicate that B = 5.5, a value commonly accepted for channel flow over smooth walls. The data of Reichardt (1940) for Re=7,700 yield a slightly smaller value of B: 5.37. The direct numerical simulation (DNS) data by Kim et a1. (1987) yield B = 4.79 whereas unpublished data by Kim (1989) yield B = 5.23. The coefficient A has also been determined from experimental data as the reciprocal — g ‘ Era" muonfisz mEoSAom 25:8» 8% 858m 3828/ :82 n: 83$ 00— o— —_~_ —_p— p 4____.__ 8%. as: coming ......... $2 a a so: ownmnom 1.... Ammo. .HEEQBMV 25$qu 23— .8351: 896qu 2mg cougar—v ocwdmuom :5 Seed 838m ___PL__. 0 U om mm Om 8 of the von Karman number it. A value for K of 0.41 has been commonly accepted even though a conceptual investigation by Tennekes and Lumley (1972) has shown that this might not entirely be justified. The experimental data by Laufer, for example, yield a value of K = 0.35 which can be seen in the larger slope of the curve. Reichardt’s data yield K = 0.41, Kim et al.’s data (1987) yield K = 0.38 and the DNS-data at a higher Reynolds number yield K = 0.41. The central region of the channel where viscosity influences are negligible is characterized by the velocity defect law which expresses the difference between the velocity and its maximum to u+ —u; =A1n(§)+B*. (1.20) This expression can be obtained by integrating the similarity expression for the velocity from the centerline (Tennekes and Lumley, 1972). The difference between Eq. (1.19) and (1.20) yields an expression for the velocity maximum to u;=A1n(5+)+B—B*. (1.21) The parameter 8+ which is defined as 6* = 81“, (1.22) V represents the ratio of the largest (i.e. 5) to the smallest scales (i.e. v / u. ). The velocity defect law as given by Eq. (1.20) is presented in Figure 1.2. This figure which shows velocity profiles for three different Reynolds numbers (see Hussain and Reynolds, 1975) indicates that similarity at the center of the channel exists if bulk flow properties are used for sealing. The fact that the integration was carried out from the center is expressed through the difference of the velocities (LHS of Eq. (1.20)) and thus removes the dependence on near wall scales (upon which the individual terms do depend). As shown in the figure, it is an increase in the Reynolds number which extends n—— —44 @505 60on 3829/ :32 ”NA oSmE L m; 3 5o Sod éy_ . _ 4 a 4 d 1— q q a a — A _ _ a _ _ _ _ _ q a o - 9% - @ODD T vaow .I N T w W 1 1 mp; Sammie . 1 v - on 0.. Sammie a - 1 2.08 8w.m_uo~_ o 1 c 1 o o r o m... 38: £83m 1 o own. 9 use 583: 1 w - on as B 3585me - I F WTA v.0 o 00 DD 1 Ofi - . 1 . . . o - B . 1 an H 8 p 1 2 D 1 D 1 %% D 1. —. D 1 1 D D 1 OD . 09 l V.— om . DD 1 on 1. N o D D 1 l MWm. H o 09> . 1 CM 0 D 1” m o 0099 I D H O 0% 1 O D 1 0 DD 7 Op 1. v o 0009 1 WM - w . 039 M U H O 0 % AU 0 1 . 1 .6 . m as 1 om . :1 += . _ _ r _ _ _ . _ _ _ _ p p _ u _ _ p u _ b . L _ NN 5140: 10 the domain of the similarity to smaller values of E. The integrated form of the momentum equation (i.e. Eq. (1.9)) with the expression for the universal logarithmic law for the velocity profile can be used to give an expression for the normalized Reynolds stress; namely, 1 1 y z _____._____ 1__ __ . u2 ( g) K5+§ t (1.23) For very high Reynolds numbers where the viscous length scale becomes vanishingly small compared to the channel half width (i.e. 8' -—> co) the Reynolds stress follows the straight line as indicated by the bracketed term on the RHS of Eq. (1.23) and indicated in Figure 1.3. For E —-> 0 and y+ (58f) finite, Eq. (1.23) expresses the fact that there is a region very close to the wall in which the Reynolds stress is approximately constant. For finite values of the Reynolds number, the influence of the second term becomes more important such that the Reynolds stress shows a maximum at finite values of g. It can be seen that for smaller Reynolds numbers (with y+ fixed) the deviation from the linear profile occurs at increasing values of 5,. The dashed line in Figure 1.3 indicates the separation of the buffer region and the inertial sublayer and is seen to proceed towards smaller values for y/S for 5‘1» oo. For channel flow the eddy viscosity as introduced in Eq. (1.2) can formally be expressed as — < u;u; > V1 =————_, (1.25) d dy which is valid for the entire domain. For the inertial sublayer (i.e. y+ > 30) the integrated form of the momentum equation (i.e. Eq.(1.9)) in combination with the prevailing logarithmic velocity distribution as given by Eq. (1.19) can be rearranged to yield _. .__.._._____ _..._..‘ ._.__.__ .._._.__ -— —"' -...___,_-— 11 $3852 mEoiom msotm> new 38% Beam mEOSAom ”mg 253m s . o._ w w c We Yo N o o H4 — q _ _ fi . u — _ q _ fl . — _ _ _ 1 id H .m o 1 o. 1 o 1 o O 1. ® 0 l T 01 I O . I. 1 nwoo o o 1 I 0a I 1 @a o o 1 I I rll v. I! I 1 8mg ELSE—:3 use 626 m_m._n+w ® OM11...» b8 :84 11 l 81% Sm Amm. 3 .cm Eases—32 Ill 23— doused ommuk @ awe EU: mmmnb o Ammonia 6 EU: ow~n+_w o — . od Nd v.0 N +A. 5. 0; 5V1 12 -\:’—t=K8+(l—§)§—1. _ (126) Figure 1.4 shows the behavior of the normalized eddy viscosity for channel flow at different Reynolds numbers. The different magnitudes for this ratio arises through 5+ in Eq. (1.26) which depends on the Reynolds number through the presence of the friction velocity u... Its definition can also be interpreted as a Reynolds number based on u. rather than the bulk average velocity ub. The form of the Reynolds number dependence of 5' can be estimated from experimental data as follows. With the definition of the Fanning friction factor f as f=Ap/(15—%u§), (1.27) and the expression for the pressure drop along the channel as given by Eq. (1.8) the friction factor can be expressed in terms of the velocity ratio ub / u. as _= , (1.28) which in turn can be evaluated by integrating the gradient of the mean velocity over the entire domain to yield 2 l g . v . £=8+J(J(1-§)/(1+—‘—1d§)d&. (129) 0 0 V For practical purposes it is, however, easier to estimate u. from measurements of the wall shear stress (i.e. from the velocity gradient at the wall) since most experimental data are commonly available in terms of rather than the eddy viscosity ratio. With given channel dimensions, the corresponding values for 5+ can readily be calculated. Computationally, the friction factor can be related to the mean field and the turbulent field using the following equation .'—~“ 1.- A-v-r dlb.» ’ 41h». 13 23m 3608; Exam 2: E aossfibma 38mm x: oSmE I in G we we v.9 oo . _ 1 u _ . . _ _ . . _ . 0 00000000000 GOOD 0 o o o o o o o..o..o..mb:o:o:o.b n no 0000000000 cococoocooooooooooooooooooooooooww m. ............ m ...... - m. E E . 11 CW ................ u I I 1 I .. . I I 1 On: .210 u v: 82? .. r @8338 05 03:38:06 ...... Cme £836 oomguk - - +w 8m 32? BE. .8333. 2 we £331: ommuk a . 3:2: 2: é an: em . 1 E8058 3:: 3326 BE. 6mg EU: mom +m . - >\.> _ _ _ A33. .3 we Eflvc pwfiutqw o _ ow“ 14 3_ +1 (111+ 2 + 111 -5 {fl—dy.) +1: 11%, (1.301 which is obtained from the integrated form of the equation for the kinetic energy of the mean field (see Appendix B). The integral over the production term which appears explicitly can be represented in terms of the integral over the dissipation rate in order to obtain Eq. (1.30). Figure 1.5 presents the ratio of ub / u. for various values of 5+ from experimental data. The Reynolds stress tensor as introduced by Eq. (1.1 l) is commonly normalized with the turbulent kinetic energy to - R 1.31 = 2k ( ) and decomposed into an isotropic and anisotropic part to 1 R=§I+R. (1.32) From Eq. (1.31) it is clear that the trace of R equals unity which in turn - with the decomposition according to Eq. (1.32) - renders R traceless. The fact that the trace of R equals unity together with the property that the individual diagonal components of R are positive - since they constitute energy contributions to R - can be represented graphically in terms of a triangle representing an energy distribution plane (Figure 1.6). The corners of this triangle represent one-component turbulent energy states with all the energy transferred to the corresponding component and lines connecting the comers are therefore two-component turbulent energy states. The center of the triangle represents energy equipartition. The perpendicular lines on the baselines, indicated by the clashes, represent transition curves between one-component and axisymmetric two-component energy states. It should be noted, however, that the shear stress is not represented. Thus, the % ———-4 15 k 5829» 8:05 8 owfiog 6:3 E 23% ”m; oSmE ll Efiofim .mmméflé. 1 mikvczmgps ”HE 0230 .II 1 3mm; SQ 1 $3: ._a 8 SE. 1 AmmoC 883* 1 Aowm C ESE—:3 was 63 23:51:54 1 OO®<1 0m ow On .39: 16 - Kim et al. (1987; Re=3,250) ° Kim (1989; Re=7,890) o Laufer (1951 ; Re=30,800) ° Kreplin/Eckelmann (1979; Re=7,700) RZZ C I I I 5:5. I . ? I ° E ;. 1 32° : o: 0 A -\ .. 1 $0 3: °| "o. Q g: o ‘0 3 ° .1 <6 S 10' a c .10 Q 0 o A 90 Q \ o / ¢ 3 \ ‘31. / ’ 0° \ l / (P o’ \ \ / / $0 (at 3 \\Bl / <33. / / I \ \ / 1 \ / I \ / / I \ \ / l \ / l \ / \ / 1 Energy. \ / / l Equrpartrtron \ \ / | \ / l \ R1111 Two-Component Energy R Figure 1.6: Energy Distribution Plane of the Turbulent Field in Channel Flow 17 center is generally a turbulent state in which the kinetic energy is equipartitioned and only in special cases truly isotrOpic. Measurements of the diagonal components of the Reynolds stress in channel flow are shown in Figure 1.6 to illustrate the path along which the energy is distributed among its components in the spatial domain. The z-direction denotes the downstream coordinatewhereas the y-direction is normal to the solid wall. The point on the baseline between Ryy and R22 is a two-component state which prevails at the wall. It can be seen that neither the DNS-data (Kim et al., 1987; Kim, 1989) nor the experimental data by Laufer (1951) or Kreplin and Eckelmann (1979) show energy equipartition at the channel center. The energy distribution path lies within the sub-triangle formed by the center (i.e. equipartition of energy), the l-component turbulent energy state (i.e. Ru) and the two- component axisymmetric energy state as by the line A—C. The initial tendency - if viewed from the center of the channel - to adhere to the transition line connecting the equipartition and l—component state indicates an energy state dominated by the downstream energy component. A supplementary view of the turbulence states within channel flow is given by analyzing the second and third invariants of the anisotropy tensor R which are defined as IIB = tr(R-R) (1.33) and IIIB = tr( "or: new ). (1.34) Note, that the first invariant vanishes identically by definition (i.e. IB = tr(R) = 0). Since R is a symmetric tensor, the second invariant is always positive. The third invariant can assume positive and negative values. For fully developed channel flow, these invariants are given explicitly by 18 IIB = Bi, + B; + B; + 213; (1.35) and 3 3 3 2 IIIB =BXX+BW+BZZ+3BYZ(BW+BZZ). (1.36) This figure translates the existence of all turbulent states into points within the triangle bounded by the lines given (see Figure 1.7). Lumley and Newman (1977) have established the boundaries of this triangle from a conceptual analysis of all possible states for the anisotropy tensor R from which it can be concluded that turbulence must occur within this triangle. The analysis to obtain these boundaries has been with respect to a transformation of the Reynolds stress tensor to its principal axis. The comers of this triangle correspond to single—component (1C), two-component (2C) and three-component (3C) turbulent states as indicated in the graph. The 2C turbulent state is characterized by the fact that both components contain equal amounts of energy and the 3C turbulent state is truly isotropic since both invariants vanish. The line connecting points 2C and 1C represent two-component turbulence for which the energy is not equally distributed among its two components. The lines connecting the 3C-state with the 1C-state and the 2C-state represent axisymmetric turbulence. The line to the left constitutes an energy state in which two components are equal and larger than the third component whereas the right line represents a state with the two equal energy components being smaller in magnitude than the third one. The functional forms of the boundaries and the comers of the domain are listed in the graph. The symbols in Figure 1.7 indicate the trajectories of turbulence as given by DNS- data (Kim et al., 1987; Kim, 1989) and by the experimental data (Laufer, 1951). As noted earlier in the energy distribution plane it can be seen here that the center of the channel is not isotropic. The path along which turbulence proceeds towards the wall is close to the right boundary indicating a strong tendency towards two-component axisymmetric II 19 0'7 f - Kim etal.(1987;Re=3,250) _. 0 Kim (1989; Re=7,890) - . Laufer(l951; Re=30,800) 1C Needle 0.6 — 0.5 — . - 2C AnisOtroprc 0.4 — q °. y = : Boussinesq . 0.3 — for 011:0'09 . o- _ 0.0. \ Prolate 0'2 :_ 0.; Axisymmetry : 2C Isotropic 1 G2' I t 80 '° 60 _ y = :8 80 0.1 — 250 100 — , 120 — Oblate / . . .;: 3C Isotroprc r Axrsymmetry 1:/ 0.0 1 LL] 1 l l l l : l l I I l l l L 1 LI 1 I l i I I 1 I 1 1 I l l i -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 IIIB Oblate Axisymmetry: 113:6(‘1113/612/3 1 Component: (2/9, 2/3) Prolate Axisymmetry: IIB=6(IIIB/6)2/3 2 Component Isotropic: (-1/36, 1/6) 2 Component Anisotropic: IIB=2l9+ZIIIB 3 Component Isotropic: (0,0) Figure 1.7: Realizability Diagram 20 turbulence. Whereas the data by Kim et a1. (1987) proceed along a monotone path within the region of y+=60 and y+=120, the DNS data for the higher Reynolds number indicate a region over which the turbulence seems to adhere in ‘one state’. The occurrence of this ‘hook’ in the invariant plane is of yet unexplained. The tendency to assume a one- component turbulent state is clearly visible in Figure 1.7. Very close to the solid boundary (proximity is indicated) this tendency is changed towards a two-component turbulent state as indicated by a deflection point. Within the invariant plane it seems that for increasing Reynolds numbers this deflection point occurs at smaller values for I13 and THE. However, within the spatial domain no characteristics can be found which may relate to the Reynolds number. The turbulent state at the wall is a two-component state. Unfortunately, no measurements at large Reynolds numbers have been made deep within the viscous sublayer in order to supplement the tendency of the wall location towards a 1C- or 2C- turbulence state with increasing Reynolds number. It is thus unclear whether there might exist a critical point along the transition line between 1C and 2C which is approached for high Reynolds numbers. From Eqs. (1.35) and (1.36) it can readily be seen that for simple shear flows a Reynolds stress model according to Boussinesq’s approximation with an energy equipartition among its components does not have a third invariant of the anisotropy tensor R and follows therefore a straight vertical line as indicated in Figure 1.7. As noted earlier this approximation yields an unrealizable turbulence model if C11 is constant and kS/e assumes values larger than 2/(3cu). This model does not reflect the correct turbulent state while approaching the upper boundary representing the two-dimensional state. The use of the Boussinesq approximation for the Reynolds stress therefore requires a functional form for the parameter C11 as introduced in the eddy viscosity representation in Eq. (1.3) in order to provide a realizable turbulence model (see Patel et al., 1985 for a review). 21 The distribution of the turbulent kinetic energy is shown in Figure 1.8 in terms of the wall coordinate y”. Several experimental and numerical results for the spatial distribution of k” are shown. Kreplin and Eckelmann (1979) used an open channel with oil as working fluid and were thus able to measure turbulence quantities very close to the wall (i.e. y+ = 1.5). However, since all three components were not measured this deep within the viscous sublayer, values for k+ are only available for y+23. These measurements were probably the only detailed measurements available for the kinetic energy this close to a solid wall in a channel. The data of Laufer do not extend into the near wall region. His measurements reflect the Reynolds number influence on the magnitude of k+ in the outer region which can be seen to decrease with increasing Reynolds number. The numerical simulations by Kim et a1. (1987) and Kim (1989), however, indicate an increase with increasing Reynolds number. The extent to which this apparent discrepancy is based in fundamental aspects has not been discussed yet. It can be seen that also in the near wall region (y+ < 30)some differences occur in the magnitudes of k+. It seems unclear whether the scaling of k”r with near wall values yields the proper similarity as was pointed out by Wei and Willmarth (1989). The budget for the kinetic energy and dissipation for the higher Reynolds number DNS-calculation show a slight imbalance in the region y+ < 10 (see Chapter 4) which might explain in part the differences between the two numerical simulations. Differences in the peak values for the normalized streamwise RMS-value can be found also in Laufer’s data which seem to decrease from a value of 2.76 for Re=12,300 to 2.54 for Re=30,800 and 2.04 for Re=61,600. These values occur at locations between y+=16 and 54 and need to be considered in combination with the fact that the hot wire used for the measurements extended over a spatial domain of Ay+=3 to 13 which has certainly some averaging effect on the values measured. Comte-Bellot (see Kreplin and Eckelmann, 1979) claimed that the scaling of fluctuating quantities with near wall values is adequate 22 :Bwom :85 :32 2: 5 $55.. 285M 222:3. ”M: ouswfl coo— BoE 353:0 we somwom 830 05 E can oo— o_ _ ‘4—q__— _ __4—__ bpbb. _ q r — 10.221. . fl . Soc. 5qu ”ammo C .5qu Sowdmuom Mama 3 £33 33.2qu €33 $qu Aooafiuom ”a3 _ v ccmEEV—om ecu EEBM Sowfiuom Mama—v EU— 8m~.mnom meC .E E SE —_~_-—-— _ —.__LP—— _ 23 for values of y+ up to 100 even though her measurements indicate a similar decrease from 2.85 to 2.65 and 2.5 for Re=57,000, 120,000 and 230,000, respectively. However, the probes lengths in these experiments were between Ay+=l3 and 36 which may cause similar effects as in Laufer’s experiments. Experimental data for the dissipation rate are usually unavailable due to the inherent difficulty in obtaining them. The DNS-database provides therefore the best source for this turbulent quantity. Figure 1.9 shows the computed profiles for 8+ for the two Reynolds numbers indicated. The only measurements of this quantity have been made by Laufer (1954) in turbulent pipe flow at a Reynolds number of Re=50,000. However, in obtaining these data several assumptions about isotropy of the small scales have been made to facilitate the evaluation of 8+. These data are included in Figure 1.9 in order to get a general view of the Reynolds number dependence of this quantity. With the assumption that the outer region of channel flow (y+>30 but §<<1) is in equilibrium (i.e. production of turbulent kinetic energy is balanced by its dissipation) the functional behavior of 8+ in the inertial sublayer where the stress is approximately constant and the logarithmic velocity profile prevails can be deduced from Eq. (1.9) to +_ 1 e _1<8+§ (1.37) This behavior is indicated by the dashed lines in the graph and seems to hold fairly well for the DNS-profiles as well as the experimental data by Laufer (1954). The shear parameter kS/e which compares the turbulent time scale k/e to the mean field time scale US is often used in the development of algebraic Reynolds stress models inasmuch as it provides some characteristic of how the turbulent field compares to the mean field. The spatial distribution of this parameter is given in Figure 1.10 for both DNS-data sets as well as the experimental pipe flow data by Laufer (1954). It can be seen that in a wide region of the channel (pipe) this parameter remains fairly constant. The 24 358E030 Bot 28 2833586 3258:: Eat. Bed 8336me no; 2sz 56 cod Gmd .vm 9 dog +w l l. . . . . . wave 03m 2: 38:2 w . 88.8132 ”>3: 85 ”43: £23 . . . . “ 83?an gm: SQ . smumnom ”5a: .a 3 SE . _ P I — h — _ 1 mod . Dd L_‘._... _ . . "9;, .-—_._... me $88.85 32m ”9 A oSwE o._ wd 9o v.0 \ Nd od IOI h M _ _ q d _ _ _ a _ _ _ _ _ — _ q _ o I n n I u no. n .03 u no. ” Ion H cos n Ion u I.- ” no. u... - 'uuolnlnlnnnlllllllul ¢ I a” In 0 o 0 III 4 I m 5 I o I (on: 2 u . . 1 u . . .. o. I . I 2 H soodmnom 36: 2E «a: £23 . . . 1 sowsnom 53: 8E . ... swamnom 53 ...w 8 EE . - 8%. p b p _ _ _ P r _ _ p _ _ _ _ _ _ _ _ ON 26 DNS-data approach values of around kS/e = 4:03 whereas the experimental data gradually approach a value of 3.3 after which the values increase. The constant behavior suggests that turbulence and mean field are in some sort of ‘temporal balance’ without any significant ‘influences’ by other physical phenomena. The magnitudes of the shear parameter indicate that the ‘average’ life span (i.e. turnover time) of a turbulent eddy is by a factor of 3-4 larger than the time in which momentum transport by the mean field takes place. Indirectly, this leads to the conclusion that the presence of turbulence enhances the transport of momentum through the comparable ‘large’ life cycle of its eddies. Towards the centerline the momentum transport caused by the mean field vanishes, indicated by a zero shear parameter. The more gradual approach apparent in Laufer’s data may be attributed to the fact that those experiments were conducted in pipe flow with significant curvature rather than channel flow. At the wall the shear parameter vanishes as well. The reasons for this are, however, different. In the viscous sublayer there exists a region in which viscous effects become dominant. The kinetic energy of the eddies decreases and ultimately vanishes. The effect of turbulent eddies decreases and does not contribute to momentum transport which - in this region - is entirely dominated by the mean field (i.e. large shear —> small time scale). It can also be seen that the shear parameter assumes similar magnitudes in both flow fields. which may be in accordance with the diminishing influence of the curvature (i.e. the geometry). A region of increased magnitudes in the shear parameter can be found close to the wall. This region occurs within the buffer region where the majority of the turbulence is produced. This effect is in accordance with the interpretation of highly agitated turbulent momentum transport. The influences of molecular transport are of no major importance within this domain. 27 1.2 Objectives and Methodology Objectives This research pursues the development of an algebraic Reynolds stress model using fully developed turbulent channel flow as the basis for this development. The amount of available experimental and numerical data support the development inasmuch as in depth comparison between the turbulence model developed herein and the data for turbulent statistics can be made. Several key issues are subject to a detailed investigation. Anisotropies among the normal components of the Reynolds stress as observed experimentally and presented graphically in form of the energy distribution plane (see Figure 1.6) form a major aspect of this research. Realizability as a fundamental aspect to which every turbulence model ought to be subjected restricts the class of applicable turbulence models to the invariant plane by Lumley and Newman (see Figure 1.7). This research attempts to provide a turbulence model which is a priori realizable and does not depend on the solution of the associated transport equations for the turbulent kinetic energy and the dissipation rate. Two further aspects are investigated in this research. The frame dependence of the Reynolds stress and the correct asymptotic behavior of the individual components of the Reynolds stress as a solid boundary is approached are addressed and incorporated into the turbulence model developed herein. The ubiquitously employed Boussinesq approximation - commonly used in combination with transport equations for the kinetic energy and dissipation rate — shows some deficiencies with respect to the key points mentioned above and as pointed out earlier. The associated transport equations are modified in this research to comply with realizability (Chapter 4, 5). 28 Methodology The development of the turbulence model is based on the equation for the fluctuating velocity 3'. A formal representation of the Reynolds stress is developed from this equation through a Green’s function technique. A smoothing approximation is introduced to reduce the non-local structure of the model to a local structure. This representation gives an algebraic approximation for the Reynolds stress in terms of the mean velocity gradient due to the coupling of the mean velocity field with the fluctuating field and a self-correlation of an effective fluctuating force induced by pressure fluctuations and fluctuations in the instantaneous Reynolds stress as well as a relaxation time relating to the temporal structure of the turbulence. The coupling of the Reynolds stress to the velocity gradient V rather than the strain rate dyadic <§> provides a means to insure the frame-dependence of the Reynolds stress (Chapter 2). This self-correlation constitutes a prestress for the actual Reynolds stress and is decomposed into an isotropic and an anisotropic part. The representation of the Reynolds stress through the isotropic prestress (preclosure) guarantees a priori realizability inasmuch as the self—correlation is a positive semi-definite tensor with non-negative eigenvalues. Some degree of anisotropy among the normal components of the Reynolds stress is introduced at this level. A zero secondary normal stress difference does, however, not support the level of anisotropy as shown in Figure 1.6. Closure of the anisotropic prestress is achieved through phenomenological modeling using a frame-independent relaxation/retardation model (Chapter 3). Adequate modeling of the associated parameters assures realizability and an extended degree of anisotropy. The correct asymptotic behavior of the Reynolds stress as the solid wall of the channel is approached is achieved through modeling of the relaxation time introduced through the preclosure. The universality in this preclosure accounts for the correct asymptotic behavior for all components of the Reynolds stress. Functional extensions of 29 the parameter cu in combination with the Boussinesq approximation as pointed out in the introduction usually provide only an adequate representation for the shear stress component as the wall is approached. The correct asymptotic behavior for other components, for example, the component normal to the wall, is usually disregarded. 1.3 Summary Outline of this Research This manuscript is divided into seven chapters plus an appendix. Chapter 1 reviews the background for this research and presents relevant data for channel flow with which the numerical predictions will be compared. The preclosure theory stemming from the Green’s function technique is developed in Chapter 2. Some of its underlying properties are discussed here. Chapter 3 presents the closure hypothesis for the anisotropic part of the prestress in conjunction with the underlying motivation for its use. The associated transport equations for the turbulent kinetic energy k and the dissipation rate 8 are presented in Chapter 4. The new formulations which are introduced into these equations are derived. Chapter 5 deals with the extensive calibration of the turbulence model proposed using different flow fields to determine the various model parameters needed. Special emphasis is given to the DNS-data for channel flow since they constitute a major source for the development. Chapter 6 presents the numerical solution of channel flow for the outer region and discusses the implementation of the new model developed in this research. Analogies to existing non—linear algebraic Reynolds stress models are given. The near wall region is discussed with regard to the DNS-generated database. Conclusions and recommendations for further research based upon the findings in this research, as well as suggestions for further research as an outgrowth of this research, are presented in Chapter 7. The appendices list data used for the development and presents a listing of the computer program used in Chapter 6. CHAPTER 2 REYNOLDS STRESS PRECLOSURE IN INERTIAL FRAMES 2.1 Preclosure Formulation for the Reynolds Stress ' The equation of motion for the fluctuating velocity can be obtained by subtracting the Reynolds equation Eq. (1.1) from the equation of motion for the instantaneous velocity to iu'=-u'-V -V-(y_'u'-+£1), (2.1) W p ‘- ’ / \ coupling of fluctuating and Reynolds stress fluctuating mean field fluctuations pressure where a 2 £=5f+V—VV (2.2) denotes the linear convective—diffusive operator associated with the mean velocity field. The first term on the right hand side (RHS) of Eq. (2.1) couples the fluctuating velocity with the mean field. This term constitutes the net transport of mean momentum through the fluctuating velocity field. The second term can be interpreted as an acceleration due to the divergence of the instantaneous Reynolds stresses and the gradient of the fluctuating pressure. This term will be denoted as f for brevity. Thus, the term pf represents an effective force per unit of volume of fluid and acts together with the first 30 31 term on the RHS as a source term for velocity fluctuations. A formal solution of Eq. (2.1) can be written as y x,t =—J jG(x, a: am < g > +£'}dvclt, (2.3) EV where G(~IA) denotes the Green’s function associated with the linear operator if. Through a formal dyadic multiplication of Eq. (2.3) with u’ and subsequent ensemble averaging a representation for the Reynolds stress of the form: —I jG(x, tlx, t)F (x, tlx, t)dth (24) tv can be derived. The dyadic function :1 is given by <2 ">i (2.5) |t\:/> Mali) <">g§ €7< The Green’s function spreads over a domain with a length scale comparable to (v1)°'5 where I is a characteristic time scale for the relaxation of the Green’s function. For large Reynolds numbers where the influence of the viscosity vanishes, the Green’s function remains sharply peaked relative to a frame moving with the mean velocity for a period of time during which the structure of the turbulence decays. Thus, the Green’s function only ‘allows’ the autocorrelation of the terms in brackets to be important (i.e. x = jg ). With this ‘spatial smoothing approximation’ the non-local character of the correlation is reduced to a local structure. The feasibility of this approximation can be examined by analyzing the autocorrelation function of velocity fluctuations in combination with measured integral sizes for two-point correlations which indicate the size of the energy containing large eddies. For example, Zaric (1972) used hot-wire anemometry to obtain the autocorrelation function in a channel at Re = 36,000 using air as fluid medium. His measurements which were made at a wall distance of y+ = 2.2 constitute measurements of 32 the fluctuations of the velocity vector such that the autocorrelation function represents some sort of average measurements composed of all three fluctuating velocity components. His measurements show that within At=4 ms the autocorrelation has decreased to about 10 %. The same time span was observed by Favre et al. (1957) in a turbulent boundary layer (i.e. y/5 = 0.25, Rex=Umxlv = 592,000) for a decrease of the autocorrelation (for longitudinal velocity fluctuations) to 10% (see also Schlichting, 1951). Given those numbers for the relevant temporal decay, the spatial decay of the convective-diffusive Green’s function can be estimated to 5'5025 mm. With 5:40 mm (channel flow; Zaric, 1972) and 8:25.] mm (boundary layer; Favre et al., 1957) the relative spatial decay lies between [78 = 0.0063 (channel) and 875 = 0.01 (boundary layer). With average sizes for the energy containing large eddies measured in channel flow at Re = 30,800 (Laufer, 1951) to 89:15 mm their relative magnitude (i.e. U5) is 0.236. Thus, it can be stated that by the time the turbulence has decreased sufficiently, the ratio of the relative magnitudes for the size of energy containing eddies and the decay of the Green’s function is still large enough for the Green’s function to be adequately represented by a delta function. In an isotropic turbulent flow field behind a grid (x/M = 48, = 7.7 m/s, ReM =25,000) a decay time of At = 10 ms could be observed (see Hinze, 1959). Under the same conditions an integral scale for the large eddy size of approximately 20 mm was measured. This, in turn, renders 8*50.4mm such that a direct comparison of L, and 8‘ for this flow field yields M" = 50 which remains in the same order as estimated for wall bounded flows. Therefore, the applicability of this ‘smoothing’ approximation is considered feasible and extended to the unknown statistical correlation appearing in Eq. (2.5). Therefore, the quantity :1 becomes independent of x and can therefore be taken outside the volume integral in Eq. ((2.4) to yield 33 < u'u'>= —jl=:](x,tlx_,i) jG(_x,tI3,E)dVdf. (2.6) t The autocorrelation of the two terms in :1 as given by Eq. (2.5) is modeled through the postulation of an empirical memory function 4). This function reduces the time dependence of 1:31 (x, tl x, f) and incorporates therefore intrinsically the character of the autocorrelation function. Formally this gives goals?) =¢<.>s,t—E)I==,= —TRFI(_x,tlx,t). (2.9) Since the dyadic quantity ;. contains the unknown correlation an analogous derivation is needed to find an expression for this term. Through a formal dyadic multiplication of Eq. (2.3) with f and subsequent ensemble averaging, a representation for the transpose of this unknown correlation is obtained: <£’.u'>= —J 16d,tls.i>_13,(x.tls.’f)dvdi. (2.10) Ev ’ The same reasoning to obtain Eq. (2.9) is applied to yield 34 = —’CRF2(x,tlx,t). (2.11) Eq. (2.9), which formally consists of two following terms, =—TR-V—Tk (2.12) can be rearranged to read -(1+1:RV)=—‘cR . (2.13) The same procedure can be applied to Eq. (2.11) to yield <£'_q'>-(l+TRV)=—TR. (2.14) Upon transposing Eq. (2.14) and inserting into Eq. (2.13) a formal preclosure formulation of the Reynolds stress is obtained: =fi§g (2.15) where the operator A is defined by é=(l+TRv)_l. (2.16) For small dimensionless relaxation times, i.e. TRIIV < u >|| << 1, (2.17) the operator A reduces to the unit dyadic. Thus, the prestress, which is formally expressed as I: < ff> , reduces to the Reynolds stress. 2.2 Turbulent Relaxation Time The turbulent relaxation time TR as introduced by Eq. (2.8) depends on a local turbulent time scale as well as the invariants of the mean velocity gradient: 1,, = 1:,(1, ,||V < u >||). (2.18) 35 IR is modeled as being proportional to the turbulent time scale “C[ as TR = cR(De,)’t,. (2.19) The explicit dependence of TR on the invariant of the mean velocity gradient is incorporated into the scaling function cR which may depend on the turbulent Deborah number Det defined as Det = t,IIV < y_ >|| = turbulent timescale . (2.20) mean field timescale The norm of the mean velocity gradient is taken to be: V= VzVT, (2.21) H and reduces to the absolute value of the mean velocity gradient for channel flow. Note that for channel flow this term is identical with the characteristic strain rate as defined in Eq. (1.4) The turbulent time scale 13 depends on local turbulent parameters such as the kinetic energy and the dissipation rate a and on the kinematic viscosity v to t, = t,(k,e,v) (2.22) k = —f R . e w< e.) fw represents a wall function which empirically incorporates effects as a solid boundary is approached. The local turbulent Reynolds number Re, is defined as k2 Re, = 72’ (2.23) and decreases as the wall is approached. In the inertial sublayer (i.e. y+ > 30) and the fully turbulent core region where Ret » 1 (see Tennekes and Lumley, 1972), this wall function fw assumes unity, i.e. fW(ReI >> 1) —> 1 . (2.24) 36 The turbulent Deborah number Det reduces therefore to the shear parameter kS/e as introduced in Chapter 1. In the near wall region where Ret —> 0, this wall function can be represented as fW (Ret —> 0) = c?” Re:1 . (2.25) where 0;, and n are empirical parameters to be determined from experimental and/or numerical data. In this region, the turbulent Deborah number can thus be represented as kS 0 if n>—% De, =ch:Ref = cfw/v/SS if n =-%. (2.26) °° 1f n<—% The relevant time scale to be used here will be determined from an asymptotic analysis for y+ —) 0 (i.e. Re, ——> 0). In previous research (Shih and Lumley, 1993) it has been argued that a turbulent time scale which uses n=-1/2 should be used for proper scaling. However, this research focuses on the introduction of a time scale which uses n<-1/2 in order to comply with the asymptotic behavior for the Reynolds stress as solid boundaries are approached. The direct numerical simulation data by Kim et a]. (1987) and Kim (1989) provide the resource for the specific determination of fw. 2.3 Properties of the Preclosure Several aspects have been incorporated into the development of the preclosure for the Reynolds stress and the turbulent relaxation time: 1) The non-local structure of the turbulent correlations represented by F1 and F2 is represented by a local structure through a spatial smoothing approximation which relates the duration for spatial relaxation of the associated Green’s function to the 37 relaxation character of the turbulence itself. This approach is based on the assumption that a high Reynolds number prevails. 2) The introduction of the empirical relaxation time IR is based on the assumption that the decay of the autocorrelation functions of all associated fluctuating quantities (Eq. (2.5)) can be represented effectively through a single parameter. This relaxation time TR incorporates the effects of the memory function and the Green’s function. 3) The consolidation of the divergence of the Reynolds stress fluctuations and the gradient of the pressure fluctuations into an effective force per unit mass constitutes an a priori measure. The adequate prediction of information contained within the individual terms of this quantity 1" is thus shifted to the closure hypothesis which models the unknown‘ correlation in terms of local kinematic properties of the flow rather than deriving explicit relations for those quantities. It should be noted that an analysis of the instantaneous fields from the DNS database may provide insight in the suitability of this approach and its benefits and/or limitations. If the prestress is represented as an isotropic tensor to tfi=2a —I, 2.27 -- 3 = ( ) some interesting features of this prestress can be extracted. The quantity 26L is given through the trace of the prestress and can be arranged to read 26L 2 2k +2tR(V < u >T:< 1_.1_' u’ >)+ tthr(V < p >T -< E' g' > -V < u>). (2.28) Since 26L can be expressed in terms of known parameters, Eq. (2.27) constitutes a basic closure for the Reynolds stress. For fully developed channel flow Eq. (2.27) reduces to 2d = 2k + ers < u;u; > +(TRS)2 < u;u; >, (2.29) where S is given by the mean velocity gradient to 38 _d dy S (2.30) The term 6t constitutes the kinetic energy for the prestress very much like k represents the kinetic energy for the Reynolds stress. It is composed of the kinetic energy k plus two additional terms. The first term represents the production rate of kinetic energy k multiplied by the local relaxation time whereas the second term represents the energy contribution from the normal component of the kinetic energy k. For channel flow the individual components of the Reynolds stress tensor as defined by Eq. ( 1.32) can be developed to - -——l———— (2 31) W 3 +(1RS)2 ’ ' —Ry, = IRSRyy , (2.32) Rn — Ryy = 0. (2.33) The streamwise component R22 is given by R22 = 1— 2Ryy , (2.34) because the trace of 3 assumes unity. Eq. (2.33) indicates that the isotropic prestress closure does not account for a secondary normal stress difference. Despite the fact that experimental and numerical observations show the existence of this secondary normal stress difference it is valuable to describe the potential of this approach since it provides the foundation for a more elaborate theory. Substituting Eq. (2.31) into Eq. (2.29) the kinetic energy of the prestress (i.e. (it) can be rewritten to read 26L = 2k(1— IRS(-—Ryz )) p (2.35) = 2k(1 — 2cR E" As Eq. (2.35) indicates it is the ratio of production of kinetic energy of the Reynolds stress to its dissipation which determines the magnitude of the kinetic energy of the 39 prestress relative to that of the Reynolds stress. The shear component can be written in terms of an eddy viscosity representation as given by Eqs. (1.2) and (1.3). The eddy viscosity coefficient 0,1 is given by c : 2fw(Re,)cR(Det) ” 3+[De,cR(De,)]2 ' (2.36) The limiting case for small values of IRS (i.e. in the channel center) yields the following expression for Cu: 2 Cu = 3CR(O)' (2.37) The algebraic form of -Ryz allows the identification of two distinct physical regions. For small values of IRS (i.e. IRS « 1) the shear component is proportional to the mean velocity gradient thus indicating a gradient type transport of momentum. For IRS » 1, Eq. (2.32) renders -Ryz to be inversely proportional to S from which an equilibrium type behavior can be inferred (see Chapter 5). The near wall analysis of the individual Reynolds stress components (see Appendix F) shows that the behavior of < u;u; > approaches the wall as O(y4) whereas — < u’yu; > behaves as O(y3). The two remaining components behave as O(y2). From Eqs. (2.31) and (2.32), it can be seen that IR is required to be inversely proportional to y to comply with the near wall behavior. With cR assumed to be constant without explicit dependence on the turbulent Deborah number Det, the exponent n in Eq. (2.26) assumes the value -3/4. The turbulent time scale It can thus be written as k -3/4 I, = -€-(l +cw Ret ). (2.38) The development of the wall function cw is subject to two constraints. This function controls when the near-wall time scale becomes dominant and thus becomes active in the region where molecular viscosity is important. On the other hand it needs to control the monotonic behavior of IRS in the near wall region such that Ryy and -Ryz reflect 40 monotonic behaviors, facts which have been observed experimentally and numerically (see Eqs. (2.31) and (2.32)). With the general acceptance that influences of molecular viscosity vanish beyond y+ > 30 (see Tennekes and Lumley, 1972), the following proposal is made for cw: ReT 2 Re:) )- (2.39) cW = ci’v exp(—( The turbulent Reynolds number Rel as defined by Eq. (2.23) is given in Figure 2.1. In order to ensure the vanishing influence of viscosity beyond y+=30, a reference value of Re1=50 is chosen which effectively reduces the value for cw to about 0.01% of its initial value for Ret=150. The monotonic behavior for IRS in the near wall region is attained by selecting a proper value for cfv. The DNS-data for 5+=395 are used to present Ryy in a log-log graph as a function of Its for different values of C3,. A value of ca, =15 is sufficient to obtain a monotonic single-valued functional relation for Ryy in the near wall region. However, to avoid very steep gradients (occurring around 1.3 E 3.3) a larger value of C3,; = 25 is chosen. The dashed line in Figure 2.2 indicates the asymptotic behavior of Ryy for large values of 1,8 (i.e. as the wall is approached). A similar behavior (i.e. monotony in the near wall region) is observed for -Ryz. Cross-references to the spatial distribution are indicated for Ryy and -Ryz. The explicit appearance of a single parameter CR does not likely render a satisfactory calibration of both the normal component Ryy and the shear component -Ryz as given by Eqs. (2.31) and (2.32). Since -Ryz directly couples to the eddy viscosity which in turn determines the mean velocity profile and thus the bulk average velocity (i.e. a quantity endowed with major significance) the shear component is used to determine an estimate of CR. Due to the quadratic appearance of CR in R”, a least square analysis rendered two values for cR, namely cR1=0.146 and CR2=2.534. These values thus provide a first quantitative assessment of the preclosure in combination with an isotropic representation of the prestress. From Eq. (2.31) a maximum value for -Ryz can be determined to 0.289. 41 536% =w>> :32 2: E 6% Mo 8:59:35 33mm Aim ocsmfi x Om + mm ON m_ O_ m 0 a _ d _ _ _ d _ _ _ _ 2 _ _ _ _ 4 _ _ _ d 0.0 OAVOBVO H o H 1 me - (IIIIIIIIIIIIIIIILWIIIIa H . H I O l r O .. 1| 0 1 CO— I O r H o o H r O O r lLllQldlllollullllllulll|lllllllllllllll Om— : O r .I o o o i l o 825/ ooceflom I II It com - rank 0 - - ofinb . - am I A _ - p p _ — — h — — — — p b — p _ — p — b — _ — OWN 14 — _ ’ o a ln(Ryy) : /2°V 12 — : /J _ o o v 2 10 _— f _. o n fi/v O C 0:5 _ w 8 -- o a A/v o : o a v 0 CW =15 ._ O D V 6 __ 0 :1I: vv A CW0=25 _ nove— y+= .. O 4 _ v CW =35 b y+=10 : — —- Asymptotes R =l/3 2 ~— _yy_ _ __ L _n I L I I 1 l I l l I l I I I I I I -5 0 5 10 ln(‘t‘S) 15 a) Normal Stress Component Ryy 1 ... -ln(-Ryz) : " O D/Bv 8 " / _ /—Jl " O D [‘V 6 ... O D/k/v o Cwo=5 - o 0 Av o _ o (:1 ,Av 0 CW =15 a - : Dav 04 — 0‘ AAVV + A Cwo=25 0 ~ Av<—— = c y o _ + v Cw =35 2 _ y =10 _ ———- Asymptotes I I I I I F I I I I I I I I I I I I I I l -5 0 5 10 mats) 15 b) Shear Stress Component -RyZ Figure 2.2: Dependence of the Normal Stress Component Ryy and Shear Component -Ryz on cwO (based on DNS-data for 82395) 43 ° Kim(1989;Re=7,890) - Kim et al. (1987; Re=3,250) R 22 C Isotropic Prestress Theory (Energy Path from Center E: to Wall indicated) 2: 3. 3 = 2 0 °: Q0 c °o 5.0 ° 03, c? o .00 0Q A $0 $4; \ / a. 0° \ \ ° / / ‘6‘ 0' \ \ e / / 90 g \\B / 332 )1 / \ / | \ / \ / I \ / \ / / I \ \ / | \ / I \ / \ / I Energy. \ / / l Equrpartition \ \ / | \ / I \ Rxx Two-Component Energy Ryy Figure 2.3: Energy Distribution Plane of the Turbulent Field in Channel Flow 44 for which TRS assume a value of 1.732. the occurrence of this maximum is intrinsic to the IPS theory and is not influenced by any particular value for CR. The two individual values (mat) for CR thus shift the occurrence of -Ryz within the spatial domain since they couple TR to It which is determined by the spatial distribution of k and 8. With the definition for the time scale used in this research, the turbulent Deborah number as extends over a semi-infinite domain (0 S ”as S 00) comprising the entire flow domain (i.e. 5+ 2 y+ Z 0). Independent of the choice for CR it is possible to illustrate the energy distribution path using a triangular plane as presented in Figure 2.3. The isotropic prestress theory is represented by the solid line originating at the energy equipartition point which for channel flow also represents an isotropic state since the shear component vanishes at the center. Due to the intrinsic feature of the IPS-theory that Rxx = Ryy the energy distribution path proceeds along the line representing axisymmetric turbulence as indicated by the arrow (i.e. prolate energy state). All the energy is ultimately transferred into the streamwise component Ru. The turbulence dose not achieve a two-component turbulence state at the wall as indicated by the experimental and numerical data. The deviation from the isotropic state at the center can be expressed in terms of an energy distribution parameter 613 which indicates the degree of anisotropy among the energy components. Thus, the normal components of the Reynolds stress are expressed as Ru=RW=%U—£L aAm Rn =%(1+208). (2.41) The parameter 08 therefore reads 61? = TRS(—Ryz ) cR < u;u; > S (2.42) — 2e ’ 45 mo 5% 8885mm— fion ES momnk H8 «ED mZQ wfims % coHoEmSm we 8:3ng Guam 3mm 8sz 84 k 0% com 8. o . _1..O.Q.O1CQOOO_OOO_OOO— _ d _ fl _ _ _ A q _ _ . _ _ OO oooooooooooooooooooooooooooooooooonH flow 1 .2. o 1 O o l Nd o 1 O I . ................ O l ............................................... o 1 Yo am 35 :m 0:: 88.52. 38cm o 1 2 . . o aze Mm-_uazoé .802 o 1 o 1 9o 0 . o 1 o . Gomnbv 3am- mzo ......... w- . Hm O 1 . $2 20 .. 01 mo .00 03 o I o o 1 coo. o one...oooooooooooooooooooooo 46 and can be interpreted as the ratio of production to dissipation. Figure 2.4 illustrates the spatial distribution of this energy distribution parameter 6? based on the DNS-data for 32395 for both estimates of CR. It can be seen from the figure that for the smaller value for cR the energy transfer into the streamwise component occurs very close to the wall whereas the larger value indicates the energy transfer occurring close to the center. The DNS-data included in the figure indicate an energy transfer of approximately 40% from the Ryy component to both remaining components. The entire energy transfer from the normal component into Rxx and R22 is initiated around y+=100. Note that due to an existing secondary normal stress difference, the energy is not just transferred into the R22 component as for the IPS-theory. The transition to a one-component turbulence state can also be seen in the realizability diagram which shows the path in terms of the second (113) and the third (IIIB) invariant of the normalized anisotropy tensor ; as a function of the turbulent Deborah number 1,8. The curve given by the values for 113 and H13 as presented in Figure 2.5 represent the entire channel domain here. A computation of the channel flow with the isotropic prestress formulation is thus redundant and serves merely as to cross-reference the value of as with a particular location in the physical domain (i.e. 1(S='ttS(y+)) and does therefore m); prohibit to infer the qualitative features presented here. The transfer of the entire kinetic energy solely into the longitudinal component is a direct consequence of the intrinsic incapability of this formulation to express secondary normal stress differences. However, even though this energy transfer seems unrealistic if compared with experimental or numerical evidence it does provide a closer representation of the turbulence states in the interior region of the channel. The most dominant feature is that this formulation provides a fully realizable algebraic turbulence model for the Reynolds stress in channel flow provided k and 8 are both positive. 47 1C, 2C and 3C indicate 1, 2 or 3-component turbulence state (Note: 3C indicates isotropy) DNS (6“:180) DNS (5“:395) Isotropic Prestress O + y: All turbulence 0.6 ~ state must lie 1C ' within the triangle 0.4L ()2 _ Isotropic Prestress IIB _ 0.0 '1 1 1 . . 1 0.5 :. 0.0 0.1 0.2 _ p 0.4 - Invariants of anisotropy tensor : IIthr(§-§) 0.3 ”_ urging-g; Z F 0.2 0.00 0.20 - IIB C 0.15 - 0.10 L 05 '_ DNS (5+=180) — 0 DNS (5“:395) _ Isotropic Prestress LIIII Illglilllllllll -0.01 0.00 0.01 0.02 0.03 111B Figure 2.5: Realizability Diagram for isotropic Prestress CHAPTER 3 PHENOMENOLOGICAL RELAXATION/RETARDATION THEORY The preclosure derived in Chapter 2 couples the Reynolds stress tensor with the gradient of the mean velocity through the operator A as defined by Eq. (2.16). The isotropic closure for the prestress rendered a first assessment of the preclosure formulation. However, qualitative properties such as the vanishing secondary normal stress difference requires an alternative closure. This chapter deals with a more elaborate closure approach and discusses a rational development of it. The calibration of the anisotropic prestress closure as presented herein is done in Chapter 5. 3.1 Motivation A closure for this anisotropic prestress formulation is presented here in the form of a frame-invariant constitutive equation which links the anisotropic prestress to the mean strain rate dyadic of the flow field. Constitutive equations have been used to link transport phenomena to the physical quantities to which they are related. Newton’s law of viscosity is an example for a simple constitutive equation relating the molecular stress in a fluid to the strain rate using the molecular viscosity as a proportionality factor. The Boussinesq approximation - given by Eq. (1.2) - has been used to relate the Reynolds stress tensor in a similar way to the strain rate in order to provide a closure for the Reynolds equation. 48 49 Here, the proportionality factor consists of an eddy viscosity which is related to the (turbulent) flow properties rather than fluid properties. However, this approach - which can be considered to be a linear expansion of the Reynolds stress in terms of the mean strain rate - has been found to be inadequate for various reasons (see Chapter 1). This research developed a turbulent model which extends the isotropic prestress theory (IPS- theory) in the following way: 1) Memory effects are included through the incorporation of relaxation and retardation effects (see also Chapter 5). 2) The explicit use of terms representing retardation effects entails nonlinear expressions in terms of the strain rate dyadic <§> through frame-invariant derivatives. One of the key ideas of using a constitutive equation of the form above stems from the observation of similarities between the turbulent flow of Newtonian fluids and the laminar flow of non-Newtonian fluids. Rivlin (1957) pointed out that the velocity profile of a turbulent (Newtonian) pipe flow has a similar shape as its laminar (non-Newtonian) counterpart. He further noted that the flow of a Newtonian fluid in a pipe of non-circular cross-section under fully turbulent conditions shows secondary flow patterns in planes perpendicular to the main flow direction. The origin for those secondary flow patterns - whose maximum velocity is approximately 1% of the magnitude of the downstream mean velocity - lies in the existence of nonvanishing normal stress differences (see also Speziale, 1982). Since non-Newtonian fluids also exhibit secondary flow patterns for flow through pipes of non-circular cross-section it was suggested (Rivlin, 1957) that the Newtonian fluid in a turbulent state may be regarded as a hypothetical non-Newtonian fluid. For viscoelastic fluids (i.e. non-Newtonian fluids) these secondary flow patterns in flows through non—circular cross-sections arise through the presence of molecular normal stress differences. In turbulent flow of Newtonian fluids it is the normal stress differences 50 of the apparent (or Reynolds) stresses which cause secondary flows to occur. However, differences in the occurrence of those secondary flow patterns do exist (see Tanner, 1985; Leonov and Prokunin, 1994). In square-duct flow, for example, the secondary flow of non-Newtonian fluids within the cross—section occurs along the diagonal from the comer towards the center (i.e. inwards) whereas turbulent flows show the opposite behavior (i.e. flow is outwards). Apparent normal stress differences have also been observed for simple turbulent shear flows (Laufer, 1951; Kreplin and Eckelmann, 1979). In analogy to this effect, the molecular normal stress differences are observed for laminar flows of viscoelastic fluids. Apparent differences in these simple flows are observed in the sign of the normal stress difierences. The primary normal stress difference in viscoelastic fluids (i.e. ’Lu-tyy)-is positive whereas the second normal stress difference is negative (Tanner, 1985; Leonov and Prokunin, 1994). For turbulent simple shear flows of Newtonian fluids the first normal stress difference tin—”Cyy =(—p)—(—p) (3.1) is negative, whereas the second normal stress difference is positive. From the general occurrence of the described normal stress difierences as well as the appearance of secondary flow patterns the use of constitutive equations for the description of turbulent, Newtonian flow seems compelling. Rivlin (1957) mentions that the 'stress tensor in general can be represented as polynomials in terms of the velocity gradients, acceleration and higher order derivatives. Linear constitutive equations arising through this idea were seemingly first proposed by Jeffreys (1929, see Oldroyd, 1958) for dilute suspensions. These models were characterized by three parameters: the viscosity; the relaxation time; and, the retardation time. Nonlinear approaches based on the above model were implemented by Oldroyd (see Leonov and Prokunin, 1994, Chapter 2). The necessity to incorporate relaxation times into a turbulence model based on constitutive equations has also been shown through the 51 experiments by Tucker and Reynolds (1968) and Choi and Lumley (1984). Their experiments in turbulent, Newtonian flows clearly show that the Reynolds stress possesses memory. Their observations indicate a finite relaxation of the Reynolds stress towards an isotropic state after a sudden removal of the strain rate (see also Chapter 5). 3.2 Closure Hypothesis for the Prestress The general formulation as set forth by Rivlin (1957) expresses the stress components of viscoelastic fluids as polynomials in the gradients of velocity, acceleration, second acceleration, etc. as 2 2 .I. — "Pl'l' (112“)‘1' (ll—fie) + 0‘32“) + (142(2) ’ (3.2) where the or, are scalars which may depend on the invariants of the kinematic tensors A“). A“) represents the mean strain rate dyadic <§> and the derivatives of these kinematic tensors are given by 8A _ =0) _ T . T . 211.11- at +2 Vg(,,+gm (V2) +Vy_ gm. (3.3) Thus, the kinematic tensors _A_(i) describe deviations of the stress tensor ; from its isotropic state. This strategy is employed here in a similar way to describe the anisotropic part of the prestress. The general form for the prestress is thus given by I + g , (3.4) isotropic anisotropic where the anisotropic part 2 will be expressed in terms of the mean strain rate dyadic. The tensor h is modeled as a traceless tensor. Thus, the mathematical expression for 52 the parameter 6L which was derived through the trace of the prestress (see Eq.(2.28)) does not change. The model for Q deviates slightly from the general form given by Eq. (3.2) inasmuch as memory effects for h are included, explicit quadratic forms of the mean strain rate dyadic are not considered and the general form of the derivatives employed are of a more general form. The constitutive equation for h; is hence set forth to 5. Edi-x1057 5_. .1. _ Q E. l] g tr(stg)3I)—B[<§>+Az(5t<§> tr(6t<§>))3l. (3.5) The trace of <§> has been omitted from Eq. (3.5) since this research only deals with incompressible fluids for which this quantity identically vanishes. Relaxation effects - defined as the response of the system (i.e. the prestress) to a change in the cause (i.e. the strain rate) are incorporated through the time scale 21 whereas retardation effects are implemented through the time scale 22. The parameter [3 acts as a viscosity coefficient for the prestress. The time derivatives in Eq. (3.5) are frame-invariant derivatives (Bird et al., 1977; Denn, 1990) They can be expressed for an arbitrary tensor to A to: SCA 8A T f=f+.Vé— -é—é-+p{<§>é+é<§>}. (3.6) <fl> represents the antisymmetric part of the velocity gradient. The left hand side (LHS) of Eq. (3.5) for which A is replaced by h leaves the parameter ‘a’ in place of ‘c’, whereas the RHS contains the parameter ‘b’ (A is replaced by <§>). Those derivatives are convected frame-invariant derivatives in reference frames embedded in the fluid which undergo the deformation of the flow. The choice of the parameter ‘a’ and ‘b’ depends on imposed restrictions on the general form of the resulting algebraic expression. The particular form of the derivatives in Eq. (3.5) and formally expressed by Eq. (3.6) is given when the algebraic expressions for channel flow are developed. It should be noted here that Eq. (3.3) constitutes the covariant derivative (i.e. c=+1). In flows of viscoelastic 53 fluids, for example, Oldroyd (1958) showed that the use of the contravariant material derivative (i.e. a=b=-1) successfully predicted the Weissenberg rod-climbing effect (this type of fluid is also denoted as Oldroyd-B fluid). This effect is caused by secondary normal stress efiects. It might therefore seem compelling to employ a covariant derivative (i.e.a=b=+1) for turbulence models in contrast to the findings for viscoelastic fluids. However, to date no conclusive evidence exists to promote either type of derivative for the description of turbulent flow of Newtonian fluids. The underlying ideas and assumptions for the choice of the derivatives (i.e. the choice of the parameters a and b) is presented at the end of this section. The parameters (or material functions, if a hypothetical non-Newtonian fluid is considered) M, B, and 2»; can be scalar functions of the invariants of the mean strain rate. The characteristic strain rate S, as given by Eq. (1.4), represents the second invariant of <§> and is the first nonvanishing invariant. Note that the first invariant is the trace of <§> and identically vanishes for simple, incompressible shear flow. For dimensional reasons the material functions need to be scaled properly. The scaling factors at hand are the characteristic time scale 11 given by Eq. (2.38), the turbulent kinetic energy k, and the dissipation rate 8. The scaling yields 7», = 01.1111 7L2 = c121,, B = 2kCB‘E, . (3.7) Through this scaling the dependence on the invariants of <§> is shifted to the dimensionless coefficients CM, c5 and cm. For the case of a simple, fully developed, steady shear flow (e.g. channel flow, pipe flow) the closure assumptions and the scaling for the material functions made in this chapter can be evaluated to yield the following algebraic relations for the normalized prestress g to 54 2 (ac,Ll - bcA2 ) + (a - b)ci1cm(t,S)2 H = c r s , 3.8 "" “( ‘ ) 3+(3— a12)(c,,r,S)2 ( ) H = 39. (1 sf (3+ 1))6,2 - (3+ a)cM + (b— a)(1+ a)c§,c,,(r,S)2 (3.9) ’y 2 ‘ 3+(3— a12)(c,,t,S)2 ’ Ha =—Hxx —H,,, (3.10) C _ _ 2 Hy, = ins? (ab 3)2C“C“(T‘SZ) (3.11) 2 3+(3—a )(c,,1.-,S) where 11 g = i° (3.12) Bearing in mind that 0=g) the isotropic prestress contribution reduces to one third of the unit dyadic. Thus, the contribution to the Reynolds stress stemming from the isotropic prestress formulation can be interpreted as consisting of an isotropic part and an anisotropic part (i.e. g) whereas the remaining terms including H“) and He) both represent anisotropic contributions to 5, thus i=5.“ +2‘”+2"’- (3.26) 1 l, 5; + g (3.27) 57 However, due to the formulation of the isotropic prestress, it is not possible to express this contribution in tensorial form explicitly in terms of Eq. (3.27). The closure hypothesis for the anisotropic prestress is implicit in the mean strain rate dyadic due to the implementation of frame-invariant derivatives. An algebraic formulation relating the prestress explicitly to components influenced by either the coefficient B or B12 is thus given here for the case of fully developed channel flow. The Table 3.1 shows how the individual Reynolds stress components can be decomposed into the various terms. It can be seen that the term containing B influences all components of the Reynolds stress. The relaxation parameter CR stemming from the preclosure and the dynamic relaxation parameter On both moderate the contributions of 05 on the Reynolds stress due to the explicit occurrence of Its in the denominator of the appropriate terms. By construction, the retardation parameter cm does not influence the component normal to the wall and the shear component. Provided B > 0 the magnitude of the shear stress -Ryz is reduced and constitutes therefore an ‘anti-viscosity’ coefficient for the Reynolds stress. This is in contrast to the shear component of the prestress where B acts as a viscosity. 3.3.1 Approximation for Small Time Scale Ratios An approximation to the earlier derived equations for the Reynolds stress can be made by assuming that the characteristic relaxation time as derived from the preclosure formulation is small compared with the mean field time scale (see also Chapter 2). This can be interpreted as the region of the flow domain where changes in the mean field happen very slowly compared to changes in the turbulence. Mathematically this is expressed as 58 5:233:50 3282: 058:0me 2: 0:: 08:0:BEE 2m mun 05 Tum £80883 RE. ”802 - ”8.23 + m 1. «8.3% + m 261m. “$.23 +m am 26 1 23:08.18 3 m mio 1 ”8.52626 1 uquew + m “$.16: m ”Queue + m a: ”$522 :6 1 me o 1 Na 9 3 i - «8.12% + m “8.923 + m his + m a: 052251 28 m _ “$52080 NAmfizovm + m ”$.28 + m N853 + m :M 0.32251 28 m _ 6,3 2303000 n 0585000 8:003:00 80:00:80 60:0 :003350 80:.“ 50:05:00 3050*: 030:8: mmohm mEOEAom 0:00:03 00 wash. 5 35:09:00 mmobm mEoiom 30:35:00 30505 0000:8844 0:“ ”am 638:. 59 tJVkxo. 62& For channel flow this region can be identified to be near the centerline where the mean velocity gradient vanishes and the mean field time scale approaches infinity. The mathematical description of this approximation neglects all terms in IRS which are quadratic or of higher order. For cR = 0(1) this restriction can be transferred to 118 being small. The components of the Reynolds stress as presented in Table 3.1 can therefore readily be reduced to the following subset. RH=RW=RH=%, 82% C C R =-i——& s. ' 330 fl (2 3XI.) ( ) Whereas the contribution to the isotropic prestress formulation as presented in Table 3.1 is moderated by IRS and the anisotropic contribution by IRS and CMItS, the approximation given here shows isotropic behavior. It is interesting to note that the shear component might serve as a first source to estimate a value for C3 in this flow domain. If one were to compare the expression for the shear components given by Eq. (3.30) with the common approach using the Boussinesq’ approximation as given by Eq. (1.2) with its definition for c”, the term in parenthesis can be expressed as C c C gwgc—ga wan This incorporates that 113 reduces to kS/e within the flow domain of consideration. This relation is, however, strictly valid in the close neighborhood of the centerline for which as « 1. The Boussinesq approximation is often employed with a value of c,Ll = 0.09. The DNS-data by Kim et al. (1987) yield cu = 0.1266 whereas the higher Reynolds number data render 011 = 0.0946 (data obtained from least squares analysis for 0.24 < ’CtS < 0.71) from which a fairly good agreement can be inferred for larger Reynolds numbers. 60 3.3.2 Other Phenomenological Models It is of particular interest to examine other turbulence models which are based on a phenomenological modeling ground. Very few models have been developed with a direct expression for the Reynolds stress in terms of a series expansion of the mean strain rate dyadic. Most models which derived some sort of algebraic prescription for the Reynolds stress used the exact transport equations for the Reynolds stress as a basis and introduced different kinds of physical assumptions to reduce the equation set to an algebraic prescription. However, some general ideas have been developed - even though they have not been implemented in actual computations - from a more general standpoint of phenomenological modeling. This section briefly describes some issues of these various modeling approaches. It should be borne in mind that every turbulence model presented here is associated with the transport equations for the turbulent kinetic energy and the dissipation rate. These equations, which are generally also subjected to modeling assumptions, constitute an integral part of the complete model. However, this section deals exclusively with the algebraic expression for the Reynolds stress. The associated equations for k and e are presented in Chapter 4. Thus, the conclusions drawn here stem solely from the algebraic Reynolds stress expression. Nonlinear Turbulence Model Using F rame-Invariant Modeling Speziale (1987) used an asymptotic expansion to obtain a nonlinear Reynolds stress model which constitutes an extension of Boussinesq’s approximation. His motivation to extend the linear Reynolds stress model was - besides the prediction of the anisotropy among the normal stress components in channel flow - the prediction of secondary flow 61 patterns in square ducts and an improved prediction of the reattachment point in turbulent flows over a backward facing step. The general form of his model is given by Eq. (3.32) to 2 1 gz—Epkl+p\/kl<§>+CDplz(<§>~<§>—-§<§>z<§>l) 8<§> _1_t 5<=S=>I (3'32) 51 3“ 51 )=)' + CEplz( The length scale appearing explicitly in Eq. (3.32) has to be solved for using the transport equations for the length scale or for the dissipation rate in conjunction with the definition of the length scale which is given as k% 5: C:— . (3.33) Speziale used the frame-invariant upper convected Oldroyd derivative in Eq. (3.32), which is given as a special case of Eq. (3.6) with c=-1. This derivative was included for consistency reasons since terms involving the derivatives of the mean strain rate dyadic consist of the same dimensions as the nonlinear terms in <§>. The calibration of this nonlinear turbulence model was done using experimental data by Laufer (1951). With the assumption that the parameters CD and CE introduced in Eq. (3.32) are constants, a one- point evaluation in the middle of the flow field yielded cD = CE = 1.68. (3.34) With this choice for the model parameters CD and CE, both the nonlinear term as well as the term containing the derivative of the mean strain rate dyadic can be consolidated to the following form 5c<§> _ _2_ 2L. _1- _ ;_—3pk;+pJEe<§>+ch1( 61 —3tr( 51 )1), (3.35) where the parameter ‘c’ determining the particular form of the derivative can be evaluated 62 to 02-1/2. The true nature of the form of the frame-invariant derivative and its basis are therefore changed. Thus, a more thorough investigation needs to be performed in order to evaluate the consequences of this form for the derivative and its use in the near wall region of channel flows. In order to evaluate the general performance of the model, Speziale (1987) did not solve any computational flow problem but took empirical data from Laufer to determine the length and velocity scales which in turn were used to algebraically predict the normal components of the Reynolds stress using Eq. (3.35). The results of this ‘performance preview’ will be presented in combination with the results from this research in Chapter 6. Shih and Lumley (1993) developed a more general form of an algebraic Reynolds stress model in terms of different orders of the mean velocity gradient following principles of invariant theory. The model, thus derived, contained eleven undetermined coefficients which - in general - can be taken as invariants of the tensors involved. These coefficients needed to be determined using further constraints and experimental data. Shih et a1. (1995) deduced an algebraic turbulence model from this more general form by truncating the series and retaining only quadratic terms in the velocity gradient V. Thus, the model for an incompressible fluid can be given to 2 k2 k3 =—kI—2c ———2c2-7{—-+-}. (3.36) 3 = “ e = e = = = = The tensor <§> denotes the symmetric part of the velocity gradient whereas the tensor <fl> constitutes the antisymmetric part. The parameters c“ and c2 are given by the following formulas: _ 1 J1—9cfi<8'%)2 “ ’ A, +A;(U‘%>’ C~ " C. +60% 10% > ’ c (3.37) where 63 .. 1 .. AS = J6cos¢ , and (b = garccos(\/6W ). (3.38) The parameters 8*, W* and U* constitute characteristic values which are defined as . ., (.): S=,/:, W=——E L 5—, 3.39 = = (<§>:<§>)% ( ) and U'=\/<§>z<§>+:. (3.40) The values of 6.5 and 1.0 are assigned to the parameters A0 and Co, respectively. With S as a characteristic strain rate according to Eq.(1.4), the individual components of the Reynolds stress are given to RM =%, (3.41) Ryy = %—%(583—)2, (3.42) R. = %+-C—22-<%>2, (3.43) _ Ry, 235:3 (3.44) From Eq. (3.37) it can be seen that for a vanishing velocity gradient (i.e. U*=0) clLl assumes a value of 0.154 (e.g. at the center of a channel) and represents a larger value as commonly used with the Boussinesq approximation. This formulation, however, does predict primary and secondary normal stress differences as can be seen from Eqs. (3.41)- (3.43). The model given by Eq. (3.36) is not intended to be valid in the near wall region where molecular influences become dominant. A comparison of this algebraic model by Shih et al. (1995) with other models mentioned in this section will be given in Chapter 6 together with the results for the theory developed in this research. 64 Algebraic Formulations from Reynolds Stress Transport Equations Some researchers have used the transport equation for the Reynolds stress as a basis for the development of an algebraic formulation. The exact equation of the Reynolds stress can be written as —D—<—E—E—>= —{< u'u'> -V < u> +(V )T' < E'E'>}+ < %{(VE')T +VE'} > Dt . r V . Production

k V Convection " Pressure/Strain — 2v < (Vu')T -V_t_1_ >— V - {< u’u'u'> -vV < u‘u'>} (3.45) Dissi‘pration Diffusion <2; (turb.+mol.) — {(V < 3'3 >)T + V < Q'E >}. \ p p l Pressure Diffusion <2>p Rodi (1976) noted that the eddy viscosity as introduced in Eq. (1.2) is direction sensitive in complex flows like, for example, swirling flows, and stated that the viscosity coefficient 011 as given through by Eq. (1.3) can not be a constant. He derived an explicit expression for the Reynolds stress by invoking an approximation for the convection and turbulent diffusive terms in the Reynolds stress transport equation. The molecular diffusion term has been neglected in his analysis due to the consideration of high Reynolds number flows. The modeling for the pressure/strain term is taken from the analysis of Launder et al. (1975) who modeled this term as ' 2 < %{(Vu' )T + Vu’} >= —cl -:—{< E!» —§kl} — y{< :13 > —%Pl}. (3.46) The term designated as <§> constitutes the production of the Reynolds stress due to the velocity gradient, as can be seen in Eq. (3.45). The symbol P denotes the trace of the “production” dyadic. The major assumption made by Rodi involves an approximation of the convection and turbulent diffusion as D Dk Dt _: k Dt D" . (3.47) = "‘ P—e, k ( ) where Dk denotes the diffusion of turbulent kinetic energy. With this approximation his expression for the Reynolds stress assumes the following form -k EI+1_Y ='3‘kl—2V1<§>+21m <§m>—§tr<§m>; . (3.50) m=l The first two terms on the right hand side (RHS) of Eq. (3.50) make up the Boussinesq approximation while the additional terms arise from Yoshizawa’s (1984) analysis. They I can be expressed as follows <§‘>=<§>-<§>+<fl>—-, (3.51) <§2>=<§>-<§>+<1>-, (3.52) <§,>=<§>-<§>—<§_fl>--, (3.53) where <§w_>s<§>.+(<§>-)T. (3.54) Tm = Cum—2. (3-55) The eddy viscosity v, is given by Eq. (1.3). For channel flow this set of equations yields the following expressions for the normal components of the Reynolds stress 1 1 k8 Rxx = 3 - g(Ct, + C13 )(—€—')2 , (3.56) l 1 kS R”, = 5 - E(Cfl — 2Ct3 )(‘:')2 , (3.57) 1 1 k8 R =—+— 2C —C —2. 3.58 22 3 6( 1:] t3)( 8 ) ( ) It can be seen that the coefficient C12 does not appear explicitly. The coefficients Cu and C13 have been given the values of 0.07 and -0.015, respectively. It is interesting to note that this nonlinear formulation predicts primary and secondary normal stress differences. For channel flow the shear component is not affected by the formulation given by Eq. 67 (3.50) and constitutes the same formulation as Boussinesq approximation. Thus, (—)- (3.59) Even though the general Reynolds stress formulation as given by Eq. (3.50) does not include further parameters which need to be adjusted, several components of the governing equations will be extended empirically using a wall damping function similar to van Driest’s (1956) damping function to extend the applicability of this Reynolds stress formulation to the solid wall. The shear component as given by Eq. (3.59) will be extended in the following way c f kS —R = “ d -— , 3.60 ,. 2 <8 ) < ) with fd given by rd =1—exp(—-):-) and A=5.2. (3.61) The results of this formulation will be presented in combination with the results of this research in Chapter 6. The modifications made to the transport equations for k and 8 are presented in Chapter 4. CHAPTER 4 TRANSPORT EQUATIONS FOR TURBULENT KINETIC ENERGY AND DISSIPATION RATE Transport equations for the turbulent kinetic energy k and the dissipation rate a are necessary to obtain proper scaling factors for the variables introduced through the closure for the Reynolds stress. Modeling assumptions for higher order terms in the exact equations required to obtain a closed form of the equations are elicited. The final transport equations for k and 8 associated with the phenomenological models presented in Chapter 3 are given here. Some additional approaches which have been discussed in the literature are introduced here for reasons of comparison. 4.1 The Exact Equations The exact equation for the turbulent kinetic energy k (see Appendix D) is given by aa—l:+< u> oVk—VVZk = — < u’u'>:V —v < Vu’:(Vu')T > I II III IV V —Vl-(< 3'3 > + < gig >). p 2 VI VH (4.1) The three terms on the LHS of Eq. (4.1) are the basic part of the structure of any 68 69 transport equation and represent the partial derivative of k (I), convective transport by the mean velocity (H), and the viscous diffusion (III) (see Hinze, 1987). Those terms do not need any modeling. The terms denoted as IV and V represent the production of turbulent kinetic energy through the work of the Reynolds stress against the gradient of the mean velocity field and the dissipation rate a, respectively. They both reflect the effect of a source and a sink term. Since both terms do not contain any new unknowns they do not need to be modeled either (it is assumed that an adequate closure for exists). The terms VI and VII do contain the unknown pressure/velocity and energy/velocity- correlation and therefore need to be modeled. The exact equation for the dissipation rate 8 (see Appendix E) is given by 985+ < u > -Ve— sze = -2v <(V_11')T ~Vu'>IV < u > -2V < V9."(VE')T >:(V < u>)T I II III IV V —2v < u'(Vg')T >EV(V < E >)T —2v < {(Vu')T . V3‘}:Vu'> VI VII — VV- < u’(Vu' )TzVu'> -2VV- < V%- Vu'> —2v2 < V(Vu‘ )T§V(Vu')T >. VIII D( X (4.2) The first three terms represent the substantial derivative of e (I and H) relative to the mean velocity and viscous diffusion (III). Like the k-equation, they do not need to be modeled. The first four terms on the RHS are commonly referred to as production terms. They are, in detail, mixed production (IV), production by mean velocity gradient (V), gradient production (VI) and turbulent production (VII). The next term designates the turbulent transport (VIII), followed by the pressure transport (DC) and the destruction of dissipation (X). All terms on the RHS need to modeled (Rodi and Mansour, 1992). 70 4.2 The Closure Hypothesis The closure for the transport equations for the turbulent kinetic energy and the dissipation rate is divided into different sections according to their respective modeled terms (i.e. transport, production and destruction). A comparison to existing modeling approaches is done where applicable. In order to obtain a model expression for the transport terms (i.e. terms VI and VII of the k-equation; terms VIII and IX of the 8- equation) a general formulation is given here which will be subsequently applied to the transport terms in the equations of consideration. This general expression stems from the equation for the fluctuating quantity analogous to that one for the fluctuating velocity u’ derived in Chapter 2. An approximation leads to the final model for the transport terms. The production and destruction terms for the kinetic energy equation do not need modeling. Thus, they will be presented only in combination with the final model equations. 4.2.1 Turbulent Transport Terms The general equation for the transport of a scalar 4) in an inertial frame of reference can be written in the following form a p%+pu-V¢-V-B,V¢=S,f”. (4.3) The first two terms on the left hand side (LHS) represent the substantial derivative of 11>. The third term is the transport due to molecular viscosity with 13¢ as the molecular diffusion coefficient and S¢M represents the net effect of all source and sink terms for 111 per unit volume. Upon decomposition of (1) and E into their mean and fluctuating parts, division through the density p and subsequent ensemble-averaging the transport equation 71 for the mean value of q) (i.e. <¢>) is obtained (Eq. (4.4): 8<¢> a1 +-V<11)>+—V-D¢V<¢>=. (4.4) D.» represents the ratio of the diffusion coefficient to the density and can be regarded as the molecular diffusivity (i.e. D¢=B¢Jp). Analogously, the value Sq, now constitutes the specific net source (i.e. per unit mass). If the mean field equation is subtracted from the transport equation for the instantaneous field, the evolution equation for 41’ is obtained to %+ < y- > .V¢'—D¢V2¢': —E'-V < (I) > —E'-V¢'+ < E'VCD > +8; . (45) Using the continuity equation (i.e. V ' g’ = 0), the second and third term on the RHS of Eq. (4.5) can be rearranged to represent the divergence of the fluctuating component of (3311’), i.e. - 2"V¢'+ < 2'-V¢ >= -V - (MY- < 2'¢'>). (4-6) The divergence of the fluctuating correlation between u’ and (j) in combination with the specific source term 8; will be denoted f; . Thus, Eq.(4.5) can be rewritten as in: —g-V < 111 > —f’, (4.7) with 58 a 2 I 1 1 1 1 I =§+-V—D¢V , and f¢=V°{u¢—}-S¢. (4.8) A formal solution to Eq. (4.7) can be obtained in the same way previously done for the fluctuating velocity (see Eqs. (2.l)-(2.3)) in terms of a Green’s function for the fluctuating scalar value (11’ to ¢'(x, 1) = —H G(x, 113, may < (11> 1.3151111. (4.9) In order to obtain a statistical expression for the transport of this quantity 111’ due to 72 velocity fluctuations, Eq. (4.9) is multiplied by u’ and ensemble-averaged. Using the same spatial smoothing approximation for the non-local statistical correlations and substituting the explicit time dependence into a composite time TR (see also Chapter 2) the transport term can be written as <_u'1’p'>=—'tR -V<¢>—IR . (4.10) The applicability of the spatial smoothing approximation made previously was based on a comparison of the spatial spread of the Green’s function during times for which the autocorrelation of the relevant turbulent quantities decays (see Section 2.1). Two observations were made for this comparison to justify the spatial smoothing: a) the diffusive character of the linear operator was given by the kinematic viscosity, and, b) the relevant decay of the autocorrelation function was taken to be same for all the turbulent quantities in question. The formal procedure to obtain a transport term for ¢’ (i.e. Eq. (4.10)) is applied to the transport equations for the turbulent kinetic energy and the dissipation rate (see Appendix C and D). For those equations the diffusivity, expressed in Eq. (4.5) as D», can be replaced by the kinematic viscosity v. It is postulated that the same composite time scale TR can be used to incorporate the temporal decay for the autocorrelations involved (cf. Section 2.1). To obtain an expression for the unknown correlation < _1_1_' f; > a similar procedure is performed which in essence parallels the development of Eq. (2.14). The starting point is Eq. (2.3) for the fluctuating velocity. A formal multiplication with f; and subsequent application of the spatial smoothing approximation yields vé_' =4, , (4.11) where _A_ is given by Eq. (2.16) and f’ denotes the divergence of the instantaneous Reynolds stresses and the gradient of the fluctuating pressure as they appear in Eq. (2.2). 73 After multiplication of Eq. (4.11) by A and insertion into Eq. (4.10), a representation for the turbulent transport of the fluctuating quantity 111’ by the fluctuating velocity is given to =—TR -V<¢>+rf, -é. (4.12) Thus, a representation for the turbulent transport of a fluctuating scalar quantity d)’ can be expressed as being proportional to the gradient of its mean value <¢> with TRl|ll representing the effective transport coefficient and the unknown correlation of the fluctuating term f with the quantity f; . As a first approach of implementing the concept developed and expressed in Eq. (4.12) the consequences of setting I: s(_) are explored. For isotropic flow fields, for example, it has been shown that any correlation between a vector-valued quantity and a scalar vanish (see esp. Hinze, 1987, p. 178 ff.). With 1,2, < f;_f_'>a 9, Eq. (4.12) reduces to a gradient transport hypothesis relating the flux of a fluctuating quantity to the gradient of its mean value to =—TR -V<¢>. (4.13) In order to apply the concepts developed herein to obtain transport expressions for the turbulent kinetic energy and the dissipation rate it is necessary to derive an equation which expresses the transport of the instantaneous quantity 4) rather than 111’. With the identity, =< E'{¢-<¢>}>=< u'¢>, (4.14) Eq. (4.13) can be rewritten to =—TR -V<¢>, (4.15) which relates the flux of an instantaneous quantity to the gradient of its mean value. 74 Turbulent Kinetic Energy As mentioned earlier, it is the pressure diffusion and the (instantaneous) energy diffusion terms which need to be modeled, (i.e., terms VI and VII). The approach made in this research is the representation of both terms using the gradient transport hypothesis (i.e. Eq. (4.15)) developed in the previous section, i.e. "'<1~1-'-2--'>=—ck'tR -Vk (4.16) The earlier used variable <¢> is replaced by k. The parameter ck is a dimensionless coefficient which has been introduced for adjustments of the transport term to experimental data and thus needs to be calibrated (see Chapters 5 and 6). This approach of combining both contributions to the turbulent transport follows common approaches applied by researchers using k-e-type modeling approaches. It is interesting to note Harlow and Nakayama’s (1967) approach to obtain a gradient type expression for the turbulent transport is similar in appearance to the one expressed by Eq. (4.16) (see Table 4.1). However, instead of including the pressure diffusion into the general gradient hypothesis they used the same concept of expressing the transport of an instantaneous quantity as being proportional to its mean gradient to the modeling of the pressure diffusion. Thus, they proposed the following expression, =—-—— , (4.17) p where 0' denotes the kinematic eddy viscosity. 0 and 7 denote dimensionless adjustable parameters. Besides the explicit incorporation of Eq. (4.17) into the modeling of the turbulent transport, it is the use of an isotropic eddy viscosity-type representation for the diffusivity which distinguishes their derivation from the one in Eq. (4.16). The approach as expressed by Eq. (4.17) was seemingly the first one to offer an explicit model proposal for the pressure diffusion term. Applications of this proposal were not pursued, though. 75 However, it seems possible to follow this approach as long as the pressure Poisson equation which can be derived by taking the divergence of the Navier-Stokes equation, i.e. —%V2p=Vu:Vu, (4.18) is satisfied. The model developed by Hanjalic and Launder (1972b) for boundary-layer flows uses the transport equation for the shear stress in combination with the transport equations for k and e. The k-equation in this context is derived from the contraction (i.e. the trace of ) of the general Reynolds stress transport equation which has been laid out in their derivation for more general flow problems. In the Reynolds stress equation only the triple velocity correlations have been modeled (see Eq. (3.45) - turbulent diffusion term, i.e. ). The explicit modeling of the pressure diffusion has not been considered. The , modeling they applied for the turbulent diffusion led to the following form for the turbulent transport of kinetic energy, I. ' k 9. I < u'-Ez—E>= —cs;{< u'u'>:V < u'u’> +< u'u'> -V <% >}. (4.19) From Eq. (4.19) only the second term in the brackets on the RHS compares with the approach presented in Eq. (4.16). However, the time scale here is taken to be the turbulent time k/e whereas in Eq. (4.16) it was TR (see also Chapter 2). The parameter cs is an adjustable parameter and was assigned a value of 0.08 obtained from computer optimization. The additional term in Eq. (4.19) arises from the modeling of the triple correlation in the Reynolds stress equation. For the transport of k normal to the wall (i.e. only the y-component is relevant) Eq. (4. 19) can be reduced to u'-u' k a < u; -“—2:- >= —cs ;{< u;u; > 53,—(k+ < u’yu; >)+ < uQu’y > B—y < u;u; >}. (4.20) 76 Experimental data were taken from homogeneous shear flows (see Hanjalic and Launder; 1972b) in order to develop an algebraic relation between < u'yu’y > ,< ugu’y > and k. With relations thus obtained and incorporated into Eq. (4.20), they arrived at: , u'-u_' k2 Bk = —O.8c —E——, (4.21) an approach commonly taken for the modeling of turbulent transport terms (see also Jones and Launder, 1972). This approach (Eq. (4.21)) has also been used by Speziale (1987), Nisizima and Yoshizawa (1987) and Shih et al. (1995) in combination with their respective nonlinear Reynolds stress models (see Chapter 3). In all of their models, however, the explicit modeling of the pressure diffusion term has been omitted. Nisizima and Yoshizawa (1987), however, introduced an additional variable parameter function fd as given by Eq. (3.61) besides the general adjustable parameter cs. This was done in order to extend the applicability of their model to solid.walls (see also Appendix F). Durbin’s (1991) expression for the transport term also uses the Reynolds stress in combination with a turbulent time scale as an anisotropic viscosity coefficient. The general form of his suggestion can be written as v =—:—°.V<¢>, (4.22) k where the diffusivity is represented by a tensorial quantity. The adjustable parameter 6k is introduced for adjustments of the transport terms to experimental data and given the value of 1.3. The only component of the diffusivity applicable to channel flow is the yy- component. This component applies to the transport equations of k and e as well as the momentum equation (see Eq. (1.7)). The diffusivity is expressed as vm = C11 < u’yu’y > 1,. (4.23) The turbulent time scale I. is taken to be the larger value of the two time scales k/e and 77 (v/ta)0'5 and the parameter cp was assigned a value of 0.2 in accordance with the analysis of the DNS data by Kim et al. (1987) and Kim (1989). This approach has also been used in Rodi’s (1976) derivation for the transport of kinetic energy. Dissipation Rate In the e-equation the terms denoted as VIII and D( represent the turbulent transport and need to be modeled. The strategy as outlined earlier in this section and represented through Eq. (4.15) is applied to yield a representation of these terms. The turbulent transport term VIII can be expressed as V=< ye», (4.24) where 8’ denotes the fluctuating value of the dissipation. Thus, application of Eq. (4.15), with <¢> replaced by 8, leads to v < u'{(Vu')T:Vu'} >= —c,1:R < u'u'> Vs. (4.25) The parameter cE is - analogue to ck in Eq. (4.16) - a dimensionless parameter which can be adjusted to conform with experimental data. The diffusional transport of e by pressure fluctuations is neglected in this modeling approach. This approach follows the reasoning of Hanjalic and Launder (1972b) who pointed out that this term leads subsequently to higher order derivatives which are presumed to be small compared with the remaining diffusional term. They derived a comparable form for the turbulent transport from an analysis of thin shear flows (i.e. boundary layers). Like before, they took experimental data from homogeneous shear flow to establish an algebraic relation between < u’yu’y > and k. Thus, Hanjalic and Launder (1972b) proposed a gradient-type model for the transport of dissipation by velocity fluctuations: 78 2 v < u'{(Vy_' )T:Vu'} >= —0.5cc L92. (4.26) 8 8y Daly and Harlow (1970) used a similar approach by setting the turbulent transport of 8 proportional to the gradient of its mean value. However, they proposed an explicit modeling expression for the contribution to the transport by the fluctuating pressure to < V% - Vg'>= G'V- < u'u'> (4.27) where 0* has the unit of s'1 and incorporates adjustable parameter. This approach follows in analogy to the derivation of Donaldson (1969) (see Daly and Harlow, 1970) who developed an expression for the pressure diffusion in the Reynolds stress equation as being proportional to the divergence of the Reynolds stress. Durbin (1991) used the analogue modeling of the turbulent transport as he proposed for the kinetic energy (i.e. Eq. (4.22)). However, (11’ is replaced by 8’ and the adjustable parameter - now denoted 0'g - assumes a value of 1.6. The explicit algebraic Reynolds stress models by Nisizima and Yoshizawa (1987) and Shih et al. (1987) both use a gradient transport hypothesis for the turbulent transport term which has commonly been applied in conjunction with the standard e-equation as given by Eq. (4.35). Nisizima and Yoshizawa (1987) used in their modeled transport term the same damping function as was used in the equivalent term for the k-equation. Rodi (1976) does not specify the exact form of the dissipation equation. The model of Speziale (1987) uses a transport equation for the length scale rather than the dissipation rate. The approach taken, however, resembles an analogue gradient transport hypothesis similar to the ones used for the dissipation rate. 4.2.2 Production and Destruction Terms The production and destruction terms of turbulent kinetic energy as denoted in Eq. 79 (4.1) as IV and V do not need to be modeled. The equivalent terms in the e-equation, however, need to be modeled and those approaches are described. The reason for the presentation of both terms here stems from the approaches taken in the literature (see for example, Hanjalic and Launder, 1972b; Mansour et al., 1989) where often the difference between those terms is modeled rather than their individual influences. The production of the dissipation in the exact transport equation for the dissipation rate consists of four different terms which are given by 0 IV: PEI = —2v < (Vu')T - Vu'>:V < q > (mixed production) 0 V: P2 = —2v < Vu'~(Vu' )T >:(V < u >)T (production by mean velocity gradient) 0 VI: P: = —2v < u'(Vu')T >EV(V < u >)T (gradient production) 0 VII: P: = —2v < {(Vu')T -Vu'}:Vu'> (turbulent production) whereas the destruction is given by the term denoted as X to o X: —y=—2v2 The modeling of those terms has undergone a variety of approaches. Some of the ideas represented in the literature are presented here. The model adopted for this research will be given in conjunction with the remaining modeled terms (i.e. transport and destruction terms) in the final equation presented in section 4.3. Tennekes and Lumley (1972) inferred that at high Reynolds numbers the turbulent production (VII) due to stretching of vortex filaments and the destruction of dissipation (X) due to viscosity tending to reduce instantaneous velocity gradients outweigh the other terms. Their difference, however, is in the same order of magnitude as the turbulent transport terms (VIII and IX). Launder et al. (1975) modeled the net effect of the terms VII and X as 80 _ _ 1 _ 8 P: — 'y = —cEl I — c£2 ? (4.28) with the turbulent time scale 1.} given by k/e. However, Hanjalic and Launder (1972b) model this net effect as being proportional to the second term on the RHS of Eq. (4.28). The mixed production of 8 (IV) as well as the production of e by mean velocity gradient (V) are closely related to the production of turbulent kinetic energy (Rodi and Mansour, 1992) and were therefore modeled as being proportional to it. Hanjalic and Launder (1972b) used this idea earlier and modeled those terms as P,‘ + P: = —c,, E < 11' y_' >:V < u > (4.29) whereas the destruction term 7 was modeled in conjunction with P84 to £2 13:- y = —c,, T (4.30) A revision by Hanjalic and Launder (1976) of this earlier developed model yielded an approach similar to that used by Launder et al. (1975) as presented in Eq. (4.28). However, the revised model contained an additional function fE as a multiplicative factor for the second term on the RHS of Eq. (4.28). This additional function incorporated the effect of the local turbulence Reynolds number Ret as given by Eq. (2.30) and assumed the following form 0.4 1 f =1—— - —R 2 4.31 and was introduced through an analysis of experimental data for isotropic decay behind a grid. Tennekes and Lumley (1972) argued that the terms denoted as P51 and P82 are smaller than the other terms by a factor of Reto'5 and PE3 is smaller by a factor of Re, and could therefore be omitted for high Reynolds number flows. Therefore Hanjalic and Launder (1976) did not use an explicit representation for PE1 and PE2 in their revised 81 model version. They did, however, model PE3 as P: = Gin/E < u'u'>:{V(V < u >)T:V(V < E >)T} (4.32) which was suggested by Taylor’s (1915, see Hanjalic and Launder, 1976) vorticity— transport theory. The models for the production and destruction terms by Durbin (1991) and Shih et a]. (1995) are of the form as given by the RHS of Eq. (4.28) whereas the model by Nisizima and Yoshizawa (1987) is given by the following expression 'Me . e P,‘ —y=2ce,k<§>-<§>—c£2fd21— (4.33) i=1 t The turbulent time scale used in Durbin’s model is the same as employed in Eq. (4.23) whereas Shih et al. and Nisizima and Yoshizawa use k/e. The additional function fd in Eq. (4.32) is the same as used for the eddy viscosity in their algebraic Reynolds stress model (see Eq. (3.61)). The introduction of this function has its origin in the otherwise divergent behavior of this term when a solid wall is approached (see also Section 5.6). 4.3 The Final Model Equations The final form of the transport equations for the kinetic energy and the dissipation as they are used within this research are presented here. The kinetic energy equation can thus be written in its final form as dk The unknown parameter ck explicitly appearing in this equation is subjected to calibration against experimental data and its determination will be explained in Chapter 5. The equation for the dissipation rate is given to 82 §+ < u > -Vs — sze = V (ch < g'u'> Vs) —c,r,, "1' - —chD;-. (4.35) The time scale It is given by k/e. The parameters cp and cD are used for calibration against experimental data (cf. Launder et al., 1975; Hanjalic and Launder, 1976). The functions fp and f1) are introduced to comply with the asymptotic behavior of those terms when a solid wall is approached. Details for their development are given in Chapter 5. 4.4 Discussion For a better understanding of the variations in the individual modeling terms of the transport equations the different model approaches have been summarized in Tables 4.1, 4.2 and 4.3. The tables include the so-called standard k-e model because it provides the basis for most variations developed by different researchers. It can be seen that the turbulent transport terms in the standard k-e model are expressed as being proportional to the gradient of their respective mean values. The effective diffusivities in those gradient- type transport terms have commonly been represented by an isotropic eddy viscosity which does not take a directional dependence into account. This approach, however, has been successfully used for modeling purposes and proved to be adequate for high Reynolds number flows and is frequently used in present Reynolds stress models using transport equations for the turbulent kinetic energy and the dissipation rate. The first modifications to better incorporate near wall effects can be found in the model of Jones and Launder (1972) some of which were introduced because of computational advantages. The function f“ (see Table 4.1 and 4.2), for example, was used to decrease the magnitude of the (predicted) eddy viscosity in the near wall region and was therefore ‘characterized’ as a damping function. The function fa served the same purpose as flLl 83 Aw _ .0 .cm 0: .000 :9 I I 8:088:00 80:: 0.0 £80003: 0:; 05. A.:.: v 093.0 I Amwa— .._0 :0 .000: 003 a: 0008 0 0. 800:0:m VS 0 3.0 31.6 M8.0.1.0 N: 1:10 £108 ém 0: .000 10 . 0 3.0 _H0_.O mmm_ :0 :0 :Em xDM—imml 4.0100 x. :16 0.01.0 1 1 . ..0 AW\C0.&¥08 H .0. .03 :58D 05. A.=.: v P11101I 20.9 .0: 9 .80 a $0. 351.01 a; 1 oodnxo .030N_:mo> 0:0 085$ Z 00. 0 02:00:: :0: ”0 .30 $2 03003 VSIxosol . N 0 0806000 :0: ”...0 03— ._00.~_ VS. A.m..m v M :0 I .L :03— II II . II 0. mo o1 0 0008004 0:0 020.83: #05. A.:.= v + A.0.: v >.A.=.= v m. 01 360003 >20 ”.0 $0— .90 0 > 0050000 :0: 8 .> 6 0800002 0:0 30:03 05 I A a v blml 8000883: 00::0m 5600:0me 0002 N 1 1 Q . - 3; 0:0 _> m8:0.C A film: v + A. :M v 020008 00 0: 0800p. aw:0:m 0:0:E 80:58:. 0:: 00.: 8000:0010 80E€0E 05 :0 00:20:02 ”:0 030% 84 :0008 010— 0:00:05 0:: 00 5:000:30 80:80:: 0800 0:: 0m: 3003 88:02 0:0 :00: .390: 80:00.: 0:0 0:083: :0 .0008 0:: 00: $3: ._0 :0 80:03 .0:0 0: :0:0: 0:: 00 :0: 0:0: 5:08:20 0:: :0: :000200 80:80:: 00:90 :0 00: :0: 00 Ammo: 2080mm 0:0 33: :00: .3100: 0800002 0:0 32:03 :0 £008 0: :1 ”0:02 Aw _ .NV .3: 0: .000 :9 I I . :0:w0: .0000 80:: .0 :0:00m0~: 03:1 05. A.:.: v 0:. 0I $00: .._.. .0 .20: 8.: .: 0 .0 . . .0 08 010. 0:00:0:m 0>Il I m._n b .oodnao 0 N0. ..:10 2.09 .Um— 0: .000 ..0 w> w 0b . 1. 82 1.30 85 11|1 m Tb N: ..0 4.01:0 a. :1... 0.010 1 1 . .0 $0.003... 1 .0 as: 55...: ”$1.... V fl..- .. 20.8 40m 0: .000 0: $0: 0>1m1 0:201 $0.01.... 0.2.02.3 05. 0:00:02 .._ . 1. £3: 1 1 .1. . 2 CI 0 .0080: 0:0 0:083: 0>. A0,: v M 0 I 880880: 00::0m 8:30:90”: :0002 1 a 1 . 1 1 . On: 0:0 :::> 08:0,: A.=> . N> v >m+ A.:>.:A.:>:.= v > 00000:: 0: 0: m8:0:. 80% 800885 0:: :0: 80:88:. 80058:. 0:: :0 w::0002 ”m0. 030 :1 85 .QN. ..v ..:m 000: 8508:8850 3080003: :3 00800080: 0.0.0:; 00.053. ...... .0.: :0 w::0008 0:: 350050: 3:0 E .1: ..0 o .> 8:0: :0: 8:30:98 20:98 :0 000 $3.: 88:02 0:0 .00: ..0008 0 0. 0:00:06 0:: 00 0:230:98 080m 0:: 0.8 @003 88:02 0:0 .00.: 0:0 $00.: ..0 :0 :_:m .0:0 0: :0:0: :0:: 00 :0: 0:0: 83093.0 0:: :0: :0::0070 80:80:: 5:98 :0 0m: :0: 00 Ammo: 0.0.025 0:0 63.: .03: A50: 08080.02 0:0 30.:0I ”0:02 80:05:00 0033..— 50 ...: ...o nemomom mEH : : Han—DUI II P .III Aka—OI. 0 A=v>uA.:.=v Ammo. ...0 :0 .80: 00w: .0008 0-0. 0:00:0:m 0.0 “a: .:~.0I A m v >“A.m.m vl .0.:0 ..:. ...: .8. :1... 0.... :1... .... .0 £1de ..0m 0: .000 10 I I 0 0 1 . I . . . I . mam. 5680.). 0:0 .00: A s v NDAA : v > . V.>:|>€0+N AAI : v ~>v>|1 080m 00 TN 0.3115 ... V_ a: V. Our—.0 ..O Nflmuo Km —H.uo :P N: :P .0 Ali, 0.4108”; .00. 8:85 I0- 0IAmv>uA.m.mv 0I . . . V II II 2me .Um— 0:. 00m : hwa 0. lw...:~.0lAWV.AWV0.:0N a .n_ 0 m. on: 0 0300.50? 0:0 08.0.0.7. 00 .I.. .I. :30. WSI 0. ..I o NI0 0 .m: .L 0 £00834 0:0 3.0.302 00 0 Amv DUAEEV w 0 80:08:80: 00::0m 8.305%: .0002 CA .:> .:> :0:0:: A :..m>:>m._:.m>:> v .31 Am? ....S . :93. v >0 1 :AA 0 v BRA :95: v >0 1 9 .>: 08:09: 0AA m v DVHA 1.35.3 > v >NI A m v DAMS. .0036 v >N I ”00.0008 0: 0: 08:0... 80% 8008005 0:: :0: :0.:0:::m0Q 0:0 :0::0000:n. 0:: :0 w::0002 ”m... .33. 86 whereas fez also served the purpose that the modeled dissipation term could be integrated to solid walls (see also Patel et al., 1985, for a review of near wall modifications). The representation for the eddy viscosity developed by Hanjalic and Launder (1972b) - even though its derivation was strictly speaking done for thin shear flows - proved to be a better choice for the diffusivity. The tensorial character of the transport coefficient accounts for directional dependencies of transport effects. This was also demonstrated by Durbin (1991) who used the recent direct numerical simulation data for simple shear flows to verify the adequacy of this approach. The transport terms derived in this research distinguish themselves in two aspects from the ones presented by Hanjalic and Launder (1972b) and Durbin (1991): 1) The derivation was done from a more general standpoint rather than for a specific flow field. The particular form used in this research stems from the a priori omission of an additional term resulting from this derivation (see Eqs. (4.12) and (4.13)). It is therefore surmised that their range of applicability extends the one previously presented. 2) The relevant time scale for the transport coefficient was derived from the Green’s function technique applied. This time scale allows to control the explicit near wall behavior and distinguishes itself from the time scales used in earlier approaches (see also Chapter 2). The pressure diffusion in the k-equation has commonly been incorporated into the gradient transport hypothesis. Approaches like the one by Harlow and Nakayama (1967) have gone unnoticed and were not pursued. However, if the explicit expression of Eq. (4. 17) does allow for the pressure Poisson equation to be satisfied, this approach might be a viable choice. The modeling assumptions made for the transport equation of the dissipation rate 8 has to most of its parts been accepted from the form in which it was originally given in Hanjalic and Launder (1972b). The assessment of the single terms appearing in the exact 87 equation for 8 has not been available until the recent direct numerical simulations (DNS) by Kim et al. (1987) and Kim (1989) for fully developed turbulent channel flow. Their simulations were done at Reynolds numbers of Re=3,250 and Re=7,890 based on the centerline velocity and channel half width. Mansour et a1. (1989) used the low Reynolds number calculations for the modification of the eddy viscosity for the Reynolds stresses. Rodi and Mansour (1992) used the DNS-data obtained at the higher Reynolds number for a systematic analysis of the modeling assumptions made in Eq. (4.20). A comparison of the sum of the mixed production (IV), production by mean velocity gradient (V), turbulent production (VH) and destruction (X) with their modeling in the standard k-e mode] as given by Eq. (4.28) revealed differences in the near wall region of the channel (5 < y+ < 20). The prediction of the net effect of turbulent production and destruction which certainly depends on the numerical values of cp and CD slightly overpredicts the equivalent terms evaluated using DNS-data. Therefore, the functions fp and f1) introduced into the standard k-e model overtake - besides the already mentioned purpose - the role of damping functions for those terms. Suggestions for fp and fD are made in Chapter 5. CHAPTER 5 MODEL CALIBRATION The model parameters appearing in the turbulence model developed in this research have to be evaluated using different flow fields. These flow fields are chosen as to isolate specific physical phenomena of the flow which are associated with a single parameter or a subset of those model parameters. Thus, their numerical values are determined by comparing the reduced set of equations resulting from the turbulence model with the corresponding experiments. The fact that the calibrated constants should be able to predict the class of flows from which they were calibrated is self explanatory. However, since the experiments used for the calibration are subjected to initial and/or boundary conditions it is clear that the representation of this class of flows through a single constant will not be able to recover all possible (i.e. subjected to all different possible initial and/or boundary conditions) experiments within this class with exact numerical predictions. The respective flow fields for the parameter estimates with a cross-reference to the equations in which they appear are listed below. The determination of the transport parameter ck is part of the numerical calculation of channel flow and will be done by matching the centerline value for k+ using experimental data by Laufer for Re=30,800. o Isotropic Decay cD Eq. (4.35) - Return-to-Isotropy c111 Eqs. (3.5), (3.7) 0 Channel Flow cR Eqs. (2.8), (2.18) CB, Cu Eqs. (3.5), (3.7) 88 89 o Equilibrium Region in Channel Flow c5 Eq. (4.35) o Asymptotic Homogeneous Shear cp Eq. (4.35) . Near Wall Analysis fp, f0") Eq. (4.35) 0 Channel Flow Calculation ck Eq. (4.34) The calibration of the function fD as it appears in the transport equation for the dissipation rate is done in two separate parts. fl) is composed of f}, , a function which can entirely be developed from the flow field of isotropic decay, whereas fl? constitutes a near wall function which becomes important for shear flows. 5.1 Isotropic Decay Isotropic turbulence represents a flow field where no mean strain rate is present and the turbulence - however generated - is homogeneous. It has been studied extensively since it reflects the simplest type of turbulence and a minimum number of quantities and relations are required to describe its structure (Hinze, 1987). However, it is also a hypothetical type of turbulence, because no actual turbulent flow shows true isotropy. The equations which represent a description of this flow field can be reduced from the general model equations (i.e. Eqs. (4.34) and (4.35)). The spatial homogeneity of this flow field requires it to be nonstationary. This can be seen from the transport equation for the kinetic energy in which - after the omission of the transport terms due to the spatial homogeneity - the only term left to balance the dissipation is given by the nonstationary term on the LHS of Eq. (4.37). A physical flow field which resembles very closely an isotropic decaying flow field is the downstream development of grid-generated turbulence. Batchelor and Townsend (1948a) measured the initial decay of the streamwise velocity fluctuations downstream of “)fD = flgfgv (see Sections 5.1 and 5.6) 90 grids with different mesh sizes. The isotropic nature of their flow field was implied. No measurements of lateral velocity fluctuations were made to verify whether the assumption of true isotropy holds. From the analysis of the experimental data, they concluded that the intensity of the streamwise velocity fluctuations (and thus k) decays asymptotically as < u? >~ ti, with n=1. (5.1) The experimental data did not allow to deduce a universal upper limit of the (temporal) extent of this initial decay period. It seemed that different initial correlation functions are likewise influencing this extent as does the Reynolds number ReM defined as M ReM — v (5-2) M denotes the mesh size used to generate the turbulence and denotes the mean stream velocity. Mohamed and LaRue (1990), who also tried to determine whether and how the decay exponent n can be related to the Reynolds number, mesh size and solidity, found no systematic dependence on these parameters. Measurements by Comte-Bellot and Corrsin (1966) revealed that a contraction behind the turbulence generating grid improved the degree of isotropy. They found that the exponent for the decay should be n = 1.25 rather than 1.0. The mathematical description of this flow field can be obtained from the modeled transport equations for k (Eq. (4.34)) and 8 (Eq. (4.35)) to 52% = —e ' (5.3) and de 82 a? : —CD ? , (54) The function fD as introduced in Eq. (4.35) assumes unity for high Reynolds numbers 91 such that its influence becomes negligible for large Re. An analytical solution for the decay of turbulent kinetic energy can be obtained to l k:ko{1+(cD—1)§k—°—(1—10)}‘-°v. (5.5) O From Eq. (5.5) in combination with Eq. (5.1) Batchelor and Townsend (1948a) obtained a value of cD =2. The proposed value of n=1 .25 by Comte-Bellot and Corrsin (1966) rather than n=1 suggests a value of cD=1.80. In view of the fact that the data obtained by Batchelor and Townsend (1948a) were most likely obtained in an anisotropic flow field, the value obtained by Comte-Bellot and Corrsin (1966) seems more reliable. Additional measurements by Comte-Bellot and Corrsin (1971) brought forth a value of cD=1.66. However, the fact that very few data were available to obtain this low value decreases the reliability of this low value. Therefore, a value of cD =1.80 is adopted for this research. Hanjalic and Launder (1976) reexamined the experimental data by Comte-Bellot and Corrsin ( 1966) and likewise concluded that the value of cD =1.80 is a reasonable choice. The final decay period of isotropic turbulence (i.e. small Ret) is characterized according to Batchelor and Townsend (1948b) by the law of energy decay as < u;2 >~ 1i , with 11:25. (5.6) Several other researchers (see Mansour and Wray, 1993) tried to determine the decay exponent for this low Reynolds number regime either experimentally or theoretically and values ranging from n=3/2 to n=5/2 have been assigned. Mansour and Wray (1993) computed the decay of isotropic turbulence using direct numerical simulation (DNS) in order to obtain a general understanding of the behavior of n for different Reynolds numbers Rex with Re}, being defined as 12 Rel : ‘1/< uV > A. (5.7) 92 A denotes the Taylor microscale and is related to the gradients of the velocity fluctuations such that it can be interpreted as the dissipative length scale for isotropic turbulence. Re}, can be related to the turbulent Reynolds number Rel as given in Eq. (2.23) in the following way 3 Re, = Eng. (5.8) The DNS-data revealed that for Re}, —> O (i.e. Re, —> 0), an asymptotic value of cD =1.4 is obtained. This compares with the derivation of n=2.5 in Eq. (5.6) by Batchelor and Townsend (1948b). However, for a different initial energy spectrum a value of CD =1.67 resulted. This value compares with Saffman’s (1967, see Mansour and Wray (1994)) analysis of n=1.5. The equivalent high Reynolds number results (i.e. Re}, —> oo) could be obtained to CD =1.83 and cD =1.7, respectively. These high Reynolds number values are within the range which was determined experimentally by Comte-Bellot and Corrsin ( 1966 and 1971). The variation of the parameter cD with the turbulence Reynolds number was usually incorporated into the transport equation for e by adding a multiplicative factor as expressed in Eq. (4.35) by fly. The proposal made by Hanjalic and Launder (1976) for the function fD was given as 0.4 1 2 fD =1- ‘1—.8—6Xp[—(g R6,) ], (5.9) which — in combination with their adopted high Reynolds number value of cD=1.8 - renders cD=1.4 for a vanishing Reynolds number. An important consequence of their model to describe the low Reynolds number decay (i.e. the final decay process) is the fact that it provides a realizable model, i.e. the turbulent kinetic energy and the dissipation rate vanish both simultaneously. This can be seen from a combination of Eqs. (5.3) and (5.4). Both equations can be rewritten to yield the following form 93 {ls-F = 3. (5.10) ko 8o The parameters kO and £0 denote the initial values for the kinetic energy and dissipation rate for a given decay process, respectively. For any finite positive value for cD the kinetic energy becomes zero when the dissipation rate vanishes, and vice versa. For the limit of vanishingly small Reynolds numbers (i.e. Rel—)0) in the final decay period the exponent in Eq. (5.10) is replaced by chD which approaches a value of 1.4. Thus, for the final decay (i.e. chp z constant), the differential equation describing the behavior for the turbulence Reynolds number as defined by Eq. (2.23) can be expressed as dRe k t=— —2. 5.11 dt V(cD ) ( ) From Eq. (5.11) it can be seen that for value of cD < 2 the Reynolds number decays until the kinetic energy becomes zero. This research adopts the high Reynolds number value of cD=1.80. The form of f9 as given by Eq. (5.9) will be taken as the part of f1) developed from isotropic decay, i.e. ff). The necessity of introducing an additional function f1? becomes clear when examining low Reynolds number asymptotes (i.e. Ref->0) for simple shear flows. Whereas the product of CD and fl; assumes a finite value for Ret—>O (here: cD fé-—>l.4), the entire dissipation term in Eq. (4.35) becomes unbounded because k—>O. The necessary additional aspects for fD (i.e. f3“ ) are developed in Section 5.6. 5.2 Return-to-Isotropy The Retum-to-Isotropy of homogeneous turbulence constitutes a flow field in which an initially strained flow field is subjected to a sudden removal of the strain rate. Experiments of this kind were conducted in distorting ducts (Choi and Lumley, 1984; 94 Tucker and Reynolds, 1968). Through the straining of the flow field in the ducts an anisotropic turbulence state is achieved. After the removal of the strain rate the development of the flow is entirely left to itself. Thus, this flow field demonstrates the self-interaction of a turbulent field since no shear is present and therefore tends towards an isotrOpic state. As noted by Lumley and Newman (1977) this decay towards an isotropic state takes place on the same time scale as the decay of energy. The return of the individual Reynolds stress components is therefore accompanied by a simultaneous decay of kinetic energy (and thus dissipation). For this flow field the Reynolds stress transport equation as given by Eq. (3.45) can be written as d < I l> ! ———§tg =<%{(Vu')T+Vu'}>—2v<(Vu_')T-V_u'>. (5.12) The pressure/strain - whose trace vanishes - is solely responsible for redistributing the energy among its components and is therefore a crucial part for the return towards an isotropic state. Lumley (1978) developed the following form of the Reynolds stress equation to d<§g> 2 =— —— I, . dt 82 38: (513) in which Q represents a dimensionless and traceless tensor which is defined as — sq =< %{(Vu' )T + Vu'} > ~{2v < (Vg' )T - Vu'> —%el}. (5.14) The term in the brackets on the RHS of Eq. (5.14) can be interpreted as the deviatoric part of the dissipation tensor. Thus, if the turbulence assumes an isotropic state the tensor g vanishes. It seems therefore reasonable to assume that 0 only depends on the anisotropy tensor Q (see Eq. (1.32)) and its invariants. Thus, a first order representation of the above yields 95 (b = C2 (5.15) Thus, Eq. (5.13) can be rewritten in terms of the anisotropy tensor B as 2 - (2 C)B (5 16) ch: — =' ' in which the dimensionless time dt has been introduced (dt=E/2kd‘t). Rotta (1951) suggested that the pressure/strain correlation can be represented as being directly proportional to the Reynolds stress. His model proposition rendered the parameter C to assume a value of three. Gence and Mathieu (1980) have found that the dynamics of the retum-to-isotropy is influenced by the sign of the third invariant of the anisotropy tensor B (see Eq. (1.34)). The comparison of their experiments with the ones conducted by Tucker and Reynolds (1968) revealed a slower return if the third invariant is positive (i.e. for the experiments by Tucker and Reynolds) which corresponds with a tendency towards an axisymmetric, prolate turbulence state (see Figure 1.7b). Lumley and Newman (1977) analyzed the experimental data by Comte-Bellot and Corrsin (1966) and found that for large Reynolds numbers (i.e. when viscous effects become negligible) the dynamic process of return-to-isotropy is entirely nonlinear. With this conclusion they argued for the apparent slow return rate for small anisotropic initial states. The operator 4 developed in the preclosure and given by Eq. (2.16) reduces to the unit dyadic. The quantity 26c (Eq. (2.27)) represents twice the kinetic energy. Thus, a combination of the preclosure with the closure as given by Eq. (3.3) shows that the deviatoric part of the prestress (i.e. H) equals the anisotropy tensor B. Therefore, the closure hypothesis as set forth by Eq. (3.4) can be evaluated to yield dB B(l—cM)+cM—E=O. (5.17) The high Reynolds number data of Choi (1984) are used to evaluate the model constant P 96 cM since they constitute experiments for which [[IB > 0. Both the experimental data sets for plane distortion and axisymmetric expansion yielded a value of on = 0.644. An evaluation of the experiments by Tucker and Reynolds (1968) for which H13 < 0 showed that CM assumes the lower value of cm = 0.295. This confirms the slower retum-to- isotropy for [H3 > 0. However, since Figure 1.8 shows that simple shear flows assume turbulent states for which IIIB > 0 only, the analysis of the experiments of Choi and Lumley (1984) are considered to be applicable here. 5.3 Channel Flow The following section presents the calibration of the algebraic form the Reynolds stress model. The particular form of the frame-invariant derivatives as introduced through Eq. (3.6) is controlled through the parameters ‘a’ and ‘b’, which are set forth to a=—1 and b=-3 (see Section 3.2). The algebraic form of the normal component Ryy is given to 1+ 2(1,S)H,, +3Hyy W ’ 3+(1RS)2 (5.18) With the specific form of the prestress Hyz as given by Eq. (3.19) and the equation relating Hyz to ~Ryz as given by Eq. (3.14) the shear stress can be written as c l —R =1:st "Ens . ’2 R' W 2 l+%(CMT,S)2 (5.19) The secondary normal stress difference, which is related to the prestress via Eq. (3.15), can be expressed as R” — Ryy = c,c,,(r,3)2. (5.20) For the evaluation of the parameters 0;; and c5 a least square minimization is performed which stems from the comparison of predicted values for Ryy and -Ryz and the numerical l I“ 97 data for 8+=395 (Kim, 1989). The quantity chosen for the minimization is expressed as: E: 2m?” — Rx”)? +2911?” —(—R;';“’S’))2 . (5.21) The data in the range of 30 S y+ S 395 were used for the estimate. Figure 5.1 graphically shows a contour plot for the range of cR and CB investigated. Two distinct areas are obtained for which E assumes a minimum. Included in the figure is the relation between CR and c5 as given by Eq. (3.31). This relation is valid near the center of the channel (i.e. 138 small) and for c,Ll = 0.09. With the idea that constant values for the parameters under investigation adequately describe the behavior of Ryy and -Ryz within the channel and the fact that one of both areas shows close proximity to the indicated relation a direct connection as indicated in the figure rendered CR = 0.428 and 03 = 0.195. The determination of cm is subjected to a more quantitative investigation. This becomes clear if one examines experimental and computational data both of which indicate that the stress difference Rxx-Ryy assume a finite value at the wall. For a constant op, Eq. (5.20) shows that CA2 should be inversely proportional to (1:18)2 for large values of 1:8. Thus, it becomes clear that the parameter cm serves as to retard the unbounded grth of the secondary normal stress difference as the wall is approached. A graphical representation of Rxx-Ryy in terms of 15 indicates a region for Its for which the secondary normal stress difference remains constant. This region occurs in close proximity to the wall and is indicated in Figure 5.2 as ‘Plateau’. In order to represent the functional behavior for on a two parameter base model with a functional extension for this ‘Plateau’ region is suggested. 0 Cu _1+b,,_(z,S)’ ”' (5.22) cm The functional extension is introduced in order to adequately describe the occurrence of the ‘Plateau’ region as indicted and is hence set forth as 2.0 98 1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4 lrlTl—TIIIITITfiIIITIIlrlIIll[lllrlll Ill IIIITTIIIIIIVIIIIITIII The solid line indicates the relation between CR and CB for small values of Its as given by Eq. (3.31) o. (53: W___,;ggb:o; //0t tr ‘ a I. '_ ‘. ‘.- ". - ~_' 0 'I ~- ‘I " -‘ -\. ‘.'- '-_‘ u i. |. ... l-. 1..“ ..H‘."' V... _... ....-""- »v-'-‘.V _.-._...~'-- ‘ ‘- I‘ l- I" AN“.- ---... . ....... - "H "' "-I'-.--' ~ ..i _... .<-" I - ‘ ' ' . _ H __....~ . .-__-“'__..-“ .. V, ‘-.. '.-"" .‘~""VV_."...- - .V’.'."'.. lllJlLLllJlLliJllllllill IIIIIIIITIIIII ll llJlilJLllLJLlllllllllJJLllJlllJllJ lllllLllllJLll 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 Figure 5 . 1: Least Squares Analysis for CR and CB CR 2.0 99 b2 T S — f,2 = l — aexp(—(—‘—:——) ), (5.23) with a=0.5, b=20 and c=40. This functional extension decreases the magnitude of cm by a factor of 2 for ”as = 20 such that a better representation of the indicated region is achieved. The function flz does not influence the asymptotic behavior for large values of as. For this asymptotic approach the secondary normal stress difference as given by Eq. (5.20) can be expressed as l2 0 c° —1n(R,, —Ryy)=—ln{ EM}, (5.24) from which the ratio of c312 to bu can be obtained. This analysis requires though the knowledge of a representative finite value for the secondary normal stress difference at the wall. From both DNS-data sets (i.e. Kim et al., 1987; Kim, 1989) an average value of (Rxx-Ryy)”=0.25 has been taken. The analysis for small values of 13S (i.e. Its—)0) yields the following form for Eq. (5.20): —ln(Rxx —Ryy)=—ln(cflc§:2ff2)—2ln('c,S). (5.25) With a value of ffz =0.61, the analysis yields 9:2 =0.083. Therefore, bu can be estimated to be 0.082. The functional behavior of the secondary normal stress difference as given by Eq. (5.20) is indicated in Figure 5.2 by the solid line. With these parameters the algebraic structure of the model is given. As for the isotropic prestress theory (see Chapter 2), the entire channel domain is characterized for 15 extending over the semi-infinite domain 0 .<_ 15 S oo. The energy distribution in terms of this parameter as is graphically illustrated in the energy distribution plane as given by Figure 5.3 and the invariant diagram as given by Figure 5.4. The arrow in Figure 5.3 indicates the energy distribution path for increasing values of 13S (i.e. from the channel center to the wall). The incorporation of the second normal stress difference into the 100 »» m-:m monocota $25 3852 ”N6 2%?“ D G x AsTmS 2.82m-.. x .3835. Ammo E 9 3288 so we .85 . awe .5 2 3288 fuzz 1 2 8ng 8an348. Ill . >52; 30:35 ofioboflcxx ll 1 S awe .Emvm ”momnkv oawfinom o - Ewe ...,N a 52 sink 838m . . Exam: 101 0 Kim (1989; Re=7,890) ° Kim et al. (1987; Re=3,250) Anisotropic Prestress Theory Isotropic Prestress Theory (Energy Path from Center to Wall indicated) Energy Equipartition \ \ R Two-Component Energy R Figure 5.3: Energy Distribution Plane of the Turbulent Field in Channel Flow with the Ansiotropic Closure for the Prestress 102 All turbulence 0.6 * state must lie ’ within the triangle 1C 1C, 2C and 3C indicate 1, 2 or 3-component turbulence state a) (Note: 3C indicates isotropy) 0.4 ~ [ Isotropic 02 _ Prestress . DNS (5+=180) . o +— . Anisotropic IIB _ DNS (5 —395) Prestress ‘ Isotropic Prestress ’ ' 0.5 - . . 0.0 A A . A A . A A A A A A _ Anisotropic Prestress 0.0 0.1 0.2 _ 0.4 — Invariants of anisotropy tensor : HB=tr(§-§=) : IIIB=tr(§'§'I=3__) 0.3 A: 0.2 - i I 0.00 0.05 0 10 0 15 111B 0.20 — a" . IIB - ”géi‘ ' ’ _ x 0 e _ (38’ go o“ ' \9 0.15 — \9 , 8 C) ” Q6 )— g. ‘0 : o. ' . $09 0.10 - O h w + 6g — - DNS (5 =180) 6» I , 0 DNS (5“:395) 965¢ ‘ 05 __ ,‘ Isotropic Prestress {‘5 j" Anisotropic Prestress .’:‘ 3C-I Llllllllllllllll -0.01 0.00 0.01 0.02 0.03 IIIB Figure 5.4: Realizability Diagram for Anisotropic Prestress 103 anisotropic formulation of the closure (of. Chapter 2) causes the energy distribution to deviate from the path for the isotropic prestress theory. Initial deviations (i.e. for small values of 138) proceed closer to the prolate energy state and are therefore in better agreement with experimental data. The deviations from the isotropic prestress theory become significant very close to the wall where a two-component turbulent state dominates. A finite secondary normal stress difference at the wall emphasizes this two- component character. This effect, i.e. the seemingly sudden break-away from the original tendency towards a one-component turbulent state to a two-dimensional state at the wall, can also be observed in the realizability diagram (Figure 5.4). The modeling of the second normal stress difference through the retardation parameter 0),; qualitatively captures this effect. The arrow indicates the path of the invariants for increasing values of {(8. It is noteworthy that for small values of Its the anisotropic prestress theory is in better agreement with the numerical data than the isotropic theory in terms of the path of the invariants. For clarity, cross-references to the spatial locations for numerical data are omitted from Figure 5.4. Those are included in Figure 1.7 (see Chapter 1). 5.4 Asymptotic Homogeneous Shear A homogeneous shear flow is characterized by a constant lateral velocity gradient and a homogeneous distribution of turbulent statistical quantities. This flow field is widely used to calibrate constants for turbulence models since the homogeneity eliminates the transport terms (Tavoularis, 1985). Through the presence of the (constant) velocity gradient the direct influence of the production of turbulent kinetic energy on the magnitude of C1) can be determined (see Eq. (4.38)). For the calibration of this parameter only the asymptotic behavior of this flow field is considered. The transient behavior may 104 be used to verify the existence of the asymptotic state and to extract additional information about the dynamic behavior of the theory developed in this research. This goal, however, is not pursued here. The governing equation for the turbulent kinetic energy for this flow field can be written as 3_k_ —:V—e. (5.26) The implementation of homogeneous shear is done in a channel with a shear flow generator. The constant velocity gradient is generated by a ‘honeycomb’ generator whose resistance varies over the height of the channel. This experimental setup used to implement this flow field results in a developing kinetic energy profile along the downstream direction of the channel. Taylor’s hypothesis (see Hinze, 1987) is used to relate the time derivative in the LHS of Eq. (5.25) to a convective derivative. Thus, Eq. (5.26) can be written as ~Vk=—:V—e. (5.27) The spatial homogeneity in the lateral direction which has been verified experimentally (Tavoularis and Karnik, 1989) to hold reasonably well yields the fact that the transport terms for this direction could be neglected in Eq. (5.27). Experimental data show that the omission of the transport terms in the downstream direction can be justified a posteriori. As Gibson and Kanellopoulos (1987) pointed out, the implementation of statistically stationary, homogeneous shear flow is strictly speaking not possible. The reason becomes clear as one analyzes Eq. (5.27). Whereas the RHS of Eq. (5.27) is independent of the lateral coordinate (the Reynolds stress and the dissipation rate have been shown to be constant, and the velocity gradient is constant through the experimental setup) the LHS depends explicitly on the lateral coordinate through the presence of the convective term (the velocity varies linearly with the lateral position). It is for this reason that this flow field is commonly termed a ‘nearly homogeneous’ shear field. 105 The turbulent kinetic energy k and the dissipation rate 8 have been found to grow exponentially at the same rate in homogeneous shear flow (T avoularis, 1985). Thus, the ratio of k/e develops an asymptotic state. The asymptotic behavior of k/e can therefore formally be expressed as E. £)=lgl(.__kz_d—EEO. (5.28) dt 8 edt 8 dt The two contributions to Eq. (5.28) are evaluated from the transport equations for k and 8 which can be written as dk d —=—< ’ ’> 2 S—8 5.29 dt uyuz dy ( ) and d8 8 ,, d £2 Ez—CP— dy —CD?' (5.30) The existence of an asymptotic state characterized by a single time scale k/e implies an asymptotic state for the Reynolds stress and the dimensionless prestress since both tensors are function of kS/e only (S=constant). From Eqs. (5.29) and (5.30) it can be seen that both contributions in Eq. (5.28) assume constant values in the limit as the turbulence structure assumes an asymptotic state. They furthermore provide an explicit algebraic expression for the asymptotic value of the shear component of the Reynolds stress to 1 l—cD -R = . ’2 2(kS/e)1-—cp (5.31) In order to obtain the possible asymptotic states for homogeneous shear, the equations for the preclosure (Eqs. (2.15) and (2.16)) and the closure hypothesis (Eqs. (3.4), (3.5) and (3.6)) have to be solved. The parameters ‘a’ and ‘b’ in Eq. (3.5) are chosen to a=~1 and b=-3 (see Section 3.2). The evaluation of the closure hypothesis yields the following expression for the deviatoric part of the prestress: 106 1 H =—c c 1S 2 , 5.32 W " "( ‘ ) 3(1+p)2+2(c,,1:,3)2 ( ) 3 1+p H =—c rs , 5.33 ’7 2 ”( ‘ )3(1+p)2+2(c,,r,3)’- ( ) C C Hxx — Hyy = —"——*l(¢,3)2. (5.34) 1+p The factor p appears because the partial derivative with respect to time in Eq. (3.5) does n_ot vanish as this was the case in fully developed, stationary channel flow. The factor p appearing in Eqs. (5.32) - (5.34) can be expressed as p : Cu __, (5.35) The algebraic form for the Reynolds stress components for homogeneous shear thus changes from the ones developed for channel flow (see Chapter 3) and can be expressed in the following way: 1 + 3cB(1:,S)2(cR(1+ p)—c,,) = , 5.36 y’ 3 +(1:RS)2 (3+(TRS)2)(3(1+ p)2 + 2(c,,t,S)2) ( ) 3 1+ p —R =1 SR ——ctS , 5.37 ’2 R W 2 B ‘ 3(1+ p)2 + 2(c,,«t-,S)2 ( ) c c Rxx -Ryy = 1‘ ’2 (1,5)? (5.38) 1+ p An alternative expression for the factor p can be developed using the k-equation to 1 dk kS =c ——=c 2—R ——l . 5.39 P M 8 dt 11{ ( yz) 8 } ( ) It can be seen that Eqs. (5.35) - (5.37) reduce to the algebraic Reynolds stress structure as developed for channel flow for p=0. The solution of Eqs. (5.35) - (5.38) is iterative because Eq. (5.30), which links the e-equation into the algebraic structure and thus the 107 shear stress to the parameters cp and CD, also needs to be satisfied. The value for cD used in this iteration process has been determined from isotropic decay to be cD=1.8. The result of this evaluation is the determination of the parameter Cp as a function of kS/e. If Cp is known, the theory predicts the existence of exactly one asymptotic state (see also Parks, 1997). Asymptotic states for homogeneous shear for values of kS/e as high as 6.0 have been reported (Tavoularis and Corrsin, 1981). However, more recent experimental results by Gibson and Kanellopoulos (1987) and Tavoularis and Kamik (1989) indicate asymptotic states with kS/8=4.27 and 4.18, respectively. Figure 5.5 shows c}: for the range 3.0wr—t wowd 02.0 0mm.0 mad mnmd 00.0 thoC .:o:mv~\:0m£0 N0; mm_.0 mmvd mmmd M200 5”? A00 Amwev fivaQSQF 00A m0_.0 0000 00— .0 0NN.0 2 .v Amv awa: ozEavCSNF mood m£0 00m.0 00N.0 0vm.0 0N.v A= uf. (5.40) The wall shear stress in Eq. (1.9) has been substituted using Eq. (1.13). From this asymptotic analysis it can be seen that the equilibrium region is often referred to as a constant stress region. The premise that the production and dissipation of kinetic energy are balanced can therefore be expressed as e . 1 "RV: = 2(kS/8)“" ’ (5'41) which provides an analytical expression for the magnitude of the Reynolds shear stress for a given value of kS/e. With the parameters cR=0.428, Cp=0.195 and cM=0.644 the expression for the shear component (i.e. Eq. (5.19)) is used to determine the equilibrium value for kS/e to 3.05. The fact that the production and dissipation of kinetic energy balance can be used in combination with Eq. (5.40) to relate the dissipation to the gradient of the mean velocity in the following way + _—_ 2*. (5.42) With a logarithmic profile for the mean velocity according to u+ zi—ln y+ +B. (5.43) The functional behavior of 8 can be given as 8’: 1 . (5.44) + icy 113 Eqs. (5.41) - (5.44) can be used to derive an expression for the eddy viscosity as defined by Eq. (1.25) to V _:= +, 5.45 v Ky ( ) which indicates a linear behavior for the eddy viscosity ratio in the inertial sublayer. The transport equation for the dissipation rate 8 (Eq. (4.39)) can be written as d k2 d8 8 , , d $2 _d—y(2CRCE?RWd—y):—Cpk dy —cD?. (5.46) The functions fp and fl) in Eq. (4.37) are assumed to have a negligible influence in the inertial sublayer where the turbulence Reynolds number is large. Therefore, in this analysis both these functions assume unity. With a constant lead term for k+ (i.e. k+=k+eq) in the inertial sublayer the behavior of Ryy is subsequently constant and can thus be removed from the derivative. Eq. (5.46) can therefore be rearranged to read d_(1d€) d c 82 (547) ——_—' _ —C _ _. . “2“ Ryydy Edy Pk dy Dk —2cRck With S+ derived from Eq. (5.43) and 8+ from Eqs. (5.44), the following analytical expression for cE is obtained (cD —c,,) =—_, 5.48 c5 2KZCRRWk+3 ( ) The values for Ryy and k+ have to be evaluated at (kS/s:)eq.=k+cq which was determined to k+eq=3.05. With cD=l.8, ep=l.51, cR=0.428 and K=0.4l, an estimate for cE=0.374 is obtained. With the eddy viscosity represented by Eq. (5.45) and the velocity profile by Eq. (5.43) the momentum equation as given by Eq. (1.7) can be expressed as d dy d , , d d—y(—)ZEW' )=0 (5.49) 114 This asymptotic analysis resembles the O‘h-order approximation. It indicates that the net change in the turbulent momentum is not balanced by the pressure gradient as illustrated by Eq. (1.7). Note that the molecular viscosity does not influence this result (y+ » 1). The ls‘-order approximation does not constitute an asymptotic analysis in the strictest sense since it allows linear terms in y/5 to enter the analysis while retaining the condition that y+ » 1 (see Appendix E). If the linear term is retained in the representation for the velocity profile (Eq. (E.4)), the momentum equation for channel flow, = —_ (5.50) is exactly satisfied. This analysis can be used to derive an analytical expression for the transport parameter ck to: 2 + a1 + cl 2b,x2ch;f_R;;A ' Ck— (5.51) The parameters a1, b1 and c1 denote the coefficients for the 15t-order expansion of the eddy viscosity Vt, the kinetic energy k and the dissipation rate a (see Appendix E). The appearance of c] requires knowledge of the behavior of e, a statistical turbulent quantity which is hardly assessable experimentally. Thus, Eq. (5.51) emphasizes the importance of direct numerical simulation data for modeling purposes inasmuch as estimates of parameters difficult to assess can be given. 5.6 Near Wall Behavior The algebraic structure of the Reynolds stress model is developed to comply with the correct near wall behavior as outlined in Appendix F. The functions fin and additional 115 provisions for fD introduced into the associated transport equations for the dissipation rate will be developed here in the context of channel flow. The algebraic structure for the normal component of the Reynolds stress is given by Eq. (3.13) to 1+ 2(Ac,,S)Hyz +3Hyy ’7 _ 3+(TRS)2 (5.52) The isotropic prestress formulation was used to develop the functional form of ”CR to extend the applicability of this theory for the near wall region. For large values of 1.8 the shear component of the anisotropic prestress, i.e. Hyz, behaves as (1.8)] whereas its normal component Hyy is on the order of one (see Eqs. (3.17) and (3.19)). Thus, the anisotropic prestress formulation renders Ryy ~ (1.8)2 or Ryy ~ y2 which complies with the correct asymptotic near wall behavior (see Appendix F). The correct near wall behavior for the shear component can readily be deduced from - Ry, = 1,511,, - Hy, (5.53) which behaves as being proportional to ('1:.S)'l or -Ryz = 0(y). For large values of as, Eqs. (5.52) and (5.53) yield restrictions as to the admissibility of combinations of cR and CB. For r.S—) oo, Eqs. (5.52) and (5.53) can be evaluated subject to Ryy > 0 and -Ryz > 0 (experimental data for Ryy and -Ryz indicate that both Reynolds stress components are positive in the near wall region), respectively, to yield 2 c. < 3 C“ (5.54) 3 cM —cR and 4 2 cB < ———9£—, (5.55) 3 2c“ —cR of which the latter includes the former restriction. With the selection for cR and c5 made based on the least square analysis, the compliance of these restrictions is given. 116 The pressure diffusion in the transport equation for the kinetic energy behaves as O(y) whereas the turbulent energy transport behaves as O(y3) (see Appendix F). The modeling approach taken to represent the turbulent diffusion in this research has been developed (see Eq. (4.16)) to be V-(< u'%>+)=—V~(CRTR -Vk). (5.56) The RHS of Eq. (5.55) behaves as O(y3). With the provision made in Chapter 4 that both terms are represented by a gradient type model it can be seen that the modeling approach taken here strictly accounts only for the correct behavior of the (instantaneous) energy diffusion term and not the pressure diffusion term. The modeling of the transport terms in the e—equation has been done in analogy to the modeling of the equivalent terms in the k- equation. The model expression for the transport term is given to vV- < u'{(Vu')T:Vu_'} >= —V-(c€'cR < u'u'> V8). (5.57) The RHS of Eq. (5.57) can be determined as O(yz). The LHS of Eq. (5.57) is, however, of O(y). A consequence of the model expression as it is developed in Chapter 4 for the transport term in the 8—equation is thus not capable of correctly describing the near wall behavior if cE is assumed to be constant. The modeling of the diffusional transport of e by pressure fluctuations which is - analogue to its counterpart in the k-equation - the dominant term as the solid wall is approached is omitted in this approach. The behavior of all production terms (see Eq. (4.2)) can be evaluated from a Taylor series expansion as 4 2?; = O(y). (5.58) In accordance with Eq. (4.35) those terms are represented by 117 2 P; = —c,f, . (5.59) 1:] TI The time scale I. - taken to be k/e - behaves as O(yz). For channel flow the only Reynolds stress component contribution to the production term is the shear component. Thus, for epfp of order unity, the RHS of Eq. (5.59) behaves as O(y). Therefore, the function fp is not used for adjustments towards the correct asymptotic near wall behavior. If the dominant terms contributing to all the production terms (Mansour et al., 1989; Rodi and Mansour, 1992) is examined, the function fp can be crafted to comply with numerical values. The high Reynolds number value for Cp (i.e. fp=1) has been as determined from homogeneous shear to be 1.51. To determine the equivalent value for small turbulence Reynolds numbers Re. (i.e. near wall approach in a channel) the numerical results from the DNS-calculations by Mansour et. al (1989) and Rodi and Mansour (1992) are examined. Eq. (5.59) can be rearranged to read d dy ' 4 2 P; = 2cPfPe(—Ryz) (5.60) 1:] The asymptotic behavior of -Ryz for large values of IS (i.e. small turbulence Reynolds number in channel flow) can be obtained using Eqs. (3.13), (3.14), (3.17) and (3.19). The numerical data available are reported in a normalized form (i.e. normalized with (u?/ v2 )) which renders the normalization of Eq. (5.60) necessary. For small Reynolds numbers an evaluation of Eq. (5.60) yields 4 . J5 7 + 2F,“ V00 : {EC‘")cPfPo(e:v )A}y , (5.61) i=1 with C(°°) given by c ct”) = 1+ 2—‘-’-(C—R— 2). (5.62) CM CA! The parameter fpo represents fp for small Reynolds numbers (Re.—>0). From Figure 5.7 118 mwm+Nwm+mm wagon. :28:on ”Nam 233m » om .. 2 o. o .1 _ _ d _ _ _ _ _ _ _ q . _ ACCOO 1 $28.3 5:5 cofiEEco: 8a ‘. H 9E“: 223355 2: c8 3% 2:. \W 11 0l1 wood 1 \ 1 - oz: 4Q - 1 o m- \ \ 1 Bed ... o o \ \ H O - o . . \ o\ - 1 o . . . x \\ 1 mac - o . . \ i - W O O O C C C \ \\ U 1 \ 1 owed I O O L - \ . H o o o \ H - a8. .5852 25 somvmamub o \ 1 Roe H 8%. .._a 339...:an 8.1% . ~22: \ H _ p _ _ _ b _ _ _ — _ _ h _ — _ _ — — omo.o 119 which presents the behavior of the production terms P; +P,2 +P: for the DNS-data obtained by Kim et al. (1987) and Kim (1989) in channel flow the coefficient for the linear term (i.e. the bracketed term on the RHS of Eq. (5.60)) is obtained. With Cp=1.51, an average value of fpo=2.66 was obtained. A proposal for the functional dependence of fp on the turbulent Reynolds number is made to fl, =1+(2.66—1)exp{—B—&}. (5.63) a The value for ‘a’ is chosen such that the function fp deviates by no more than 2% from unity in the equilibrium region which is taken to occur at Ret=150 (see Section 2.3). This criterion renders a to a value of 30. The function fl) is decomposed into two independent terms, ft!) and f9,” , respectively. The function f 11) has been determined from an isotropic flow field (see Section 5.1). f: is necessary in order to be able to adequately represent channel flow. The reason becomes clear by examining the structure of the destruction term in the transport equation for the dissipation rate (Eq. (4.35)). Since the kinetic energy vanishes at a solid wall, the destruction term becomes unbounded and can therefore not be integrated towards the wall. Furthermore, the DNS-data by Mansour et al. (1989) and Rodi and Mansour (1992) show that this term assumes a finite value at the wall. Therefore, the functions ff, and fgv are used such that both flow fields (isotropic decay for small turbulent Reynolds numbers and channel flow in wall proximity) are adequately represented. A formal representation of the destruction term in the 8-equation is given to 2 w§_ 'y = —ch[‘)fD k . (5.64) The function fl; is taken from the proposal made by Hanjalic and Launder (1976). The function f9)” constitutes the necessary extension in order to properly describe the near wall behavior in channel flow and insures that the destruction term can be integrated 120 towards the solid wall. A formal proposal for fl? can be given as 1 fw =———. 5.65 D 1+b(1:,S)2 ( ) The reason for the introduction of the turbulent Deborah number 1.8 stems from the necessity that ff)” has to assume unity in an isotropic decaying flow field and needs to provide a means to compensate for the kinetic energy appearing explicitly in Eq. (5.64). The parameter ‘b’ is used as an adjustable parameter to comply with the DNS-data for channel flow. With respect to the limiting behavior as a wall is approached (i.e. 1.8—>00; Re.—>0) the destruction term can be rewritten as 2% 7w = -CDf§; W. (5.66) The subscript ‘0’ indicates the use of values for which Ret—>0. S represents the mean velocity gradient and c3, arises through the time scale 1:.. Available numerical simulation data present the destruction term in normalized form similar to the production terms (i.e. normalized with (uf/v2)). Thus, after a normalization of the RHS of Eq. (5.84), the parameter b can be expressed as +2.5 8 b = — f‘” —‘1—— . 5.67 CD Do 7:110:38” ( ) With available DNS-data (7;, =-0.02 for 52180; 7;, =-0.03 for 52395) and S+ assuming unity at the wall, an average value of b=1.48-10’3 was determined. With the equilibrium value of 1:.S=3.05 the function 13 assumes unity within 2% such that b can be taken to be constant. 121 5.7 Summary and Discussion The parameters of the turbulence model developed in this research are summarized in Table 5.2 with respect to their origin. The determination of the transport parameter ck is integral part of the numerical predictions and is described in Chapter 6. The flow fields have been analyzed to the degree that was necessary for the determination of the model parameters used in this research. The algebraic structure of the theory with the time scale as introduced in Chapter 2 for compliance with the asymptotic near wall region has led to a realizable Reynolds stress model provided the kinetic energy and dissipation rate are both positive. The decay of isotropic turbulence has been decomposed into two temporal regions for the purpose of analysis, i.e. the initial and final decay period. The high Reynolds number experimental data (i.e. the initial decay period) seem to be fairly reliable since those data could be reproduced well through direct numerical simulations experimental data. The final decay period (i.e. Re.—->0) has been analyzed theoretically and numerically. The limiting behavior for Rep->0 developed theoretically to yield the parameter CD to be 1.4 was confirmed through the DNS-data for very small Reynolds numbers. The initial energy spectrum used as an input to the simulations by Mansour and Wray (1994) and the increasing lengthscale for very small Reynolds numbers endow the exact determination of CD with some ambiguity. The numerical proximity to the theoretical predictions, however, amplify the reliability of those simulations. The functional behavior for f1; was adopted from the proposal for fD made by Hanjalic and Launder (1976) and constitutes only one of many proposals available in the literature. Improvements of this function are primarily based on empirical fits to DNS- or experimental data (if available). The proposal adopted satisfies the limiting behaviors and provides the important feature that it renders a realizable decay process inasmuch as the turbulent kinetic energy k and the dissipation rate 8 vanish simultaneously (see Section 5.1). 122 48:3. 28:08:: 5:85 BEE—Bow no Buofifiam tommcfib ”E 2626632 .2499 Ea 8.9 .360 ....E 2 .ooe 2 use ea ..c 2885 ”e Amway .wm 8 9.6808 2.. cone—Sm ”5 A _ m.mv .cm 9 wEEoan no 25 go 5253 cog—om ”S Amwi .Emvm ”$2 5:32:30 32m .2520 so- no 8M3: EE A23: .3 8 .EM 30E 35:20 Qmm do Ammo: fiEavEmtho§H $25 mzoocoonom 600m; .6 Gem: EmonHB—omotcou >803 oEobofl 503.— no e3 Ba .6 .8 a6 .6 co mozaEEm sewed 635555 42d .6 6%; ea 5&2 . . €23 0 «Mo :95 So: osanxmm .._a go E50 San—-mZQ was Bo: .occwzo Emfl .o no E2213 5 cewe SHED EwNvd ~.o Ammo C 650 Eobofléfifiaom 340.0 20 cousom «SQ mafiwomhxofi 30E 0383mm “80:8qu $82 socmccwaoméoumxflom oEobomE< .50 £225ch ”NW 033. 123 The model developed for the flow field initially anisotropic through the presence of a shearfield returning to an isotropic turbulence state constitutes the precursor to the isotropic decay process. The time scales for the return process are of the same order as the mean field time scale such that the return to isotropy occurs simultaneously to the decay of kinetic energy. This flow field is modeled through the closure hypothesis set forth in this research as a linear model as given by Eq. (5.16). Lumley and Newman (1977) concluded that the rate of return for high Reynolds numbers is entirely nonlinear. The model set forth in this research to describe retum-to-isotropy does not account for those nonlinearities. However, the experimental data of Choi (1984) for plane distortion and axisymmetric expansion can be fairly well represented with the proposed form of the model (see Section 5.2). The calibration of the algebraic form of the Reynolds stress model has been done in strong compliance with DNS-data acquired for turbulent channel flow. The correct limiting behavior of the turbulence model for 1.8—>00 is achieved through the empirical extension of the time scale as presented in Chapter 2. An analytical expression relating cR and c5 is obtained from the asymptotic behavior of the shear component for small values of 1.8 (i.e. in the core region of the channel). The normal component Ryy does not enter the analysis for ItS—>0 since the structure of the model predicts an isotropic turbulence state whereas an anisotropic state with a somewhat smaller magnitude for Ryy is observed experimentally. The particular form of the time scale chosen provides a realizable theory for simple shear flow (provided k>0 and e>0). The retardation term introduced in the phenomenological closure prohibits the unbounded growth of the secondary normal stress difference. An important issue arises through the question of whether the general form of the new time scale I. used for the algebraic formulation of the model restricts the theory in its applicability in wall proximity. The asymptotic near wall behavior as outlined indicates a unique behavior for the Reynolds stress components. In a simple shear field 124 such as channel flow the components are destined to proceed in the range of 0 —0.936, 8":395 —> -1.088) and indicates a more ‘rapid’ energy transfer into the longitudinal component (i.e. Ru) than the DNS- data show. This effect may reflect the fact that the numerical results show an anisotropic turbulent state at the channel center whereas the anisotropic as well as the isotropic theory reflect isotropy. This can be seen in Figure 6.1 in the enlarged section for the center region. A transfer of energy to equal amounts from Ryy into Rxx and Ru would be attained both for cR=0.188 and cR=1.377. It appears that in the outer region (see cross-reference to the spatial locations) larger parameters for CR agree qualitatively better with the energy 128 255 53:55me Edam 55:5 beam Shoshana ”fie oEmE eewom =55 Pmnuao 111 consume || 33".; 1 1| 33m 0525ch 8% ucowoq 91mg KSJoug sword :omwom 3:50 mmmnb Lon. mop—2330-880 38% D ow _u+w .8 cocoBEMémocU Exam 4 Sednao was 1 .5 $0585 93.8254 '1 $31.8 mzn . 821.3 mzn . a _/ / //\%A.% _... / / 46¢ _.../ I // WONW ... / // \90 n... / , Oahu. / .. .w n L . _ ., a _ , Mm 81.4 %o _ an a. _ L L 90 _ ,, .H 10 " Tm 21..» ....w? _ zlq ..QOJ _ n. \. ow. _ F. pd _ o/..1\. «W. _ . .. _ mt. _ . Sn.» _ _ _ 129 states as illustrated by the DNS-data whereas within the equilibrium region smaller values show better agreement. In the near wall region, for which the influence of the second normal stress difference plays a dominant role (see also Chapter 5), it is the smaller magnitude for cR which provides more influence of the normal component Ryy. For small values of Ryy, it is the secondary normal stress difference which causes the two-dimensional turbulence energy state near the wall. The almost parallel energy distribution path of the DNS-data in the near wall region indicate that once Ryy is small enough the energy distribution takes place entirely among Rxx and R22 both of which are controlled by the secondary normal stress difference in this region (see enlarged region for the wall region). The value of cR=O.428 developed for the APS-theory adequately describes this behavior. The turbulence states associated with the parametric study are further illustrated in the realizability diagram as given by Figure 6.2. The tendency that for small values of 1.8 (i.e. close to the center) larger values for CR would represent the invariant path more accurately, as the energy distribution plane indicated, can not be confirmed. Small values for CR do not show good agreement with data. A value of cR=O.188 even causes the flowfleld to follow an initial path towards an oblate turbulence state for which IIIB=rw(1——g-). (6.2) The molecular viscosity was neglected in Eq. (6.2) since the outer region of the channel is under consideration. Similarly, under omission of molecular transport, the transport equations for k and 8 can be reduced to 0 d( k< ’ ’>dk) < ’u’>d e (63) =—cc— uu —- u — . d "Re ’y dy ’2 dy and O ( k‘< ,,>de) e< ,,>d 82 (64) =—cc— uu ——c— uu —c —. . dy‘Re de Pk H dy Dk The functions fp and fD (see Chapter 5) assume unity. Explicit influences of molecular viscosity enter through the boundary conditions (esp. in the equilibrium region). However, for the numerical implementation Eqs. (6.2) through (6.4) are normalized using the kinematic viscosity v, the friction velocity 11* and the channel half width 5 as given below. This normalization follows the practice often encountered for this type of flowfleld (see Yang and Shih, 1993; Demuren and Sarkar, 1993). 5n. nzi, 11*: 1‘1, 8+: u+:, 5 p V u. k 8V - < u’u’ > [(+2 M 8+: , —R =___Y_Z__, R :__Y__Z__. (,_5 u: u? ’2 2k ” 2k ( ) Eqs. (6.2) through (6.4) can thus be rewritten to read __1__du+ _c*(l-n) 6+ dn — gk+2 ’ (6.6) 134 d k+2 1dk" 1c+2 1 du" —— 2 -—R — = — 2 —5+ +, 6.7 dT]( Cch 8+ yy 5+ dT] ) g 8+ 8+ (11] ) 8 ( ) and d k+2 1 de“ , 1 du‘“ ,e” __(ZC£CR FR” 57-31?)=Cpgk g:( d“ )2 -C05 k+ ' (6'8) The normalized shear stress has been incorporated into Eqs. (6.3) and (6.4) using the following definition - R = ———'+—— , (6.9) where g is given by CB g=2c R — . R ” 1+—§-(ch8/8)2 (6.10) The velocity gradient S is given by Eq. (2.30). Eq. (6.10) resembles the analogue of the coefficient cu for the eddy viscosity as defined by Eq. (1.3). The boundary conditions in the equilibrium region at y+ = 30 are derived from the Om-order analysis as given in Chapter 5 whereas the boundary conditions on the centerline of the channel are given through the symmetry condition which states that the flux of k and 8 must vanish. Thus, 11 = neq = (%),q.: k‘“ = 3.225, 3* = 0.081, u+ = 13.8, (6.11) and = O and =0. (6.12) The value for neg. is found by substituting expressions for y and 8 in terms of y+ and 5+ (i.e. Eqs. (1.12) and (6.5)). 135 6.2.2 Solution Strategy The numerical scheme uses a finite difference scheme to approximate the differential equations for k and e. The formal implementation for the finite difference scheme includes the convective term in these equations which can thus be given as 316 a k+2 1 316 k+2 1 du‘“ , + + u BE, =an(2ck kc—R 8+ R”8—+ an)+g8+ 8+(dn) —88 (6.13) and +855 3 k+2 1 38+ + 1 du 2 +5” u BE =a—n-(2c cR— 8+ RW— 5* all —-)+cP g—-k 8—; cDS k” , - (6.14) where g represents the coordinate normalized by 5. The normal velocity component v+ =< uy > /u. is set to zero such that convection is only represented through streamwise terms. The finite difference scheme resulting from this formulation is implicit in the streamwise direction (i.e. z-direction) inasmuch as the diffusive terms are evaluated at the location of interest whereas the coefficients of the convective terms are evaluated at the previous streamwise location. The momentum equation is used in its integrated form to update the velocity gradient and the velocity at each new location along the channel. The source terms are linearized such that numerical stability is ensured. This approach follows “common practice” and can be found explained in detail in Patankar (1980). Details of this procedure and the complete implementation of the finite difference scheme is given in Appendix G. Convergence of the numerical solution is given when the dependent variables do not change in the streamwise direction, thus representing a fully developed state. To verify this, several local dependent variables as well as some integral properties are monitored. The monitored local variables are the kinetic energy, the dissipation rate, and the velocity at the centerline. The global properties chosen for monitoring are an average velocity (for 136 neq.(n)dn, (6.15) 8V. and values for the global production and dissipation given by 1 ‘ k” 1 du” , = — d 6.16 01. 1_neq‘ neq‘g 8+ 8+ dn ) Tl ( ) and 1 1 8+2 801.: j cD5+FdTL (6.17) 1- n... n, Two different computational grids have been investigated, one resembling the (cosine- distribution) mesh used by Kim et al. (1987) given through j-l M— n . n=(1-neq.)(1—cos( 13))+n,q, J=1,2,...,M (6.18) while the other is given using a power law expression according to j—l M—l n=(1-n,q.)( )" +n,., j=1,2,...,M (6.19) The difference between the global production and dissipation (i.e. the difference between Eqs. (6.16) and (6.17)) can be interpreted as the flux of kinetic energy at 11:11“). which can be verified by integrating Eq. (6.7) and applying the boundary conditions (i.e. Eq.(6. 12)). Thus, a connection between a local property of the flow to a global property is given. This fact was also used to investigate the influence of the different computational grids. It was found that a mesh size consisting of 150 grid points is sufficient for both meshes under consideration. 137 6.2.3 Results The computations for the outer region of fully developed channel flow are presented here. The unknown transport coefficient ck in the k-equation is evaluated by matching the value of k+ on the centerline with the one given by experimental measurements (Laufer, 1951). The experimental value was found to be k+CL=0.77. The fact that the level of k+ is directly influenced by its transport term led to the choice of this quantity for the determination of ck. The value thus obtained for ck was 0.558. The provisions made to a priori estimate this transport coefficient using Eq. (5.51) (see Section 5.6) is employed to investigate the feasibility of this approach. A least squares analysis for the DNS—data by Kim et al. (1987) and Kim (1989) for 30 :82 0350800 “0.0 oSmE 000— + 00_ 0_ _ .____H— u _ _____——_ _ ________ q o 10— 88.1% 33%: ......... 1 2 SENS boo; 10. 503230 I I 1 sank $2250 .. m Ginkgo: .36 EE 0 11 cm Sofinbv 002:1 cogxfiom I I 1 Gamubv 385 852.1% 000.3Hom o 00%?“qu n. O 4 1mm 1 00m.m_uom 00N.N_no.m _ _ _ _ _ _ _ ~ _ _ _ _ . _ _ _ _ _ _ _ — p p _ _ _ _ om 142 indicated (Re based on 8 and U0). Included in the figure are the DNS-data. The data are presented in wall-coordinates. All profiles show the anticipated logarithmic profile in the equilibrium region. The arrow shows where deviations from the log-law can be seen (821,200), indicating that the equilibrium region is indeed a small region of the flow domain. The calculated velocity profile for 5+=180 is in accordance with the DNS-data. The underprediction in comparison with Laufer’s data has been attributed to the small value for the Karman constant of K = 0.35 (see Chapter 1). A slight overprediction can be seen for 82395. The region outside the log-layer is the core region (see Figure 1.2). This region is characterized by the velocity defect law which resembles the similarity law for the mean velocity profile at and near the centerline (see Chapter 1). Figure 6.7 shows the velocity defect profile for flows at different Reynolds numbers (Hussain and Reynolds, 1975) and computed profiles at two different Reynolds numbers. The high Reynolds number DNS-data are also included. The computed velocity profiles show good agreement with the experimental data by Hussain and Reynolds (1975) as presented. The shear stress profile (normalized by 2k) is given in Figure 6.8. Experimental data by Laufer (1951), direct numerical simulation as well as computations by Nisizima and Yoshizawa (1987) are included. The predictions done with the relaxation/retardation model show a slight Reynolds number dependence at the boundary in the equilibrium region which diminishes once the Reynolds number is large enough. The fact that Reynolds number influences only enter through the boundary condition at ncq causes all curves to collapse for n—>1. This behavior is moreover expected since for small values of 118 all profiles satisfy Eq. (3.31). Thus, the influence of the Reynolds number is only given for predictions with 10w Reynolds numbers. The largest deviations (up to 50%) of the predictions occur near the equilibrium region if compared with the DNS-data. Comparable results are obtained with respect to the experimental data. 143 200.5 Semen 5328/ :82 ”N10 oesmmm _0.0 _00.0 —~ou« q — Amomnkv Eg «30075 AO0N._n+mv hoof. cocaxflom I II Ammmubv boo; 1r cosmos—om .III 00m.mmnom 00N.mmnom 0006 _ Mom - D O O 3%: 5E E 83 325652 .33: fioium Ba Emmmsm 3 San .mucoEtomxm —PPP-p _ p l 1N 100 1N0 13 100 1m: 10m NN 5-40: + + 144 358:0 05 mo newom $30 2: $88 coca—Ema mmubm Seam ”We oSmE w} o_ oo mo no we We we mo me _o oc Jul—__—__—__——d____u_—q__——uq___.__d—__q_—dq_~_qfifiq 8o 8 co 09 O 88.7%.3233355 Ba 25522 4 .. m- I acqfibv 382 8:333 -II . %1 mod , $3me 382 852$ I | .. w - - swims 382 832% I! N Am 1 I ASN._n+m._8:§=3 a .. % 1 r Gomwmdwmcei o . - .. . 1 3o 8w_n+w .333 ..a 6 EM 2d r_L_—__—__—_—_FL___—p_____p_~—____—~.__—_-_—_~_P ON.O 145 Turbulent Shear Parameter Figure 6.9 shows the computed profile for the turbulent shear parameter 1:8. The dependence on the local turbulent Reynolds number does not affect the calculation since it only appears in the near wall function cw which has been omitted from this calculation in the outer region. Thus, with cR, CB and on specified (see Chapter 5), it is the spatial distribution of the turbulent shear parameter which uniquely determines the spatial distribution of the Reynolds stress and the prestress. Therefore it is crucial to examine the behavior of this shear parameter in order to draw upon the behavior of the Reynolds stress. In the equilibrium region the shear parameter assumes a finite value of ItS=3.05 (for Re large). as is determined from the asymptotic solution of the equilibrium region for which the dissipation rate 8+ equals the velocity gradient 8*. Thus, the value for as equals k+ and can be determined from the algebraic structure of the normalized shear component (see Chapter 5). The value for Its increases slightly to a value of ItS=3.10 which prevails to about n=0.4 and decreases monotonically thereafter towards zero at the centerline. The almost constant behavior of this parameter over a wide region of the channel is a phenomenon which has also been observed in the DNS-calculations by Kim et al. (1987), Kim (1989) and the channel flow computation by Demuren and Sarkar (1993). Since the predictions of the mean velocity profile compare well with experimental (Laufer, 1951) and numerical data (Kim et al., 1987; Kim, 1989) it is presumed that the dynamics of the k and e-equation determine the behavior of Its in this region. The experimental data by Laufer (1954) appear to agree in magnitude near the equilibrium region. Budgets 0f Kinetic Energy and Dissipation Rate In order to understand the behavior of the shear parameter it is necessary to examine the budgets for the transport equation as given by Eqs. (6.7) and (6.8). Their respective 146 CBoESmm Seam 00235.5 23% 23m 083. 22m :moEbcoEnSH 05 mo cogngma :3an no.0 8305 0.0 0.0 md Yo md Nd ..o o.o _d___.._«—__dq___fiu——~__—d4——__—__—qfi4—o 35c 2E «me .5230 swank 33 .50: wank - Ahwmfi .._0 8 EU: OwHHer o l k mo .82?» got? 08 1 2003028 003808 8:: 23. 0 Banana? - o mmmltb r . I‘k‘ll‘ m o o owN.N\OON._H m n O O + . L 0 o 0 o l 0 o o o o 0 0 o o 000 0 o o o o o o o 000 00 . l 0 O 000 O . O 0 DD 000 0 00 00000 0 I] V 0 0000000 0 0 00 000 l 0 000 o 0 l o C o I. __b_—__p——L_—p—r_—__—_P_—_—_——____—D___m 8T+w .08 l N we 147 budgets are given in Figures 6.10a and 6.10b for 521,200. The dissipation terms in both figures is accounted for as a loss and therefore presented as negative values. The first pair of the small graphs inside the figures show the relative magnitude of the individual terms. The balance equation for k states that the production of (turbulent) kinetic energy is balanced by the net transport of this energy and its dissipation. The small graph for the k- budget shows that the balance of production and dissipation is given within ilO % for a very large region (llcq.»».. N.O cod r +0 83 88 com. 82 com o 1 1121..................._._. a .1 ..........................o...o.&8xzm.¢ 12:- .1 n oil!!! . ,1 I-.1¢/ 1. cc. 1 m / m mod r 1] // I1. wd- m / m 3 00002005 I 1| / .11. 1 x 4 1 V H w; :50 T L.. _ .. _»pr._.ti_LF11 0.0- a 1 . . com 711% In - . p p _ . p L p — _ _ _ _ — . _ _ _ O~.O 149 8nd 023935 05 08 835m 53 .0 8:me . w} . . . . . o _ w o c o v o N o o o _ _ q 4 T d _ _ a _ _ _ _ _ _ _ . _ V8.0: 1 a H 1 3 + 3 ed to 3 ed 1 1 1.1.21.1...3112o... .1 I . I moo o1 1 1 .EQfiEm 1 w o- 1 H 1 1 co- 853329 H I . I mood- H - - a- H H 1 352:0 1 3- H I . I 50.0. .I p b r 1p — > _ F— > > CO I. ............. lg 80.0 I +@ L 1 83 88 com. 82 Sr. o W h .... _ I. 1.. ..__- 02255 H I m. ........................................ [W O._1 I“ —8 o 1 1 1. ad- H m1 1 T IIIIII 111/ 1. we- 1“ . 1 n // n 0002605 1 moo o H m // 1” he- .1 fl / . 1 1 .1 1. ed- 1 . I h V a? I Moo 0 .l 1 p _ p (F _ r F _ _ . p » W131 a U 1 CON _H+w .III 1 T _ p _ _ _ _ _ _ _ — _ p _ — _ _ _ _ l wood .0.qu END 150 30250800 mmobm mEoEAom RES Z 05 @503 0009mm 380m 2:. ”2d Rama 3 o._ ad wd do 0.0 We to Go Nd _ — _ _ _ Adda. 7+3 m3a§§o> Us“ «EEEZ o 8801.8 essay. 1 ___ ________—____________—_____—1__u__. sfiafle .252 ceaéam .11 1 sawmubv 5:3 ._ 880118 382 8:.333 .l 1 88.718 4.33 n_e ___— _______r___—__ ~__.__________—___p———__ 0 0 O O Q I I O m 4 S<</2k remains fairly constant. Since this is not observed experimentally, a rigorous application of Eq. (3.49) for the entire channel region may prove inadequate. Friction Factor The friction factor f can be estimated using Eq. (1.29). Due to the fact that computations are only available in the outer region, exact (computational) estimates are not given. The following analysis shows, however, how (upper and lower) bounds are established. Neglecting the viscosity ratio in Eq. (1.29) an exact expression for the friction factor in laminar flow is obtained. Extrapolating this result into the turbulent region (i.e. for higher Reynolds numbers) this expression serves as a true lower bound for f (or upper bound for (2/00‘5). A closer estimate of the lower end magnitude for f, though not true upper and lower bounds, can be obtained using the following approach. The bulk average velocity ub+=ub/u* is given by | u; = Ju+dn, (6.20) O and can be decomposed into 2 parts, a contribution near the wall given by the equilibrium value for which u+ remains constant (see Eq.(6.11)) and the contribution in the outer region (i.e. neanSI) to u+‘n,1q.u+d l{+ tfifl "} 11- cq_n+_1 ueq_+81 dn cm, (6.21) where F represents the eddy viscosity ratio. Eq. (6.21) employs the fact that molecular influences are unimportant for neq_0m 2: :0 30:85 2: :0 8:255 ”2 0 Eswi w: 0.0 M3 5.0 0.0 m0 v.0 m0 m0 _ .0 od IIIIIIIIIIIIIIIIIIIIIIIIIITIIIIIIIIIIII .11 dq—fiJdfid—__d——_—__——_—__-—a_—_u___——_____Jq—__—u xmaocm 00 :0::80300 _P__——___—r_pL_~pb_—.___———________~__»L-_pb_—__~ lLIlLlLllllllllllllllllllllLlIlllllllll cod mod 2.0 2.0 ONO mmd Omd mmd 0:6 .3 164 38% mw_0:>0m 05 00 35:00:80 00:3:0mm 05 0:: 0238085 05 :0 30:00E 05 m0 00:0:EE ”E .o 0::wE e: o._ 0.0 m0 5.: 0.0 m0 v.0 m0 No _.o 0.0 Jqq___qfifi——_____fi_—___.dl—qqq4__.—fi——q._—d_.dfi-_d —.o Nd md v.0 m6 0.0 N: .xx llllllilJlllllllllllJJll‘lllL llllllIIITTITT‘IIIIIITIITIIIII ____F—~_P__—p_—_F___—___—__—_———_p—-_—___.__—p___ NHO 165 0.— im- 30:09:00 30.5 0:: :0 $0505 05 00 00:0:EE H S 0 0::mE 3 ad wd Nd 0d md vd Md Nd _ d dd — d u — d _ — — q u n J - u 4 — d — Ono u d a — _ _ — — W l— — u u u — _ q — — I— d u q — _ — d — — _ _ > > I a. .3 .00 2 05288 mm- 25 NE- a .. In om: l —.O 1 / 1 / r / 1 / / - / x I Ink - - / N m- > a > L / .m- i . 1| // N» Em: l N O / m- I / : 1 / / I , ./ / f / rt 9 < \\ 1 x % [fillliriflllT ||||| I. - ”may- - am- I -ffiill < luk\\ - — p — p — — — p — — P — — — — — — — — — — — p b — b p p p — h — — _ _ — b p — — — — — — — — _ p — md 166 contribution). In general there are two distinct effects contributing to -Ryz from the isotropic prestress contribution. Even though the first term appearing on the RHS of Eq. (5.19) is formally the same as for the IPS-theory developed in Chapter 2, the additional terms included in Ryy (see Eq. (5.18)) impose additional effects on this term. These additional effects can be attributed to the parameter CB and thus can be consolidated with Hyz as Table 6.2 indicates. However, the figure indicates the decomposition into both parts appearing explicitly in Eq. (5.19) as well as a decomposition according to Table 6.2. It can be seen that the isotropic prestress contribution overpredicts the shear stress by a factor of almost two. It is interesting to note that the first term on the RHS of Eq. (5.19) which contains both the isotropic prestress contribution as well as the additional effects according to Eq. (5.18) does already diminish the magnitude and thus serves as a corrective contribution. The (positive) parameter c5 therefore acts as a reducing agent, a ‘pre-viscosity’ for the Reynolds stress, but as a viscosity-type coefficient for the anisotropic part of the prestress. Summary The energy transfer among the components of the prestress is entirely controlled by the three parameters CM, c5 and cm. A linear closure hypothesis (i.e. no retardation) does result in a zero secondary stress difference for the prestress. The contributions of the linear term (i.e. cB-term) provides an energy transfer from Hu and Hyy to H22 in equal parts. The nonlinear interaction of the relaxation term moderates the rate of this transfer. Retardation provides a means for nonzero secondary normal stress differences through its energy-retaining effect on the HXX-component. Relaxation influences are more pronounced for, the shear component Hyz than for the normal components of g. This moderation effect on Hyz through CM counteracts the linear increase through c5. 167 The interactions of the prestress with the Reynolds stress through the implicit algebraic preclosure adds another level of complexity to the assessment of the energy transfer. Retardation effects on Hxx are directly passed on to R“, and constitutes the mechanism by which a non-zero secondary normal stress difference Rxx-Ryy is created. Both isotropic prestress contributions as well as viscosity effects introduced through c5 do not render Rxx-RyyatO. The energy retaining mechanism of the retardation increases the energy level for H“. For very large values of the shear parameter, the effect of retardation is not only to distribute energy by retaining it but also to prohibit unbounded growth of Rxx-Ryy. The magnitude of the shear component is lowered due to the presence of the prestress. A nontrivial part of this adjustment is caused by the influence of the parameter CB on the contributions from the isotropic prestress formulation. CHAPTER 7 CONCLUSIONS AND RECOMMENDATIONS 7.1 Conclusions A new representation for the Reynolds stress has been developed based on a Green’s function technique which reduces the non-local structure of turbulent two-point correlations to a local structure through the application of a spatial smoothing approximation. This smoothing approximation makes use of the fact that the decay of turbulent two-point correlations relax faster than the associated Green’s function. The validity of this assumption has been shown to hold in regions of large turbulent Reynolds numbers (Chapter 2). The introduction of the phenomenological time scale 1R (introduced through the presumed existence of a memory function common to all correlations appearing explicitly in the formulation) ‘allows’ the extension of this structure into the near wall region for which the local turbulent Reynolds number is small and ultimately vanishes. The algebraic structure obtained through this preclosure relates the Reynolds stress to IR, the mean velocity gradient and a correlation which was subject to further modeling. This correlation can be considered as a prestress to the Reynolds stress. A first assessment of this new model is invoked by closing the theory with an isotropic representation of the prestress. The unrestricted realizability of this closure for all values of as can be guaranteed a priori and therefore constitutes its most valuable property. The symmetry of the prestress and the way the Operator A acts on the prestress 168 169 (see Eqs. (2.15) and (2.16)) does not change the non-negative character of the Reynolds stress. Thus, non-negative eigenvalues of the prestress translate into non-negative eigenvalues for the Reynolds stress. For simple shear flows, the algebraic structure of this closure predicts a primary normal stress difference (Eqs. (2.31) and (2.34)) through the coupling with the operator A. Vanishing velocity gradients, as observed on the symmetry line of channel flow, cause the prediction of an isotropic turbulent state. An important key feature of this closure lies in the fact that the operator A, through the explicit appearance of the velocity gradient, preserves the frame-dependent character of the Reynolds stress. The time scale In is used to bridge the gap between the high (turbulent) Reynolds number region in which the smoothing approximation was applicable and the near wall region for which viscous effects become dominant. Its form is specifically chosen to guarantee a correct asymptotic behavior for the Reynolds stress component Ryy and the shear stress -Ryz. For wall- bounded simple shear flows, such as channel flow, a direct consequence of this closure and the representation of the time scale IR is the transfer of the entire energy from the normal components, i.e. Rxx and R ,, into the streamwise Reynolds stress component Rzz as the wall is approached. Even with this one-component energy state at the wall, a condition which is attributed to the intrinsic incapability of predicting secondary normal stress differences (see Eq. (2.33)) and which is not observed experimentally, this closure bears its own potential as it constitutes a significant improvement over the Boussinesq approximation commonly employed (Eq. (1.2)). Two distinct regimes are identified, a gradient type transport region for small values of 1,8, and an equilibrium type regime for large “as, consistent with the analysis of the equilibrium region (Section 5.5). The magnitude of the shear stress exceeds those observed experimentally, with a maximal value of -Ryz=0.289 for 1:,S=1.73. Through this change of transport regimes, the eddy viscosity coefficient assumes a shear thinning 170 behavior. This behavior arises naturally through the theory whereas various turbulence models introduce artificial damping to the eddy viscosity. Nontrivia] representations for the anisotropic part of the prestress H are implemented using a phenomenological modeling approach with frame-invariant derivatives which include both relaxation and retardation terms. The observation of secondary flow patterns in non-circular ducts for laminar flow of viscoelastic fluids as well as for turbulent flows of Newtonian fluids, both of which are caused by normal stress differences, lead to the idea of using the mathematical structure of constitutive equations as they are used in the field of viscoelasticity. Terms to describe relaxation effects became necessary to describe flow fields where explicit memory effects of the Reynolds stress are present, such as retum-to-isotropy. The associated parameter cl] obtained from the calibration of the turbulence model against this flow field was done for a positive third invariant of the anisotropy tensor (i.e. IIIB>O) and estimated to 0.644. Parameter estimates for which IIIBO. It is therefore concluded that for flow fields with IIIB<0 and IIIB>O a more general formalism may be necessary to adequately represent cl]. It was found that retardation is necessary to reproduce two effects, both of which are directly related to the second normal stress difference: a) retardation causes the retention of energy in the spanwise component of the Reynolds stress Rxx for moderate values of Its (i.e. O S 1.33 S 3) and is thus directly responsible for the existence of a second normal stress difference within the channel domain (see Figures 6.9 and 6.11). b) For ItS —-) oo (i.e. as the solid wall is approached), the unbounded growth of the secondary normal stress difference is prohibited and a unique two-component turbulent energy state is attained (Section 5.3). The latter of these effects, i.e. the partition of energy among two Reynolds stress components, is paralleled by similar observations in simple shear flows 171 with high shear but without the presence of walls and lead to the conclusion that the unique existence of this asymptotic state seems justified (Lee et al., 1990). The extension of the isotropic prestress theory (IPS) with an anisotropic prestress (APS) representation added a level of complexity which made it impossible to a priori guarantee realizability. However, for channel flow as a simple shear flow, and with the parameter estimates made (Chapter 5), realizability for all values of 1,8 was shown a posteriori (Figure 5.4). The choice of the parameters a=-l and b=-3, describing the particular form of the convected derivatives, was done for practical purposes and enabled an independent analysis of influences of the scaling parameter CR, as introduced through the time scale TR, on the behavior of the energy partition and the turbulence states. From paths in the energy distribution plane (Figure 6.1), it may be concluded that a variation of cR with local turbulent parameters, such as the turbulent kinetic energy and the dissipation rate, yields a better representation of energy partition among the normal components of the Reynolds stress. For example, larger values for cR might better represent conditions at or near the channel center (i.e. for small 135), whereas smaller values seem to better represent the equilibrium region (i.e. “CIS z 3.5). A mapping into the invariant plane (Figure 6.2) indicate similar trends and moreover show all case studies to be entirely realizable. The application of the Green’s function technique used to develop the algebraic coupling of the Reynolds stress to the prestress was also employed to develop the associated transport equations for the turbulent kinetic energy and the dissipation rate. This approach rendered a new expression for the representation of the turbulent transport of a scalar quantity (Eq. (4.12)). The additional term appearing explicitly in these expressions and the fact that the time scale IR is an integral part indicate a range of validity which is surmised to have more universality as common approaches (e.g. Hanjalic and Launder (1972b)). Common approaches have been followed inasmuch as 172 the pressure diffusion term in the kinetic energy equation was incorporated into the resulting gradient type representation (Eq. (4.16) through an adjustable transport parameter, the equivalent term in the dissipation rate equation was omitted. The numerical procedure chosen to solve the nonlinear two-point boundary value problem proved to be adequate and robust for the calculation of the outer region of fully developed channel flow. The solution of the associated transport equations was done sequentially using a standard tri—diagonal matrix solver. The sequential solution of the transport equations in combination with initial profiles for the velocity, kinetic energy and dissipation rate indicate an iterative nature of this procedure. A solution is attained by calculating a developing flow towards a fully developed state, implemented through convective terms for the streamwise direction. Computational time was negligible. 7.2 Recommendations The recommendations presented here reflect two aspects. On the one hand, they stem from the findings obtained through this research and are therefore a direct consequence of them. The second aspect aims at the outgrowth of the theory set forth at the onset of this research. This outgrowth comprises extended research from fundamental aspects as well as the application of the new theory to different flow fields. The latter aspect of the recommendations therefore consists of a more qualitative nature. Preclosure and Closure The spatial smoothing approximation is based on the assumption that the Green’s function associated with the operator given by Eq. (2.2) peaks spatially for time intervals 173 long enough for the turbulence structure to relax. The local character of the two-point correlation thus obtained reflects, however, a high Reynolds number approach not through the operator .8 but through the space-time structure of the turbulence. For wall- bounded flows where viscosity becomes dominant it is the spatial and temporal decay of the Green’s function which picks up non-local effects from the two-point correlation of the terms given by Eq. (2.4). The approach taken in this research employs the phenomenological extension of the time scale ‘L'R to incorporate low Reynolds number regions. A more direct approach is to incorporate effects of the wall on the behavior of the Green’s function, thus employing some sort of Greens function for a semi-infinite domain. This can be done using the method of images (Morse and Feshbach, 1953). Along with a modified Green’s function, an adequate representation of the two-point space-time correlation can be incorporated as solid walls are approached. Two issues need to be addressed with respect to this suggestion: a) It was shown that the spatial smoothing can hold up to y+=2.2 (Chapter 2). With a postulation that this approximation can be extended towards the wall, the research here used the empirical time scale IR to describe near wall phenomena. A better description of the autocorrelation of the associated fluctuating quantities might be a viable choice to avoid the introduction of the empirical time scale TR and yield a better description of the turbulent statistics. b) The implementation of the Green’s function for a semi-infinite domain does not only ‘allow’ values for which x = 3 to be of importance. Therefore, its use entails a complete description for the two-point space-time correlations of the associated turbulent quantities. Analytical expressions for the form of the two-point correlations have been given for homogeneous flow field (see Hinze, 1987). For inhomogeneous flow fields and for correlations involving higher order moments very little is available. The assumption of a 174 general exponential decay of the space-time correlation, a practice employed by various researchers, can therefore be understood as a first step towards a more detailed analysis. The applicability of a suitable representation for the space-time correlation of the prestress needs further research. The detailed investigation of the instantaneous fields obtained through the direct numerical simulation might serve as a guideline for the feasibility of this aspect. The assumptions that the gradient of the fluctuating Reynolds stress and the gradient of the fluctuating pressure can be consolidated to the single quantity 1” constitutes a facilitation inasmuch as the physical effects of the individual terms might be different. The description of the influences of these terms was shifted to the relaxation and retardation terms in the closure approximation. The analysis of the (instantaneous) DNS- data might serve as a guideline and direct the thinking of whether this consolidation is justified. The choice of the parameters a and b determining the particular form of the frame- invariant derivatives in the closure enabled an adequate representation of the problem at hand. From a more general point of view it might prove necessary to subject those parameters to a similar optimization routine as employed for the determination of the parameters cR and c5. Relaxation effects have been observed experimentally such that a mathematical representation capable of describing this effect is compelling and some sort of time derivative is justified (LHS of Eq. (3.5)). The use of retardation for turbulent channel flow was justified a posteriori as it could be related to the secondary normal stress difference. Whether or not the same level of description of secondary normal stress differences could be obtained using nonlinear strain rate terms rather than retardation effects needs to await further research and is thus a strong recommendation. 175 Transport Equations The modeling of the transport equation for the kinetic energy is considered adequate through the new formulation developed in this research. The direct application of this formalism to the pressure diffusion terms in both associated transport equations has not been pursued in this research. However, the open question remains whether this violates any basic principles. From a fundamental point of view it seems generally possible as long as, for example, the pressure Poisson equation is satisfied at all times. The general idea of doing so has been employed by other researchers (Harlow and Nakayama, 1967) who modeled this transport term as being proportional to the gradient of the mean pressure (Eq. (4.17)), a result which could be obtained in a similar way using the formalism used in this research. The modeling of the transport terms in the equation for the dissipation rate needs further evaluation. Modeling approaches of the omitted correlation explicitly appearing through the derivation here might seem a suitable way to proceed in order to obtain an improved representation for this term. The representations of the production and destruction terms in the e-equation proved an adequate description of channel flow and, moreover, provide a realizable representation for the decay of isotropic turbulence. However, recent direct numerical simulation data for this flow field indicate differences in the decay process depending on the initial energy spectrum supplied to the simulation process. Thus, a more thorough investigation of this term in particular might prove useful. Numerical Aspects The numerical method for the solution of the two-point boundary value problem employed in this research was adequate. The explicit calculation of the near wall region was, however, not possible. Preliminary calculations in this region showed that the 176 implicit appearance of the shear parameter (i.e. the inverse proportionality) in the description of the shear stress seem to be the reason for the failure of the method at hand. A further test was initiated to determine which of the equations to be solved pose problems with the implementation of the new theory as laid out for the near wall region. An existing turbulence model based on a Boussinesq representation was used as a basis. The substitution of the turbulent transport terms with the descriptions developed in this research did not show any problems related to the numerical solution. Moreover, it could be confirmed that the differences between the different implementations of the transport terms in general did not show any notable differences in the mean velocity profile. The further substitution of the eddy viscosity representation in the momentum equation did not result either in any problem. Only the substitution of the eddy viscosity into the production terms of both the k- and e-equation showed an unstable behavior of the numerical method with the result that no convergence was obtained. From this preliminary study it was concluded that the steep gradients of both the kinetic energy and the dissipation rate in the near wall region in combination with the earlier mentioned implicit appearance of the shear parameter contribute to the observed numerical instability. A recommendation which results from this observation is the search for more robust numerical methods to compute regions where the mentioned effects are important. APPENDICES 177 A. Direct Numerical Simulation Data During this research extensive use has been made of the direct numerical simulation data (DNS-data) which were acquired during computations of turbulent channel flow by Kim et a1. (1987) and Kim (1989). Those data have been provided by Dr. Kim from the Turbulence Research Center (TRC) at Stanford University. For purposes of reproduction of any modeling steps made during this research the available data are presented here. 1) Direct Numerical Simulation Data for 8+=180 (Kim et al., 1987) y+ < u1u; >+ < u;u; >+ < uQu; >1 — < u;u; >+ 8+ u+ 0.000 0.000 0.000 0.000 0.000 0.165 0.000 0.054 0.000 0.000 0.000 0.000 0.164 0.053 0.217 0.002 0.000 0.006 0.000 0.161 0.213 0.488 0.008 0.000 0.030 0.000 0.157 0.479 0.867 0.023 0.000 0.095 0.000 0.152 0.850 1.354 0.051 0.000 0.231 0.002 0.148 1.325 1.948 0.092 0.000 0.475 0.006 0.143 1.902 2.650 0.148 0.002 0.864 0.014 0.138 2.576 3.459 0.215 0.004 1.430 0.030 0.131 3.340 4.374 0.290 0.009 2.173 0.057 0.124 4.180 5.394 0.370 0.017 3.058 0.098 0.117 5.077 6.520 0.452 0.028 4.008 0.152 0.1 13 6.007 7.751 0.534 0.045 4.925 0.218 0.111 6.940 9.085 0.615 0.066 5.716 0.291 0.1 12 7.850 10.522 0.694 0.093 6.321 0.365 0.113 8.713 12.061 0.769 0.126 6.713 0.436 0.114 9.513 13.702 0.840 0.164 6.900 0.500 0.1 13 10.240 15.442 0.905 0.206 6.908 0.555 0.1 10 10.892 17.282 0.964 0.252 6.774 0.600 0.106 11.470 19.220 1.014 0.300 6.536 0.636 0.101 11.980 21.254 1.056 0.349 6.229 0.664 0.096 12.430 23.384 1.089 0.397 5.883 0.683 0.090 12.825 25.609 1.1 16 0.445 5.519 0.697 0.084 13.174 27.926 1.137 0.490 5.156 0.705 0.078 13.484 30.335 1.153 0.531 4.805 0.709 0.072 13.759 32.835 1.165 0.569 4.471 0.708 0.066 14.007 35.423 1.173 0.602 4.159 0.705 0.061 14.232 38.098 1.177 0.630 3.871 0.699 0.056 14.438 40.858 1.176 0.653 3.608 0.692 0.051 14.628 43.702 1.170 0.671 3.369 0.682 0.047 14.806 46.629 1.160 0.685 3.154 0.671 0.043 14.974 49.636 1.146 0.693 2.959 0.658 0.040 15.133 52.721 1.129 0.698 2.783 0.645 0.037 15.286 55.883 1.108 0.699 2.625 0.631 0.034 15.432 59.119 1.083 0.696 2.481 0.616 0.031 15.574 178 62.429 1.056 0.691 2.352 0.600 0.029 15.71 1 65.809 1.028 0.683 2.236 0.584 0.027 15.844 69.258 1.000 0.672 2.130 0.567 0.025 15.974 72.774 0.971 0.660 2.032 0.549 0.023 16.101 76.355 0.940 0.646 1.942 0.531 0.021 16.226 79.997 0.908 0.631 1.858 0.512 0.020 16.348 83.700 0.876 0.615 1.779 0.492 0.018 16.467 87.462 0.843 0.599 1.705 0.473 0.017 16.583 91.278 0.809 0.582 1.634 0.453 0.016 16.696 95.149 0.776 0.565 1.566 0.434 0.015 16.805 99.070 0.744 0.549 1.500 0.414 0.014 16.911 103.040 0.712 0.532 1.434 0.393 0.013 17.014 107.060 0.681 0.516 1.369 0.373 0.012 17.113 111.120 0.651 0.501 1.305 0.352 0.011 17.209 115.220 0.622 0.487 1.241 0.331 0.010 17.301 119.360 - 0.594 0.473 1.178 0.309 0.010 17.389 123.540 0.567 0.460 1.117 0.287 0.009 17.474 127.750 0.542 0.448 1.057 0.265 0.009 17.555 131.990 0.518 0.436 0.999 0.243 0.008 17.632 136.260 0.496 0.424 0.943 0.221 0.008 17.705 140.560 0.476 0.414 0.890 0.199 0.007 17.771 144.880 0.458 0.404 0.839 0.177 0.007 17.832 149.230 0.441 0.395 0.793 0.155 0.006 17.887 153.590 0.426 0.388 0.750 0.132 0.006 17.936 157.970 0.412 0.382 0.713 0.110 0.006 17.978 162.360 0.400 0.377 0.682 0.088 0.006 18.013 166.760 0.390 0.373 0.657 0.066 0.006 18.040 171.170 0.382 0.371 0.640 0.044 0.005 18.060 175.580 0.377 0.369 0.629 0.022 0.005 18.072 180.000 0.375 0.369 0.625 0.000 0.005 18.076 2) Direct Numerical Simulation Data for 8+=395 (Kim, 1989) y+ < u;u; >’ < u;u; >" < u;u’z >+ — < u;u; >+ 3+ u+ 0.000 0.000 0.000 0.000 0.000 0.221 0.000 0.053 0.000 0.000 0.000 0.000 0.219 0.053 0.211 0.003 0.000 0.007 0.000 0.213 0.210 0.476 0.013 0.000 0.035 0.000 0.206 0.473 0.846 0.037 0.000 0.1 1 1 0.000 0.198 0.841 1.321 0.080 0.000 0.270 0.002 0.190 1.312 1.902 0.146 0.000 0.553 0.007 0.181 1.886 2.588 0.233 0.003 1.005 0.018 0.172 2.557 3.379 0.337 0.007 1.651 0.038 0.161 3.317 4.275 0.452 0.014 2.489 0.071 0.149 4.152 5.275 0.572 0.025 3.465 0.120 0.139 5.041 6.380 0.692 0.042 4.487 0.184 0.132 5.958 7.588 0.809 0.065 5.450 0.259 0.130 6.873 8.901 0.921 0.094 6.264 0.340 0.130 7.760 10.317 1.029 0.131 6.874 0.422 0.130 8.596 179 11.835 1.132 0.175 7.267 0.499 0.129 9.366 13.457 1.228 0.225 7.458 0.567 0.126 10.065 15.180 1.315 0.280 7.481 0.627 0.122 10.689 17.005 1.394 0.340 7.371 0.676 0.117 11.244 18.932 1.462 0.401 7.166 0.717 0.110 11.733 20.959 1.520 0.464 6.896 0.749 0.104 12.165 23 .086 1.567 0.527 6.589 0.774 0.097 12.547 25.312 1.605 0.588 6.263 0.794 0.091 12.885 27.638 1.636 0.646 5.931 0.808 0.084 13.187 30.062 1.659 0.701 5.602 0.818 0.078 13.456 32.583 1.677 0.753 5.284 0.825 0.072 13.699 35.202 1.691 0.799 4.983 0.829 0.067 13.921 37.917 1.702 0.841 4.706 0.830 0.062 14.123 40.727 1.710 0.877 4.453 0.830 0.057 14.311 43.633 1.715 0.909 4.227 0.828 0.05 3 14.485 46.632 1.715 0.936 4.024 0.826 0.049 14.649 49.725 1.712 0.958 3.843 0.821 0.045 14.803 52.910 1.705 0.976 3.681 0.816 0.042 14.948 56.186 1.697 0.990 3.535 0.809 0.039 15.088 59.554 1.686 1.000 3.403 0.802 0.037 15.221 63.01 1 1.673 1.007 3.279 0.793 0.034 15.352 66.557 1.660 1.01 1 3.163 0.783 0.032 15.479 70.191 1.645 1.012 3.056 0.773 0.030 15.604 73.91 1 1.627 1.010 2.958 0.764 0.028 15.727 77.718 1.606 1.006 2.871 0.754 0.027 15.848 81.610 1.582 1.000 2.793 0.745 0.025 15.966 85.585 1.555 0.993 2.722 0.736 0.024 16.082 89.644 1.526 0.985 2.657 0.726 0.022 16.195 93.784 1.496 0.977 2.596 0.717 0.021 16.306 98.004 1.468 0.969 2.537 0.707 0.020 16.416 102.300 1.441 0.959 2.482 0.697 0.019 16.524 106.680 1.414 0.949 2.431 0.688 0.018 16.632 1 1 1.140 1.388 0.938 2.382 0.679 0.017 16.739 1 15.670 1.365 0.926 2.337 0.670 0.016 16.844 120.280 1.344 0.913 2.294 0.662 0.015 16.948 124.960 1.323 0.899 2.251 0.653 0.015 17.049 129.710 1.299 0.886 2.211 0.643 0.014 17.148 134.530 1.272 0.872 2.171 0.632 0.013 17.245 139.430 1.242 0.858 2.132 0.620 0.013 17.342 144.390 1.212 0.845 2.096 0.607 0.012 17.437 149.420 1.182 0.830 2.063 0.594 0.011 17.531 154.510 1.153 0.815 2.032 0.582 0.011 17.622 159.670 1.125 0.800 2.001 0.570 0.010 17.711 164.890 1.099 0.784 1.967 0.557 0.010 17.799 170.170 1.072 0.769 1.930 0.544 0.009 17.885 175.520 1.046 0.755 1.891 0.531 0.009 17.971 180.920 1.020 0.741 1.851 0.518 0.009 18.056 186.380 0.994 0.728 1.81 1 0.504 0.008 18.141 191.890 0.968 0.715 1.770 0.491 0.008 18.226 197.460 0.942 0.702 1.730 0.477 0.008 18.31 1 203.080 0.915 0.689 1.690 0.463 0.007 18.395 180 208.760 0.888 0.675 1.652 0.449 0.007 18.478 214.480 0.862 0.661 1.613 0.436 0.007 18.558 220.250 0.838 0.647 1.574 0.422 0.006 18.638 226.070 0.815 0.633 1.535 0.408 0.006 18.715 231.940 0.793 0.619 1.494 0.394 0.006 18.790 237.840 0.768 0.605 1.454 0.380 0.005 18.863 243.790 0.743 0.591 1.415 0.366 0.005 18.934 249.780 0.718 0.578 1.378 0.352 0.005 19.003 255.810 0.695 0.564 1.342 0.338 0.005 19.071 261.880 0.673 0.550 1.306 0.324 0.005 19.137 267.980 0.653 0.537 1.269 0.310 0.004 19.200 274.120 0.635 0.524 1.232 0.295 0.004 19.262 280.280 0.618 0.513 1.194 0.281 0.004 19.321 286.480 0.602 0.502 1.156 0.266 0.004 19.378 292.710 0.586 0.492 1.117 0.251 0.004 19.433 298.970 0.569 0.484 1.076 0.235 0.004 19.486 305.250 0.553 0.477 1.037 0.220 0.003 19.538 31 1.550 0.540 0.472 0.998 0.204 0.003 19.587 317.880 0.530 0.468 0.958 0.188 0.003 19.635 324.230 0.521 0.464 0.917 0.172 0.003 19.681 330.590 0.513 0.461 0.876 0.157 0.003 19.725 336.980 0.506 0.458 0.836 0.141 0.003 19.766 343.380 0.498 0.455 0.799 0.125 0.003 19.804 349.790 0.492 0.454 0.765 0.109 0.003 19.838 356.210 0.486 0.453 0.735 0.094 0.003 19.869 362.650 0.481 0.452 0.710 0.078 0.003 19.896 369.090 0.476 0.452 0.691 0.062 0.003 19.918 375.550 0.472 0.452 0.677 0.046 0.003 19.935 382.000 0.469 0.452 0.667 0.031 0.003 19.948 388.460 0.467 0.452 0.662 0.015 0.003 19.956 394.920 0.466 0.452 0.660 0.000 0.003 19.959 181 B. Friction Factor Analysis The equation for the mechanical energy - derived from the Navier-Stokes equation - and the equation for the turbulent kinetic energy form the basis for the derivation of the relation given by Eq. (1.30) which relates the friction factor to the mean velocity gradient and the dissipation rate for channel flow. The mechanical energy balance is obtained by multiplying the equation of motion (Eq. (1.7)) with the mean velocity . Minor rearrangement yields: 1a< > d2 d — p =v 2’ ——. (B.1) p 02 dy dy ’ i E E The individual terms (i.e. I, II and [11) can be integrated over half the channel domain (i.e. 0 10 1=jp az dyzg (3:102 >dy ] 1:w 2 = -E(F)(ub0) == —ubu.. (B2) 5 d2 8 d d d H=Jv—Tz—dy=IV{—( z )—( )2}d 0 dy o dy dy d55d 5d =v< . >—( ‘ )Zdy}=—v( 2 >20. (13.3) dy 0 i dy 1 dy L V0 8 8 d d d 111: ——dy= —{-——(< uz>< u’u;>)———--—z }dy 1 dy ’ 1 dy ’ ’ dy d dy 5 d< > dy =J-—-ul—dy. (3.4) 0 dy =v —()| +£ Icy 182 Upon arranging the individual contributions: d dy 5 d 5 —ubu3 2 —VJ( z )2dy+J< u’yu; > 0 dy (B5) 0 dy is obtained. In the second part, an expression is obtained relating the integral of the production of kinetic energy to the integral of the dissipation. The transport equation for the turbulent kinetic energy as given by Eq. (C3) is evaluated for statistically stationary, fully developed channel flow to 1d dzk , , d d ,u'-y_' = — y z ——dy+V'([':—y—dy— Zfi-Zdy , , dd u'u' —J yif——dy, (B.7) 0 with the result I ’ 25 dk 25 28 0=—— +v-—- —Iedy P y 0 dyo 0 g T 25 1 1 28 d< > — I u, dy—< u’ E—E-> (B8) Y Z d y 0 y 0 Thus, an expression relating the integral of the production to the integral of the dissipation has been developed to 28 28 d 0=—J8dy—J dy. (B9) 0 o dy 183 With the fact that channel flow is symmetric across the center it readily follows that d dy 5 5 0=~jedy—j dy. (13.10) 0 0 Substituting Eq. (B.10) into Eq. (B.4) one obtains d dy 5 5 —ubu3=—vj( )Zdy—Jedy, (13.11) 0 0 which -with the normalization given by Eq. (6.5) and the definition in Eq. (1.28)- yields 2 .‘ du*. . \Eza {[(df) +8 116, (3.12) 184 C. Derivation of the Exact Turbulent Kinetic Energy Equation The derivation of the exact transport equation for the turbulent kinetic energy k stems from the transport equation for the fluctuating velocity which is given by ' 1 a—Et+ < u > ~Vu': —3Vp'+szu'—u'oV < u > —u'-Vu’+ < u'-Vu'>. (CD I H IH IV V VI VII Through a scalar multiplication of Eq. (CD with u‘ and subsequent ensemble-averaging the exact transport equation for k is obtained. The derivation is devided into 2 parts, one of which is the scalar multiplication to obtain the instantaneous value for the kinetic energy of the fluctuating field. This intermediate equation is used in the derivation of expressions for the transport terms. For purposes of clarity the symbol denotes the mean kinetic energy of the fluctuating velocity field whereas k represents the instantaneous value. This convention does only apply within this appendix. The subsequent ensemble average of the instantaneous equation finally renders the transport equation for the mean turbulent kinetic energy. The terms in Eq. (C.1) are treated separately. Constant properties (i.e. density and viscosity) are assumed. 0.32.; ziafl'! .. 3_k ‘ at 2 81 —0t H: u'-(< u > -Vu_')=< 11> -(u'-(Vu')T)=< u > ~V—E—é—E— =< u > -Vk 1H: —lg'-Vp'=—lV-(u'p')+lp'V‘u'=-1V°(u'p') p p 9 :0 9 IV: VE'.VZE': VE'°(V ' V2) = V{V ’ (VE.'E. ) _ VE':(VE' )T} = V{V%V(u'u')— Vu':(V2’)T} 185 = V{V2k — Vu’:(Vu’ )T} V: — u'-(u’-V < u >) = —(u"-V < u >)- '= u E'ZV < u >2 —u'u':V < 11 > V1: - 2'12”!) >= —(2"Vu' ) - 2'= Q'E'IV! VH: u'-(< u'-Vu'>) = u" < u’-Vu'> The transport equation for k (instantaneous) is thus given by 3512+ .Vk z—lV -(p'E')+(VV2k —€)-_1_1'_l}_'2V < E> t P I H HI IV V —V-(u'-E—;;)+u_'.. (C.2) VI VH The designation with roman numbers used here refers to the individual term in the equation for the fluctuating velocity (i.e. Eq. (C.1)) such that their origin is apparent. This designation differs from the one used in Chapter 4. Upon ensemble-averaging of Eq. (C2) the transport equation for the mean kinetic energy of the fluctuating velocity field is obtained. a+-V=-lV-+(VV2 —<£>) I H HI IV —2V—V-. (C.3) V VI It is noteworthy that the Reynolds stress itself - appearing in Eq. (C.1) as term VH - contributes to Eq. (C.2) but not to the mean kinetic energy equation. The terms in the equation for are designated as follows 186 1: Time Derivative IV: Diffusion and Dissipation H: Convection V: Production HI: Pressure Transport by u' VI: Energy Transport by u' The terms HI and VI are combined to the turbulent diffusion since the pressure fluctuations and the instantaneous turbulent kinetic energy are simultaneously transported by velocity fluctuations. Their explicit form need to be modeled in order to obtain a closed form of the transport equation for the turbulent kinetic energy. After applying the same formalism as used in Chapter 2 for the fluctuating velocity field (i.e. Eqs. (2.1)-(2.3)) to the equation for the fluctuating kinetic energy (Eq. (C2)), the following form for k' is obtained my) 2 —j jog, 13,1){629 <12> +fig}dvdE, (C.4) EV where f; is given by f; : lV-{p'11_'— < p'u'>}+(e— < e >)+ {gy- < u'g'>}:V < u> P + {u'-Vk'— < u'-Vk'>}— {u’-(V- < u'u’>)— < u'-(V~ < u'u'>) >}. (C5) A formal representation for the average transport of k' due to velocity fluctuations can therefore be obtained to < u'k'>= —j jG(x,tlg,1){< 30> V < 12> + < u'f; >}dVd1. (C.6) Ev With the same reasoning as applied before (see Chapter 2) that the turbulence structure decays faster than the Green’s function such that only the autocorrelation dominates the RHS of Eq. (C6) and that the temporal structure can effectively be consolidated into the time scale TR (i.e. it is assumed the same common memory function (1) exists) the transport term can be written as 187 < u'k'>= —‘CR .V < k > —‘CR < u'fk' >. (C7) The second term in brackets on the RHS of Eq. (C.6) can analogously be developed using the equation for the fluctuating velocity (i.e. Eq. (2.3)) as a basis to : —TR ~V—1R , (C.8) where f‘ represents the consolidation of the divergence of the fluctuating Reynolds stress and the gradient of the fluctuating pressure. From Eqs. (C7) and (C8) the following form for the transport of the fluctuating kinetic energy due to velocity fluctuations can be obtained. = —‘ER -V7--'cfz -A. (C9) The operator A is defined by Eq,. (2.16). 188 D. Derivation of the Exact Turbulent Dissipation Rate Equation The derivation of the exact transport equation for the dissipation rate 8 stems from the transport equation for the fluctuating velocity which is given by 1 I a—E+ < 11 > Vu': —V£ + VVzu—Q'V < u > —u_'-Vu_'+ < u'-Vu_'>. (D.1) 1- P I H [H IV V VI VH To obtain the transport equation for e the following steps are performed. First, the gradient operator V is applied to Eq. (D.1). Second, a scalar multiplication with (Vti')T is performed. After multiplication with the kinematic viscosity the scalar equation for the instantaneous dissipation rate is obtained. This equation is used in the derivation of the transport term. The subsequent ensemble averaging yields the transport equation for the mean value of the dissipation rate. For purposes of clarity the symbol <8> denotes the mean dissipation rate of the fluctuating velocity field whereas 8 represents the instantaneous value. This convention does only apply within this appendix. The terms in Eq. (D1) are treated separately. Constant properties (i.e. density and viscosity) are assumed. 0(Vu') 1 B(Vu')T:Vu' 182 1; V 'T:——-— =— - -=—— W E) 01 2V 8: 201 II: V(Vu')T:V(< u > Vu') = V(V_u')T:(V < u > -V_11'+ < u > ~VVu') = V(Vu')T:(V < u > ~Vu') + V(Vu')T:(< u_ > -VVu‘) = V(Vu')T:{(Vu')T -(V < .11 >)T}T Jrév < u > -V((Vu')TIVu') = V(Vu’):{(Vu')T -(V < u >)T}+%< u > Vs 189 = v{(Vu')-(Vu’)T}:(V f +—;—< 11> -Ve . T 1 . V . r . V . . 111: v(Vu) :V(——Vp)=-—(Vu) :V(Vp)=——V(Vp.):)72 p p p ' =_%{V-(Vp'-Vu')—Vp'°(V-(V.11')T)}=-%V-(VP"VE') x 1v; v(Vy_')T:V(7V22') = VZV ~ {(Vu')T:V(Vu')T}" V’V(Vu')TEV(Vu')T . 2%V2V2((Vy_')TI(VQ'))-V2V(VE.)TEV(VE')T =%VV28— V2V(Vu')TEV(Vu')T V: v(Vu’)T:V(—u'-V < 2 >) = —v(Vu')T:(Vu'-V < u > +g'-VV < u >) = —v(Vu')T:(V2"V )—V(Vu’)T:(u'-VV <2>) = —v{(Vu')T -Vu'}:V < E > —v(u'(Vu')T)EV(V < u >)T VI: V(V1_1_' )T:V(—u’-V_u') -.: —V(Vu_’ )T:(Vu_'-V_u'+u'-VVE') = —v(Vu' )T:(Vy_'-Vg' ) — v(Vu')T:(g'-VV9') o 1 _ r 1 1 1 . 1 : —v{(Vu )T .170 }.Vu—Ev_u-V{(Vg )T.Vu} VII: V(V_l_l_')T2V(< u'-Vu'>)=v(V_u')T:V The complete transport equation for the instantaneous dissipation rate is obtained after multiplication by a factor of 2 to EJF {2VVQ'.(V_11')T:(V < u >)T+ < 11 > 'V8} 01 1 11a Hb 190 = -2XV . (Vp'-V_u')+{VV28 - 2V2V(Vu')TEV(VE')T} p IH IVa 1V6 — 2v((V.11')T -V1_1.’):V < u > -ZV(2'(V2')T )3V(V < 11 >)T V, VI, - 2v{ (V1.13 )T ‘Vu' } :V2'-vu'-V{(Vu' )T:V2'} VIa VIb +2v(Vu_')T:V. (D.2) VII The designation with roman numbers used here refers to the individual term in the equation for the fluctuating velocity (i.e. Eq. (D.1)) such that their origin is apparent. The designation in Chapter 4 is done with respect to the individual terms appearing in the 8- equation and differs therefore from the one used here. Through ensemble-averaging the equation for the mean dissipation rate is obtained. 8 + {2v < Vg-(Vu'f >:(V < E >)T+ < E > -V < e >} 1 Ha Hb = —2XV- < Vp'-Vu'> +{VV28 - 2V2 < V(VE' )TEV(VH' )T >} 9 HI IVa IVb — 2v < (Vu')T 'Vu'>:V < u > —2v < u'(Vu')T >EV(V < E >)T Va Vb — 2v < {(Vu’ )T -Vu' } :Vu'> —v < u'-V{(Vu’ )T:Vu'} >. (D.3) V 1a VIb 191 The Reynolds stress itself - appearing in Eq. (D. 1) as term VH - does not contribute to the transport equation for the dissipation rate . Through the application of the gradient Operator on Eq. (D.1) additional terms arose in the equation for which are designated as follows (Rodi and Mansour, 1992): 1: Time Derivative IVb: Dissipation Ha: Convection Va: Mixed Production Hb: Production by Mean Velocity Gradient Vb: Gradient Production 1H: Pressure Transport V13: Turbulent Production , 1V3: Viscous Diffusion VIb: Turbulent Transport The terms Hb, Va, Vb and VIa are combined to the turbulent production and terms [H and VIb are combined to the turbulent diffusion Both turbulent production and turbulent diffusion as well as dissipation need to be modeled in order to obtain a closed form of the transport equation for the dissipation rate. After applying the same formalism as used in Chapter 2 for the fluctuating velocity field (i.e. Eqs. (2.1)-(2.3)) to the equation for the fluctuating dissipation rate (Eq. (D2)) the following form for 8' is obtained 8'(x, 1) = —j 1 (3(5, 113,1){0'17 <12> +i,’}dvc1f, (0.4) W where fe’ is given by f,’ : 2vV . {VR- Vg— < Vg- V05} + 2(8— < e >):V < g > +2{82Vu'— < e:Vu'>} P P = = " = = + 2{_§__:Vy_'— < g:Vu'>} + 2v2{V(Vu')TEV(Vg' )T— < V(Vu')TEV(Vu')T >} (D5) + 2v{(_11'(Vu')T )EV(V < 11 >)T— < g'(V_1_1_')T >EV(V < E >)T}+ V - {3'8— < u'e'>}. A formal representation for the average transport of 8' due to velocity fluctuations can therefore be obtained to 192 < u'e‘>=—j jG(x,tlg,1){ V < 8 > + <3"; >}d\?di . (D6) \7 H) With the same reasoning as applied before (see Chapter 2) that the turbulence structure decays faster than the Green’s function such that only the autocorrelation dominates the RHS of Eq. (D6) and that the temporal structure can effectively be consolidated into the time scale TR (i.e. it is assumed the same common memory function (1) exists), the transport term can be written as =—1:R -V—IR - (D7) The second term in brackets on the RHS of Eq. (D.6) can analogously be developed using the equation for the fluctuating velocity (i.e. Eq. (2.3)) as a basis to = —tR -V—’tR , (D.8) where f’ represents the consolidation of the divergence of the fluctuating Reynolds stress and the gradient of the fluctuating pressure. From Eqs. (D7) and (D8) the following form for the transport of the fluctuating dissipation rate due to velocity fluctuations can be obtained: = —TR -V—IZR ~é. (D9) The operator A is defined by Eq,. (2.16). 193 E. First Order Approximation of the Equilibrium Region The basis for the development of a first-order solution are the expansions of k, e and vt around the equilibrium region. These expansions are expressed as vt =v1cy’(1+a,1), (E.1) k“ = kg, + b, %, (E2) 1’ 2+ = Iii? (E.3) Ky The momentum equation as given by Eq. (5.49) in combination with the definition for the eddy viscosity (i.e. Eq. (1.25)) can be integrated twice with respect to y in order to obtain the corresponding velocity profile. It is expressed in terms of a series representation (Takemitsu, 1990) to + u+ zilny’+B— a ____l_y —+-—'(a +00%)-g;(a1+1)(1)3+0(y+4)- (E4) K K 8 5 With the velocity profile extended by the linear term only and the eddy viscosity according to Eq. (5.50) the evaluation of the LHS of Eq. (5.50) reads d d uf' d—y-(Ve dy )z— 8 , (E-S) which satisfies the momentum equation exactly. The transport equation for k and e are given to d dk’ du’ , _._ v .... ___ V _g+, E.6 dy+ (Ck D dy+) t( d +) ( ) d de’ 8+ du’ 2 01,8"2 8 D dy+) k+ Cp ((der) k+ ( ) 194 The term vD arises through the consolidation of various terms given by V k+2 .._D : 2CR V R (E8) + yy ’ E and constitutes the ratio of the effective transport viscosity to the molecular viscosity for this particular coordinate direction. Ryy may be expanded around the equilibrium region 10 _ . y Ryy _ R3 +83, (E9) so that VD can be expressed as VD y y —=N1c—1+d—, 13.10 V 5( 1 5) ( ) with N = 25*cif-kijgj- . (E.11) Upon inserting ve, k1, 8+, u+ and VD into Eq. (B.7) and comparing the coefficients which have the same power in y one obtains 0(1); 0:5——8—, (E12) y 1c 1c 8+ 0(1): —ckNxb,=—?(2+a,+c,). (E13) Eq. (E. 12) shows that the Olh-order solution is equivalent to the balance of production and dissipation. Eq. (E. 13) provides an a priori expression for the transport coefficient ck to 2+a,+cl Ck :— 2 +2 eq. ' 2b11< cheqRyy (E. 14) However, the requirement of c. which stems from experimental data for 8+ may render this expression unreliable for an a priori estimate (cf. Chapter 6). The evaluation of the 8- 195 equation yields the following result 04,) —EE—N£fc—:‘:—C—E—fl3-, y 8 K“ K“ 1 c N + c c O(;): — 5+2 (blcR +f1keq.):-K—:(Cl_al—2)—2K—3Cl' (E.15) (E. 16) Eq. (E.15) yields the analytical expression for the transport coefficient ca (see Eq.(5.48)). The evaluation of Eq. (E.16) can provide information about the behavior of cR, which may, however, render unreliable because of the earlier mentioned requirement of c1. The results of the lS’-order analysis are used for an a posteriori evaluation of the transport parameter ck. The general concept of expanding the kinetic energy profile in the inertial sublayer proves helpful in the discussion of the results of the numerical study of the channel flow. 196 F. Near Wall Analysis The general derivation of the near wall behavior of the individual terms appearing explicitly in the transport equations for the turbulent kinetic energy k and the dissipation rate 8 are given here. The basis for this derivation is given through the Taylor series expansion of the fluctuating velocity and pressure. The near wall behavior for the fluctuating velocity can be written as u’ = aly + by2 + higher order terms 11 = bzy2 + higher order terms (F1) )1 i u; = a3y + by2 + higher order terms The coefficients a,, b, and the coefficients for the higher order terms depend in general on x, z and t. Due to the property of the fluctuating velocities which must vanish when ensemble averaged, the ensemble averaged values of the coefficients must likewise vanish. The fact that u; behaves as O(yz) stems from the continuity equation for the fluctuating velocity which can be written as V - u“: 0 (F2) Eq. (F2) - evaluated at the wall - renders the coefficient of the linear term for u; zero. The near wall behavior of the fluctuating pressure can be obtained from the Navier— Stokes equation for the fluctuating velocity. This equation can be determined by subtracting the mean field equation as given by Eq. (1.1) from the instantaneous equation and is given by D!‘ 1 l 2 I O l I l 1 —DT:_5VP+VV u—u-V—u~Vu+ (F3) This equation - when evaluated at the wall - yields a boundary condition for the pressure gradient to 197 —Vp'= VVzu' (1:4) From Eq. (F.4) the pressure gradient with respect to y can be estimated to behave as O( l ). Thus, the fluctuating pressure can be written in a general form as p’ 2 p(,, + 2ubzy + higher order terms (F.5) Using Eq. (F.1) the behavior of the individual components of the Reynolds stress can be written as < u u >=< alal > y2+ < blb, > y4 + higher order terms >= < bzb2 > y4 + higher order terms 2 4 . (F.6) < u u >=< a3a3 > y + < b3b3 > y + h1gher order terms I I x x I I Y y I I z z I I y z < u u >=< bza3 > y3 + higher order terms From the summation of the normal components - which sum up to 2k - it can be seen that the leading term of the kinetic energy behaves as O(yz). The next order term for k behaves as O(y4). The behavior of the dissipation rate can be developed using its definition (see Appendix C) and Eq. (F. 1 ). Thus, the near wall behavior of e is given to e = 8w + asy + b5y2 + higher order terms (F7) The coefficients in this expansion are — like the ones in Eq. (F.1) - functions of x, z and t with vanishing ensemble averages. The value for the dissipation rate at the wall can be written in terms of the coefficients for the fluctuating velocities as 8w 2 v(< a12> + < a§ >) (ES) The mean velocity profile may similarly be expressed in terms of a series expansion to < uz >= ay + by2 + cy3 + higher order terms (F.9) However, the coefficients in the series expansion for the mean velocity are not functions of x, z and t neither do they vanish when ensemble averaged. The individual terms of the transport equations for the kinetic energy and the dissipation rate can be obtained using 198 these basic relations as given by Eqs. (F.1) — (F9). Their evaluation is given here with respect to channel flow. Only the leading term is presented. Transport Term for k d p_’ —-—= —2b_ w —y (F.10a) dy p p u' -'u —’=2b 2133(+)y (F.10b) dy 2 . Transport Term for e d — vd—y < u;(Vu')T:Vu'>= —2v2(< bzalal > + < b2a3a3 >)y (F.11a) d p. I 7 7 —2v——=—8v‘ (F.11b) dy 9 Production Terms for 8 d PE' =—2v d Z =—4vay (F.12a) y P 2vd 2v{< 8‘" >+< aa>}:.1 (F12b) =— — —— =— a — a . 5 8y 82 dy ' 8 3 8—2 y 3 ,du; d2 P, :-—2v—2=—8vby (F.12c) y dy dy ' ' P: =—2v< {(Vu‘)T ‘Vu'}:Vu'> aa, (1212.1) 8a 8 =(< a,(2a]b2+a —+a3—(_;l—l)>+)y 82 2 8a ' 8x 199 Destruction Terms for e — Y = -—2V2 < V(Vu' )T§V(Vu' )T > : —4v2{< (% 2 8a] 2 8a. 7 )->+<(—) >+<(—’)‘> (EB) 82 8x 8a3 2 2 2 2 +<(E) >+2 The modeling approaches for the transport terms in the k— and e-equation as developed in Chapter 4 (i.e. Eqs. (4.16) and (4.25)) are likewise given here in terms of their near wall behavior. Thus, the order of modeled transport term in the k—equation given by d , , dk aim, <11qu >E}=O(y3) (FM) and the equivalent term in the e-equation is given by d , , d8 3 EiCETR < uyuy >E} = O(y) (F.15) 200 G. Computer Program For Fully Developed Turbulent Channel Flow The general methodology applied to obtain the numerical solution is explained in Chapter 6. The finite difference scheme (Yang and Shih, 1993; Patankar, 1980) is illustrated here and the resulting matrix formulation is given here. The governing equations for the numerical implementation are presented in Chapter 6. The momentum (equation is used in its integrated form as given by Eq. (6.2). The treatment and numerical implementation of the transport equation for the kinetic energy is explained here. The equation for the dissipation rate parallels this procedure. The finite difference scheme is developed for interior grid points. ' z ———-> For ease of reading, the k-equation is rewritten in terms i 1+1 of the following dimensionless variables (i.e. y=n, k=k+, 1+1 . + + + - - - 8:8 , u=u , 5:5 ). The dimensronal vanable z y j introduced through the convective term is normalized with 5. j-I 0 Thus, the (normalized) k-equation is given by: u£_i(v_d_k)+8 V_ S243,3 (G1) dz dy 5 dy ' The viscosities v... and vI are given by +2 k+2 vk = 2cch .E—R” vl = g 8+ (G2) and S is defined as _E_ _ 2 (G 3) — dn ”'1 7 5 ' The finite difference scheme for the nodes as indicated by the sketch above thus reads 201 u km.) _ km _ 1 V ki+l.j+l "' ki+l.j _ V ki+l.j — ki+l,j—l H AZ 5(Yj+.5 — yJ-__5) “1+5 Yj+1 — Yj kid-'5 Y; " Yj-I (6.4) VI 2 + (:8 )i'j - 555“- With Vi<+) = Vki,j+.5 vi.) = Vki,j—.5 Ay = Yj+.5 _ yj-.5 = yJ'+l " Yj = Y) — yj-l ((35) and VI 2 8 skl : (33 )i'j sk2 = —6(E)“jk”"j (G6) the difference equation can be expressed as V(,_) uj'. V(,') +v(+) V(+) ui.‘ _ “1:302 ki+l.j-| +{ All + k5(Ay); —Sk2 km.) —_5(Aky)2 km.“l :Skl +7A—ziki'j ((3.7) which can be written in matrix notation as {bi c, 0 0 o )(ka, W” (d, \‘j’ a2 b2 c2 0 0 k2 d2 0 .. .. .. O .. = .. (6.8) O 0 a b c k d n-l r00 0 a bn)\kn} id") in which denotes a continuation in this tri-diagonal matrix scheme and n is the number of grid points. The coefficients a, b and c as well as the parameter d can be readily taken from Eq. (G7). The boundary condition at y=y(°q’) leads to the following result for the matrix coefficients at this points: a1=c1= O bI =1 (11:19.1 (G9) Eq. (6.7) is also applied at the nth grid point. The zero derivative condition (i.e. dk/dy=0) at the centerline is used to eliminate the coefficient cm. Thus, the transport coefficient vk”) is given by w”. Whence, 2v‘k" ui n 2v? ui n an=————, "2 ‘ + ,—sk0 C =0 (1 = ' k. (G.10) 5(Ay)“ Az 5(Ay)‘ “ AZ ‘ The expression for dn incorporates the fact that no turbulent production exists at the centerline. The solution of this is obtained by applying the ‘Thomas’-algorithm, i.e. a readily available tri-diagonal matrix solver (Press et al., 1992). The computer program used to calculate the boundary value problem as well as input files are given here. The variable names are consistent with the notation used in the general body of the text. Exemptions are explicitly stated here. For purposes of readability, the font to display the source code has been changed. The first file serves as input. The numerical values can be taken as an example. Input File &input1 re=395., 8‘ ymax=1.0, Outer Boundary dz=0.80, Step Size in z ystrch=1.03, Stretch Factor for Grid 1 /&end &input2 np=150, Number of Points nst=800, Number of Steps inp=O, 0 —> Top Hat Initial Profile, 1 —> Specified Input Profiles Ipr=10, Monitor Control Parameter ng=1, Grid Indicator apow=1.5, Power Exponent for Grid 2 /&end &input3 yw=30.0, Inner Boundary uw=13.8, BC for u‘ hw=3.225, BC for k’ ew=0.0813, BC for 8‘ (Note: eps=re*e) vs=0.0, Control parameter for molecular viscosity /&end &input4 u0=1.0, Top Hat Profile for u’ akO=1.0e-O, Top Hat Profile for k’ e0=1.0e-1, Top Hat Profile for e’ tiny=1 .Oe-ZO, Small Parameter to control overflow /&end &input5 Model Parameter cro=0.6667, br=0.237, ttSf=20.0, ecr=0.5, ef1 =2.0, cbo=0.355, bb=0.1 1 9, ttsb=20.0, ecb=0.5, ef2=2.0, Afb=0.30. be=3.3, Cfb=2.0, cgo=0.0068, bg=0.027, ecg=2.0, Afg=0.50, Bfg=1 2.0, Cfg=8.0, /&end &input6 203 cp=1.56, cd=1.80, ck=0.400, ce=0.240, Production Coefficient Destruction Coefficient Transport Coefficient for k‘ Transport Coefficient for 8’ /&end C Program Source Code program relax parameter (n=401) common/velo/ u(n),dudy(n) common/turb1/ak(n),eps(n),anutao(n),anutak(n),anutae(n) common/para/ re,np,ist,nst,tiny,lpr,ng,vs,apow c call init call mesh ist=0 istpr=1 c --- Loop Start 10 ist=ist+1 if((lpr.gt.0).and.(ist.eq.istpr*lpr)) then write(',110) ist,ak(np),eps(np)/re,u(np) 1 10 format(1x,i4,2f8.4,f8.2) istpr=istpr+1 endif call stress call velocity 204 call model call kequat call eequat if(ist.lt.nst) goto 10 c Loop End call stress call result stop end subroutine init parameter(n=401) common/velo/u(n),dudy(n) common/grid/x dx ,,ymax ystrch,y(n)l,dy(n) common/turb1/ak(n) pnqln) ‘ ' ’ \, ‘ f : common/para/ re ,,np ist, nst ,tiny, lpr, ng, vs ,apow common/ketra/cp, cd, ck ce common/msu1/ryy(n), ryz(n ), g(n), rxxyy(n), tts(n ), rey(n) "’ :‘ h(n) nnfn) hyxyy(n) hyy(n ) hyz(n) common/relax1/cro,br,ttsr,ecr,ef1,cbo bb, ttsb ,ecb,ef2 common/relaxZafb,bfb,cfb,cgo,bg,ecg,afg,bfg,cfg common/bound/yw,uw,hw,ew namelist/input1/re ymax dx ,ystrch r nct inn Inr ng apnw namellst/inputS/yw, uw, hw, ew, vs namelist/input4/u0, akO, e0 ,tiny "" hr ttsr err ef1,cbo,bb,ttsb,ecb ef2, + Afb, Bib, Cfb ,c,go bg, ecg Afg, Bfg, Cfg namelist/inputG/cp, cd, ck ce open(3,fl|e='Relax.inp') read(3,input1) read(3,input2) read(3,input3) read(3,input4) read(3,input5) read(3,input6) close(3) c c Convert First Point yw=yw/re c Initial Profiles if(inp.eq.0) then do10j=1,np u(j)=u0 ak(j)=ak0 eps(j)=e0 10 continue elseit(inp.eq.1) then open(10,fi|e='Relax.ini') do 20j=1,np read(10,100) y(j),ak(j),eps(j),u(j) eps(j)=eps(j)'re 100 format(7f10.4) 20 continue close(10) 205 end” return end D subroutine kequat parameter (n=401) real a(n),b(n),c(n),d(n),phy(n) common/velo/u(n),dudy(n) common/grid/x,dx,ymax,ystrch,y(n),dy(n) common/turb1/ak(n),eps(n),anutao(n),anutak(n),anutae(n) common/turb2/pk(n),sk1(n),sk2(n),pe(n),se1(n),seZ(n) common/para/ re,np,ist,nst,tiny,lpr,ng,vs,apow common/bound/yw,uw,akw,ew do 10 j:1,np phy(i)=Pk(l)/dy(i) 10 confinue c Coefficients for Matrix A do 20 j=2,np-1 phyp=(phy(i)+phy(i+1))/2- phym=(phy(i)+phy(i-1))/2- a(j)=-phym/dy(i) b0)=u(j)IdX+(phyp+phym)/dy(i)-sk2(i) c0)=-phyp/dy(i) 20 d(j)=u(j)'ak(j)/dx+sk1(j) c BC (Equil. Region) ak(1)=akw b(1)=1.0 C(1)=0.0 d(1)=ak(1) c --- BC (Center Line) a(np)=-2-'phy(np)/dy(np) b(np)=u(np)/dX+2.‘phy(np)/dy/(np-1-))"aPOW+yw 45 continue elseif(ng.eq.3) then pi = 2.‘acos(0.0) do 55 j=1,np y(j)=(ymax-yw)'(1.-cos(pi'(j-1.0)/(np-1.)))/2.+yw 55 continue endif dy(1)=(-y(3)+4.0‘y(2)-3.0'y(1))/2. dy(np)=(y(np-2)-4-0‘y(np-1)+3.0'y(np))/2~ do 40 j=2,np-1 dYO)=(Y(J+1)-Y(J-1))/2- 40 continue return end subroutine model parameter(n=401) common/velo/u(n),dudy(n) common/turb1/ak(n)eps(n), ‘ , ,, 'f ), ‘ common/turb2/pk(n),sk1(n),sk2(n), e(n),se1(n),se2(n) common/para] re,np,ist,nst,tiny,|pr,ng,vs,apow common/ketra/cp,cd.ck,ce common/msu1/ryy(n),ryz(n),g(n),rxxyy(n),tts(n),rey(n) ’ "’ 1Irho”)Pg(n).hxxyy(n).hyy(n).hy2(n) I\ l\ \ I do 10 j=1,np anuta00)=g(j)'rey(i) anutak(j)=2.‘ck‘cr(j)'rey(j)'ryy(j) anutae(j)=2.'ce’cr(j)'rey(j)'ryy(j) pk(j)=(vs+anutak(j))/re pe(j)=(vs+anutae(j))/re 10 continue do 20 j=2,np-1 sk1(j)=anutao(j)/re‘dudy(j)“2 sk2(j)=-eps(j)/(ak(j)+tiny) 208 tt=re'ak(j)/(eps(j)+tiny) fp=1,+(3.24—1.)'exp(-(rey(j)/30.)“1.0) se1(j)=cp'fp/(tt+tiny)'re'sk1 (j) fd=1.-O.4/1.8'exp(-(rey(j)/6.)"2.) fdw=1./(1.+0.0837'tts(j)'exp(-(rey(j)/50.))) se2(j)=—cd'fd'fdw/(tt+tiny)'re 20 continue I. v return end write(1,')’ y',‘ U <-uv> 1 . subroutine result parameter (n=401) common/velo/u(n ), dudy(n) common/grid/x, dx,ymax, ystrch y(n), dy(n ). common/turb1/ak(n) pncfn) ‘ . common/para/ re, np, ist, nst ,tiny, lpr, ng, vs ,apow dimension uvbar(n) open(1,fi|e=‘Relax.dat‘) k‘,' eps ' do 10 j=1,np uvbar(j)=anutao(j)'dudy(j)/re write(1,100)y(j)‘re,u(j),uvbar(j),ak(j),eps(j)/re 10 continue 100 format(1x,6e12.4) close(1) return end subroutine stress parameter (n=401) common/velo/u(n ), dudy(n ) common/para/ re ,np ist, nst ,tiny, lpr, ng, vs ,apow common/turb1/ak(n) pn<(n) common/msu1/ryy(n ), ryz,(n) g(n), rxxyy(n ),tt sn( ),rey(n) "’ :‘ ,r‘hf‘n) ng(n) hxxyy(n )hyy(n),hy2(n) common/relax1/cro, br, ttsr, ecr, ef1, cbo, bb ,ttsb,ecb,ef2 common/relax2/afb,bfb,cfb,cgo,bg,ecg,afg,bfg,cfg do 10 j=1,np ak(j)=abs(ak(j)) ePSU)=ab5(ePS(i)) dudy(j)=abs(dudy(j)) rey(j): ak(j) 2."/(eps(j)+tiny) re tts(j): ak(j)/ e(ps(j) +tiny)'dudy(j) 11 =exp(- (tts(j)/ttsr)"ef1) f2=exp(-(tts(j)/ttsb)"ef2) cr(j)=cro/(f1+br‘sqrt(tts(j))) fb=1 .-Afb'exp(-((tts(j)-be)/Cfb)“2.) cb(j)=cbo/(f2+bb'sqrt(tts(j)))‘fb fg=1.-Afg'exp(-((tts(j)-Bfg)/Cfg)"2.) cg(j)=cgo/(1.+bg'(tts(j)"2.))'fg ryy(j)=1./(3.+(cr(j)'tts(j))"2.) g(j)=2.‘cr(j)'ryy(j)-cb(j)/(1.+2./3.‘(cr(j)'tts(j))"2.) 10 209 anutao(j)=g(j)'rey(j) rxxyy(j)=cg(i)‘(1180)"?-) WZG)=9(J')/2-‘tt8(l) hyy(j)=-cb(j)‘cr(j)/3.'(tts(j)“2.) + /(1.+2./3.*(cr(j)‘tts(j))"2.) hxxwfll=rxxyv01 hyz(j)=cb(j)/2./(1.+2./3."(cr(j)*tts(j))“2.)*tts(j) confinue return end 11 12 subroutine tridag(a,b,c,d,u,n) parameter (nmax=401) dimension gam(nmax),a(n),b(n),c(n),d(n),u(n’) if(b(1).eq.0.)pause bet=b(1) u(1)=d(1)/bet do 11 j=2,n gam(j)=c(j-1)/bet bet=b(i)-a(j)*gam(j) if(bet.eq.0.)pause U(j)=(d(j)-a(i)*U(i-1))/bet conflnue do 12 j=n-1,1,-1 U0)=U(j)-gam(i+1)’U(i+1) confinue return end LIST OF REFERENCES LIST OF REFERENCES Abid, R. and C. G. Speziale, 1993, “Predicting equilibrium states with Reynolds stress closures in channel flow and homogeneous shear,” Phys.Fluids A, 5(7), 1776—1782. Batchelor, G. K. and A. A. Townsend, 1948a, "Decay of isotropic turbulence in the initial period", Royal Society of London Series A, Cambridge Press. Batchelor, G. K. and A. A. Townsend, 1948b, "Decay of isotropic turbulence in the final period", Royal Society of London Series A, Cambridge Press. Bird, R. B., R. C. Armstrong and O. Hassager, 1987, Dynamics of Polymeric Liguids, Vol. I, New York, John Wiley & Sons. Bird, R. B., R. C. Armstrong and O. Hassager, 1987, Dynamics of Polymeric Liquids, Vol. H New York, John Wiley & Sons. Boussinesq, J ., 1877, “Theorie de l'ecoulement tourbillant,” Mem. Pres. Div. Sav., Acad. Sci. Inst. France, 23, 46-50. Choi, 1993, “A study of the return to isotropy of homogeneous turbulence,” PhD. Thesis, Cornell University. Choi, K.-S. and J. L. Lumley, 1984, "Return To Isotropy Of Homogeneous Turbulence Revisited", Turbulence And Chaotic Phenomena In Fluids, Kyoto, Japan, North-Holland. Chou, P. Y., 1945, “On Velocity Correlations and the Solutions of the Equations of Turbulent Fluctuation,” Quart.Appl.Math., 3(1), 38—54. Comte-Bellot, G., 1965, “Ecoulement Turbulent Entre Deux Parois Paralléles,” Report No. 419, Ministere de l'Air. Comte-Bellot, G. and S. Corrsin, 1966, “The use of a contraction to improve the isotropy of grid-generated turbulence,” J.Fluid Mech., 25(4), 657-682. 210 211 Comte-Bellot, G. and S. Corrsin, 1971, “Simple Eulerian time correlation of full- and narrow-band velocity signals in grid-generated, 'isotropic' turbulence,” J.Fluid Mech., 48(2), 273-337. Daly, B. J. and F. H. Harlow, 1970, “Transport Equations in turbulence,” Phys.Fluids, 13(1 1), 2634-2649. Dean, R. B., 1978, “Reynolds Number Dependence of Skin Friction and Other Bulk Flow Variables in Two-Dimensional Rectangular Duct Flow,” J.Fluids Eng., 100, 215-223. Demuren, A. O. and S. Sarkar, 1993, “Perspective: Systematic Study of Reynolds Stress Closure Models in the Computations of Plane Channel Flows,” J.Fluids Eng., 115, 5-12. Denn, M. M., 1990, “Issues in Viscoelastic Fluid Mechanics,” Annu. Rev. Fluid Mech., 22, 13-34. ‘Durbin, P. A., 1991, “Near-Wall Turbulence Closure Modeling Without "Damping Functions",” Theoret.Comput.Fluid Dynamics, 3, 1-13. Durst, F., J. Jovanovic and J. Sender, 1993, "Detailed Measurement of the Near-Wall Region of a Turbulent Pipe Flow", 9th Int'1.Symp.Turb.Shear Flows, Kyoto, Japan, August 16-18, 1993, Berlin: Springer-Verlag. Favre, A. J ., J. J. Gaviglio and R. Dumas, 1957, “Space-time double correlations and spectra in a turbulent boundary layer,” J.Fluid Mech., 2, 313-342. Gatski, T. B. and C. G. Speziale, 1992, “On Explicit Algebraic Stressmodels for Complex Turbulent Flows,” Rep. No. 92-58, NASA Langley Research Center. Gence, J. N. and J. Mathieu, 1980, “The return to isotropy of an homogeneous turbulence having been submitted to two successive plane strains,” J.Fluid Mech., 101(3), 555-566. Gibson, M. M. and V. E. Kanellopoulos, 1987, "Turbulence Measurements in a Nearly Homogeneous Shear Flow", Transport Phenomena in Turbulent Flows: Theory, Experiment, and Numerical Simulation, Hemisphere, Tokyo. Hanjalic, K., 1970, “Two-dimensional asymmetric turbulent flow in ducts,” PhD. Thesis, University of London. Hanjalic, K. and B. E. Launder, 1972a, “Fully developed asymmetric flow in a plane channel,” J. Fluid Mech., 51(2), 301-335. Hanjalic, K. and B. E. Launder, 1972b, “A Reynolds stress model of turbulence and its application to thin shear flows,” J. Fluid Mech., 52(part 4), 609-638. 212 Hanjalic, K. and B. E. Launder, 1976, “Contribution towards a Reynolds-stress closure for low-Reynolds-number turbulence,” J. Fluid Mech., 74(part 4), 593-610. Harlow, F. H. and P. I. Nakayama, 1967, “Turbulence Transport Equations,” Phys.Fluids, 10(11), 2323-2332. Hinze, J. 0., 1959, Turbulence, New York, McGraw-Hill. Hoffman, G. H., 1975, “Improved Form Of The Low Reynolds Number k-e Turbulence Model,” Physics of Fluids, 18(3), 309—312. Hussain, A. K. M. F. and W. C. Reynolds, 1975, “Measurements in Fully Developed Turbulent Channel FLow,” J.Fluids Eng, 97(4), 568-580. Jeffreys, H., 1929, The Earth, New York, Cambridge University Press. Johansson, A. V. and P. H. Alfredsson, 1986, “Structure of turbulent channel flow,” in Flow Phenomena and Measurement, Vol. 1, N. P. Cheremisinoff (editor), London, Gulf Publishing Company, 825-869. Jones, W. P. and B. E. Launder, 1972, “The Prediction of Laminarization with a Two- Equation Model of Turbulence,” Int.J.Heat Mass Transfer, 15, 301-314. Jones, W. P. and B. E. Launder, 1973, “The Calculation of Low-Reynolds-Number Phenomena with a Two-Equation Model of Turbulence,” Int.J.Heat Mass Transfer, 16, 1119-1130. Kim, J ., P. Moin and R. Moser, 1987, “Turbulence statistics in fully developed channel flow at low Reynolds number,” J.Fluid Mech., 177, 133-166. Kim, J ., 1989, “Personnal Communications (1994).” Kim, J. and R. A. Antonia, 1993, “Isotropy of the small scales of turbulence at low reynolds number,” J.Fluid Mech., 251, 219-238. Kreplin, H.-P. and H. Eckelmann, 1979, “Behavior of the three fluctuating velocity components in the wall region of a turbulent channel flow,” Phys.Fluids, 22(7), 1233- 1239. Lam, C. K. G. and K. Bremhorst, 1981, “A Modified Form of the k-E Model for Predicting Wall Turbulence,” Journal of Fluids Engineering, 103, 456-460. Laufer, J., 1951, “Investigation of Turbulent Flow in a Two-Dimensional Channel,” NACA Report No. TR 1053. 213 Laufer, J ., 1954, “The Structure of Turbulence in fully Developed Pipe Flow,” NACA Report No. 1174. Launder, B. E. and D. B. Spalding, 1972, Mathematical Models of Turbulence, Academic Press. Launder, B. E., A. Morse, W. Rodi and D. B. Spalding, 1972, "Prediction Of Free Shear Flows - A Comparison Of The Performance Of Six Turbulence Models ", Free Turbulent Shear Flows, Langley Research Center, Hampton, Virginia, National Aeronautics and Space Administration. Launder, B. E., G. J. Reece and W. Rodi, 1975, “Progress in the development of a Reynolds-stress turbulence closure,” J. Fluid MEch., 68(part 3), 537-566. Lee, M. J ., J. Kim and P. Moin, 1990, “Structure of turbulence at high shear rate,” J.Fluid Mech., 216, 561-583. Leonov, A. I. and A. N. Prokunin, 1994, Nonlinear Phenomena in Flows of Viscoelastic Polvfimer Fluids, London, Chapman & Hall. Lumley, J. L. and G. R. Newman, 1977, “The return to isotropy of homogeneous turbulence,” J.Fluid Mech., 82(1), 161-178. Lumley, J. L., 1978, “Computational Modeling of Turbulent Flows,” in Advances in Applied Mechanics. C.-S. Yih (editor). New York, Academic Press. 18: 123-176. Mansour, N. N., J. Kim and P. Moin, 1988, “Reynolds-stress and dissipation-rate budgets in a turbulent channel flow,” J.Fluid Mech., 194, 15-44. Mansour, N. N., J. Kim and P. Moin, 1989, “Near-Wall k-e Turbulence Modeling,” AIAA, 27(8), 1068-1073. Mansour, N. N. and A. A. Wray, 1994, “Decay of isotropic turbulence at low Reynolds number,” Phys.Fluids, 6(2), 808-814. Mitsoulis, E., 1986, “Numerical Simulation of Viscoelastic Fluids,” in Polymer Flow Engineering, N. P. Cheremisinoff (editor), London, Gulf Publishing Company. 9: 649- 704. Mohamed, M. S. and J. LaRue, 1990, “The decay power law in grid-generated turbulence,” J.Fluid Mech., 219, 195-214. Morse, P. and H. Feshbach, 1953, Methods of theoretical physics, part I, McGraw-Hill. 214 Nisizima, S. and A. Yoshizawa, 1987, “Turbulence Channel and Couette Flows Using an Anisotropic k—e model,” AIAA, 25(3), 414—420. Oldroyd, J. G., 1958, “Non-Newtonian effects in steady motion of some idealized elastico-viscous liquids,” Proc.Roy.Soc.Lon.Ser.A, 245, 278-297. Parks, S. M., 1997, “Relaxation mode] for homogeneous turbulent flows,” PhD. Thesis, Michigan State University. Patankar, S. V., 1980, Numerical Heat Transfer and Fluid Flow Washington, Hemisphere. Patel, V. C., W. Rodi, et al., 1985, “Turbulence Models for Near-Wall and Low Reynolds Number Flows: A Review,” AIAA, 23(9), 1308-1319. Pope, S. B., 1975, “A more general effective—viscosity hypothesis,” J.Fluid Mech., 72(p.2), 331-340. Potter, M. C. and J. F. Foss, 1982, Fluid Mechanics East Lansing, Great Lakes Press, Inc. Prandtl, L., 1925, “Bericht fiber Untersuchungen zur ausgebildeten Turbulenz,” ZAMM, 5(2), 136—139. Press, W. H., S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, 1992, Numerical Recipes in FORTRAN, New York, Cambridge University Press. Reichardt, H., 1938, “Messungen turbulenter Schwankungen,” Die Narurwissenschaften, 26(24/25), 404-408. Reichardt, H., 1940, “Die Warmeiibertragung in turbulenten Reibungsschichten,” Z.A.M.M., Bd. 20(6), 297-328. Rivlin, R. S., 1957, “The relation between the flow of non-Newtonian fluids and turbulent Newtonian fluids,” Q.Appl.Maths, 15(2), 212—215. Rodi, W., 1976, “A New Algebraic Relation for Calculating the Reynolds Stresses,” ZAMM, 56, T219-T221. Rodi, W. and N. N. Mansour, 1992, “Modeling The Dissipation-Rate Equation With The Aid Of Direct Simulation Data,” in Studies In Turbulence. T. B. Gatski, S. Sarkar and C. G. Speziale (editors), Springer-Verlag: 17-38. Rotta, J ., 1951, “Statistische Theorie nichthomogener Turbulenz. 2. Mitteilung,” Z.Phys., 131, 51—77. 215 Rotta, J ., 1951, “Statistische Theorie nichthomogener Turbulenz. 1. Mitteilung,” Z.Phys., 129, 547-572. Schlichting, H., 1951, Boundag-Layer Theory, McGraw-Hill, New York. Schumann, U., 1977, “Realizability of Reynolds-stress turbulence models,” Phys.Fluids, 20(5), 721-725. Shih, T.—H. and J. L. Lumley, 1993, “Remarks on Turbulent Constitutive Relations,” Math. Comput.Modelling. 18(2), 9-16. Shih, T.-H., J. Zhu and J. L. Lumley, 1995, “A new Reynolds stress algebraic equation model,” Comput.Methods Appl.Mech.Engrg., 125, 287—302. Speziale, C. G., 1982, “On turbulent secondary flows in pipes of non-circular cross- section,” Intl J.Engng Sci, 20(7), 863—872. Speziale, C. G., 1987, “On nonlinear K—l and K-e models of turbulence,” J.Fluid Mech., 178, 459—475. Speziale, C. G., 1991, “Analytical Methods for the development of Reynolds-stress Closures in Turbulence,” Annu. Rev. Fluid Mech., 23, 107—157. Takemitsu, N., 1990, “An Analytical Study of the Standard k-e Model,” J.Fluids Eng, 112, 192-198. Tanner, R. I., 1985, Engineering Rheology, Clarendon Press, Oxford. Tavoularis, S. and S. Corrsin, 1981, “Experiments in nearly homogeneous turbulent shear flow with a uniform mean temperature gradient. Part 1,” J.Fluid Mech., 104, 311-347. Tavoularis, S., 1985, “Asymptotic laws for transversely homogeneous shear flows,” Phys. Fluids, 28(3), 999-1001. Tavoularis, S. and U. Karnik, 1989, “Further experiments on the evolution of turbulent stresses and scales in uniformely sheared turbulence,” J.Fluid Mech., 204, 457—478. Tennekes and Lumley, 1972, A first course in turbulence, MIT Press, Cambridge, Massachusetts. Townsend, A. A., 1976, The structure of turbulent shear flow, Cambridge, Massachusetts. Tucker, H. J. and A. J. Reynolds, 1968, “The distortion of turbulence by irrotational plane strain,” J.Fluid Mech., 32(4), 657-673. I, ’ 'Fil'i ' I, b I! ’Lll’llli Ill 1 216 van Driest, E. R., 1956, “On Turbulent Flow Near A Wall,” J. Aero. Sci., 23, 1007—1011. von Karman, T., 1930, “Mechanische Ahnlichkeit und Turbulenz,” Nachr. Ges. Wiss. Gottingen, Math.-Phys.Kl., 58-76. Wei, T. and W. W. Willmarth, 1989, “Reynolds-number effects on the structure of a turbulent channel flow,” J. Fluid Mech., 204, 57-95. Yang, Z. and T. H. Shih, 1993, “New Time Scale Based k-e model for Near—Wall Turbulence,” AIAA Journal, 31(No. 7), 1191-1198. Yoshizawa, A., 1984, “Statistical analysis of the deviation of the Reynolds stress from its eddy viscosity representation,” Phys.Fluids, 27(6), 1377-1387. Zaric, 2., 1972, “Wall Turbulence Studies,” Adv.Heat Transfer, 8, 285-350. 14311.31... . . 2.1.. .1.... .. \..-4i.(.uf.l.l . .1.. . I11.tt.. .t .11......l.1.t.1.1\vl. .1. wixtill. 1131.11... 11.11.1111! 1......) . J . .v 1111.11.11.32... . 1.11). 1.111.111.1111... . . .. «1.3.11. I ... . . .. ..t....1o.1..l.xcf11 . . . 11.11151. . .t....ax(.rr to .tttrr45101k 10.111111 .\.-4.111.315.1111 . 1:)! 1911111512012 ‘11 111.151. . 11.111. 111.. 1.111. . 114.11.51.11; :1.) .1111 .1141“? 1. “481g! Illicit. 5). .i . .1 ...: . 1.1.1.0. 11.11211 11.1.1114}... ...tlsti. 1.1.3.1111} .11.. 11.3 1.31) 541111... ..‘11111111. cilia-0111.21.10.11! \cti.l?zh\.1\rif\.\\af\rinl ....7.~vt.1.tittr..tt..f.iii.-.1.§1111 1111.. . . 1.. . r.~..§‘.\.7%‘l‘n§“lt‘nv‘i§ nit-52.115 E . 111.511.111.111. .1....1‘3911‘61 . 511151...) .75 .1. ...4.. . 3.1.1.1..in tits}: ...... . . . . . ."11..\Ufl1..s~111.|\1.1.1¢11¢. _. . . . . ......~.1sltl.111~1. .\ 11...... 1‘01 Iiitiv‘l‘ ‘4 .\ittvlgltttftoii 1.1. 111.1 I 1141911111113; 4 11131111111 . 1:1... 111.041.151.11... ..1 .. . (.\.-)1... 11.1. I #2:. ‘3 :1... 1.21.1111... 1:... H1111; (1111011111111!!! 15.414111111101111. 114 1.111 5.1.1.10. 1.111131% v9.1.1 1% 11%).! a. 1“ filli.§.1..i.tlt.n11¢iofii§ 111.citllbtllltiut14‘“ 1111.1 .. . 41...». it... 1.2111. I111. 111.111.). 1.1.1.3. . . l {1151.11.11.13 .931}... . .. .11....‘3... . stilt. ...u .!\1 11.91.14: 11.11.411.355. .... .11 ll .... . .11 flukflflvllnf..1111.ht . u .04.) Jig - 3.1.1.111 91.11.10.016 115.194.1141.. . 3.11.1111 91» 7130.1. 15113.05!!! . . 1‘4; . .1181?! 0.31.11“: . 9‘ 1.10000”... . . innvlfluufloflhioihiiiitifil .1.»... ......tif. ...-.14.. 11:151.... 5. Lip... .\..11...H.)|t...\.1.i...1.1. 1.5.1.. . .... .izi..117\t..rll.roaf. . 17.4.11... II... 111.11.11.5512196Liliit '1 .s. 5.7.1.11 Et.kt§n.rwii agilzlinist. 1.1-1.111.411.3311} . . .a \...I._(.1£Sr.\.\11p~.fi\£r)3u Eihyslisxt . tall}. .1 I J(. (.(f. _ G. .t .1 .... I... .....i: .1.. 1!. . 1..-}...(23.--) .... .. 3.... ....Iiltt. 3.3.. :5. . 75.32”: i.e..-. -111. 3 11412.72... ... . 11131911111 . . ..rl... I 411.1111. 1 11 ....I L . 11.11.11.» n 14.11.31. 1. t flail! .. ...: 7ft}... . ‘ 1!. . ... V. 4... 11:. 137211. 111.24.)... . .rrllrill1\.>\ 7... ......ltttitio .11.. su ltirfltiot.lnstt..: . . . . 1.111.111.1311 Ellisflttiid 5.7.5117. 1: ... 3.1.1.1.. 1.11.1: 4.1.11 ». it . 11‘s....- \1ttllf-)..liyllh§xr.\t¢xilr11.1 c241 . . .\Ktt‘tsitslwtf‘ilt. 1 .1... .3... . 2411711111251»; 1.511.. 1 it? 1 11112111.! . liftinkiftfvts 1.1.1 1.1.71.4 . . $1111.14); 1.11.2711! .. ‘ 1 )1! .v 38%.: lehrsLUnP. fisttdvutltiflnfilfliiifi t... . . I .\.. I. . . .... Irritant n. ..1 . r .21.. If... t itwu stilt... . .11! t 1103‘. \ I t .... .. .....i...h.t.finn.. . fl. i5) 15111.1. trivial... 7.1.1.31 7771511.)... . \ €3.11... 1.9.1.5 .5. A iii‘ktriffstx. 1.77:) 1 .. . . . .1245}; .\I0. .15.... .tJfi.‘Mth—LHNN....~JH.13A! . v.... 1.11.! 141 13.1.1. .31....Itttl.zu Let I...“ lt{.-‘nld..1.\p~l‘.|§ 1.1.071“ ‘ o - lili‘.‘ 0., .- tin. tat-'1‘ 11x31 HICSfiauhi. 11 1.41.713...“ IthrtOssZtttswtslvfxli its ..n.......I..flfiliflahsift.ziivla.llzliititi1l2fizillv I. . . ...“..DLD... 1.. . . 6.1.5.!!!11317‘...11 1 $1.11 kl.» . . . 3.11.3. .....r..~ .~ $11§V¥Idfiixxh _. I...» ...4}... 12 .