The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics (1990)

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

RNG Modeling Techniques for Complex Turbulent Flows S. A. Orszag and G. E. Karniadakis Princeton University, Princeton, USA A. Yakhoti Ben-Gurion University of the Negev, Beesheva, Istael V. Yakhot Princeton University, Princeton, USA Permanent address: Ben-Gurion University of the Negev Beesheva, Israel Abstract In this paper, we combine RNG modeling tech- niques with spectral element discretization proce- dures to formulate an algorithm appropriate for simulating high Reynolds number turbulent flows in complex geometries. Three different approaches of modeling are followed based on RNG algebraic, differential k - c, and subgrid-scale models for the turbulent viscosity. Results obtained for the fully developed channel flow, and for the separated flout in a cavity and over a backwards- facing step sug- gest that all three formulations are suitable for tur- bulent flow simulations. The implementation of RNG models within the framework of a high-order discretization scheme (i.e. spectral element meth- ods) is essential, as it results in resolution of fine turbulence structures even with the simple differ- ential k - ~ model at relatively small number of degrees of freedom. 1. RENORMALIZATION GROUP FORMULATION Renormalization Group Theory (RNG) has proven to be a very useful tool in theoretical physics and statistical me- chanics to study systems with a very large number of de- grees of freedom. The most prominent use of this approach was in developing the theory of phase transitions of the second kind. Although RNG theory has found its way into fluid mechanics relatively recently, particularly in develop- ing turbulence theory [1], it has already proved to be a valuable research tool, and it has provided a series of inter- esting theoretical and numerical results. A milestone for RNG theory in fluid mechanics was the establishment of the so-called correspondence principle, stating that in the inertial range the behavior of the small-scale Navier-Stokes turbulence is statistically equivalent to the modeled Navier- Stokes equation with the additicjn of a random noise term. This principle makes it possible to use all the formalism of classical RNG theory. Recently, renormalization group methods have been devel- oped [1i, [2] to analyse a variety of turbulent flow problems. For homogeneous turbulent flows, such important quanti- ties as the Kolmogrov constant, Batchelor constant, tur- bulent Prandtl number, rate of decay to isotropy, skewnes factor, etc. have been obtained directly from this theory in good agreement with available data. Efforts in developing RNG methods for sub-grid (large-eddy simulation) model constants have also been notably successful [34. RNG methods involve systematic approximations to the full Navier-Stokes equations that are obtained by using perturbation theory to eliminate or decimate infinitesimal bands of small scale modes, iterating the perturbation pro- cedure to eliminate finite bands of modes by constructing recursion relations for the renormalized transport coeffi- cients, arid evaluating the parameters at a fixed point in the lowest order of a dimensional expansion around a certain critical dimension. The decimation procedure, when ap- plied successively to the entire wavenumber spectrum leads to the RNG equivalent of full closure of the Reynolds av- eraged Navier-Stokes equations. The resulting RNG trans- port coefficients are differential in character as opposed to ad hoc algebraic coefficients of conventional closure meth- ods. All constants and functions appearing in the RNG closures are fully determined by the RNG analysis. In essence, the RN(; method provides an analytical method to eliminate small scales from the Navier-Stokes equations, thus leading to a dynamically consistent description of the large-scales. The formal process of successive elimination of small scales together with re-scaling of the resulting equa- tions results in a calculus for the derivation of transport approximations in turbulent flows. 35

The renormalization group method has been used to derive three different types of models: an algebraic eddy viscos- ity model, a differential k-e, and a subgrid scale turbulence model for large eddy simulationst1~. The governing equa- tions for the flow motion are the Navier-Stokes equations: Bt + Vj00Vi = _ PIP + ~ ~ deli along with the incompressibility constraint dvi = 0 (2) Here z, is the total viscosity defined as z' = zoo+ Z'T (the sum of the molecular and turbulent viscosity, respectively). RNG Algebraic Model This model is the simplest of all, however it is not very general as it requires the a priori postulation of the charac- teristic integral length scale, and thus an assumption needs to be made based on physical considerations. The eddy vis- cosity ~ is obtained from the following relation: ~ = ',o~1 + H¢~2 ye 3~4 _ C'4~/3 (3) where H(x) is the Heaviside function defined by H(x) = :z: for x > 0 and H(~) = 0 otherwise, and /\ is an integral length scale of the turbulence in the inertial range. The constants a = 0.120 and C = 100 are derived in t13. The mean dissipation rate ~ can be expressed entirely through the resolvable field as L, ~ Levi + ~Vj )2 (4) The total viscosity is then obtained by solving an algebraic cubic equation at every node of the computational domain at each time step. RNG k-e Model This is a differential model based on differential recursion relations discussed in t13. It is free of adjustable parame- ters, however two additional equations are derived for the turbulent kinetic energy k, and the rate of dissipation e. Dk p ~ ok Dt Ski ( 6mi ) Dt 0zi ( B~i (5) (6) Here we define the production term P based on the turbu- lent viscosity p AT ( Levi + ~Vj)2 and the terms P. Y are defined subsequently from: P = 0.656~/2YCP (7b) y = 63/2y (7c) The RNG procedure provides the additional two differ- ential equations for Ye and ti/2 ~ which can be integrated as follows yul/2 ~ _ ~ ~ (t3 - 1 + C)~/2 (8a) ~/7 = C\k + ilk ~ (t3— 1 + C)~/2 (8b) Here the coefficients (c, and (°rk,~k) are known con- stants derived from the asymptotic behavior of the model, while the parameter ct involved in equations (5-6) is the inverse total Prandtl number defined from the RUG alge- braic relation (coo refers to molecular properties) —1.3929 ~0.632~ a+2~3929 ~0.3679_ ~0 (9) a0 - 1.3929 DO ~ 2.3929 ~ The total viscosity z' is implicitly defined from equations (8) at each node of the computational domain. In the high Reynolds number region where ~ ~ HOD we obtain k2 ~ = 0.0845— (10) This high Reynolds number approximation is similar to that commonly used in algebraic models. Thus, the al- gebraic models contain terms like A, which diverge when k ~ O (i.e., near wall regions or separation zones) thus creating immense difficulties computationally. The differ- ential relations obtained via RNG techniques, however, do not contain such singular behavior and thus have great pm tential for success with separating flows. The differential relations provide definitive interpolation formulae to con- nect low and high Reynolds regions of the flow. RNG subgrid-scale Model In large eddy simulations (LES) the velocity field vi is de- composed into large scales vi and subgrid components v,. It is the modeling of terms involved these latter compo- nents that is crucial for accurate computations as previous attempts using various subgrid models t4] have indicated. The RNG subgrid model is derived by elirrunation of modes from the interval AO A, where A = ~y(~/4 (] = 0.20) is a dissipation cut-off limit, and the wa~revector AO = lie~r can be expressed through the computational mesh size A. The general RNG derived viscosity given by (il) can then be expressed in terms of a length scale /`, which in LES represents the width of a suitably chosen Gaussian (7a) filter. Thus, equation (11) reduces to an identicalequation as in (3), where /& is given by the integral /~4 = J ; J (k2 + p2 + q237/2 (12) 36

Here Ai = or// (for i = 1,2,3) and /\, = 24. The above integral can be readily evaluated by breaking it up to three asymptotic integrals and a finite triple integral that can be computed numerically. In cases of simple geometries (i.e. plane channels) the integral can be evaluated through algebraic relations [3~. 2. SPECTRAL ELEMENT METHODOLOGY Let us denote by Q the three-dimensional computational domain, and by ~Q the computational boundary-surface; we can then rewrite the governing equations (1) in the form, dvj _0pi + [3 ~z'~Vi) + Hi in Q (13) where Hi includes all nonlinear and forcing terms (e.g. the random force to be included in LES). The separation of linear and nonlinear terms in the above equation leads naturally to a mixed (explicit-implicit) time advancement scheme. The temporal discretization proceeds by employ- ing the splitting scheme; accuracy and relative advantages of the scheme in the current applications are discussed in [5~. The resulting system of equations is a system of sepa- rately solvable elliptic equations for pressure and velocity. The discretization of these equations in space is obtained through the spectral element method. For simplicity we only present the two-dimensional equations, however our implementation is general for truly three-dimensional ge- ometry and flow. Here we follow a 'layered' approach, ac- cording to which the discretizations and solvers are con- structed on the basis of a hierarchy of nested operators proceeding from the highest to the lowest derivatives. This approach is motivated by the fact that the highest deriva- tives in an equation govern the continuity requirements, conditioning, and stability of the system. Given the brevity of the current paper we shall limit our description of spec- tral element methods to the innermost layer, the elliptic 'kernel', which represents both the pressure equation and the viscous corrections in equation (13~. A typical elliptic equation for a field variable ~ can be put in a standard Helmholtz equation form with variable coefficients as follows ~ `p 6+ ~ >2¢ = f in Q for j = 1,2 (14) In addition, let us assume homogeneous boundary condi- tions ~ = 0 on Q. Equation (14) can then be further discretized using planar spectral elements in plane x—y. If we define Ho the standard Sobolev space that contains functions which satisfy homogeneous boundary conditions, and introduce testfunctions ?,6 ~ Ho, we can then write the equivalent v~ariational~tatement of (14) as, / unit?/) ~ ~ ds + A2 J.n Adds = - J.n Fed (15) The spectral element discretization corresponds to numer- ical quadrature of the variational form (15) restricted to the space Xh C Ho. The discrete space Xh is defined in terms of the spectral element discretization parameters (`K, N:, Ad), where K is the number of "spectral elements", and N., N2 are the degrees of piecewise high-order polymo- mials in the two directions respectively that fill the space Xh. By selecting appropriate Gauss-Lobatto points (k and corresponding weights ppe = pppq, equation (15) can be re- placed by, i=1 p=0 g=0 ~ 0=j 0~j]{Pq ~ ~ ~ ~ ~ Pm art = K N1 N. ~ ~ ~ Pm]"[¢f]~[16~ k=1 p=0 q=0 Here Jkq is the Jacobian of the transformation from global to local coordinates (m, y) ~ (r, I), for the two-dimensional element k. The Jacobian is easily calculated from the par- tial derivatives of the geometry transformation r-, r`,, so, s.... The next step in implementing (16) is the selection of a basis which reflects the structure of the piecowise smooth space Xh. We choose an interpolant basis with components defined in terms of Legendre-Lagrangian interpolants, h`(rj) = 6,j. Here, rj represents local cordinate and [,j is the Kronecker- delta symbol. It was shown in [6], [7] that such a spectral element implementation converges spectrally fast to the ex- act solution for a fixed number of elements K and N ~ oo, for smooth data and solution, even in non-rectilinear ge- ometries. Having selected the basis we cam proceed in writing the local to the element k spectral element approximations for For ~k) as follows, ok = ~m"~(r~kn(~) Vm, n ~ (O. ..., Nay, (O. ..., N2), (17a) where ilk n iS the local nodal value of ¢. The geometry is also represented via similar type tonsorial products with same-order polynomial degree, i.e. 8, Y) = ~ sm"'Ymn~hm(7~)hn~s) Vm, 7: ~ (O. ..., N'.), (O. ..., N2), (17b) Here Ok n' yk n are the global physical coordinates of the node mn in the k element. This isoparametric mapping leads to a compatible pressure formulation without the presence of spurious modes t54. We now insert (17) into (16) and choose test functions Can, which are non-vanishing at only one global node to arrive at the discrete matrix system. This procedure is straight- forward and here we cite the final matrix system, 37

K Ni N2 ~~ 1 ~~ ~~ (H~!mn+H,~Y~!mn+) 2 JO Bk Bkn)<J)\ n k=1 m=0 m=0 = _ ~ / ~ ~ J,JB'mB.'nfm", (18a) k=1 m=0 m=0 where / denotes direct stiffness summation for the globe' system to insure that the ensemble is performed in space Hi. The x-component, for example, of the Helmholtz op- erator is defined as follows, IT nk=ppqJ~qvpql(rz)pqDpiDpminq+(~)pqDq9~qntmp + (r2~)p'DpiD~`mp + (r2s:)pqD~Dpm§~ (lab) Here low = u(~.` ), and Dij = £; all other parame- ters have been defined previously. The mass matrix Bski is diagonal and is defined as B,k; = Pitij. The y-component of the Helmholtz operator is defined sirrularly. The natural choice of solution algorithm for a time- and space-dependent viscosity z'(:~:,y~t) is an iterative proce- dure; to date both conjugate gradient techniques and multi- grid methods have been implemented for elliptic equations 6,8. The advandage of the formulation proposed here as compared to the formulation in [6] is that the splitting scheme results in separate,elliptic equations for the pres- sure and velocity that can be very efficiently solved using those iterative techniques. We have found that, for spec- tral element discretizations involving elements of low as- pect ratio, multigird methods converge much faster than conjugate gradient methods; their convergence rate, how- ever, is greatly deteriorated for large aspect ratio or very deformed spectral elements. A more quantitative analysis of the computational complexity as well as of the conver- gence properties of the aforementioned iterative methods in spectral element approximations, and in particular in the context of parallel implementation, is given in A. 3. NUMERICAL RESULTS In this section we first present results for the flow in a lid-driven cavity obtained using the algebraic model and two-dimensional spectral element simulations. We then validate the RNG methodology by simulating two different flows for which detailed experimental data and results from direct numerical simulations are available: the fully devel- oped turbulent channel flow, and the flow over a backwards- facing step. These results are obtained using the k—~ and the subgrid-scale RNG models. Turbulent Flow in a Lid-Driven Cavity - Here we consider a square two-dimensional cavity in 2:—y with mean flow being two-dimensional also. The exact geometry and the corresponding spectral element mesh is Fig. 1 Spectral Element Mesh employed for the simulation of £10w in a lid~riven cavity shown in figure 1. Flow is generated by driving the lid of the cavity in the x-direction at a constant velocity U = 1. To avoid the spurious effects due to the step changes of the velocity at the two upper corners a sufficiently smooth ve- locity profile that asymptotes to U in the immediate vicin- ity of the corners was imposed. Experimental evidence suggests that the flow in the cavity becomes turbulent at a Reynolds number slightly above R = 5, 000 with the ini- tiation point of transition being at the downstream lower corner [93. Fig. 2 Velocity vectors of mean turbulent flow at Reynolds number R = 50,000 38

We simulate the lid-driven cavity flow at Reynolds num- ber R = SO,OOO with the algebraic model (equation (3~) employed in the computation. As a characteristic length ~ needed for the model at each nodal point we consider the distance from the point to the nearest wall which is proportional to the characteristic size of the largest eddies of turbulence in this flow. In figure 2 we plot the velocity vectors of the mean turbulent flow and in figure 3 several ~ velocity profiles in the neighborhood of the right lower corner of the cavity. It is seen that the flow undergoes sep- aration as we move from the center of the cavity to the corner, as expected. However, the size of the recirculation zone as well as the strength of the reverse flow is weak due to the effects of the increased apparent viscosity. Indeed, direct numerical simulations on the same mesh at Reynolds number R = 10, 000 shows regions of multiple small and large size instantaneous eddies (figure 43; the time-a~reraged flow however exhibits only very small recirculation zones sirrular to the ones shown in figures 2-3. : ; i: Xa0.85 X.O.9O Xe0.95 BOTTOM WALL Fig. 3 Velocity (x—component) profiles at station x= 0.85 (a); 0.90 (b); 0.95 (c) in the neighborhood of the downstream secondary eddy To examine the spatial variation of the eddy-viscosity after a stationary state is reached we plot in figure 5a the total viscosity at a vertical station through the center of the cavity. We see that the distribution is smooth and peaks at the lid where there is enhanced mixing, while it achieves a minimum at the wall (molecular viscosity) and in the region close to the center of the cavity where the fluid is almost stagnant. Thus, the total viscosity increases by almost three orders of magnitude from its molecular to its maximum value at the moving wall. In figure Sb we plot the total viscosity at the same position as in figure 5a but at Reynolds number R = 25, 000. We see in this case that the region close to the bottom wall where the fluid motion is weaker the region of molecular viscosity is broader, whereas the outside distribution is similar to the one in figure 5a. r7 r ,. t t . ~ ~ `_ \~ ~ ~ N`_ `_ `_ . _ __ _ ~ ~ __ if l I: _ .. , ~ Fig. 4 Instantaneous velocity vectors at R = 10,000: The DOW is computed via direct spectral element simulation Lit ~ , ~ As, _~ _~` _ ~ 1 1 1 1 1 1 l J 1 /' K. ~ <.,~,~'A ~~ ~ _. _ J in 2x10-5 0.1283E-O1 (b) / 4xlO-S ~ o.sooo X o. Y l.ooo I/ 1 \/ 0.5000 X o. Y 1.000 Fig. 5 Total viscosity variation at a vertical station through the cavity center (a) R = 50,000 (b) R = 25,000 39

k—~ Modeling of Turbulent Channel Flow We simulate first flow in a plane channel at Reynolds num- ber R* = 1000. Here the Reynolds number is based on the wall-shear velocity U* and the channel half width. The computational domain extends from y+ = 8.5 to y+ = 1991.5, so that the constant shear regions next to the walls are excluded from the simulation. Two-dimensional spec- tral elements are employed for the discretization, (K = 12), of order N = 9, so that there are 97 collocation points across the channel. The boundary conditions are prescribed u, k and ~ at the y-sides of the domain, while periodicity conditions in the streamwise direction (x) en- sure one-dimensionality of the mean flow. The simulation starts from an equilibrium state and integration proceeds for a long time. Results of the mean quantities are pre- sented in figures 6 and 7. As a next test we simulated the same flow but with the computational domain extended to walls. The boundary conditions at the walls are now no-slip for velocity, zero for kinetic energy, and Neumann /: Ox v 8.5258O 5 X 50.0 I .0 Fig. 6 Mean velocity and eddy-nscosity profiles at R* = 1000 for the turbulent Dow in a channel 2.9 K / Q 14 0 ·— ~ ~ '0.0014 50.0 X 50.0 8.5 Y 1991.5 - Fig. 7 Turbulent kinetic energy (k) and dissipation rate (~) profiles at R* = 1000 for the turDtllent flow in a channel condition for the dissipation rate. In order to minirruze discretization errors and ensure accurate imposition of the zero flux condition on ~ we increased the resolution to or- der N = 15 per element. The results of this simulation were identical to results displayed in figures 6 and 7, which suggests that the zero flux condition on ~ is a meaningful one. k—~ Modeling of Turbulent Flow Over a Backward—Facing Step The geometry for this flow and the corresponding spectral element mesh are shown in figure 8; this geometry is identi- cal with the geometry used in experiments of Kim et al [10] and in the computations of Avva et al [11~. The outflow length is taken to be twenty times the step height H. A total of K = 92 elements of order N = 9 were employed in the discretization. The Reynolds number R* = It - ~ is 1870 (or R. = 45, 000~. In this simulation the computational domain includes the walls, where we impose the zero flux condition for the disipation rate e; at the outflow Neumann conditions are specified for all field variables. The inflow conditions match the measured profiles in the experiment of Kim et al [10] and employed in the computations of Awa et al t11~. Finally, the initial data are based on the results reported in t11) and some initial perturbation. Fig. 8 Spectral Element Mesh employed for the simulation of flow over a backward—facing step. In figure 9 we. first present the results of this simulation in the form of streamlines. In addition to the large re- circulation zone shown in this plot eddies of smaller size of opposite rotation appear at the step test section floor juncture consistent with the experimental findings [123. To the best of our knowledge no other simulation has resolved such fine structures previously using a single or a zonal modeling approach. Furthermore, the length of the recir- culation zone (figure 9-10) is computed to be L = 7.3H in close agreement with the experimental value t10~. Most studies todate fall short of predicting the correct value due to errors both in numerics and turbulence models. In fig- ure 10 and 11 we plot all mean variables at a station very close to the reattachment point. The eddy viscosity distri- bution attains its maximum close to that point where the turbulence intensity exhibits its extremum. 40

- - - ~ ~ - ~ ~ ~ ~ - - ~ Fig. 9 Spectral Element—RNG/k—~ Simulation of Dow over a step at R* = 1870. The curves are mean Dow streamlines while color represents turbulent kinetic energy (red=max, :Uue—min). In this calculation, no ad hoc fitted parameters or experimental input is used. The recirculation zone length is 7.3 step heights in agreement with the experiments of Kim, et al (1979~. The small vortices in the step corner are observed experimentally too. -Q 027 I/'/ . /-/ ~.v 13,800 X ~,800 O Y 1500 Fig. 10 Mean velocity and eddy—viscosity profiles at R* = 1870 for the DOW described in Figure 9 30.3:j K j) it/' ~ i / 1/ V/ t3180° o Large-Eddy Simulations of Turbulent Channel Flow In the last example we employ the RNG subgrid-scale model to simulate the turbulent channel flow at R* = 185. This simulation corresponds to identical conditions as the di- rect simulation recently reported by Kim, Main and Moser (1987) t13~. In our simulation we have used however 60 times fewer grid points. In particular, for this case we em- ploy a global spectral discretization based on Fourier ex- pansion in strearnwise and spanwise directions, and Cheby- shev expansion in the inhomogeneous direction (16 x 64 x 643. The initial conditions are based on three-dimensional Tollmien-Schlichting waves. The results of our simulation are essentially identical with the results of the direct simu- lation as shown in figures 12 and 13, and in close agreement with the experiments t14~. The agreement extends also to higher order statistics as well as to flow structure and the streak spacing; the results of our simulations are plotted in figure 14 as color contour plots of the fluctuating ve- locity component at a plane close to the wall (here red indicates low velocity; blue high velocity). The mean sep- aration between streaks is )* ~ 90 t3] in close agreement with the experimentally observed spacing; all previous LES have failed to predict the correct value of streak separation [15~. 0.027 X 13,800 Y 1500 Fig. 11 Turbulent kinetic energy (k) and dissipation rate (~) profiles at the reattachment point (R* = 1870) (DOW as in Figure 9) .' ,~ Y+ _ ~ 10' :~ Y+ Fig. 12 Mean Velocity Profile at R$ = 185 of the turbulent channel DOW computed via RNG - arge eddy simulations —: 5.0 + 2.5Iny+ 41

3~0 2.5 2.0 1.0 0.5 - O - !~\urms o o :W o Yrms - 20 40 60 Y+ Turbulent intensities for the Dow described in Firure 12. Also shown are experimental results of Kreplin and Ackerman (1979) Contours (low velocity: red; high velocity: blue) of near—wall fluctuations x—velocity in a turbulent channel £10w large eddy simulation using the renormalization grou~based subgrid scale eddy viscosity at R* = 185. Observe the alignment of the £10w in the streamwise direction and the formation of low~peed streamwise streaks. 4. CONCLUSIONS Turbulent flows are difficult to compute because they re- quire order R3 work and order R9/4 storage to resolve dy- namically significant velocity fluctuations at large Reynolds number R. These requirements become more severe in com- plex geometry flows where a broader spectrum of significant scales is present. Direct numerical simulation although can provide a very faithful description of the flow is limited severely at the present time and in the next few years by the aforementioned computational requirements to flow simu- lations at a Reynolds number of a few thousands. Progress in simulating realistic flows of engineering importance can be made if reliable and well validated modeling techniques are combined with highly-accurate numerical methods. In the current work, we demonstrated how RNG methodol- ogy con be used to represent all scales of turbulent flow using for example a differential k—~ model, or to repre- sent only the small-scale dynamics using a subgrid-scale model. This methodology is very robust and can be ap- plied to a variety of flows as a totally prognostic tool of analysis, since it requires no apriori known parameters or any experimental input, which is typically the case with the currently used turbulence modeling techniques. A compu- tationally efficient implementation of the RNG methodol- ogy is obtained if it is combined with spectral or spectral element discretization methods, which are used today pri- marily in direct computations of transitional and turbulent flows. We are currently working in further validating our RNG/spectral element methodology in unsteady turbulent flows (e.g. vortex streets) in complex geometries. 5. ACKNOWLEDGEMENTS This work was supported by AFOSR under Contract F49620- 87-C-0036, by ONR under Contract N00014-82-C-0451, and by DARPA under contract N00014-86-K-0759. Most of the computations were performed on the CRAY-X/MP48 at the Pittsburgh Supercomputing Center. References t1] Yakhot V. and Orszag S.A. Renormalization Group analysis of turbulence.I. Basic theory. J. Sc. Comp., 1:3, 1986. t2] Yakhot V. and Orssag S.A. Relation between Kol- mogorov and Batchelor constants. Phys. Fluids, 30:3, 1987. [3] Yakhot A., Orszag S.A., and Yakhot V. Renormaliza- tion Group formulation of large-eddy simulations. J. Sc. Comp., to appear, 1989. t4] Piomelli U., Ferziger J., and Moin P. Models for Large Eddy Simulations of Turbulent channel Bows including transpiration. Technical Report Report TF-32, Ther- mosciences Division, Dept. Mech. Engineering, Stan- ford University, California,94305, 1987. [5] Karniadakis G.E. Spectral element simulations of lam- inar and turbulent flows in complex geometries. Appl. Norm. Math., to appear, 1989. 42

[6] R0nquist E.M. Optimal spectral element methods for the unsteady three-dimensior~al incompressible Navier- Stokes equations. PhD thesis, Massachusetts Institute of Technology, 1988. t7] Patera A.T. A spectral element method for Fluid Dy- namics; Laminar flow in a channel expansion. J. Com- put. Phys., 54:468,1984. t8] Fischer P., Ho L.W., Karniadakis G.E., R0nquist E., and Patera A.T. Recent advances in parallel spectral element simulation of unsteady incompressible flows. Computers arid Structures, 30:217, 1988. t9] Prassad A.K. and Koseff J.R. Reynolds number and end-wall effects on a lid-driven cavity flow. Phys. of Fluids A, 1:208, 1989. t10] Kim J., Kline S.J., and Johnston J.P. Ir~vestigatior' of separation and reattachment of a turbulent shear layer: f ow over a back?vard-facir~g step. Technical Re- port Report MD-37, Dept. Mech. Engineering, Stan- ford University, California,94305, 1978. t11] Avva R.K., Kline S.J., and Fersiger J.H. Computation of turbulent flow over a backward-facing step-zonal approach. In AIAA-88-0611, 1989. t12] Abott D.E. and Kline S.J. Experimental investiga- tion of subsonic turbulent flow over single and double backward facing steps. J. Basic Engineering, 84, 1962. t134 Kim J., Main P., and Moser R. Turbulence statistics in fully developed channel flow at low Reynolds number. J. Fluid Mech., 177:133,1987. t14] Kreplin H. and EckelmAn H. Behavior of the three fluctuating velocity components in the wall region of a turbulent channel flow. Phys. Fluids, 22:1233, 1979. t15] Main P. and Kim J. Numerical investigation of tur- bulent channel flow. J. Fluid Mech., 118:341, 1982. 43

DISCUSSION by E.P. Rood For some applications a requirement is free surface flows near their intersection to predict with boundary layers or in the turbulent wake of a ship. What experimental or other physical information is needed to develop an adequate turbulence model for such free surface flows? Author's Reply In the renormalization group approach, full turbulence closure is obtained without a priori experimental input. This is true for both free surface and other flows. However, in developing turbulence models for free surface flows, it may be best to "tune" the model using experimental observations of free surface turbulence. It is hoped that many of these questions can be answered in the near future with detailed application of the renormalization group models. 44