Microsoft Word - 025.docx CHEMICAL ENGINEERING TRANSACTIONS VOL. 48, 2016 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Eddy de Rademaeker, Peter Schmelzer 2016, AIDIC Servizi S.r.l., ISBN 978-88-95608-39-6; ISSN 2283-9216 Modelling of Stably-Stratified Atmospheric Boundary Layers with Commercial CFD Software for Use in Risk Assessment Rachel Batt*a, Simon E. Gant*a, Jean-Marc Lacomeb, Benjamin Truchotb a Health and Safety Executive, Buxton, SK17 9JN, UK b INERIS, Dept. PHDS, Parc Technologique ALATA, BP 2, 60550 Verneuil-en-Halatte, France rachel.batt@hsl.gsi.gov.uk © Crown copyright (2015) One the of the significant challenges faced in simulating atmospheric boundary layers (ABLs) with standard computational fluid dynamics (CFD) models is to preserve the correct ABL profiles throughout the flow domain. Appropriate ABL profiles may be imposed at the inlet to the CFD domain but they are often progressively modified downwind by the CFD model until they are no longer representative of the correct stability class and/or wind speed. To address this issue, this paper reviews solutions proposed in the literature with respect to the choice of turbulence model and CFD boundary conditions (inlet profiles, wall conditions, outlet and top boundary conditions). The main focus is on modelling stable ABLs which often produce the largest dispersion distances in hazard analyses. CFD results are presented for a 2 km long domain with flat terrain using ANSYS-CFX with standard and modified k-ε turbulence models and two different choices of inlet profiles. The results show that improvements in the ABL profiles are obtained using a modified turbulence model with a consistent choice of inlet profile. 1. Introduction There is growing interest in the use of CFD to assess the risks posed by atmospheric releases of toxic and flammable gases, due to the potential for CFD to simulate the effects of complex terrain and obstructions. However, there are a number of challenges to be overcome in modelling atmospheric boundary layers, especially under stably-stratified conditions. A central problem is that standard CFD models are unable to preserve the correct ABL profiles throughout the flow domain. Even though appropriate profiles may be imposed at the inlet, they are progressively modified downwind by the CFD model. Eventually, further downwind, equilibrium is reached between the modified profiles and the CFD model, but these profiles are no longer representative of the correct stability class and/or wind speed. Turbulence profiles are also affected, which is important since they may have significant impact on mixing and dilution. Dispersion predictions obtained using the modified profiles may be adversely affected by the incorrect ABL profiles. Different dispersion model predictions may also be produced depending upon whether the CFD model incorporates a short or a long upstream fetch. A number of different approaches have been proposed in the literature to address this issue and are reviewed in the following section. CFD results are then presented using one of the more promising approaches. 2. Literature Review 2.1 Turbulence Model The focus of the present work is on the k-ε turbulence model that is widely used in CFD codes for risk assessment. The standard constants for the k-ε model were originally chosen to produce good predictions for a range of classical shear flows found in engineering (Jones and Launder, 1974). However, it has been demonstrated in many previous studies that these model constants do not preserve the correct profiles in ABLs. Various researchers have therefore proposed turbulence model modifications specifically aimed at improving predictions for ABLs (e.g. Richards and Hoxey, 1993, Alinot and Masson, 2005; Pontiggia et al., 2009; Vendel, 2011). Whilst these modifications may be beneficial for simple ABLs over flat, unobstructed DOI: 10.3303/CET1648011 Please cite this article as: Batt R., Gant S., Lacome J., Truchot B., 2016, Modelling of stably-stratified atmospheric boundary layers with commercial cfd software for use in risk assessment, Chemical Engineering Transactions, 48, 61-66 DOI:10.3303/CET1648011 61 terrain, it is unclear whether they will actually lead to worse predictions in more complex cases, for example involving flow around obstacles. Recently, models have therefore been proposed that take a zonal approach, where the turbulence model modifications are removed near obstructions (e.g. Balogh et al., 2012). Alinot and Masson (2005) [henceforth referred to as A&M] and Pontiggia et al. (2009) both proposed modifications to the k-ε model to simulate stable ABLs. A&M applied their modifications by changing the k-ε model constants whilst Pontiggia et al. (2009) took a different approach of introducing additional source terms. Both models have been examined but efforts have focused on A&M. 2.2 ABL Inlet Profiles In a neutral ABL, the shear stress within the near-wall layer of the ABL (i.e the surface layer) is approximately constant and the velocity profile is well-described by a log-law. There are standard corrections to this log-law to account for the surface roughness of the ground (see Stull, 1988). Under stably-stratified conditions, modifications are usually made to the log-law based on Monin-Obhukov similarity theory (MOST). One of the inputs to these corrections is the Monin-Obukhov length, L, which represents the height at which the production of turbulent kinetic energy due to buoyancy is balanced by that due to wind shear. Lacome and Truchot (2013) found that one of the issues responsible for significant variability in CFD results for atmospheric dispersion was that different CFD modellers used different approaches to translate from the Pasquill stability class to the value of L. To address this issue, they prescribed a particular approach as given by Golder (1972) and the TNO Yellow Book (Van den Bosch and Wetering, 2005). The log-law-based wind profiles are only strictly valid within the surface layer, which is less than 100 m deep in stable ABLs. For some hazard analyses, particularly those involving terrain or obstacles, it may be necessary to use a higher domain. Different approaches have been proposed to extend profiles above the surface layer (see Zilitinkevich et al., 1998). Lacome and Truchot (2013) recommended the velocity profile of Gryning et al. (2007). Temperature profiles are commonly defined in a similar manner to the velocity, based on the log-law with MOST corrections (Stull, 1988). A complication with temperature is that it is affected by the pressure and humidity within the atmosphere. In a dry atmosphere, the temperature naturally falls with height at rate of approximately 1°C per 100 m (i.e. the lapse rate). Although this temperature drop is modest, if it is not accounted for properly in the CFD model it can lead to unphysical flow behaviour, usually near the outlet boundary. To account for this effect, some CFD models solve for the “potential” temperature, θ, which in a neutral ABL remains constant with height. This approach is well-suited to CFD models using the Boussinesq approximation, where buoyancy forces are expressed in terms of a vertical temperature gradient rather than a density gradient. The inlet profiles for the turbulence quantities are usually derived from analytical solutions to the relevant turbulence transport equations, assuming that the ABL is fully developed: i.e. zero streamwise gradients, constant shear stress and zero vertical velocity. Suitable profiles for neutral and stable ABLs were derived by, for example, Richards and Hoxey (1993) and Vendel (2011). In order to maintain the correct velocity profile throughout the domain, the turbulence profiles should be consistent with the velocity and temperature profiles. 2.3 Wall boundary conditions Close to the wall in a turbulent flow, there are usually rapid changes in the velocity, temperature and turbulence characteristics. Most engineering CFD models do not attempt to resolve these sharp gradients but instead use “wall functions” within the near-wall cell. There are two significant issues with wall functions for modelling ABLs. Firstly, for CFD simulations over urban terrain it is common practice to use a large aerodynamic roughness length, z0, of the order of 0.3 m. This length represents the height where the velocity falls to zero, when extrapolated from higher elevations assuming a logarithmic profile. The majority of commercial CFD codes do not allow the input of z0 directly but instead use the equivalent sand-grain roughness, ks, where ks ≅ 30 z0, or around 9 m for urban roughness. Usually, CFD models require the node within near-wall cell to be placed at a height greater than ks. This means that for urban-scale roughness, the minimum near-wall cell height is 18 m, which is clearly inappropriate for modelling dense-gas dispersion when the entire gas cloud may be less than 18 m deep. For other CFD codes that use wall functions based on z0, a smaller cell size can be used in practice, but the physical basis of the model is unclear. Alternatively, it is possible to assume a smooth wall and accept that dispersion model predictions may be affected, or to resolve the urban obstacles, which may require a fine grid and long computing times. The second issue with wall functions is that they are based on assumed profiles for velocity, temperature and turbulence within the near-wall cell that may be inconsistent with the ABL inlet profiles. Blocken et al. (2007) discussed various remedial measures and both Hargreaves and Wright (2007) and Parente et al. (2011) presented improvements, although these required complex modifications to the CFD code. 62 2.4 Top boundary conditions Different approaches have been proposed for modelling the top boundary of the domain, which include: • Shear stress (e.g. Hargreaves and Wright, 2007, Parente et al. 2011) • Inlet condition (e.g. Alinot and Masson, 2005, Blocken et al. 2007, Pontiggia et al. 2009) • Symmetry (e.g. Franke et al. 2007) • Flux (e.g. Vendel, 2011) • Pressure (e.g. Montavon, 1998) It is generally accepted that the top boundary of the computational domain should be as far away from the flow field of interest as possible in order to minimise its effects. Franke et al. (2007) recommended that the domain height should be at least five times the height of the largest obstacle whilst Pontiggia et al. (2009) recommended that in dispersion simulations it should be at least double the maximum gas cloud height. 2.5 Outlet boundary conditions CFD models of the ABL commonly apply a constant pressure condition at the outlet that is derived from analytically solving the vertical momentum equation, assuming fully-developed flow. If the air density is constant then the relative pressure can be set to zero at the outlet. This applies to a neutral ABL where a Boussinesq approach is used. In all other cases, it is necessary to apply a pressure profile at the outlet to avoid unphysical flow behaviour. This includes cases where the atmosphere is neutral but the air is treated as an ideal gas. For stable ABLs, the pressure profile requires integration of the temperature or density profile (see Vendel, 2011). All other variables on the outlet are usually assigned a zero-gradient condition across the outlet boundary. 3. CFD Modelling 3.1 Methodology CFD simulations have been performed using the commercial CFD software, ANSYS-CFX v15, of a moderately stable ABL with Pasquill Class F and a wind speed of 2.4 m/s (i.e. F2.4). This is a standard weather condition that is used in dispersion modelling of toxic substances for land-use planning in the UK. Two different methods have been tested to specify the inlet profiles: those proposed by Lacome and Truchot (2013) and A&M. A summary of the relevant ABL inlet conditions is given in Table 1. Table 1: Input parameters for the ABL profiles Lacome and Truchot (2013) Alinot and Masson (2005) Velocity, u(H) 2.4 m/s 2.4 m/s Reference height (H) 10 m 10 m Roughness length, z0 0.005 m 0.005 m Monin-Obukhov length, L 8.73 m 8.73 m Friction velocity, u* 0.0623 m/s 0.0756 m/s Boundary layer height, zi 29.5 m 32.5 m Temperature, T(H) 5 °C 5 °C Two different turbulence models have been tested: the standard k-ε model in ANSYS-CFX and the modified A&M model. A summary of the constants in the two models is provided in Table 2. To implement the A&M model in ANSYS-CFX involves complex user-coding with Fortran, since Cε3 is a function of height. Verification tests were undertaken to demonstrate that the results presented in A&M’s paper could be reproduced. Buoyancy effects were modelled using a Boussinesq approach and the potential temperature was solved, rather than the real temperature. In terms of its implementation in ANSYS-CFX, solving for potential temperature with a Boussinesq approach simply involves removing the lapse rate from the inlet temperature profile and ensuring that consistent temperatures are applied at the other boundaries. There are appropriate default source terms in the k-ε equations in ANSYS-CFX and no additional source terms are needed, but if real temperature were used, user-defined source terms would be needed in the momentum and k-ε equations. Table 2: Summary of turbulence model constants Cε1 Cε2 Cε3 Cμ σk σε Standard 1.44 1.92 1.0 0.09 1.0 1.3 A&M 1.176 1.92 Profile 0.033 1.0 1.3 63 The top boundary condition was treated as an inlet with prescribed values of velocity, temperature, k and ε, following the approach taken by A&M. The relative pressure was prescribed on the outlet boundary using a user-defined profile that involves the integral of the vertical temperature profile, similar to the correction of Vendel (2011). The bottom boundary was set as a rough wall with ks = 29.7 z0. The computational domain was two-dimensional with a length of 2 km and height of 245 m. A Cartesian mesh was used with 502 cells in each direction, which were clustered near the ground to produce a near-wall cell height of 0.4 m. This mesh and domain size were chosen to be comparable to the work previously presented by A&M. Since the roughness height was 5 mm and ANSYS-CFX uses a wall-function based on sand-grain roughness, the minimum allowable near-wall cell size was around 0.3 m. Second-order accurate differencing was used for all of the transport equations (including both k and ε). 3.2 Results The ABL velocity profile produced by the CFD model with the standard k-ε turbulence model is shown in Figure 1. The inlet profiles for this case are prescribed for F2.4 using the Lacome and Truchot (2013) method. The results show that velocity at the outlet is significantly different from the profile prescribed at the inlet. The velocity is higher near the ground and lower above 10 m. The correct profiles for F2, F3 and E2.4 are shown for comparison. It is clear that the outlet profile is not representative of the same stability class as the inlet profile. If this model were used for a toxic dispersion calculation, the higher velocity near the ground could artificially increase dilution rates, which would produce an inappropriate reduction in the hazard range. Figure 1: Velocity profile at the inlet and outlet simulated with Lacome and Truchot (2013) inlet profiles and CFX standard k-ε model. Typical velocity profiles of relevant stability classes are shown for comparison. Figure 2 shows the profiles of velocity, turbulence and temperature at the inlet and outlet using the A&M turbulence model with the Lacome and Truchot (2013) inlet profiles. The agreement between the inlet and the outlet velocity is slightly improved near the ground, but there is still a significant difference above 6 m. Figure 3 shows the final case where the A&M turbulence model is used with the A&M inlet profiles. There is better agreement between the inlet and outlet profiles in this case. However, the agreement is still far from perfect, especially for k and ε. The reason why the results presented here are not as good as those published by A&M is principally that they only considered a domain length of 450 m, whereas a 2 km long domain is considered here. A&M also considered a higher velocity of U(H) =10 m/s in a less strongly stratified ABL with a Monin-Obukhov length of 191 m. 64 Figure 2: Profiles of turbulence kinetic energy, turbulence dissipation rate, temperature and velocity at the inlet and outlet for Pasquill Class F2.4 with Lacome and Truchot (2013) inlet profiles and the A&M modified turbulence model. Variables are non-dimensionalised using the same method as A&M. Figure 3: Profiles of turbulence kinetic energy, turbulence dissipation rate, temperature and velocity at the inlet and outlet for Pasquill Class F2.4 with A&M inlet profiles and the A&M modified turbulence model. Variables are non-dimensionalised using the same method as A&M. Note that the T profiles are obscured by the u profile. 4. Conclusions A brief review has been presented of approaches proposed in the literature for CFD modelling of stable ABLs. The central issue that has been investigated is whether CFD models can maintain the correct ABL profiles along the length of a 2 km long, flat, unobstructed domain. The review discussed a number of ongoing issues and identified the Alinot and Masson (2005) model as deserving further investigation. CFD simulations were then presented for moderately stable F2.4 weather conditions using ANSYS-CFX with two different turbulence models (the standard k-ε model and A&M variant) and two choices of inlet profiles (those proposed by Lacome and Truchot, 2013, and A&M). The results showed that the A&M model improved the predictions over the standard k-ε model and that the choice of A&M inlet profiles produced more consistent ABL profiles than those from Lacome and Truchot (2013). However, even when using A&M’s approach for both the turbulence model and inlet profiles, the ABL profiles were still modified along the length of the 2 km 65 domain. These differences probably arose from inconsistencies with the wall functions, but these are very difficult to adjust in commercial CFD codes, such as ANSYS-CFX. Two significant issues remain to be investigated: firstly, the effect of the modified ABL profiles on dispersion predictions; and secondly, the effect of modified turbulence models (like A&M) on more complex flows, such as those involving dispersion around buildings. It is proposed to investigate the first of these issues in future work using data from the Prairie Grass and Thorney Island experiments. Acknowledgments The contribution made to this paper by Rachel Batt and Simon Gant was funded by the Health and Safety Executive (HSE). The contents, including any opinions and/or conclusions expressed, are those of the authors alone and do not necessarily reflect HSE policy. References Alinot C., Masson C., 2005, k-ε model for the atmospheric boundary layer under various thermal stratifications, J. Solar Energy Engineering, 127, 438-443. Balogh M., Parente A., Benocci C., 2012, RANS simulation of ABL flow over complex terrains applying an enhanced k-ε model and wall function formulation: Implementation and comparison for FLUENT and OpenFOAM, J. Wind Eng. Ind. Aerodyn., 104-106, 360-368. Blocken, B., Stathopoulos, T and Carmeliet, J., 2007, CFD simulation of the atmospheric boundary layer: Wall function problems, Atmos. Env., 41, 238 – 252. Van den Bosch, C.J.H. and Weterings, R.A.P.M (Eds.), 2005, Methods for the calculation of physical effects due to releases of hazardous materials (liquids and gases), ‘TNO Yellow Book’, 3rd Ed. Second revised print, TNO, The Hague, The Netherlands. Franke J., Hellsten A., Schlünzen H., Carissimo B. (Eds.), 2007, Best practice guideline for the CFD simulation of flows in the urban environment, COST Action 732 “Quality assurance and improvement of microscale meteorological models”. COST Office, Brussels, Belgium. Golder, D., 1972, Relations among stability parameters in the surface layer, Boundary-Layer Meteorol., 3 (1), 47-58. Gryning S.-E., Batchvarova E., Brummer B, Jorgensen H., Larsen S., 2007, On the extension of the wind profile over homogeneous terrain beyond the surface boundary layer, Boundary-Layer Meteorol., 124, 251-268. Hargreaves D.M., Wright N.G., 2007, On the use of the k-ε model in commercial CFD software to model the neutral atmospheric boundary layer, J. Wind Eng. Ind. Aerodyn., 95, 355-369. Jones W.P. Launder B.E., 1974, The numerical computation of turbulent flows, Comput. Methods Appl. Mech. Eng., 3, 269-289. Lacome J.-M., Truchot B., 2013, Harmonization of practices for atmospheric dispersion modelling within the framework of risk assessment, 15th conference on “Harmonisation within Atmospheric Dispersion Modelling for Regulatory Purposes”, Madrid, Spain, 6-9 May 2013. Montavon C., 1998. Simulation of atmospheric flows over complex terrain for wind power potential assessment. PhD Thesis. Ecole Polytechnique Fédérale de Lausanne, Switzerland. Parente A., Gorle C., van Beeck J., Benocci C., 2011, Improved k-ε model and wall function formulation for the RANS simulation of ABL flows, J. Wind Eng. Ind. Aerodyn., 99, 267-278. Pontiggia M., Derudi M., Busini V., Rota R., 2009, Hazardous gas dispersion: A CFD model for atmospheric stability classes, J. Hazard. Mater., 171, 739-747. Richards P.J., Hoxey R.P., 1993, Appropriate boundary conditions for computational wind engineering models using the k-ε turbulence model, J. Wind. Eng. Ind. Aerodyn., 46-47, 145-153. Stull, R., 1988, An Introduction to Boundary Layer Meteorology, Kluwer Academic Publishers, London. Vendel F., 2011, Modélisation de la dispersion atmosphérique en presence d’obstacles complexes: application à l’étude de sites industriels, PhD Thesis, L’Ecole Centrale de Lyon, France. Zilitinkevich, S.S., Johansson, P.-E., Mironov, D.V., Baklanov, A., 1998, A similarity-theory model for wind profile and resistance law in stably-stratified planetary boundary layers, J. Wind Eng. and Ind. Aerodyn., 74-76, 209-218. 66