CFD ANALYSIS OF TERRAIN INDUCED TURBULENCE AT KRISTIANSAND AIRPORT, KJEVIK

. The paper presents the finding of a study conducted in order to analyse the characteristics of flow in a real complex terrain in the vicinity of the Kristiansand Airport (Kjevik), Norway. The paper starts with a brief description of the governing equations, their numerical discretization and the safety criteria used by the aviation industry to classify turbulence as mild or extreme. Directional dependence of the flow characteristics in the region is studied followed by an attempt to validate the experience of the pilots operating in the region.


Introduction
Kristiansand Airport, Kjevik ( Fig. 1) is situated 4.3 NM (8.0 km) north-east of the city Kristiansand, Vest-Agder in southern Norway, located 16 km from the city centre. The airport serves the Agder district with domestic and international flights. The airport is operated by Avinor. Surrounded by water on three sides and hills on the fourth the airport occupies an interesting location (Fig. 2). The terrain in the vicinity of the airport is highly irregular with hills stretching up to a height of 300 m above sea level (Fig. 3). The airport itself is just at sea level. Presence of the sea in the south-west direction and valley in the north-east direction make it easier for the flights to take off and land along this direction. Although the hills are not too high, the irregularities are enough to induce vigorous turbulence in the region. According to various pilot's reports the turbulence intensity in the region is strongly influenced by the terrain and wind magnitude and direction. In the case of Kristiansand Airport, maximum turbulence intensity has been reported for the north westerly wind. The main cause is associated with the high hills on the north-west side of the airport. The flight path and airport, thus, lie on the leeward side of the hill, and therefore, get exposed to strong vortices generated by the topography. Our main aim in this work is to investigate these claims and to relate them to the risk involved in flight operations.

Governing equations
The equation of motion for incompressible flow may be generalized to atmospheric flows by the use of the socalled anelastic approximation. This formulation is often applied in meteorological models, and may be written in the following conservative form (Bannon 1995;Durran 1998): (2) Here (u, p, θ , ρ ) represent velocity, pressure, potential temperature and density, respectively. Furthermore, τ is the stress tensor, f is a source term that may include rotational effects, g is the gravitational acceleration, γ is the thermal diffusivity and q is the energy source term. Subscript s indicates hydrostatic values, and subscript d the deviation between the actual value and its hydrostatic part, i.e. p = p s + p d ; where the hydrostatic part is provided by / s s p z g ∂ ∂ =-ρ . In addition, the following expression for hydrostatic density may be derived from the state equation and the definition of potential temperature: where R denotes the gas constant and Cp is the specific heat at constant pressure. Hence, once the hydrostatic (potential) temperature profile is given, the hydrostatic pressure and density may be calculated, and then substituted into Equations (1) and (2). It may be noted that the Boussinesq approximation is obtained from the system of Equations (1) and (2) by assuming constant values ( 0 ρ ; 0 θ ) instead of the hydrostatic values, and that formulation may well be used for incompressible flow and ordinary temperature.
The aim of the present study is to solve these equations for high Reynolds-number flows. For this purpose we apply an unsteady Reynolds-averaged modelling of the equation system, together with a turbulence model. Presently a standard high-Reynolds ( k -ε ) turbulence model is used for this purpose. With these assumptions, the model equations take the following form: where turbulent viscosity is obtained by 2 T k C ν ν = ε . The Reynolds stress tensor is obtained by: while the eddy diffusivity appearing in the energy equation is / , T T T T γ = ν σ σ , being the turbulent Prandtl number. The production and stratification terms in the turbulence model are given by: Conventional constants for the high-Reynolds ( k -ε ) model are given by: The value for C 3 is more uncertain. In the present study we assume that 3 max( ,0) C G G θ θ = , i.e. C 3 = 0 in stably stratified flows, elsewhere C 3 = 1 (Rodi 1987).

Algebraic formulation
The governing equations presented in the last section are discretized in space by the use of a finite element method. The time integration is performed using a semi-implicit two-level formulation. In the compressed form the discretized equation system can be written in the following way: Here M represents the mass matrix, A is the sum of diffusion and advection matrices (subscripts indicating the actual variable), C in the gradient matrix, and s (with subscripts) represents source terms. The implicit parameter a may be chosen in the interval (1/2, 1), and A* indicates the advection velocity taken at n u +a . The variables (u, p) are redefined here as nodal vectors for velocity and pressure, and i s represents nodal vectors for each of the scalar variables ( , , K θ ε).

Segregated implicit projection algorithm
In this study we use a segregated, implicit projection method that is non-iterative, with corrections within each time-step. This algorithm has several features in common with the SIMPLER-like pressure projection method described in (Haroutunian et al. 1997), but instead of iterations it applies corrections similar to the PISO method. The algorithm is achieved by completing the following steps: 1. Predict the pressure field via pseudo-velocity prediction from the system:ˆ; where L represents a discretized Laplacian operator. 2. Compute the velocity field from the semi-implicit momentum equation: where: where: 4. Correct the velocity and pressure fields by use of the projection step: The advection matrix A* indicates use of a timecentered advection velocity, which may be calculated as 1 (1 ) n a n n u u u and the mass matrix may be lumped.
More details, description and validation of results can be found in (Utnes 2007a, b).

Safety analysis
The simplest meteorological variable considered most important for aviation safety is called the F-factor or wind shear and the so referred to turbulence, represented by 1/3 ε . These quantities are given by Equations (22) and (23): 1/3 1/2 3/2 1/3 1/2 1/3 Here c is the flypath, g is the acceleration due to gravity, u is the wind component along the y path, w is the vertical wind component, ε is the turbulent dissipation, K is turbulent kinetic energy, l t turbulent length scale and l f is the minimum response distance for landing configuration and is of the order of ~500 m, which corresponds to a time interval of about t = O (7 s). An average over this distance is indicated by the overline.
Coefficient C µ is given by 0.09 C µ ≈ . A good review of this theory can be found in the paper by (Eidsvik et al. 2004). Prevalence of the two conditions F < -0.1 and 1/3 2/3 1 0.5m sε > corresponds to severe turbulence for commercial aircrafts and represents potential danger (Clark et al. 1997). These conditions are easily met when 1 3 K ms -> .

Simulation setup
As already explained in the introduction, terrain induced turbulence is known to be problematic near the airport, especially in case of the north-westerly wind.
To simulate these effects, we use the terrain data for the area of interest, and generate unstructured hexahedral mesh for it, in such a way that the vertical mesh lines are normal to the ground surface. More details regarding the domain size, terrain, mesh, boundary conditions and simulations are provided in the following sub-sections.

Domain size and mesh
Terrain data (Fig. 3) for the area in the vicinity of the airport is available at a resolution of 100 m. This put a constraint on the finest resolution that could be used for the study as any finer resolution would imply an interpolation of the terrain data. We conducted two sets of simulations to understand the impact of resolutions on the result. One set of simulation was conducted with 150×150 cells (Fig. 4(a)) and another with 300×300 cells (Fig. 4(b)), in the horizontal directions. The horizontal expanse of the domain was 30×30 km with the airport occupying almost the central position, as shown in figure 2. The mesh was intentionally refined in the vicinity of the airport (in horizontal directions giving a resolution of about 50 m for a 300×300 mesh). Since terrain data is not available at this resolution, interpolation was unavoidable. However, with the terrain close to the airport being relatively flat the error induced due to such interpolation is expected to be insignificant. In the vertical direction, level 41 with a stretching factor of 1.1 was used to discretize a vertical expanse of 3000 m. This resulted in a vertical resolution of 3 m near the ground and 300 m near the top of the domain.

Boundary conditions
In complicated mountainous terrain it is generally difficult to specify a realistic inlet profile. Therefore, a standard profile for wind speed and turbulent kinetic energy was used to specify the boundary conditions and to initialize the domain. The profiles for the wind speed u 0 (z) and turbulent kinetic energy K 0 (z) are given by: where u * , z 0 , z and D represent friction velocity, surface roughness, height above the ground surface and boundary layer thickness, respectively. The so-called wake function W is defined by the formula W(z/D) = (A -1)(z/D) -A/2(z/D) 2 such that W(1) = 1. The values of coefficients are: k = 0.42 and A = 4.0. Synoptic wind (on a mesoscale) U is given by U = u 0 (D). In the present simulations we have used (z 0 ; D; U) = (0.3 m; 1500 m; 20 m/s), where the friction velocity is u * ≈ 0.9 m/s and wind speed at 10 m above the ground is u 0 ≈ 7:5 m/s. A surface roughness value of 0.001 has been used for the sea surface. Along with the magnitude, direction of the synoptic wind is also specified. Several simulations were conducted for different wind directions. The convention to specify the wind direction a , used in this report, is demonstrated in figure  5(a). It should be noted here that the meteorology and aviation community use a slightly different convention, as shown in figure 5(b). All the simulations conducted and presented in this report correspond to neutral stratification, and hence, the results are scalable.

Results
In this section we present the effects of resolution on the flow characteristics. We also present the results for the simulations conducted for four different wind directions ( a = 260°; 290°; 320°; 340°). Since these results are also directed towards the aviation community, we have not only presented the velocity vectors and TKE ( K ) contours on horizontal and vertical planes but also on a conical section with a cone angle of 3.5° (Fig. 3). The velocity field and the TKE contours on the surface of this cone are then projected on to a 2-D plane for the purpose of presentation. In the following subsections we present our observations and their explanations.

Effects of changing the resolution
In figure 6 we present the effect of different resolutions on the K and velocity field along the surface of a cone containing the path of any plane. It is very clear that increasing the resolution leads to a greater penetration of turbulent kinetic energy in the vertical direction. The two circles represent the height above the ground. Compared to figure 6(a), in figure 6(b) an increase of K is observed. A possible cause of this may be the resolution of small flow structures contributing to higher turbulent kinetic energy.

Comparison of flow characteristics for different wind directions
In figure 7 we present the contours of turbulent kinetic energy at 10 m above the ground. The observation here is somewhat intuitive. We observe regions of high turbulent kinetic energy mostly on the leeward side of the hills. This is expected. Also because most of the high mountains are located on the north-west side of the airport and are aligned in the north-east to southwest direction, maximum turbulence is generated when the wind direction is normal to these mountains, i.e. a ~ 320°. To investigate the flow characteristics even more closely, we have plotted the velocity field in 3-D in figure  9. We can see that when the wind direction is aligned along the direction of the valleys ( a = 260°) the flow is channelled through these long valleys ( Fig. 9(a)). Flow acceleration takes place in these valleys but the wind rarely gets a vertical component. However, things changes dramatically when the wind direction is perpendicular to these valleys ( a = 320°). It becomes difficult for the flow to be channelled through the valleys, as a result of which the air ascends along the slope of the mountains resulting in a substantial increase of vertical wind component. This can be clearly seen in figure  9(b). This ascent, as it appears from figure 10, results in the alignment of regions of high turbulent kinetic energy along the valleys just above the hills. Pilots report that the wind direction at which they experience most vigorous turbulence is a = 320°. Figures  7 and 8 appear to confirm such claims. We therefore, chose this particular wind direction to carry the investigation further. As already explained in Section a, value of K > 3.0 is considered unsafe in the aviation industry. We therefore plot this particular quantity ( K ) in figure 10 to identify the zones of maximum risk. The figure consists of four sub-figures with different values of K . It is clear from the figures that for the particular case (synoptic wind speed of 20 m/s) there is hardly any unsafe zone for flight operation in the domain of investigation as K does not exceed 3 m/s. However, it should be stressed that this is true only for a synoptic wind speed of 20 m/s. Since the results are scalable (for neutral stratification), the quantity K will also scale accordingly and might cross the threshold. Figure 10, thus, identifies the high risk zones that can pose danger to the flight operations in the region. Figure 11(a) shows a prediction of turbulence ( K ) in a vertical section of the computational domain, parallel and passing through the runway at Kristiansand Airport. Figure 11(b) is corresponds to a zoomed in view near the airport.

Conclusions
In this report we have presented the findings of our study on the terrain induced turbulence in and around Kristiansand airport. Several simulations were conducted for different wind directions and the flow behaviours were explained. Below we enumerate the most important findings from the study.
1. Effects of resolution: it appears that the resolution of 100×100 m is a little coarse to resolve the flow separations around the hills. However, for conducting high resolution simulations high resolution terrain data is required. Pilots' reports state that vigorous turbulence is experienced up to a height of 1500 m.
Our simulations failed to predict high turbulence above 700 m. However, a comparison of results for two different resolutions indicated that increasing the resolution led to greater penetration of turbulent kinetic energy in the vertical direction. It would be interesting to investigate this issue with finer resolution terrain data. Such data will also facilitate the study of terrain modi-fications in the region to decrease risks associated with terrain induced turbulence.
2. Confirmation of the pilots' reports: in spite of the underestimation of the turbulent kinetic energy high up in the atmosphere, the simulations confirm the claims made in the pilot reports. Maximum turbulence intensity was indeed predicted for the north-westerly wind. This has been attributed to the vortices forming on the leeward side of the mountain.
3. Identification of the high risk zones: the simulations resulted in the identification of the zones where the turbulence may reach a maximum. Although for a free wind speed of 20 m/s there is no potential danger, it can significantly increase when the free wind speed increases significantly reaching a value of (say) 30 m/s, in which case the magnitude of K will become greater than 3 (Fig. 10(c); (d)) posing potential danger to flight operations in the region.
It can be said that the investigation has given a better insight into the terrain induced turbulence at the specific airport. However, it should be remembered that there are some unresolved issues (underprediction of turbulence at higher altitudes), which can be related to prevailing meteorological conditions (like fog, cloud formation, wind gusts, etc.) during the period when the reports were compiled. Modelling of such phenomenon is beyond the scope of our modelling tool (SIMRA) at the moment; therefore, we used the pilots' reports to confirm our prediction about the general nature of turbulence experienced at the airport. Had the aim been to reproduce the observations in the pilots' report, we would have conducted the simulations with more realistic boundary conditions, based on the prevailing weather condition during the period of interest.