- Email: [email protected]
Numerical prediction of neutral atmospheric boundary layers (RUSHIL experiment)
Journal of Wind Engineering and Industrial Aerodynamics 74—76 (1998) 263—272
Numerical prediction of neutral atmospheric boundary layers (RUSHIL experiment) S. Duranti!, G. Pezzuto", F. Pittaluga!,* ! IMSE, Institute of Thermal Machines and Energy Systems, University of Genoa, via Montallegro 1, 16145 Genova, Italy " CMIRL, Meteo-Hydrological Center of Liguria Region, via Dodecaneso 33, 16145 Genova, Italy
Abstract In the first part of the paper, code NASTENV is presented with its main characteristics and latest developments: the version here discussed is that of an advanced finite-volume 3D time-dependent Navier—Stokes solver, suitable for thermo-fluid-dynamical analyses of atmospheric, variable density, turbulent flow fields. Rooted on a conservative, arbitrary Lagrangian—Eulerian (ALE) scheme, it was very recently converted from explicit to a largely implicit formulation, and, to further increasing numerical efficiency, a multigrid accelerator has been implemented, thus reaching substantial savings in the number of iterations required to reach convergence. A specific aim of the study is to test the predictive capacity of NASTENV’s new turbulence model by comparing predictions with the wind data obtained from US EPA’s wind tunnel experiment RUSHIL, on the flow structure over isolated two-dimensional hills. To this end, the new turbulence model (TSDIA, two-scale direct-interaction approximation) recently implemented into the code in a simplified first-order version, has been extended for the present application to a higher-order level: more specifically, use is made of asymptotic expansions based on a scale parameter capable to distinguish the slow variation of mean flows from the fast variation of turbulent fluctuations. These contributions greatly help in correctly capturing significant shear stress levels in flow regions where velocity gradients are very low. The level of accuracy of the results is high, confirming the reliability of the predictive capacity attained by the solver. ( 1998 Elsevier Science Ltd. All rights reserved. Keywords: Neutral boundary layers; Rushil experiment; Navier—Stokes predictions
1. Introduction When the dominant mechanism, governing turbulence generation within atmospheric boundary layers, is mostly wind shear, and convective forces are small, the * Corresponding author. E-mail: [email protected]. 0167-6105/98/$19.00 ( 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 1 6 7 - 6 1 0 5 ( 9 8 ) 0 0 0 2 3 - 3
264
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
experimental situation can become a quite important test-case to be used as a reference for validating turbulence models in their aerodynamic, “neutral”, aspects. This is particularly true when such types of data are taken under neatly controlled conditions, such as those afforded by the long test section of a meteorological wind tunnel, equipped with proper floor-roughness elements, inlet turbulence generators and advanced aerodynamic/anemometric instrumentation. Only after a positive validation is performed within this type of context, can any turbulence modeling activity safely be extended toward considering additional mechanisms, e.g. buoyancy. The systematic set of data corresponding to the “RUSHIL experiment”, reported in Refs. [1,2], fully satisfies the above requisites, and has thus been selected as particularly suited for performing the said validation, in spite of some flow features really challenging to be correctly predicted. The numerical Navier—Stokes solver here utilized is named NASTENV, and represents the most recent “atmospheric version” of NAST, a finite-volume ALE (arbitrary Lagrangian—Eulerian) computer code, for 3D time-dependent turbulent-flows predictions, progressively developed at IMSE, University of Genova, during the last ten years [3—8]. A previous version of the code, with rather simple turbulence models [10], was already confronted with the RUSHIL data, yielding generally nice predictions for most of the flow features, other than the shear stress profiles. This outcome prompted adoption of a new turbulence-modeling strategy and a corresponding revision for the wall functions. As later on discussed, no “full Reynolds stress model” was chosen, mainly for the heavily modelized assumptions that are required, so often showing no clear connection between reliability and analytical complexity. Rather, and in line with the several SGS (sub-grid scale) turbulence schemes that, in the past, have characterized the conceptual development of NAST code [4—6], an approach was preferred capable to differentiate, in a statistical sense, the effects produced by the high-frequency turbulent fluctuations from the slow mean-flow variations. In authors’ opinion, the TSDIA (“two-scale direct-interaction approximation”) approach adequately meets the above criteria, whilst retaining a correct balance between complexity and physical interpretation, and yet showing potential for further developments [11—13]. A recent, TSDIA-compatible, internal-flow version of solver NAST was validated in connection with the rear-facing step test-case, yielding an excellent outcome [14]. The RUSHIL test-case displays flow features even more challenging, to be predicted, than the backward-step case, in particular, when attention is devoted to Reynolds stress behavior: this point is quite interesting, and, very peculiarly, its discussion, seldom touched [15], is almost avoided in theoretical—numerical investigations about this subject, no comparisons being usually given in terms of the detailed behavior of shear stress profiles [16,17]. In the present paper, a deeper analysis of this point is performed and an explanation is given of the higher difficulties encountered in predicting, for the RUSHIL case, the turbulent shear stress behavior: in this connection, a higher-order extension of the TSDIA scheme, specifically developed for this situation, is presented, which succeeds in correctly capturing all the relevant turbulent flow characteristics.
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
265
2. Analytical and numerical structure of solver NASTENV In NASTENV every time-step consists of a Lagrangian phase, with the flow velocities applied to the (moving) grid vertexes, followed by a convective phase, i.e. a rezoning of the (updated) flow field, performed by displacing the grid back to its initial position. Partial-donor-cell differencing is applied for convection, with the use of cell-face velocities in order to lower numerical diffusion [3,5,6]. In synthesis, the formerly explicit pressure-updating Lagrangian-phase is here converted into a fully implicit phase; then, the resulting Poisson equation is solved by an iterative method (Gauss—Seidel SOR) accelerated by means of a multigrid technique that utilizes the similarity between the solutions obtained on multiple grids, each with a different refinement level [14]. Grid types are staggered, boundary fitted, structured. The governing equations for turbulent flow are written, and solved, in conservation form as follows: Lo #+ ) (ou)"0, Lt
(1)
A
B
L(ou) 2 #+ ) (ouu)"!+ p# ok #+ ) r, Lt 3
(2)
L (oe)#+ ) (oue)"!p+ ) u#oe, Lt
(3)
where u is velocity vector, p is pressure, o is density, e is internal energy, r is viscous and turbulent (i.e. k-related) stress tensor. Turbulence kinetic energy k and dissipation e are related through the k-transport equation L(ok) 2 #+ ) (ouk)"! ok+ ) u#r: +u#+ ) (k+k)!oe. Lt 3
(4)
No e-transport equation is solved as such: instead, a higher-order TSDIA-compatible approach has been adopted, along the lines that will be introduced shortly after. With regard to the wall-boundary conditions, they are determined by matching to a “rough-wall” logarithmic velocity profile [9,10], to take care of the roughened floor of the EPA wind tunnel. The near-wall solution, being the regime completely rough, is obtained by imposing that the shear-stress at the wall must satisfy a standard momentum balance, applied to the near-wall cell, almost in the same terms as done for the core-flow cells [14]: at the wall, assumed “displaced”, due to the roughness equivalent-size k , a non-zero velocity is applied, namely, that given by the lower limit 4 of validity of the rough-wall logarithmic profile itself (Fig. 20.20, p. 582, in Ref. [9]). In addition, for the near-wall cells, an “eddy” viscosity is used of the same formulation as the core-flow cells: the related characteristic length which comes out in this way, is assumed as the wall-boundary condition for the turbulence characteristic length equation (see later). The corresponding boundary condition for k is likewise peculiar [14]: in the near-wall cells, the “flow” kinetic energy loss, due to wall actions, is
266
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
directly assumed as the turbulent kinetic energy production term, and introduced as such into Eq. (4). Its dissipative character will arise naturally by solving the same Eq. (4), strictly linked with Eq. (3).
3. The two-scale direct interaction approximation TSDIA is a turbulence modeling scheme that, by taking advantage of a statistical analysis, tries to overcome the known deficiencies of the classical k—e approaches, e.g. the limitations deriving from the assumptions of (standard) turbulent-viscosity formulations or of turbulence field isotropy. In this approach, an asymptotic expansion is performed, based on a scale parameter suitable to differentiate the high-frequency turbulent fluctuations from the slow mean-flow variations [11—13]. In Ref. [14], the basic formulation of a simplified, first-order, TSDIA approach is synthetically presented, together with its numerical implementation and its successful predictions, in connection with a rear-facing step test-case. The RUSHIL flow-situation, however, is different, mainly for the vanishing velocity gradients toward the upper part of the flow channel, typical of all external flow fields, but, most peculiarly, not accompanied by proportional decreases in shear stresses. This fact, seldom addressed in the published papers on this topic [15—17], is, in the authors’ opinion, extremely interesting, and points out the weakness of standard eddy-viscosity formulations. To overcome these limits, a higher-order extension to the basic TSDIA approach has been implemented and will be here, very synthetically, introduced. Whilst the corresponding first-order analysis is given in Refs. [12,14], the extended approach is, in detail, developed in Ref. [11]. In TSDIA the Reynolds stresses are defined in the following way:
A
B
2 Lº Lº i# j #q , R "!Su@u@ T"! kd #l i j ij ij e ij 3 Lx Lx j i
(5)
where “ST” means ensemble average; k"0.5Su@u@T is turbulent kinetic energy; u@, i i U are, respectively, the fluctuation and the ensemble mean of velocity; l is isotropic e eddy viscosity, whereas q is the anisotropy-related stress contribution. In a simplified ij formulation, k and l are expressed by assuming an inertial range form (in wave e number dependence) and, in addition, the fluctuating field is assumed as consisting of eddies whose scale is smaller than a characteristic turbulence scale j. The dependence of k upon e and j turns out, to first-order expansion in the turbulence-scale parameter d, of the form given as in Eq. (16) in Ref. [14]. Applying the inertial range form also to second-order (in d-expansion) contributions to k, and solving this (extended) k-dependence, with respect to j in a perturbational way [11], we get, after an adjustment of coefficients given in Ref. [11], in order to satisfy the solenoidal condition [18,19]: Dk De j"c k3@2e~1#c k3@2e~2 !c k5@2e~3 !c k7@2e~3S !c k7@2e~3S j jk je jS1 1 jS2 2 Dt Dt
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
267
with c "1.84, c "4.36, c "2.55, c "0.119, c "0.119, j jk je jS1 jS2 Lº Lº Lº Lº j, S " i i S " i 2 Lx Lx 1 Lx Lx i j j j
(6)
(summations implied).
Notice that the first term in the j-equation shows the well-known dependence, on k and e, as does ¸ , the “turbulence scale” adopted in classical (Reynolds-averaged) 5 two-equation turbulence models. We now assume also for j (i.e. the largest turbulence scale used in the inertial range approximation) the same dependence usually adopted for ¸ , i.e., ¸ "c k3@2e~1, from which the dependence is derived: 5 5 L j"c k3@2e~1. j
(7)
It is interesting to observe that implementing this hypothesis into solver NASTENV can be interpreted, at least formally, as a peculiar SGS (sub-grid scale) turbulence model, with j acting as the “filter” length. Notice, obviously, that same dependence does not mean same numerical coefficients: c Oc . In order that Eq. (7) j L could satisfy Eq. (6), all the remaining terms in the r.h.s. of Eq. (6), i.e. excluding the first term, must yield a zero net sum. This condition, using again Eq. (7), gives the final j-equation
A
Dj j2 Lº Lº j Dk Lº i i# j "!0.2010 #2.22]10~2 Dt k Dt Lx Jk Lxj Lxj i
B
(8)
which is the fundamental equation actually solved, together with Eq. (4), by the turbulence modeling routines of code NASTENV. The meaning, and treatment, of Eq. (8) are almost the same as in its first-order counterpart, given as Eq. (20) of Ref. [14], where more details are given. At this point, the higher-order analysis is applied also to Eq. (5), in order to evaluate the term q of Eq. (5), usually neglected. From Eqs. (74), (75) of Ref. [11], its TSDIA ij formulation is q "&IJ#ZII #ZIII, ij ij ij ij
(9)
where the second-order term &IJ is a summation of contributions similar to the terms ij S and S encountered in Eq. (6), and the third-order term ZIII implies contributions ij 1 2 like L3º /Lx3: at least in the present application, both above terms turn out of i j minor importance, but have been anyhow programmed and introduced. More relevant appears the other second-order term in Eq. (9), for which, Ref. [11] gives the expression
A
B
D Lº Lº i# j . ZII "0.00922 j2 ij D¹ Lx Lx j i
(10)
268
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
The operator D/D¹ is rather peculiar: it should represent a total-derivative-like information, such as produced by observing, from the “slow” reference frame of the mean local velocity, what happens in the flow field at the time-space scales of the energy containing eddies. This trans-information (aimed at capturing the “direct interaction” between the two different scales) requires performing a time—space coordinate transformation, which here is done only on the basis of physical reasoning. In the present case, but also with (likely) wide validity due to its general structure, a successful transformation has turned out to be the following: D L L L L * aR C, " #º , with "0, + , *X +*x i i i D¹ L¹ LX L¹ LX *X j i i i
(11)
*X being the respective side-lengths of the local computational grid cell, and R the i C local mean-flow radius of curvature. In Eq. (11) the assumption is evident that the “direct interaction” takes place between a mean-flow representative spatial scale (given by R ) and the scale of the energy containing eddies, given by j. The near-unity C numerical coefficient a, for only “fine tuning”, has been set here at 1.25.
4. The RUSHIL test-case: Theoretical/experimental comparison As is well known, the RUSHIL experiment was performed in the US EPA meteorological wind tunnel [1,2,15], with the aim of measuring boundary layer development, in a first case, over a flat-floor situation (case H0, to check attainment of an equilibrium profile) and then, for three different geometries of 2D-hills, one with attached flow (case H8), one with separation (case H3), and one on the verge of flow separation (case H5). A neutral turbulent boundary layer is produced, with a logarithmic “equilibrium” profile of about 1 m height, characterized by an outer-stream velocity º "4 m/s and a friction velocity u*"0.19 m/s. The equivalent surface = roughness size is k "0.178 mm. 4 A previous version of code NAST, with turbulence models based on simple prescriptions for the dissipation length-scale, was already applied to the RUSHIL experiment [10]: the results, in terms of velocity and turbulent energy k profiles, turned out, except for some (local) near-wall k-overpredictions, generally good, but not for the shear stress profiles. In fact, a major problem was, for all cases, the shear stress distinct under-prediction, invariably appearing in the flow regions well above the hills: the reason being that, since the velocity gradients fall, there, to very low levels, no shear stress formulation based on the classical eddy-viscosity concept could yield correct predictions. This point, extremely important, is seldom addressed in the published papers on the subject [16,17]; in Refs. [18,19] such a topic is discussed, attributing it to a “rapid distorsion effect”, operating on the turbulent eddies by the upstream history of the mean flow, whose effects remain active even in absence of velocity gradients. In any case, no publication has been found where a full theoretical prediction of shear stresses is presented, taking into careful account this fundamental issue.
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
269
As already said, in the previous numerical investigation [10], the RUSHIL cases H3 and H8 were considered: here, the case H5 is analyzed, characterized by a still attached downhill flow, but right on the verge of separation. For the above discussion, one of the crucial points motivating the present study is common to all above cases, i.e. the presence of significant Reynolds stress levels, in extended flow regions above the hills, for locally vanishing velocity gradients. A further issue is important: i.e. the trend to overpredict k in some near-wall situations. The answer to the first above problem was, within a previous and unsuccessful strategy, sought mainly by attributing it to possible three-dimensional effects (i.e. longitudinal vortexes): however, a fully 3D analysis gave no real difference [20]. It was at this point that, by taking advantage of the very recent, and positive, outcome obtained by a TSDIA-compatible version of code NAST [14], this evolution to the turbulence model was implemented into code NASTENV: although first-order analysis already yielded some improvements, only the higher-order extension of the scheme, specifically developed for this study along the lines discussed above, succeeded in capturing (and thus, explaining) the trends shown by the Reynolds stresses. In addition, the adoption of a non-zero velocity at the “displaced” rough-wall boundary (as already discussed) greatly improved the k behavior near the wall. All the numerical results now being presented were obtained on a 200 MHz Pentium PC (32 MB ram; 1 GB hd). Typical cpu-times were 7000 s, for a 32]76 grid and about 1500 iterations (implicit scheme, multigrid-accelerated), such as required to reach steady state. Fig. 1 shows the tunnel stations’ locations where measurements were performed (see Ref. [2] for the exact dimensions), whilst in Fig. 2 the comparisons are presented between the experimental (horizontal) mean-velocity profiles and those predicted by code NASTENV. The agreement, for all the stations, both the up-hill and the critical down-hill ones, is excellent, even for the near-wall locations, in spite of a not-refined grid. This latter reason, on the other hand, leads to the straight-line profile appearing very near the wall, simply because no other calculated point is available: however, the experimental trend is always correctly captured, even when a strong inflection point is present (station 10), a clear symptom of a local danger of separation. It is clear that for a more rigorous capture of the very-near-wall velocities a locally more refined grid should be needed. Fig. 3 presents the corresponding turbulent kinetic energy profiles, which allow a deeper interpretation of the flow processes and also of their numerical counterparts.
Fig. 1. Descriptive locations of measurement stations along EPA tunnel.
270
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
Fig. 2. Experimental and predicted horizontal mean-velocity profiles.
Fig. 3. Experimental and predicted turbulent kinetic energy profiles.
First of all, it must be remarked that the overall prediction, all the way down to the wall, is correct, both in values and in local trends. Notice that this result, which is free of usual near-wall over-predictions, is due to the above discussed adoption of an accurate rough-wall condition. In addition, even when some local k-values are not exactly predicted, their maxima and minima are always well localized, i.e. the convection of k is extremely well captured. Some near-wall under-predictions are shown for the 10, 11 stations: although a more refined grid could improve predictions, probably the assumption of a near-wall logarithmic profile, in particular if “strong” inflection points are present, would yield some intrinsic inaccuracies. The behavior of Reynolds stresses is presented in Fig. 4, where the predictions are shown for both the first-order and the higher-order extension of the TSDIA turbulence model. As previously anticipated, the TSDIA higher-order contributions succeed in improving the capture of the stress levels measured experimentally in the flow channel above the hill, where velocity gradients are very low, particularly, in correspondence
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
271
Fig. 4. Experimental and predicted turbulent shear stress profiles.
of stations 7 and 8. Similarly to k-behavior, also for shear stresses it can be seen that, even when their levels are not exactly met, their trends are always well captured. A significant under-prediction of stress levels is visible for the up-hill station 6. This fact should be seen in a wider picture, i.e. that given by the general observation that predictions already improve at hill-top station 7, become excellent at station 8 and show some local over-predictions for the down-hill station 10. In spite of the contribution from Eq. (11), the overall behavior can thus be interpreted in the light of a still too-quick adjustment of the stress predictions to the locally prevailing situation, whilst a higher level of “memory” for the upstream stress history should be active. Thus it seems likely that a slightly “upwind-biased” interpretation of the same transformation (11) should further improve stress predictions. Notice, finally, the highly accurate (and not so frequent) capture of near-wall shear stress levels, surely due to the more physical rough-wall boundary condition, with respect to a smoothwall assumption. 5. Conclusions The presented theoretical results, backed by their positive experimental comparisons, give validity to the long-term research strategy upon which they are rooted. Main qualifying points of this strategy are, in the authors’ opinion, the primary attention (conceptual and numerical) given to turbulence characteristic length-scales’ interactions more than, e.g., to an “averaged” transport of dissipation rate, and the statistical analysis applied toward their modelization. The solver NASTENV, now successfully validated, is thus becoming a powerful tool for predicting more complex atmospheric situations. In this connection, next step, which is already under way, will be to test its predictive capability within a fully-3D sub-regional territorial situation. To this end, introduction of the gravitational terms into the governing equations has already been performed, whilst a GIS-based, graphical description interface, inclusive of grid generation, is under study. Particular care is being devoted to the treatment of the fully open atmospheric boundaries.
272
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
References [1] L.H. Khurshudyan, W.H. Snyder, I.V. Nekrasov, Flow and dispersion of pollutants over twodimensional hills, US EPA Report N.EPA-600/4-81-067, 1981. [2] F. Trombetti, P. Martano, F. Tampieri, Data sets for studies of flow and dispersion in complex terrain, Rapporto Tecnico n.1, FISBAT-RT-91-1, CNR, Bologna, 1991. [3] A. Amsden et al., Reports LA-10245-MS and LA-10534-MS; Los Alamos National Laboratory, USA, 1985. [4] P. Barbucci, M. Gremigni, F. Pittaluga, An advanced algorithm for the solution of fully 3D Navier—Stokes equations, in: F. Angeli (Ed.), Supercomputing Tools for Science and Engineering, CNR, Milan, 1990. [5] F. Pittaluga, Un criterio Lagrangiano-Euleriano di soluzione delle equazioni di Navier—Stokes per flussi turbolenti comprimibili, in: 47th ATI Congr., Ass. Termotecnica It. Proc., Parma, Vol. II, 1993. [6] F. Pittaluga, NAST: An advanced 3D compressible Navier—Stokes solver, Proc. Workshop of COST Action F1, EPFL, Lausanne, 1993. [7] C. Leoncini, F. Pittaluga, Turb. Flow analysis within exp. models of boiler comb. chambers, Science and Supercomputing at CINECA, 1995 Report, ISBN 88-86037-02-3, 1995, pp. 421—424. [8] R. Luccoli, S. Traverso, F. Pittaluga, F. Beltrami, M. Mennucci, ALE 3D Navier—Stokes Simul. of Comb. Proc. In a Gas Turb. Combustor, Proc. Workshop Turbomachinery ’96, IMSE, Univ. of Genoa, 1996. [9] H. Schlichting, Boundary Layer Theory, McGraw-Hill, New York, 1968, p. 582. [10] S. Cevasco, G. Pezzuto, F. Pittaluga, Analysis of Dissip. Length in Neutral Layers Over 2D Hills, Proc. Euromech Colloq. 338 on Atm. Turbul. and Disp. in Complex Terrain; Ist. FISBAT, Bologna, 1995. [11] A. Yoshizawa, Statistical analysis of the deviation of the Reynolds stress from its eddy—viscosity representation, Phys. Fluids 27 (6) (1984). [12] A. Yoshizawa, Statistical modeling of a transport equation for the kinetic energy dissipation rate, Phys. Fluids 30 (3) (1987). [13] A. Yoshizawa, S. Nisizima, A non-equilibrium representation of the turbulent viscosity based on a two-scale turbulence theory, Phys. Fluids A 5 (12) (1993). [14] S. Duranti, F. Pittaluga, Improvement of Navier—Stokes ALE predictions by adoption of a TSDIA turb. model, Proc. Workshop COST Action F1 on 3D Navier-Stokes Modeling, Courchevel, France, 1997. [15] S. Di Sabatino, Analisi e modellistica di flussi turbolenti non-omogenei, con riferimento allo strato limite atmosferico, Graduation Thesis in Physics, Univ. of Bologna, 1995. [16] D. Pierrat, D. Delaunay, Num. Simul. of the Flow Past a 2D Model Hill, Proc. Euromech Colloquium 338 on Atmospheric Turbulence and Dispersion in Complex Terrain; Ist. FISBAT, Bologna, 1995. [17] S. Finardi et al., Boundary Layer Flow over Analyt. 2D Hills, Bound. Layer Meteorol. 63 (1993). [18] M. Okamoto, Theoretical Investigation of an eddy-viscosity-type representation of the Reynolds stress, J. Phys. Soc. Japan 63 (1994) 2102. [19] A. Yoshizawa, 1996, private communication. [20] F. Pittaluga, S. Duranti, Effetti tridimensionali in flussi di tipo atmosferico su colline di geometria 2D, Report Interno IMSE 3/97, Ist. Macchine e Sist. Energetici, Univ. Genova, 1997.
Numerical prediction of neutral atmospheric boundary layers (RUSHIL experiment) S. Duranti!, G. Pezzuto", F. Pittaluga!,* ! IMSE, Institute of Thermal Machines and Energy Systems, University of Genoa, via Montallegro 1, 16145 Genova, Italy " CMIRL, Meteo-Hydrological Center of Liguria Region, via Dodecaneso 33, 16145 Genova, Italy
Abstract In the first part of the paper, code NASTENV is presented with its main characteristics and latest developments: the version here discussed is that of an advanced finite-volume 3D time-dependent Navier—Stokes solver, suitable for thermo-fluid-dynamical analyses of atmospheric, variable density, turbulent flow fields. Rooted on a conservative, arbitrary Lagrangian—Eulerian (ALE) scheme, it was very recently converted from explicit to a largely implicit formulation, and, to further increasing numerical efficiency, a multigrid accelerator has been implemented, thus reaching substantial savings in the number of iterations required to reach convergence. A specific aim of the study is to test the predictive capacity of NASTENV’s new turbulence model by comparing predictions with the wind data obtained from US EPA’s wind tunnel experiment RUSHIL, on the flow structure over isolated two-dimensional hills. To this end, the new turbulence model (TSDIA, two-scale direct-interaction approximation) recently implemented into the code in a simplified first-order version, has been extended for the present application to a higher-order level: more specifically, use is made of asymptotic expansions based on a scale parameter capable to distinguish the slow variation of mean flows from the fast variation of turbulent fluctuations. These contributions greatly help in correctly capturing significant shear stress levels in flow regions where velocity gradients are very low. The level of accuracy of the results is high, confirming the reliability of the predictive capacity attained by the solver. ( 1998 Elsevier Science Ltd. All rights reserved. Keywords: Neutral boundary layers; Rushil experiment; Navier—Stokes predictions
1. Introduction When the dominant mechanism, governing turbulence generation within atmospheric boundary layers, is mostly wind shear, and convective forces are small, the * Corresponding author. E-mail: [email protected]. 0167-6105/98/$19.00 ( 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 1 6 7 - 6 1 0 5 ( 9 8 ) 0 0 0 2 3 - 3
264
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
experimental situation can become a quite important test-case to be used as a reference for validating turbulence models in their aerodynamic, “neutral”, aspects. This is particularly true when such types of data are taken under neatly controlled conditions, such as those afforded by the long test section of a meteorological wind tunnel, equipped with proper floor-roughness elements, inlet turbulence generators and advanced aerodynamic/anemometric instrumentation. Only after a positive validation is performed within this type of context, can any turbulence modeling activity safely be extended toward considering additional mechanisms, e.g. buoyancy. The systematic set of data corresponding to the “RUSHIL experiment”, reported in Refs. [1,2], fully satisfies the above requisites, and has thus been selected as particularly suited for performing the said validation, in spite of some flow features really challenging to be correctly predicted. The numerical Navier—Stokes solver here utilized is named NASTENV, and represents the most recent “atmospheric version” of NAST, a finite-volume ALE (arbitrary Lagrangian—Eulerian) computer code, for 3D time-dependent turbulent-flows predictions, progressively developed at IMSE, University of Genova, during the last ten years [3—8]. A previous version of the code, with rather simple turbulence models [10], was already confronted with the RUSHIL data, yielding generally nice predictions for most of the flow features, other than the shear stress profiles. This outcome prompted adoption of a new turbulence-modeling strategy and a corresponding revision for the wall functions. As later on discussed, no “full Reynolds stress model” was chosen, mainly for the heavily modelized assumptions that are required, so often showing no clear connection between reliability and analytical complexity. Rather, and in line with the several SGS (sub-grid scale) turbulence schemes that, in the past, have characterized the conceptual development of NAST code [4—6], an approach was preferred capable to differentiate, in a statistical sense, the effects produced by the high-frequency turbulent fluctuations from the slow mean-flow variations. In authors’ opinion, the TSDIA (“two-scale direct-interaction approximation”) approach adequately meets the above criteria, whilst retaining a correct balance between complexity and physical interpretation, and yet showing potential for further developments [11—13]. A recent, TSDIA-compatible, internal-flow version of solver NAST was validated in connection with the rear-facing step test-case, yielding an excellent outcome [14]. The RUSHIL test-case displays flow features even more challenging, to be predicted, than the backward-step case, in particular, when attention is devoted to Reynolds stress behavior: this point is quite interesting, and, very peculiarly, its discussion, seldom touched [15], is almost avoided in theoretical—numerical investigations about this subject, no comparisons being usually given in terms of the detailed behavior of shear stress profiles [16,17]. In the present paper, a deeper analysis of this point is performed and an explanation is given of the higher difficulties encountered in predicting, for the RUSHIL case, the turbulent shear stress behavior: in this connection, a higher-order extension of the TSDIA scheme, specifically developed for this situation, is presented, which succeeds in correctly capturing all the relevant turbulent flow characteristics.
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
265
2. Analytical and numerical structure of solver NASTENV In NASTENV every time-step consists of a Lagrangian phase, with the flow velocities applied to the (moving) grid vertexes, followed by a convective phase, i.e. a rezoning of the (updated) flow field, performed by displacing the grid back to its initial position. Partial-donor-cell differencing is applied for convection, with the use of cell-face velocities in order to lower numerical diffusion [3,5,6]. In synthesis, the formerly explicit pressure-updating Lagrangian-phase is here converted into a fully implicit phase; then, the resulting Poisson equation is solved by an iterative method (Gauss—Seidel SOR) accelerated by means of a multigrid technique that utilizes the similarity between the solutions obtained on multiple grids, each with a different refinement level [14]. Grid types are staggered, boundary fitted, structured. The governing equations for turbulent flow are written, and solved, in conservation form as follows: Lo #+ ) (ou)"0, Lt
(1)
A
B
L(ou) 2 #+ ) (ouu)"!+ p# ok #+ ) r, Lt 3
(2)
L (oe)#+ ) (oue)"!p+ ) u#oe, Lt
(3)
where u is velocity vector, p is pressure, o is density, e is internal energy, r is viscous and turbulent (i.e. k-related) stress tensor. Turbulence kinetic energy k and dissipation e are related through the k-transport equation L(ok) 2 #+ ) (ouk)"! ok+ ) u#r: +u#+ ) (k+k)!oe. Lt 3
(4)
No e-transport equation is solved as such: instead, a higher-order TSDIA-compatible approach has been adopted, along the lines that will be introduced shortly after. With regard to the wall-boundary conditions, they are determined by matching to a “rough-wall” logarithmic velocity profile [9,10], to take care of the roughened floor of the EPA wind tunnel. The near-wall solution, being the regime completely rough, is obtained by imposing that the shear-stress at the wall must satisfy a standard momentum balance, applied to the near-wall cell, almost in the same terms as done for the core-flow cells [14]: at the wall, assumed “displaced”, due to the roughness equivalent-size k , a non-zero velocity is applied, namely, that given by the lower limit 4 of validity of the rough-wall logarithmic profile itself (Fig. 20.20, p. 582, in Ref. [9]). In addition, for the near-wall cells, an “eddy” viscosity is used of the same formulation as the core-flow cells: the related characteristic length which comes out in this way, is assumed as the wall-boundary condition for the turbulence characteristic length equation (see later). The corresponding boundary condition for k is likewise peculiar [14]: in the near-wall cells, the “flow” kinetic energy loss, due to wall actions, is
266
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
directly assumed as the turbulent kinetic energy production term, and introduced as such into Eq. (4). Its dissipative character will arise naturally by solving the same Eq. (4), strictly linked with Eq. (3).
3. The two-scale direct interaction approximation TSDIA is a turbulence modeling scheme that, by taking advantage of a statistical analysis, tries to overcome the known deficiencies of the classical k—e approaches, e.g. the limitations deriving from the assumptions of (standard) turbulent-viscosity formulations or of turbulence field isotropy. In this approach, an asymptotic expansion is performed, based on a scale parameter suitable to differentiate the high-frequency turbulent fluctuations from the slow mean-flow variations [11—13]. In Ref. [14], the basic formulation of a simplified, first-order, TSDIA approach is synthetically presented, together with its numerical implementation and its successful predictions, in connection with a rear-facing step test-case. The RUSHIL flow-situation, however, is different, mainly for the vanishing velocity gradients toward the upper part of the flow channel, typical of all external flow fields, but, most peculiarly, not accompanied by proportional decreases in shear stresses. This fact, seldom addressed in the published papers on this topic [15—17], is, in the authors’ opinion, extremely interesting, and points out the weakness of standard eddy-viscosity formulations. To overcome these limits, a higher-order extension to the basic TSDIA approach has been implemented and will be here, very synthetically, introduced. Whilst the corresponding first-order analysis is given in Refs. [12,14], the extended approach is, in detail, developed in Ref. [11]. In TSDIA the Reynolds stresses are defined in the following way:
A
B
2 Lº Lº i# j #q , R "!Su@u@ T"! kd #l i j ij ij e ij 3 Lx Lx j i
(5)
where “ST” means ensemble average; k"0.5Su@u@T is turbulent kinetic energy; u@, i i U are, respectively, the fluctuation and the ensemble mean of velocity; l is isotropic e eddy viscosity, whereas q is the anisotropy-related stress contribution. In a simplified ij formulation, k and l are expressed by assuming an inertial range form (in wave e number dependence) and, in addition, the fluctuating field is assumed as consisting of eddies whose scale is smaller than a characteristic turbulence scale j. The dependence of k upon e and j turns out, to first-order expansion in the turbulence-scale parameter d, of the form given as in Eq. (16) in Ref. [14]. Applying the inertial range form also to second-order (in d-expansion) contributions to k, and solving this (extended) k-dependence, with respect to j in a perturbational way [11], we get, after an adjustment of coefficients given in Ref. [11], in order to satisfy the solenoidal condition [18,19]: Dk De j"c k3@2e~1#c k3@2e~2 !c k5@2e~3 !c k7@2e~3S !c k7@2e~3S j jk je jS1 1 jS2 2 Dt Dt
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
267
with c "1.84, c "4.36, c "2.55, c "0.119, c "0.119, j jk je jS1 jS2 Lº Lº Lº Lº j, S " i i S " i 2 Lx Lx 1 Lx Lx i j j j
(6)
(summations implied).
Notice that the first term in the j-equation shows the well-known dependence, on k and e, as does ¸ , the “turbulence scale” adopted in classical (Reynolds-averaged) 5 two-equation turbulence models. We now assume also for j (i.e. the largest turbulence scale used in the inertial range approximation) the same dependence usually adopted for ¸ , i.e., ¸ "c k3@2e~1, from which the dependence is derived: 5 5 L j"c k3@2e~1. j
(7)
It is interesting to observe that implementing this hypothesis into solver NASTENV can be interpreted, at least formally, as a peculiar SGS (sub-grid scale) turbulence model, with j acting as the “filter” length. Notice, obviously, that same dependence does not mean same numerical coefficients: c Oc . In order that Eq. (7) j L could satisfy Eq. (6), all the remaining terms in the r.h.s. of Eq. (6), i.e. excluding the first term, must yield a zero net sum. This condition, using again Eq. (7), gives the final j-equation
A
Dj j2 Lº Lº j Dk Lº i i# j "!0.2010 #2.22]10~2 Dt k Dt Lx Jk Lxj Lxj i
B
(8)
which is the fundamental equation actually solved, together with Eq. (4), by the turbulence modeling routines of code NASTENV. The meaning, and treatment, of Eq. (8) are almost the same as in its first-order counterpart, given as Eq. (20) of Ref. [14], where more details are given. At this point, the higher-order analysis is applied also to Eq. (5), in order to evaluate the term q of Eq. (5), usually neglected. From Eqs. (74), (75) of Ref. [11], its TSDIA ij formulation is q "&IJ#ZII #ZIII, ij ij ij ij
(9)
where the second-order term &IJ is a summation of contributions similar to the terms ij S and S encountered in Eq. (6), and the third-order term ZIII implies contributions ij 1 2 like L3º /Lx3: at least in the present application, both above terms turn out of i j minor importance, but have been anyhow programmed and introduced. More relevant appears the other second-order term in Eq. (9), for which, Ref. [11] gives the expression
A
B
D Lº Lº i# j . ZII "0.00922 j2 ij D¹ Lx Lx j i
(10)
268
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
The operator D/D¹ is rather peculiar: it should represent a total-derivative-like information, such as produced by observing, from the “slow” reference frame of the mean local velocity, what happens in the flow field at the time-space scales of the energy containing eddies. This trans-information (aimed at capturing the “direct interaction” between the two different scales) requires performing a time—space coordinate transformation, which here is done only on the basis of physical reasoning. In the present case, but also with (likely) wide validity due to its general structure, a successful transformation has turned out to be the following: D L L L L * aR C, " #º , with "0, + , *X +*x i i i D¹ L¹ LX L¹ LX *X j i i i
(11)
*X being the respective side-lengths of the local computational grid cell, and R the i C local mean-flow radius of curvature. In Eq. (11) the assumption is evident that the “direct interaction” takes place between a mean-flow representative spatial scale (given by R ) and the scale of the energy containing eddies, given by j. The near-unity C numerical coefficient a, for only “fine tuning”, has been set here at 1.25.
4. The RUSHIL test-case: Theoretical/experimental comparison As is well known, the RUSHIL experiment was performed in the US EPA meteorological wind tunnel [1,2,15], with the aim of measuring boundary layer development, in a first case, over a flat-floor situation (case H0, to check attainment of an equilibrium profile) and then, for three different geometries of 2D-hills, one with attached flow (case H8), one with separation (case H3), and one on the verge of flow separation (case H5). A neutral turbulent boundary layer is produced, with a logarithmic “equilibrium” profile of about 1 m height, characterized by an outer-stream velocity º "4 m/s and a friction velocity u*"0.19 m/s. The equivalent surface = roughness size is k "0.178 mm. 4 A previous version of code NAST, with turbulence models based on simple prescriptions for the dissipation length-scale, was already applied to the RUSHIL experiment [10]: the results, in terms of velocity and turbulent energy k profiles, turned out, except for some (local) near-wall k-overpredictions, generally good, but not for the shear stress profiles. In fact, a major problem was, for all cases, the shear stress distinct under-prediction, invariably appearing in the flow regions well above the hills: the reason being that, since the velocity gradients fall, there, to very low levels, no shear stress formulation based on the classical eddy-viscosity concept could yield correct predictions. This point, extremely important, is seldom addressed in the published papers on the subject [16,17]; in Refs. [18,19] such a topic is discussed, attributing it to a “rapid distorsion effect”, operating on the turbulent eddies by the upstream history of the mean flow, whose effects remain active even in absence of velocity gradients. In any case, no publication has been found where a full theoretical prediction of shear stresses is presented, taking into careful account this fundamental issue.
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
269
As already said, in the previous numerical investigation [10], the RUSHIL cases H3 and H8 were considered: here, the case H5 is analyzed, characterized by a still attached downhill flow, but right on the verge of separation. For the above discussion, one of the crucial points motivating the present study is common to all above cases, i.e. the presence of significant Reynolds stress levels, in extended flow regions above the hills, for locally vanishing velocity gradients. A further issue is important: i.e. the trend to overpredict k in some near-wall situations. The answer to the first above problem was, within a previous and unsuccessful strategy, sought mainly by attributing it to possible three-dimensional effects (i.e. longitudinal vortexes): however, a fully 3D analysis gave no real difference [20]. It was at this point that, by taking advantage of the very recent, and positive, outcome obtained by a TSDIA-compatible version of code NAST [14], this evolution to the turbulence model was implemented into code NASTENV: although first-order analysis already yielded some improvements, only the higher-order extension of the scheme, specifically developed for this study along the lines discussed above, succeeded in capturing (and thus, explaining) the trends shown by the Reynolds stresses. In addition, the adoption of a non-zero velocity at the “displaced” rough-wall boundary (as already discussed) greatly improved the k behavior near the wall. All the numerical results now being presented were obtained on a 200 MHz Pentium PC (32 MB ram; 1 GB hd). Typical cpu-times were 7000 s, for a 32]76 grid and about 1500 iterations (implicit scheme, multigrid-accelerated), such as required to reach steady state. Fig. 1 shows the tunnel stations’ locations where measurements were performed (see Ref. [2] for the exact dimensions), whilst in Fig. 2 the comparisons are presented between the experimental (horizontal) mean-velocity profiles and those predicted by code NASTENV. The agreement, for all the stations, both the up-hill and the critical down-hill ones, is excellent, even for the near-wall locations, in spite of a not-refined grid. This latter reason, on the other hand, leads to the straight-line profile appearing very near the wall, simply because no other calculated point is available: however, the experimental trend is always correctly captured, even when a strong inflection point is present (station 10), a clear symptom of a local danger of separation. It is clear that for a more rigorous capture of the very-near-wall velocities a locally more refined grid should be needed. Fig. 3 presents the corresponding turbulent kinetic energy profiles, which allow a deeper interpretation of the flow processes and also of their numerical counterparts.
Fig. 1. Descriptive locations of measurement stations along EPA tunnel.
270
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
Fig. 2. Experimental and predicted horizontal mean-velocity profiles.
Fig. 3. Experimental and predicted turbulent kinetic energy profiles.
First of all, it must be remarked that the overall prediction, all the way down to the wall, is correct, both in values and in local trends. Notice that this result, which is free of usual near-wall over-predictions, is due to the above discussed adoption of an accurate rough-wall condition. In addition, even when some local k-values are not exactly predicted, their maxima and minima are always well localized, i.e. the convection of k is extremely well captured. Some near-wall under-predictions are shown for the 10, 11 stations: although a more refined grid could improve predictions, probably the assumption of a near-wall logarithmic profile, in particular if “strong” inflection points are present, would yield some intrinsic inaccuracies. The behavior of Reynolds stresses is presented in Fig. 4, where the predictions are shown for both the first-order and the higher-order extension of the TSDIA turbulence model. As previously anticipated, the TSDIA higher-order contributions succeed in improving the capture of the stress levels measured experimentally in the flow channel above the hill, where velocity gradients are very low, particularly, in correspondence
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
271
Fig. 4. Experimental and predicted turbulent shear stress profiles.
of stations 7 and 8. Similarly to k-behavior, also for shear stresses it can be seen that, even when their levels are not exactly met, their trends are always well captured. A significant under-prediction of stress levels is visible for the up-hill station 6. This fact should be seen in a wider picture, i.e. that given by the general observation that predictions already improve at hill-top station 7, become excellent at station 8 and show some local over-predictions for the down-hill station 10. In spite of the contribution from Eq. (11), the overall behavior can thus be interpreted in the light of a still too-quick adjustment of the stress predictions to the locally prevailing situation, whilst a higher level of “memory” for the upstream stress history should be active. Thus it seems likely that a slightly “upwind-biased” interpretation of the same transformation (11) should further improve stress predictions. Notice, finally, the highly accurate (and not so frequent) capture of near-wall shear stress levels, surely due to the more physical rough-wall boundary condition, with respect to a smoothwall assumption. 5. Conclusions The presented theoretical results, backed by their positive experimental comparisons, give validity to the long-term research strategy upon which they are rooted. Main qualifying points of this strategy are, in the authors’ opinion, the primary attention (conceptual and numerical) given to turbulence characteristic length-scales’ interactions more than, e.g., to an “averaged” transport of dissipation rate, and the statistical analysis applied toward their modelization. The solver NASTENV, now successfully validated, is thus becoming a powerful tool for predicting more complex atmospheric situations. In this connection, next step, which is already under way, will be to test its predictive capability within a fully-3D sub-regional territorial situation. To this end, introduction of the gravitational terms into the governing equations has already been performed, whilst a GIS-based, graphical description interface, inclusive of grid generation, is under study. Particular care is being devoted to the treatment of the fully open atmospheric boundaries.
272
S. Duranti et al./J. Wind Eng. Ind. Aerodyn. 74–76 (1998) 263–272
References [1] L.H. Khurshudyan, W.H. Snyder, I.V. Nekrasov, Flow and dispersion of pollutants over twodimensional hills, US EPA Report N.EPA-600/4-81-067, 1981. [2] F. Trombetti, P. Martano, F. Tampieri, Data sets for studies of flow and dispersion in complex terrain, Rapporto Tecnico n.1, FISBAT-RT-91-1, CNR, Bologna, 1991. [3] A. Amsden et al., Reports LA-10245-MS and LA-10534-MS; Los Alamos National Laboratory, USA, 1985. [4] P. Barbucci, M. Gremigni, F. Pittaluga, An advanced algorithm for the solution of fully 3D Navier—Stokes equations, in: F. Angeli (Ed.), Supercomputing Tools for Science and Engineering, CNR, Milan, 1990. [5] F. Pittaluga, Un criterio Lagrangiano-Euleriano di soluzione delle equazioni di Navier—Stokes per flussi turbolenti comprimibili, in: 47th ATI Congr., Ass. Termotecnica It. Proc., Parma, Vol. II, 1993. [6] F. Pittaluga, NAST: An advanced 3D compressible Navier—Stokes solver, Proc. Workshop of COST Action F1, EPFL, Lausanne, 1993. [7] C. Leoncini, F. Pittaluga, Turb. Flow analysis within exp. models of boiler comb. chambers, Science and Supercomputing at CINECA, 1995 Report, ISBN 88-86037-02-3, 1995, pp. 421—424. [8] R. Luccoli, S. Traverso, F. Pittaluga, F. Beltrami, M. Mennucci, ALE 3D Navier—Stokes Simul. of Comb. Proc. In a Gas Turb. Combustor, Proc. Workshop Turbomachinery ’96, IMSE, Univ. of Genoa, 1996. [9] H. Schlichting, Boundary Layer Theory, McGraw-Hill, New York, 1968, p. 582. [10] S. Cevasco, G. Pezzuto, F. Pittaluga, Analysis of Dissip. Length in Neutral Layers Over 2D Hills, Proc. Euromech Colloq. 338 on Atm. Turbul. and Disp. in Complex Terrain; Ist. FISBAT, Bologna, 1995. [11] A. Yoshizawa, Statistical analysis of the deviation of the Reynolds stress from its eddy—viscosity representation, Phys. Fluids 27 (6) (1984). [12] A. Yoshizawa, Statistical modeling of a transport equation for the kinetic energy dissipation rate, Phys. Fluids 30 (3) (1987). [13] A. Yoshizawa, S. Nisizima, A non-equilibrium representation of the turbulent viscosity based on a two-scale turbulence theory, Phys. Fluids A 5 (12) (1993). [14] S. Duranti, F. Pittaluga, Improvement of Navier—Stokes ALE predictions by adoption of a TSDIA turb. model, Proc. Workshop COST Action F1 on 3D Navier-Stokes Modeling, Courchevel, France, 1997. [15] S. Di Sabatino, Analisi e modellistica di flussi turbolenti non-omogenei, con riferimento allo strato limite atmosferico, Graduation Thesis in Physics, Univ. of Bologna, 1995. [16] D. Pierrat, D. Delaunay, Num. Simul. of the Flow Past a 2D Model Hill, Proc. Euromech Colloquium 338 on Atmospheric Turbulence and Dispersion in Complex Terrain; Ist. FISBAT, Bologna, 1995. [17] S. Finardi et al., Boundary Layer Flow over Analyt. 2D Hills, Bound. Layer Meteorol. 63 (1993). [18] M. Okamoto, Theoretical Investigation of an eddy-viscosity-type representation of the Reynolds stress, J. Phys. Soc. Japan 63 (1994) 2102. [19] A. Yoshizawa, 1996, private communication. [20] F. Pittaluga, S. Duranti, Effetti tridimensionali in flussi di tipo atmosferico su colline di geometria 2D, Report Interno IMSE 3/97, Ist. Macchine e Sist. Energetici, Univ. Genova, 1997.