Large Eddy Simulation of stratified flows over structures

We tested the ability of the LES model CLMM (Charles University Large-Eddy Microscale Model) to model the stratified flow around three dimensional hills. We compared the quantities, as the height of the dividing streamline, recirculation zone length or length of the lee waves with experiments by Hunt and Snyder[3] and numerical computations by Ding, Calhoun and Street[5]. The results mostly agreed with the references, but some important differences are present.


Introduction
The atmospheric flows affected by a stable temperature stratification are known to be connected with problematic pollution dispersion and are particularly difficult to simulate. In this contribution we present a set of computations of idealised stratified flows around an isolated hill as a test of capabilities of the numerical model CLMM (Charles University Large-Eddy Microscale Model). This configuration can serve as a model of flow around an island. The model CLMM has been previously validated for the case of moderately stratified turbulent atmospheric boundary layer over a flat terrain [2].
The chosen problem was previously studied both experimentally and numerically. Hunt and Snyder [3] used a water tank with salinity induced stratification. Baines [4] presented a detailed review of his and other results for this problem. Ding, Calhoun and Street [5] used large eddy simulation (LES) to numerically reproduce results of Hunt and Snyder. Bodnár and Beneš [1] solved this problem in 2D and compared two numerical methods and showed the influence of the inflow conditions. We will [3] and [5] use these two references for comparison.
It is important to use a correct scaling for interpretation of flows. In the case of stably stratified flow over an obstacle a useful non-dimensional scale is where U is the average velocity of the incoming flow, h is the height of the obstacle and N is the Brunt-Väisälä (or buoyancy) frequency, which is defined as where g is the gravitational acceleration, ρ θ is the potential density and θ is the potential temperature assuming that density depends linearly on temperature and that the density changes due to pressure are negligible. The buoyancy frequency is a frequency of harmonic oscillations of an air a e-mail: vladimir.fuka@mff.cuni.cz parcel moved from its stable position. For a stably stratified fluid this quantity is positive. We will neglect the effects of rotation in the simulations. This is admissible for flows in low latitudes.

Numerical model
The model CLMM uses the incompressible Navier-Stokes equations in the Boussinesq approximation. The subgrid turbulence is modelled using the large eddy simulation and the eddy viscosity approximation. The set of solved equations is where ν t and κ t are the eddy viscosity and the eddy diffusivity. The eddy viscosity is computed using the Vreman subgrid model [6]. We use a constant value of the Prandtl number to compute the eddy diffusivity, i.e. κ t = 0.6/ν t . This means we did not introduce any stability correction to the subgrid model in this study. These equations are discretized separately in space and time. In space we use the finite volume method on a uniform staggered rectangular grid. The fluxes are discretized with second order central differences, with the exception of the scalar advection which is computed by means of a third order positive method by Hundsdorfer at al. [7]. The time discretization is based on the combination of the third order Wray's Runge-Kutta and the second order Crank-Nicolson method with the pressure correction in each stage [8].
Because of the rectangular grid, it is necessary to use the immersed boundary method, where additional sources of momentum and scalars are imposed inside solid bodies. We used the direct forcing variant of the immersed boundary method by Kim, Kim and Choi [9,10].

Boundary conditions
The computational domain comprised a flat terrain and a rotationally symmetric hill with the shape defined by the function: where r is the distance from the centre and h is the height of the hill. This is slightly different from the previous studies. We also simulated the flow around the original obstacle shape by Hunt and Snyder for validation: This shape has a plateau on the top and is wider. At the inlet, a uniform velocity profile was prescribed. At the outlet we applied the Neumann boundary conditions with zero derivative and a buffer layer where all fluctuations were damped [11]. At the lateral boundaries and at the top, the free slip boundary conditions were prescribed. For U/(Nh) ≤ 0.2 the damping layer was also placed at the top boundary to avoid refraction of the gravity waves.
The Reynolds number was not constant, but was different for various U/(Nh). Its value was much larger than in laboratory experiments [3] and was in the range [7 × 10 6 , 35 × 10 6 ] which is more realistic for atmospheric flows.
The flow near the walls can not be resolved, due to the limited resolution and the very high Reynolds number. Therefore we applied a wall function based wall model with a stability correction according to the Monin-Obukhov similarity theory.

Results
The flow was solved for the following values of U/(Nh): infinity (neutral), 10, 5, 2, 1.41, 1, 0.71, 0.5, 0.2 and 0.1 (the smaller the value the stronger is the stratification). The character of the flow strongly depended on the level of stratification. The slightly stratified cases (U/(Nh) = 10, 5) were similar to the neutral one, with a shorter recirculation region behind the hill and with some indication of a very long lee wave. For moderately stratified cases (U/(Nh) = 2 to 0.71) the recirculation length kept shortening with increasing stratification and the lee waves were short enough to be comparable with the hill dimensions (figures 1 and 2). For the strongly stratified cases, the air does not have enough kinetic energy to overflow the hill and tends to flow horizontally. Therefore a Kárman vortex street develops in the wake. Also, the lee waves are shorter than the hill dimension and they are propagating more upwards, but their amplitude is smaller ( figure 3).
Another important characteristic is the shape of the recirculation region behind the hill. We concentrated mainly on the length of the recirculation zone l r . For increasing  stratification l r is decreasing up to U/(Nh) ∼ 0.7 where it is no longer closed. For even stronger stratification levels the wake structure is very unsteady and a vortex street develops and the length of the region with the reversed average flow gets much larger (l r > 4). A snapshot of streamlines for an example of a vortex street is in figure 4. Comparison of these results with experimental data is in the graph in figure 5. The agreement is not very good. The reason can be in the different Reynolds number and the lack of near wall resolution and usage of a wall model. We also note, that even though [5] do not report recirculation length values, from their figure with streamlines for U/(Nh) = 1 it can be seen that their result is quite close to our result.
We investigated the streamline reaching the top of the hill as defined by [3]. The quantity used is the height of this streamline far upstream. The results are in figure 6. For high U/(Nh) (i.e. low U/(Nh)) the results fit the empirical curve: For U/(Nh) >= 0.5 the vertical velocity near the inflow boundary was not equal to 0, therefore the results are not accurate.
We also analysed the length of the lee waves. Here we will present the numbers from computation with a coarser resolution on a larger domain (51.2h×12.8h×12.8h) for the cosine shaped hill. For strong stratification the results are equal to the high resolution ones, but in addition, the longer waves in weaker stratification were possible to measure. We compare the results with the computations of Ding et al. [5] a linear theory (see [4]). The results are in figure 7.

Conclusions
We simulated the flow over a hill in a stably stratified atmosphere. The behaviour agreed with the theoretical expectations and with the laboratory measurements in many aspects. The are still unexplained differences, especially in the recirculation zone shape. One of the possible reasons can be the fact that our model does not fully resolve the 01031-p.2  wall layer nor it enhances resolution near the walls. The differences in z t might be caused by an insufficient domain length (∼ 15h in front of the hill). The work is not finished and will be followed by the flow in a fully turbulent atmospheric boundary layer with realistic boundary conditions.   2 π U/(Nh) Fig. 7. Wavelength of the lee waves for different levels of stratification.