Application of a non isotropic turbulence model to stable atmospheric flows and dispersion over 3D topography

Abstract

A non isotropic turbulence model is extended and applied to three dimensional stably stratified flows and dispersion calculations. The model is derived from the algebraic stress model (including wall proximity effects), but it retains the simplicity of the "eddy viscosity" concept of first order models. The "modified k-epsilon" is implemented in a three dimensional numerical code. Once the flow is resolved, the predicted velocity and turbulence fields are interpolated into a second grid and used to solve the concentration equation. To evaluate the model, various steady state numerical solutions are compared with small scale dispersion experiments which were conducted at the wind tunnel of Mitsubishi Heavy Industries, in Japan. Stably stratified flows and plume dispersion over three distinct idealized complex topographies (flat and hilly terrain) are studied. Vertical profiles of velocity and pollutant concentration are shown and discussed. Also, comparisons are made against the results obtained with the standard k-epsilon model.

Atmospheric dispersion; flow over hills; anisotropic k-epsilon; numerical simulation

Application of a Non Isotropic Turbulence Model to Stable Atmospheric Flows and Dispersion over 3D Topography

Fernando T. Boçon

Departamento de Engenharia Mecânica - Universidade Federal do Paraná - C.P. 19011 / 81531-990 Curitiba, PR, Brazil

[e-mail] Clóvis R. Maliska

Departamento de Engenharia Mecânica - Universidade Federal de Santa Catarina - CP 472 / 88040-000 Florianópolis - SC

[e-mail]

A non isotropic turbulence model is extended and applied to three dimensional stably stratified flows and dispersion calculations. The model is derived from the algebraic stress model (including wall proximity effects), but it retains the simplicity of the "eddy viscosity" concept of first order models. The "modified k-

" is implemented in a three dimensional numerical code. Once the flow is resolved, the predicted velocity and turbulence fields are interpolated into a second grid and used to solve the concentration equation. To evaluate the model, various steady state numerical solutions are compared with small scale dispersion experiments which were conducted at the wind tunnel of Mitsubishi Heavy Industries, in Japan. Stably stratified flows and plume dispersion over three distinct idealized complex topographies (flat and hilly terrain) are studied. Vertical profiles of velocity and pollutant concentration are shown and discussed. Also, comparisons are made against the results obtained with the standard k-

Keywords: Atmospheric dispersion, flow over hills, anisotropic k-e , numerical simulation

Introduction

Atmospheric boundary layer flows are object of intense study over the last years. A more comprehensive understanding of the complex phenomena involved in this particular type of flow is being sought, aiming the analysis of structural implications due to strong winds (neutral atmosphere), the pollutant dispersion under neutral or stable conditions and also for meteorological purposes. The phenomenal increase in computer power over the last two decades has led to the possibility of computing such flows by the integration of the (modeled, time-averaged) Navier-Stokes equations.

Raithby et al. (1987) employed the k-e model (with modification in the Cm value) to calculate the neutrally buoyant flow over the Askervein hill, and compared their numerical results with the experiment made over the real terrain in Scotland. Dawson et al. (1991) also used the k-e model (with some modification in the constants of the dissipation equation) to simulate the flow and dispersion over Steptoe Butte (Washington, USA) under neutrally and stably stratified atmosphere. Their results were favorably compared with experimental data, indicating that mathematical models using the eddy viscosity assumption in the turbulence closure could be used to predict the flow and pollutant dispersion over complex terrain. In Brazil, Dihlmann (1989) studied numerically the thermal discharge (from chimneys) into neutral and stable stratified environments. Santos et al. (1992) applied the standard k-e model to simulate the discharge of a chimney and the correspondent plume dispersion over a flat terrain. Queiroz et al. (1994) applied the standard k-e model to study (in two dimensions) the effect of heat islands in the atmospheric diffusive capacity.

Koo (1993) developed a non-isotropic modified k-e to account for different eddy diffusivities in the lateral and vertical directions in the atmosphere. His model is derived from the algebraic stress model and was applied in one dimensional problems to predict the vertical profiles of velocity, potential temperature and turbulence variables for the horizontal flow in a homogeneous atmospheric boundary layer. Also, the model was applied in two-dimensional problems to simulate the sea breeze circulation and the manipulation of the atmospheric boundary layer by a thermal fence. Koo’s model is similar to the level 2.5 model of Mellor and Yamada (1982). Recently, Castro and Apsley (1997) compared numerical (using a "dissipation modification" k-e model, as named by the authors) and laboratory data for two dimensional flow and dispersion over topography. Also, Boçon and Maliska (1997a, 1997b) extended the non isotropic k-e model of Koo (1993) to numerically simulate the flow and pollutant dispersion over complex idealized topography, under neutral stratification. Computational results were compared with experimental data obtained from a wind tunnel simulation.

More sophisticated models, like the Reynolds stress model, were also applied to predict environmental flows and pollutant dispersion, for instance the work of Andrén (1990). Sykes and Henn (1992) applied the Large Eddy Simulation technique to simulate plume dispersion. Our view is that for the time being, because of limitations in computer resources, those more complex turbulence models (like Reynolds stress and LES) are not suitable for most engineering problems, due to large CPU time and memory required.

In the present work we extend the application of Koo’s modified k-e model to predict three-dimensional stably stratified flows and pollutant dispersion over complex terrain. The prediction of the plume dispersion downwind from a pollutant source is obtained from the solution of the concentration equation. To do so, it’s necessary firstly to calculate the velocity field and eddy viscosities in the region of interest. Thus, convection and turbulent diffusion of the plume may be calculated.

Mathematical Model

The task of computing the concentration field downstream from a pollutant source is divided into two decoupled steps. Firstly we calculate the flow (velocity, temperature and turbulence variables) in the region of interest. Secondly, we use the computed velocity field and eddy diffusivities to solve the concentration equation. This separation can be done as we consider that the pollutant release does not disturb the flow. In fact, in the wind tunnel experiment, against which we compare our results, the tracer gas was released with practically no momentum nor buoyancy force.

Flow and Dispersion Modelling

The governing equations for the stratified flow are the conservation of mass, momentum and energy, written below in the usual tensor notation. Dispersion of a pollutant is computed from the concentration equation, after the flow is resolved.

where p is the pressure deviation with respect to the hydrostatic pressure. Primed variables denote turbulent fluctuations. As we are simulating wind tunnel flows, the Coriolis effect is neglected. Modelling of fluctuation terms are described in the next section.

Turbulence Modelling

In environmental flows the non-isotropic character of turbulence is notable, specially in the case of dispersion of a scalar (pollutant) in the flow. For the case of stably stratified flows, for instance, vertical fluctuations are much inhibited due to buoyancy forces (arising from the positive vertical temperature gradient), while horizontal fluctuations are not. Even neutrally stratified flows feature some anisotropy. So, it is not expected that isotropic turbulence models may well reproduce the non-isotropic turbulent diffusion. However, standard k-e is successfully applied for environmental flows calculation where horizontal gradients (of velocity, temperature and turbulence variables) are smaller than the vertical gradients. In these situations, turbulent diffusion is significant only in the vertical direction, and an isotropic model can handle it appropriately. On the contrary, in the problem of pollutant dispersion from a point source, both vertical and horizontal concentration gradients are significant, so are the corresponding turbulent diffusion. For this situations, a better description of the anisotropy in turbulent exchanges is necessary.

In his Ph.D. thesis, Koo (1993) proposed a modification on the classic k-e model, through use of algebraic stress model including wall proximity effects. The resulting model was compared to data and higher order (turbulence closure) simulations for one and two-dimensional atmospheric flows. The modified k-e reproduced well the observed behaviors.

In our work we extend the application of the Koo’s modified k-e model to three-dimensional flow and to dispersion problems. A description of the turbulence model is given below. Detailed description of derivation of the model can be seen in Koo (1993). Following the Boussinesq’s eddy viscosity concept, Reynolds stresses are related to the gradient of the velocity components as

where K_m ^j is the turbulent eddy viscosity in the j direction. Analogously, turbulent heat exchange and mass dispersion are expressed, respectively, by

where K_h ^j and K_c.^j are the eddy diffusivity in the j direction, respectively for heat and concentration. Eddy viscosities (for momentum) and eddy diffusivities (for energy and concentration) are expressed as functions of turbulent kinetic energy and its dissipation rate. For the vertical direction:

And for the horizontal directions:

C_m is the proportionality coefficient for eddy viscosity in the vertical direction. C_h and C_c are the proportionality coefficients for eddy diffusivities in the vertical direction, respectively, for heat and mass concentration. These three coefficients are defined by functions of flow structure (from the algebraic stress model). Pr_t is the turbulent Prandlt number (=0.5) and Sc_t (=0.5) is the turbulent Schmidt number.

Except for G_M, G_H and f, the coefficients in equations (14) and (15) are model constants, which can be found in Koo (1993). f is the wall function which reflects the effect of the ground proximity on the Reynolds stresses and turbulent heat and mass fluxes, l is the turbulence length scale, k_v is the von Karman constant (=0.4), z is the distance from the ground and C = 0.13.

The C_m, C_h and C_c proportionality coefficients are functions of G_M, the production of turbulent kinetic energy by mean velocity shear, and G_H, the production (or destruction) of turbulent kinetic energy by buoyancy effects

Turbulent kinetic energy and its dissipation rate are computed from their well known prognostic equations:

P is the production term due to mean velocity gradients

G is the production (or destruction) term due to buoyancy

Constants in equations (11), (19) and (20) are those from the standard k-e model, and can be seen in table 1.

Numerical Method

The finite volume method is employed to solve the governing equations, in a non orthogonal, generalized curvilinear coordinate system. Co-located arrangement is used for variables storage in the grid, and the QUICK interpolation scheme with source deferred correction term Lien (1994) is applied on the convection terms, except for turbulence variables where a hybrid scheme (WUDS of Raithby and Torrance, 1967) is adopted. Our own codes NAVIER (1991) and SMOKE (1997) are used to solve the governing equations, respectively, for the flow and concentration.

As the grid used for computing the flow is not adequate for the concentration calculation, a second grid (refined near the source) is used for the last purpose. Velocities and eddy diffusivities obtained from the flow solution are interpolated into the second grid for the concentration calculation. Also, in order to verify grid dependent errors, the computations are made in a coarse and in a fine grid. Figures 1 and 2 illustrate some of the coarse grids used for flow and concentration (inflow boundary at left). Fine grids are 95x41x41 and 128x64x64 for flow and concentration, respectively. Coarse grids have half the number of volumes in each direction, with respect to the fine grids. Only half domain is resolved, because of symmetry.

To verify the model performance, in a first step, the above described modified k-e is applied to simulate wind tunnel experiments.

The Wind Tunnel Experiment

Pollutant dispersion wind tunnel experiments were conducted at the Mitsubishi Heavy Industries, in Nagasaki, Japan, 1991. A report containing the results was obtained directly from that company. Wind tunnel test section is 2.5m wide, 1m high and 10m long. Axisymetric hills of different heights (0, 100 and 200mm), were positioned with the top located at (x, y)=(0,0). Hill shape can be seen in figures 1 and 2. Streamwise direction is x, lateral is y and vertical is z. Source of tracer gas was positioned at (x, y, z)=(-500 mm, 0, 50mm) for hill heights 0 and 100mm, and at (x, y, z) = (-500mm, 0, 100mm) for hill height 200mm. Cases of neutral (T=0, Pasquill class D) and stable atmosphere (T=20 ° C, Pasquill class E) were performed. Streamwise velocity, velocity fluctuations, temperature and concentration were measured at various locations.

Numerical Experiments and Boundary Conditions

Three different wind tunnel experiments were computationally simulated. They are designated with a letter - indicating stability class - followed by a number indicating hill height in mm. Hill heights of 0, 100 and 200mm were simulated. Neutral flows (Pasquill class D) and pollutant dispersion over these topographies were already numerically studied by Boçon and Maliska (1997a, 1997b). At this time, stably stratified flows (Pasquill class E) are considered. At the inflow boundary, velocity, temperature and turbulent kinetic energy are specified according to experimental measured values. As the dissipation rate of turbulent kinetic energy was not measured during the experiment, its inflow profile is calculated according to a prescribed turbulence length scale. For neutral boundary layer atmospheric flows, this length scale increases linearly with the distance from the wall (height above the surface, in the present problems). However, in the case of stably stratified flows the turbulence length scale does not increase linearly with the height, but it is limited to a maximum value (Castro and Apsley, 1997).

where z is the distance from the ground, k_v is the von Karman constant (=0.4) and L is the Monin-Obukhov length (=0.13m), which was calculated from the experimental values of velocity and temperature near the ground.

Outflow conditions are that of zero gradient for all variables. For velocity, lateral and upper boundaries are impermeable, with zero tangential stresses. For all other variables, lateral and upper boundary conditions are of null fluxes. Wall functions are invoked to apply boundary conditions appropriate to a rough wall (z₀ = 1.5e-4m) at the ground. Symmetry conditions are applied at the boundary coincident with the plane of symmetry (y = 0).

Treatment of Near Source Diffusivity

After applying the modified model and computing concentrations, we constated that, for all the cases studied (neutral and stable stratification), there was a large unrealistic plume spread near the source and, consequently, low concentrations everywhere in the domain (specially up to 500mm downstream the source). Taking a look at the turbulence length scale near the source, we noticed that it is larger than the plume dimensions. It means that the turbulent eddy sizes present in the flow are bigger than the plume, and could not promote such a observed diffusion in the numerical simulations. Therefore, we speculate that the length scale to be applied in the eddy diffusivities for the concentration should be appropriately reduced for the initial stages of plume spread, according to local plume dimensions. Based on a Gaussian plume distribution near the source, as a first investigation, we decide to reduce linearly the eddy diffusivities computed from the flow solution, to be applied in the concentration calculations. Using this simple procedure, the quality of the results improved considerably. Reduction of eddy diffusivity near the source is made taking its value (at the source location), obtained from the flow solution, and applying it in the Gaussian model for diffusion from a point source, to calculate how far from the emission point the plume width is about five times the local turbulence length scale. This value (five) was empirically determined from the analysis of the results. It was found, however, that the value is roughly the same for all the cases studied. At a given distance from the source, plume width is defined as the distance from the plume center line to the point where the concentration is 10% of its peak value. Indeed, further work is needed to better model the initial stages of plume spread, where its dimensions are smaller than the characteristic turbulence length scale of the flow.

Results and Discussion

In this section, some results of the flow and concentration calculations are presented, for the three stably stratified cases which were computationaly simulated. Figures 3, 4 and 5 show vertical profiles of concentration on the symmetry plane (y=0) for the cases E0, E100 and E200. The graphs in each figure refer to different positions downstream the source. In figure 3, for the case of flat terrain (E0), it can be seen a good agreement between numerical and experimental values. At the position x=200mm (third graph in fig. 3), we believe that there possibly was a mistake with respect to the report of the experimental results, which are inconsistently underestimated (as it can be noticed by comparison between the two peak concentration values corresponding to the positions x=0 and x=500mm).

Figures 4 and 5 show, respectively, the concentration profiles for the cases of hilly terrain E100 and E200 (hill heights 100 and 200 mm). From the view of a numerical analist, these are the most critical cases, due to the characteristics of stable flow and complex topography. Although the peak concentration values are fair well predicted, their locations are not. For the problem of pollutant dispersion over complex terrain, the plume path is dictated by the deviations in the mean flow caused by the irregular topography. A correct description of the plume path requires a sufficiently accurate prediction of the flow field (which "drives" the plume).

Figures 6 and 7 show vertical profiles of the streamwise component (u) of velocity on the symmetry plane (y = 0) for the cases E100 and E200. For both cases, the modified and the standard k-e model produced nearly the same velocity profiles, and the recirculation zone in the lee side of the hill was underestimated. Different inflow turbulent length scales were tested at the inflow to verify a possible influence, but it was noticed that the flow after the hill top is essentially determined by local conditions. A possible explanation for this model defect would be that the pronounced velocity gradients in this region, due to the three dimensional open recirculation zone (see figure 8), increase the production of turbulent kinetic energy and consequently enhance the eddy viscosities there, thus diminishing the size of the recirculation. However, a comparison between numerical and wind tunnel measured values did not reveal that the turbulent kinetic energy level has been overestimated by the mathematical model. Thus, regarding to the above cited problem of high eddy viscosities in the recirculation zone, the model drawback should be attributed to the dissipation equation, which is underestimating e , and not to an overestimation of the production of turbulent kinetic energy (P).

The underestimation in the size of the open recirculating three-dimensional zone by the flow model leads to an incorrect determination of the plume path. As the numerical model foresees a smaller recirculation, the flow, after passing over the hill top, brings (convectively) the plume down to the ground. Thus, in the numerical simulation, the center line of the plume (the place of the peak concentration values along the plume path) is nearer the ground than the experimental plume, and the ground level concentration numerical values result higher than the measured ones.

For the case E200, figure 5 also shows the concentration results produced by the standard k-e model (isotropic), which clearly are more diffusive than those from the modified model. In this sense, it is demonstrated that the non-isotropic character of the turbulent diffusion under stably stratified flow conditions is relevant, and that the modified k-e has better ability to predict the dispersion of a plume in such conditions.

Conclusions

A modified non-isotropic k-e model is applied to simulate three dimensional stably stratified atmospheric flows (Pasquill class E) and dispersion over an idealized complex terrain. The results for the velocity field are similar to those obtained with the standard model, because the vertical gradients of flow variables are not great (when compared with the concentration field), resulting nearly the same turbulent diffusion for both the standard (isotropic) and the modified (non-isotropic) k-e models. However, for the concentration results, the differences between the numerical solutions obtained with the standard and the modified k-e models are quite distinct. The agreement of the concentration values - against the wind tunnel results - is very satisfactory for the case of flat terrain. Over hilly terrain, the concentration peak values are fair well predicted, but their locations (height above the ground) are not. The problem is attributed to a model fault in predicting the open recirculating three-dimensional zone, which occurs in the lee side of the hills. In the recirculation zones, the eddy viscosities are overestimated and thus the size of the recirculation is underpredicted. Consequently, the plume path over the top and the lee side of the hill is not well predicted, with respect to its height above the ground, which is underestimated. Despite of these drawbacks, the modified non-isotropic k-e produces better results than the conventional isotropic model.

Acknowledgments

We are grateful for the financial support provided by CNPq and CAPES.

Article presented at the 1st Brazilian School on Transition and Turbulence, Rio de Janeiro, September 21-25, 1998.

Technical Editor: Atila P. Silva Freire.

• Andrén, A., 1990, "A Meso-Scale Plume Dispersion Model. Preliminary Evaluation in a Heterogeneous Area", Atmospheric Environment, v. 24A, n. 4, p. 883-896.
• Boçon, F. T. and Maliska, C. R., 1997a, "Numerical Modelling of Flow Over Complex Terrain", Proceedings, XVIII Iberian Latin American Congress on Computational Methods in Engineering, Brasília - DF, published in CD-ROM.
• Boçon, F. T. and Maliska, C. R., 1997b, "Numerical Modelling of Flow and Dispersion Over Complex Terrain", Proceedings, XIV Brazilian Congress of Mechanical Engineering, ABCM, Bauru-SP, p. 211-218.
• Castro, I.P. and Apsley, D.D., 1997, "Flow and Dispersion Over Topography: A Comparison Between Numerical and Laboratory Data for Two-Dimensional Flows", Atmospheric Environment, vol. 31, no 6, pp. 839-850.
• Dawson, D., Stock, D.E. and Lamb, B., 1991, "The Numerical Simulation of Airflow and Dispersion in Three-Dimensional Atmospheric Recirculation Zones", J. Applied Meteorology, vol. 30, pp. 1005-1024.
• Dihlmann, A., Maliska, C. R., Silva, A. F. C., 1989, "Solução Numérica da Descarga de Jatos Poluentes em Meio Estratificado", Anais do X Congresso Brasileiro de Engenharia Mecânica, ABCM, Rio de Janeiro-RJ, p. 101-104.
• Koo, Y.S., 1993, "Pollutant Transport in Buoyancy Driven Atmospheric Flows", Ph.D. Thesis, The Louisiana State University and Agricultural and Mechanical Col.
• Lien, F.S. and Leschziner, M.A., 1993, "Upstream Monotonic Interpolation for Scalar Transport With Application to Complex Turbulent Flows", Int. J. For Numerical Methods in Fluids, vol. 19, pp. 527-548.
• Mellor, G.L. and Yamada, T., 1982, "Development of a Turbulence Closure Model for Geophysical Fluid Problems", Reviews of Geophysics and Space Physics, vol. 20, no 4, pp. 851-875.
• NAVIER, 1991, "Desenvolvimento de Códigos Computacionais para Solução de Problemas de Escoamentos de Alta Velocidade", Relatório preparado para o Instituto de Atividades Espaciais do Centro Técnico Aeroespacial - UFSC - Dep. Eng. Mecânica - Parte VII.
• Panofsky, H. A., Dutton, J. A., 1984, "Atmospheric Turbulence - Models and Methods for Engineering Applications", John Wiley & Sons, New York.
• Raithby, G.D. and Torrance, K.E., 1967, "Upstream-Weighted Differencing Schemes and Their Application to Elliptic Problems Involving Fluid Flow", Computer and Fluids, vol. 2, pp. 12-26.
• Raithby, G. D., Stubley, G. D., and Taylor, P. A., 1987, "The Askervein Hill Project: A Finite Control Volume Prediction of Three-Dimensional Flows over the Hill", Boundary-Layer Meteorology, v. 39.
• SMOKE, 1997, "Código Computacional para a Solução da Dispersão de Escalares em Geometrias Tridimensionais", Laboratório de Simulação Numérica em Mecânica dos Fluidos e Transferência de Calor do Departamento de Engenharia Mecânica da Universidade Federal de Santa Catarina - SINMEC.

Publication Dates

• Publication in this collection
  18 Jan 2001
• Date of issue
  2000