Effect of different source terms and inflow direction in atmospheric boundary modeling over the complex terrain site of Perdigão
- 1Environmental and Applied Fluid Dynamics Department, TS1von Karman Institute for Fluid Dynamics, Waterloosesteenweg 72, 1640 Sint-Genesius-Rode, Belgium
- 2Department of Mechanical Engineering, Université de Sherbrooke, 2500 Boulevard de l’Université, Sherbrooke, QC J1K 2R1, Canada
- 3Department of Mechanical and Structural Engineering and Materials Science, University of Stavanger, PB 8600, Stavanger 4036, Norway
- These authors contributed equally to this work.
Correspondence: Kartik Venkatraman (email@example.com)
Assessing wind conditions in complex terrain requires computational fluid dynamics (CFD) simulations incorporating an accurate parameterization of forest canopy effects and Coriolis effects. This study investigates how incorporating source terms such as the presence of trees and the Coriolis force can improve flow predictions. Furthermore, the study examines the impact of using different sets of atmospheric boundary layer inflow profiles, including idealized profiles with a logarithmic velocity profile, and a set of fully developed profiles from a pressure-driven precursor simulation. A three-dimensional steady Reynolds-averaged Navier–Stokes (RANS) equations model is set up using OpenFOAM to simulate the flow over a complex terrain site comprising two parallel ridges near Perdigão, Portugal. A 7.5 km×7.5 km terrain of the Perdigão site is constructed from elevation data centered around a 100 m met-mast located on the southwest ridge. A 30 min averaged stationary period is simulated, which corresponds to near-neutral conditions at met-mast Tower 20 located at the southwest ridge. The period corresponds to the wind coming from southwest at 231∘ at 100 m height above ground at Tower 20. Five case setups are simulated using a combination of different source terms, turbulence models and inflow profiles. The prediction capability of these models is analyzed for different groups of towers on the southwest ridge and, on the towers further downstream inside the valley, on the northeast ridge. Including a canopy model improves predictions close to the ground for most of the towers on the southwest ridge and inside the valley. Large uncertainties are seen in field measurement data inside the valley, which is a recirculation zone, and large prediction errors are seen in the wind velocity, wind direction and turbulent kinetic profiles for most of the models. The predictions on the northeast ridge are dependent on the extent of recirculation predicted inside the valley. The inflow wind direction plays an important role in wind profile predictions.
Lack of terrain availability in flat terrain pushes wind farm developers to look for alternatives along complex terrain sites. Flat terrain availability is becoming scarce and 70 % of the Earth’s surface is complex terrain, which presents a large potential for wind energy harvesting. Winds in complex terrains are governed by the surface properties of the flow (land class/roughness length) and the local elevation such as hills, ridges and mountains (Emeis, 2018). Local features such as ridges or canyons can also be advantageous for wind energy harvesting, due to the creation of local flow accelerations on top of ridges and flow channeling through valleys. However, the wind fields depend on either the local pressure or temperature gradient. Such flows are also dominated by strong thermal stratification effects and are influenced by the presence of forest canopies. Complex terrain sites remain very challenging areas to consider for wind farm siting and require extensive validation of modeling tools in representative environments. Wind farm modeling for complex terrains requires a more advanced approach than commonly used cost-effective linearized models such as WAsP (Jørgensen et al., 2005), which cannot handle complex phenomena (i.e., flow separation) to ensure reliable and accurate results. Computational fluid dynamics (CFD) tools are increasingly used to predict flows over complex terrain sites to account for such phenomena and provide more accurate wind resource predictions (Blocken, 2014). However, improving numerical modeling tools for complex terrains demands accounting for various microscale phenomena such as flow recirculation and the forest canopy effect. Such microscale flow features significantly impact the local wind resource assessment and wind turbine loading.
Initial studies in literature accounting for the effects of forest canopy were performed over forested flat terrains. Brower (2012) showed that the presence of a forest canopy increases the modeling uncertainty by up to a factor of five irrespective of the chosen modeling approach. Finnigan (2000) parameterized the effect of foliage and forest canopy by accounting for the drag force in the momentum equation. The effect of the canopy on turbulence was further taken into account by Sogachev and Panfyorov (2006) and Sogachev (2009) by adding additional source terms for turbulence kinetic energy and turbulence dissipation rate. Desmond et al. (2017) modeled both forest canopy and buoyancy effects by modeling them using source and sink terms and showed that thermal stratification plays an important role in determining the flow over canopies. The canopy model is typically implemented by specifying the tree height and leaf area density (LAD), which represents the area of the leaf and branches of trees. These parameters are typically defined from field measurements or LiDAR point clouds generated from field measurement campaigns and aerial surveys Queck et al. (2012).
Validation of numerical models for complex terrains requires field measurements. One of the earliest field campaigns was at Askervein hill Salmon et al. (1988), a smooth isolated hill. Similarly, the Bolund hill campaign was another campaign over an isolated hill in Denmark as detailed by Bechmann et al. (2009), with data under conditions that could be classified as neutral. However, no canopy is seen on these sites. The Alaiz field campaign was performed over a large-scale homogeneous mountain valley topography in Spain with forested regions Rodrigo et al. (2021). Chavez et al. (2014) showed that using a canopy model improved the accuracy of velocity and turbulence predictions for simulations over Alaiz.
More recently, a field experiment campaign has been carried out over the complex terrain site at Perdigão, Portugal, by Fernando et al. (2018), which is a double ridge with a forested canopy experiencing different thermal stratification conditions. Significant field data are available for further studies, and these are elaborated in Sect. 2. Several numerical studies have been performed for validation with field data. Laginha Palma et al. (2020) studied the choice of using the appropriate grid size for numerical modeling. Quimbayo et al. (2022) implemented a forest parameterization in the weather research and forecasting model (WRF) with an average tree height of 30 m, resulting in an improvement in the prediction of near-surface wind speeds. This tree height was chosen for practical reasons related to the model setup and is not representative of actual tree height. There is a gap in the existing literature on microscale simulations using a canopy model on a heterogeneous forested complex terrain with a rich dataset such as Perdigão.
Objectives and outline
This paper aims to evaluate the impact of different physical source terms and turbulence models when performing CFD simulations in complex terrain and compare simulation prediction with field measurement profiles at the different groups of towers located at the Perdigão test site. The influence of the following parameters are investigated: canopy effects, Coriolis force and the effect of using two different sets of inflow profiles. One is an idealized set with a log-law velocity profile, and one is a set of fully developed profiles based on a precursor simulation. Furthermore, the influence of wind direction on the prediction of wind profiles is studied for the different case setups.
An improved understanding of the importance of these phenomena will enable the development of more efficient and reliable tools to perform wind simulations in complex terrain. This paper is structured as follows: a detailed description of the methodology is provided in Sect. 2, covering the computational domain and meshing, the numerical settings, and different modeling capabilities that have been added as source terms to the conservation equations. The results are presented and discussed in Sect. 3, in terms of different groups of towers. Section 4 is dedicated to concluding remarks and future perspectives. Finally, additional details on the grid spacing, inflow profiles and wind direction study are provided in the Appendix.
An intensive observation period was carried out from 1 May to 15 June 2017 at the double ridge site of Perdigão, Portugal, by a consortium of American & European universities (Fernando et al., 2018). The location of the different met-masts of heights 60 and 100 m used in this study is shown in Fig. 2, which could be grouped by their location on top of the southwest ridge, inside the valley or at top of the northeast ridge. From the analysis of the wind rose, the predominant directions at Tower 20 and Tower 29 are perpendicular to the ridges, i.e., northeast and southwest directions, as shown in Fig. 1, while at Tower 25, flow turning is seen and the dominant wind direction is from the south.
A stationary period was found on the date of 4 May 2017 for the 30 min averaged time interval of 22:00–22:30 tilt corrected high-frequency dataset from NCAR-EOL, based on the conditions at Tower 20 on top of the ridge, which corresponded to a bulk Richardson number of approximately −0.03 which qualifies as near-neutral condition. The bulk Richardson number as defined by Kaimal and Finnigan (1994) is calculated using Eq. (1). However, the flow conditions are non-neutral at other met-mast locations.
where θ is the potential temperature, g is the acceleration due to gravity and U and V are the wind velocity components at the reference height. The same period was simulated by Laginha Palma et al. (2020) based on the stationary periods predicted by Carvalho (2019).
The terrain is a 7.5 km×7.5 km square centralized around a 100 m met-mast located on the southwest ridge. Figure 3 presents the computational domain and the dimensions are listed on the right side of the figure.
A cylindrical computation domain was developed, which provides the flexibility to simulate wind from any wind direction. The authors have successfully applied this approach in previous studies pertaining to urban areas (Hågbo et al., 2021; Hågbo and Giljarhus, 2022). A smoothing region from complex to flat terrain was applied towards the outer boundaries, with a minimum radial distance of 15 Δht. Several best practice guidelines have been formulated for grid generation for simulating complex terrain sites, such as by Sørensen et al. (2012) and Laginha Palma et al. (2020), and have been closely followed. The height of the domain is set to 10 times the difference in the elevation height of the terrain, Δht, as recommended by Sørensen et al. (2012) when simulating wind flow over complex terrain sites. It consists of approximately 88 million cells (“Fine” case as detailed in the Appendix) and is produced with
terrainBlockMesher developed by Schmidt et al. (2012), capable of generating structured meshes over complex terrain exclusively consisting of hexahedra cells. The
terrainBlockMesher tool uses a blending function to smooth the transition from the terrain patch to the outer cylindrical block. Around 50 radial block cells are defined and a radial grading factor is used to enable a stretching in the horizontal direction to cluster cells across the center of the domain and expand towards the boundaries. In terms of the number of cells per main direction () the mesh comprises of across the terrain patch. The vertical mesh resolution is 24.5 m with uniform stretching applied across the entire domain. The minimal mesh height Δz next to the ground is close to 3 m. The average value of the wall y+ is around 32 000.
The mesh also follows the recommendation for having at least three cells from the ground to the height at the first sampling point for comparison with field measurements. The horizontal mesh resolution over the terrain was set close to 12.5 m. It satisfies the minimum resolution of 40 m recommended specifically for the Perdigão site by Laginha Palma et al. (2020). Figure 4 illustrates the grid structure used inside the domain from the side. Terrestrial data were obtained from the shuttle radar topography mission (SRTM) database of the Perdigão field experiment (Fernando et al., 2018). A grid sensitivity study was performed and presented in the Appendix. The results obtained with five different grids of increasing mesh resolution show significant differences in the wind profiles especially close to the ground and inside the valley. Hence, a further grid spacing study has to be done as a part of future work and the results presented could change if a finer resolution is used.
All simulations have been conducted using the OpenFOAM (version 2012) toolbox. The simulations are steady state and performed by solving the incompressible, three-dimensional steady Reynolds-averaged Navier–Stokes (RANS) equations with the finite volume method. Second-order discretization schemes were used for spatial discretization. The initial iterative convergence criteria were that the scaled residuals should drop four orders of magnitude for all flow variables as per the best practice guidelines (BPGs).
Two steady-state solvers for turbulent flow of incompressible fluids have been used:
buoyantBoussinesqSimpleFoam (BBSF) and
simpleFoam (SF). All thermal effects are neglected in the simulations using both the solvers, such that the atmospheric stability of the simulated atmosphere is always neutral. In these simulations, the air density is assumed constant, and the gravitational force is neglected. The BBSF solver is capable of simulating the effect of buoyancy forces, but these source terms are set to zero for the present neutral case. However, in addition to solving the continuity equation and the momentum equation, which is solved using SF, the energy equation is also included, allowing for the modeling of non-neutral atmospheric conditions.
3.1 Inlet profiles and boundary conditions
Two sets of fixed inlet profiles have been used, both forming a homogeneous atmospheric boundary layer (ABL), either idealized with a logarithmic velocity profile or fully developed profiles obtained from a precursor simulation.
The idealized profiles provide the turbulence quantities and a log-law type ground-normal wind velocity based on the generalization and modification of the well-used set of equations from Richards and Hoaxey (1993). This modification has been implemented by Yang et al. (2009) and uses a more mathematically consistent formulation allowing the turbulent kinetic energy to vary with height. Also, using a log-law wind velocity profile and the associated turbulent inflow conditions are considered only to be valid in the atmospheric surface layer, which can be roughly estimated as 10 % of the atmospheric boundary layer Temel et al. (2018). The idealized velocity profile has no wind veer, meaning it does not have a changing wind direction with height.
An alternative to using the idealized profiles is using fully developed profiles obtained from a one-dimensional precursor simulation following the strategies of Koblitz (2013), Alletto et al. (2018). Cyclic conditions are applied on the sides, and the simulations are run for sufficient iterations to allow the development of the profile. The setups are identical to the setups used in the final (successor) simulations except for the computational domain and mesh, the cyclic conditions applied and the number of iterations. It enables the precursor simulations to produce inlet profiles valid for various conditions, including non-neutral atmospheric conditions. Adjusting the parameters of these source terms is an efficient way to obtain the desired inflow. The developed velocity profile has wind veer, meaning the wind turns clockwise with height caused by the balance between friction, the Coriolis force and the pressure gradient force, and that Perdigão is located in the Northern Hemisphere. Also, the roughness length, z0, was adjusted in the precursor simulations to obtain the right turbulent kinetic energy (TKE) at reference height. The inflow profiles are added in the Appendix.
An overview of the boundary conditions (BCs) used is presented in Table 1.
Robin boundary conditions that act as an inlet or outlet are used on the sides of the domain to simulate wind from any horizontal direction. The direction of the flux automatically determines the inlet and outlet regions. The inlet profiles are fixed to form a homogeneous atmospheric boundary layer (ABL), either idealized profiles representing neutral atmospheric conditions or fully developed profiles obtained from precursor simulation based on the work of van der Laan et al. (2020). Outlet conditions were used with constant pressure in the simulations using the SF solver, or equally, constant (p−(ρgh)) in the simulations using the BBSF solver. Zero gradient was set for the remaining variables. The direction of the flux automatically determines the inlet and outlet regions. No-slip conditions with wall functions were used for the terrain. Specifically, a rough wall condition was applied (Hargreaves and Wright, 2007). The roughness length was set to z0=0.02 m, based on average values of a roughness length map provided in the database of the Perdigão field experiment (Fernando et al., 2018).
3.2 Source terms
Three different source terms have been added to the momentum equation for some simulation cases: tree canopy as a porous medium, Coriolis force and pressure gradient force. Additionally, a limiter in the turbulence dissipation equation was added to one of the cases.
The canopy is modeled as a porous medium based on the work of Costa (2007). A drag term is added to the momentum equation. The canopy source term was utilized with the
simpleFoam (SF) flow solver. The porosity model is based on the
powerLawLopesdaCosta model implemented in OpenFOAM. It is a variant of the power law porosity model with a spatially varying drag coefficient. This source term is applied to the momentum equation to reproduce the momentum dissipation that the trees and their foliage should produce in the flow. The following parameters are used: porosity surface area per unit volume or the leaf area density (Σ=1.0), drag coefficient (Cd=0.25) and the power law model exponent coefficient (C1=2.0) for the Eq. (2). A mean tree height of 3 m is chosen based on the mean canopy height at different locations provided by Vasiljevic et al. (2017). This was applied all over the entire flow domain. A cell set was used to select a volume of cells 3 m above the ground (which forms a canopy zone where the source terms are activated). An overall tree height of 3 m is also supported by analyzing a LiDAR point cloud acquired from the database of the Perdigão field experiment (Fernando et al., 2018) as shown in Fig. 5. The source term in the canopy model for velocity is as shown below:
where Σ is the leaf area density and Cd is the plant canopy drag coefficient. As highlighted in Lalic and Mihailovic (2004), forests have a higher foliage density at the half top than at the bottom zones. However, for the present case, simulations are performed with a uniform leaf area density.
The Coriolis force source term is included as a momentum source term (U and V momentum equations). The Coriolis force is calculated based on the planetary rotational period (Ω=24 h) and the latitude (λ) for Perdigão (39.68∘ N), where fc=2Ωsin (λ). Similar to Koblitz (2013), the Coriolis force in the vertical direction is neglected since it is negligible compared to the gravitational acceleration.
A pressure gradient forcing term is also included in the simulation cases using the BBSF solver. In setups using this solver, the wind is driven by the balance of the lateral pressure gradient force and the Coriolis force, which are added in the momentum equation as shown in equations provided in the Appendix. This balance is called geostrophic balance, and the resulting wind is a geostrophic wind, which in the Northern Hemisphere goes counterclockwise around low-pressure systems.
The level of turbulent kinetic energy is controlled in the domain by the use of a parameter called the maximum turbulent mixing length scale, used to limit the production and dissipation terms of the turbulence model. Ambient source terms were added to avoid zero turbulence values above the ABL similar to those mentioned in van der Laan et al. (2021a) and applied to the entire domain. The values were set to kAmb=0.001 and as summarized in Table 3. The geostrophic wind is denoted by G and the Coriolis frequency parameter is denoted by fc. The governing equations for the BBSF solver and the different source terms are further detailed in Appendix.
3.3 Simulation cases
The setup for all the simulations was identical except for certain parameters. The solver used is either
simpleFoam (SF) or
The inlet profiles are either idealized with a logarithmic velocity profile representing neutral atmospheric conditions or fully developed profiles obtained from precursor simulation. Source terms are added, including canopy (tree height set to 3 m) and the Coriolis force.
Three turbulence models have been used: the standard k−ϵ turbulence model (SKE), a modified k−ϵ turbulence model (MKE) and the k−ω turbulence model (KO). The MKE model is applied in the simulation case where trees/vegetation (canopy) has been added as a porous medium with the SF solver. The turbulent constants are modified with Cμ=0.033. This model was tuned using experimental data and LES (large eddy simulation) for a homogeneous forest (Costa, 2007). The boundary conditions are overall the same but had to be adjusted to account for the different solvers and inlet profiles used. An overview of all the simulation cases is provided in Table 2. The set of chosen turbulence constants is shown in Table 3.
Numerical convergence difficulties can be encountered when using the global turbulence length scale limiter (van der Laan et al., 2021b) when applied to the present case for simulations over the complex terrain of Perdigão. These convergence difficulties can be particularly challenging when setting the values of the maximum limiting length scale to low values. Here, the model is trying to restrict the turbulence scale, but physically the large length scales of turbulence are produced from hills or features that are in the order of the maximum value that is set by the turbulence model.
The results are discussed based on three groups of towers of interest: on top of the southwest ridge, inside the valley and on top of the northeast ridge. The inlet profiles for all the simulations are calibrated to match the measured velocity magnitude and direction at 100 m at tower 20, which corresponds to an elevation of 573 m. Different metrics are analyzed for all the simulation models at different towers and are further explained. The root mean square error (RMSE) between the averaged profiles of the measured data and the simulated results is computed and presented in Table 4. A hit rate metric for a given model is defined as the number of predictions within one standard deviation over each height of the field measurement. The relative error for the turbulent kinetic energy and wind direction is shown in Tables 5 and 6 respectively. The relative error is given in percentage and is the difference between the simulated and field measurement divided by the field measurement value. The error bars represent the standard deviation of the mean measurements.
4.1 Model prediction on SW ridge
The towers on the SW ridge as shown in Fig. 2 are Tower 20, 34 and 37. This is a region of flow acceleration. The inflow profile is calibrated to match the velocity at 100 m height for Tower 20 as shown in the wind velocity profiles shown in Fig. 6a at 100 m. All models except the canopy model (SF3) show an overprediction of the velocity profiles close to the ground at Tower 20. The canopy model (SF3) underpredicts the velocity close to the ground but accurately captures the shape of the profile. This prediction could be improved by considering a non-uniform tree height and removing the canopy source terms on top of the ridges. On the other hand, at Tower 34, shown in Fig. 6b, the canopy model (SF3) shows a good match with the field measurement and on Tower 37, seen in Fig. 6c, shows a slight underprediction. In models that do not account for a canopy (SF1, SF2, SF4 and BBSF1), a speed-up is seen close to the ground with large RMSE errors as shown in Table 4. The results suggest that the inclusion of the canopy is a good choice for predicting of velocity profiles. However, the canopy parameters need to be adapted based on the forest point cloud data. BBSF1 model (k−ϵ Lim) shows a higher speed-up close to the ground at most of the towers. The turbulent kinetic energy profiles are shown in Fig. 7, showing a good prediction for all the models except the canopy model (S3), which produces excessive turbulence levels compared to field measurements. Figure 8 shows the predicted wind direction profiles at the three towers showing a good prediction for most of the models within one standard deviation of the measurements. As highlighted in the methodology section, there is a large uncertainty in wind direction seen in the field experiment that varies with location and height, which is further investigated in Sect. 4.4 for two additional wind directions. This uncertainty highlights a limitation in the models' abilities to represent actual conditions for such a complex site. In terms of the turbulence model, the SF1 model (k−ϵ) is seen to produce a slightly greater speed-up compared to the SF2 model (k−ω). Still, overall the prediction capability for the wind speed, turbulent kinetic energy and wind speed is similar.
4.2 Model prediction inside the valley
The towers of interest inside the valley are Towers 25, 7, 27 and 22 as seen earlier in Fig. 2. The valley is a region of strong flow separation and flow recirculation. The predicted velocity profiles are shown in Figs. 9 and 10 with different predictions of the extent of recirculation zones, resulting in various shapes of the wind profiles. At Towers 25, 22 and 27 the canopy model significantly underpredicts the velocity profiles. The SF2 (canopy), SF1 (k−ϵ) and SF2 (k−ω) show the best prediction at Tower 7, which is located at the lowest altitude among the other towers. Overall the SF1 (k−ϵ) model shows the best prediction at all the towers inside the valley.
Figures 11 and 12 show a comparison of turbulent kinetic energy profiles between the present model predictions and the field measurements. The measured values of turbulent kinetic energy are high inside the valley at all met-mast locations, indicative of flow mixing and high turbulence. However, all the models seem to strongly underpredict the turbulent kinetic energy with large relative errors as seen in Table 5, except for Tower 22, which is located the closest to the southwest ridge.
The wind direction profiles are shown in Figs. 13 and 14 respectively. A sudden shift in wind direction, indicative of flow separation and recirculation, is seen in Fig. 13. Table 6 provides the hit rate for different models at all the towers. Figure 15 shows the velocity contour for the predicted recirculation zone inside the valley for the SF3 model (canopy) and the BBSF1 model (k−ϵ Lim). The SF3 model (canopy) produces a single large recirculation zone compared to the other models, and a smaller and double recirculation zone is seen for the BBSF1 model (k−ϵ Lim). The large recirculation zone results in a significant underprediction of the velocity profiles seen earlier in Figs. 9 and 10 for the canopy model. Menke et al. (2019) also investigated the flow recirculation zones around the Perdigão valley under different atmospheric stability conditions.
4.3 Model prediction on NE ridge
The wind velocity profiles for Tower 29 and Tower 10 on the northeast ridge are shown in Fig. 16. This area on the top of a ridge is a region of flow acceleration. It is just downstream from the recirculation zone inside the valley, which makes the predictions quite challenging. The SF3 model (canopy) is seen to underpredict the velocity profiles significantly on top of this ridge. All other models appear to provide a good prediction of the velocity profile around one standard deviation of the field measurements. The region near Tower 10 comprises a mountain gap (Vassallo et al., 2020), where the local flow features could play an important role. The SF3 (Coriolis) and the BBSF1 models predict the profiles close to the field measurement. Excessive turbulence levels are seen close to the ground in the field measurements shown in Fig. 17. All models appear to underpredict the field measurements of turbulence close to the ground, and the SF3 (canopy) shows the closest match. The BBSF1 model (k−ϵ Lim) underpredicts the turbulence levels due to the use of the maximum length scale limiter. The wind direction profiles at the towers are shown in Fig. 18. As seen for the TKE profiles, a large uncertainty is seen in the wind direction on this ridge; this is expected as the wind is from the southwest direction. Most of the model predictions fall within one standard deviation of the measurements at Tower 29. As seen in Fig. 15, the prediction of wind direction at Tower 10 is dependent on the extent of recirculation zone and flow reattachment location on top of the ridge.
4.4 Influence of wind direction
The inflow direction plays an essential role in the wind predictions for complex terrains. A wind direction standard deviation of around 7∘ is seen in the field measurements at Tower 20 (tse04). Simulations have been performed for two additional wind directions for each model to look at the differences observed between the three incoming wind directions.
The results for the model predictions using the SF1 (k−ϵ) model are presented in this section, while the results obtained using the other models of this study are shown in the Appendix. For each model, the wind profiles at three different met-masts, corresponding to Tower 20 (tse04) on the southwest ridge, Tower 25 (tse09) inside the valley and Tower 29 (tse13) on the northeast ridge, are shown.
All simulations are calibrated for Tower 20 as seen in Fig. 19 with the same inflow velocity at the reference height of 100 m on the tower and only changing the wind direction. The predicted profiles inside the valley and on the northeast ridge appear to vary quite significantly with the inflow wind direction as seen in Figs. 20 and 21. Streamline trajectories for the inflow with three different wind directions at the Towers are shown in Fig. 22. Wind passing through Tower 29 with inlet wind from 234.5∘ exhibits a very different trajectory compared to the other inlet wind directions and the set of trajectories passing through the mast upstream, Tower 20. The wind here was initially deflected off from the NE ridge and led into the valley by a channeling effect. As a result, the wind speed decreases significantly before passing Tower 29, before re-gaining in intensity downstream of the ridge as the wind accelerates downhill.
A Reynolds-averaged Navier–Stokes (RANS) model is set up using OpenFOAM (version 2012) to simulate a 30 min averaged stationary period corresponding to near-neutral conditions at met-mast Tower 20 located at the southwest ridge. In that period and for that tower, the wind comes from the southwest at 231∘ at 100 m height above ground. Five different models are simulated comprising of different source terms to account for the effects of the canopy, the Coriolis force and pressure gradient force, and two different inflow profiles. One is an idealized set with a log-law velocity profile, and one is a set of fully developed profiles based on a precursor simulation. Based on the flow topology, the predicted profiles are analyzed in terms of the different groups of towers on top of the ridges and inside the valley. The complex terrain site of Perdigão represents a large spatial variability of forest canopy and surface elevation, which contribute to variable flow topology at different met-masts. The key conclusions for different groups of towers are summarized as follows:
For the towers on the southwest ridge. The region at the southwest ridge is a zone of flow acceleration at the first oncoming ridge downstream of the inlet for wind coming from the southwest. The inflow profiles are calibrated to closely match the wind speed and direction at Tower 20. Using a canopy model (SF3) decreases the velocities near the surface and is a closer match with field data at Tower 34 and 37. Other models overpredict the velocity profile close to the ground. However, the canopy parameters need to be tuned as the surface heterogeneity is not considered, as the prediction accuracy varies at the different locations along the ridge.
For the towers inside the valley. The valley is a zone of flow recirculation and comprises lower velocities and higher variability, which remains challenging for prediction models. Moreover, large uncertainties are seen in the wind velocity, wind direction and turbulent kinetic energy profiles for the field measurements. The prediction capabilities of the models vary with the location of the tower inside the valley. At all towers inside the valley, the SF1 model (k−ϵ) provides the best prediction for wind velocity. Most models show large relative errors in wind speed and turbulent kinetic energy profiles, especially close to the ground.
For the towers on the northeast ridge. The region at the northeast ridge is a zone of flow acceleration downstream of the recirculation zone from the valley. The canopy model (SF3) provides a strong underprediction, while all other models provide a prediction within one standard deviation of the field measurements. Predicting the extent of the recirculation inside the valley and the re-attachment location plays a key role in the prediction profiles on the northeast ridge. Significant turbulence is seen close to the ground in the field measurements and is underpredicted by most models.
Influence of wind direction. A significant difference in the wind profiles is seen using different inflow directions. The extent of the recirculation zone and the re-attachment downstream of the valley is different due to different trajectories taken by the inflow wind profiles coming from the southwest ridge. These uncertainties also depend on the turbulence model and source terms utilized.
The choice of the best-performing turbulence model is inconclusive in terms of overall prediction capability for different parts of the terrain. In the future, the surface heterogeneity of the canopy could be modeled based on the surface roughness length map at Perdigão. It would be more correct to use different patches with different heights, as seen in the forest point cloud data. Simulations using non-uniform leaf area density could also be performed. Finally, an additional grid spacing study has to be done with a finer resolution as a part of future work.
The influence of grid spacing has been studied using five different grids as shown in Table A1. The “Very fine” grid mesh is the maximum number of cells that could be generated due to computational constraints. The structured grids are generated using the
terrainBlockMesher tool Schmidt et al. (2012), which interpolates the SRTM terrain data and creates the terrain patch, which is blended into a cylindrical domain. All simulations are performed using the simpleFoam (SF1) solver with the k−ϵ turbulence model setup. The simulation inflow profiles for all the meshes are calibrated to match the field experiment at 100 m of Tower 20 (tse04) on the southwest ridge.
Significant differences are seen between the “Very coarse” and “Medium” grids especially close to the ground near Tower 20 (Calibration tower) on the southwest ridge seen in Fig. A1, inside the valley as seen in Fig. A2, and at Tower 29 close to 100 m on the northeast ridge seen in Fig. A3. Larger differences are seen close to the ground as the topology of the terrain is further resolved near the surface upon grid refinement. Increasing the number of cells refines and slightly modifies the surface mesh close to the ground. Furthermore, a small change in the prediction of the extent of the recirculation zone could have a significant change and uncertainty in the predictions inside the valley and on top of the northeast ridge. The order of these differences is similar to the order of the differences between the different model setups shown in the paper. This suggests that grid spacing is an important parameter for flow prediction around a complex terrain, and a grid-independent solution could be challenging to achieve based on the complexity of the flow topology. Further studies have to be performed with a finer grid spacing.
As shown in Fig. A4, two sets of inflow profiles are utilized: an idealized set of inlet profiles, including a logarithmic velocity profile, and a fully developed profile using a precursor driven by a Coriolis force and a pressure gradient force. The inflow velocity profiles are calibrated to reach the desired inflow conditions at the met-mast Tower 20 at a height of 100 m, for a time period identified as neutral based on the bulk Richardson number. The wind velocity magnitude profiles are close to logarithmic. For the wind direction, the idealized profile fixes a uniform wind direction, but the precursor has a source term to account for the Coriolis effect, so a wind veer is seen over the entire height. Finally, for the TKE, a profile is set in the idealized case while for the precursor developed profile it is limited by the maximum mixing length scale, increasing the roughness length z0 of the precursor to 6 (m), while also tuning pressure gradient magnitude and direction to get the right wind speed, wind direction, and now TKE at the reference height to match the meet the average field experiment value at the calibration point.
The equations solved for the
buoyantBoussinesqSimpleFoam solver are the continuity equation (Eq. C1), momentum equation (Eq. C2) and turbulence transport (Eqs. C4, C5), and temperature (Eq. C3) as described by Alletto et al. (2018). The pressure gradient (πi) drives the momentum equation, and the source terms for Coriolis and canopy effects are included in the momentum equation. The buoyancy terms are only added in the turbulence transport equations but are set to zero for the present neutral case. The equation for turbulent dissipation rate contains the maximum length turbulence scale limiter (lmax) to modify the mixing-length scale estimations for setting different atmospheric stabilities. This description has been added to the Appendix of the paper.
D1 SF2 k−ω
The influence of wind direction change using the SF2 (k−ω) model at different towers is shown in Figs. D1, D2 and D3 respectively. The results appear to be similar to the SF1 (k−ϵ) model. However, in comparison a lesser difference is seen in the profiles between the 231∘ and the 234.5∘ degree cases at Towers 25 and 29.
D2 SF3 canopy
The influence of wind direction change using the SF3 (canopy) model at different towers is shown in Figs. D4, D5 and D6 respectively. A much larger difference in wind profiles is seen at Tower 29 on the northeast ridge compared to the SF1 (k−ϵ) model, indicating larger uncertainties with wind direction. Presently, the canopy is modeled using a uniform tree height, and perhaps there would be higher uncertainties with a non-uniform canopy across the domain; the path taken by the wind for different inflows would vary significantly. A much different recirculation zone is developed, as seen in the wind direction profiles.
D3 SF4 Coriolis
The influence of wind direction change using the SF3 (Coriolis) source term model at different towers is shown in Figs. D7, D8 and D9 respectively. Interestingly, for this case, the flow difference between the 231∘ and the 234.5∘ cases at Towers 25 and 29 is relatively small compared to the SF1 k−ϵ case. With the Coriolis source term, the flow turning due to the Coriolis effect is accounted for, and hence the inflow wind for calibration is set accordingly to obtain the required wind direction at the calibration mast. Consequently, the path taken by the wind across the terrain is different.
Figures D10, D11 and D12 show the wind profiles for the BBSF simulation model with the precursor inflow. Here the inflow from wind direction 234.5∘ shows a larger difference between the bother profiles even at the calibration mast and a decrease in wind velocity at Tower 29 on the northeast Ridge, similar to the other models. The wind direction standard deviation bounds in Fig. D12 appear to be more closely predicted by this model.
KV conceived, coordinated and was responsible for both the work, the paper writing and the review process. TOH carried out the numerical simulations, paper writing and preparation of figures, and attended in the review process. SB helped with conception, paper writing and review. KEG helped with model setup and review.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors acknowledge the European Commission for its financial support through the project H2020-MSCA-ITN-2019 zEPHYR (Grant Agreement No. 860101). Perdigão data were provided by NCAR/EOL under the sponsorship of the National Science Foundation (https://data.eol.ucar.edu/, last access: 9 January 2023). The authors would like to thank Paul van der Laan for his useful insights, review and comments which helped improve the quality of the paper.
This research has been supported by the European Commission, Horizon 2020 Framework Programme (zEPHYR (grant no. 860101)).
This paper was edited by Johan Meyers and reviewed by Paul van der Laan and one anonymous referee.
Alletto, M., Radi, A., Adib, J., Langner, J., Peralta, C., Altmikus, A., and Letzel, M.: E-Wind: Steady state CFD approach for stratified flows used for site assessment at Enercon, J. Phys.-Conf. Ser., 1037, 072020, https://doi.org/10.1088/1742-6596/1037/7/072020, 2018. a, b
Bechmann, A., Berg, J., Courtney, M., Jørgensen, H., Mann, J., and Sørensen, N.: The Bolund Experiment: Blind Comparison of Models for Wind in Complex Terrain (Invited), American Geophysical Union (AGU) Fall Meeting Abstracts, https://ui.adsabs.harvard.edu/abs/2009AGUFM.A33H..08B/abstract (last access: 10 January 2023), 2009. a
Brower, M.: Findings of investigations into under-performing sites, in: Proceedings of EWEA Technology Workshop, Lyon, France, European Wind Energy Association (EWEA), http://www.ewea.org/events/workshops/past-workshops/analysis-of-operating-wind-farms/programme-proceedings/ (last access: 10 January 2023), 2012. a
Carvalho, J. P. D. B.: Analysis of stationary periods during the Perdigão 2017 campaign, Master's thesis, University of Porto, Portugal, https://repositorio-aberto.up.pt/bitstream/10216/122036/2/348346.pdf (last access: 10 January 2023), 2019. a
Chavez, R., Rodrigo, J., and Gancarski, P.: Modelling of atmospheric boundary-layer flow in complex terrain with different forest parameterizations, J. Phys.-Conf. Ser., 524, 012119, https://doi.org/10.1088/1742-6596/524/1/012119, 2014. a
Costa, J. C. L.: Atmospheric flow over forested and non-forested complex terrain, PhD thesis, Faculdade de Engenharia da Universidade do Porto, https://repositorio-aberto.up.pt/bitstream/10216/11375/2/Texto integral.pdf (last access: 10 January 2023), 2007. a, b
Desmond, C., Watson, S., and Hancock, P.: Modelling the wind energy resources in complex terrain and atmospheres. Numerical simulation and wind tunnel investigation of non-neutral forest canopy flows, J. Wind Eng. Ind. Aerod., 166, 48–60, https://doi.org/10.1016/j.jweia.2017.03.014, 2017. a
Emeis, S.: Wind Energy Meteorology, vol. 2 of Atmospheric Physics for Wind Power Generation, Springer International Publishing, Belgium, https://doi.org/10.1007/978-3-642-30523-8, ISBN 978-3-642-30523-8, 2018. a
Fernando, H., Mann, J., Laginha Palma, J., Lundquist, J., Barthelmie, R., Belo-Pereira, M., Brown, W., Chow, F., Gerz, T., Hocut, C., Klein, P., Leo, L., Matos, J., Oncley, S., Pryor, S., Bariteau, L., Bell, T., Bodini, N., Carney, M., and Wang, Y.: The Perdigão: Peering into Microscale Details of Mountain Winds, B. Am. Meteorol. Soc., 100, 799–819, https://doi.org/10.1175/BAMS-D-17-0227.1, 2018. a, b, c, d, e
Hågbo, T.-O. and Giljarhus, K. E. T.: Pedestrian wind comfort assessment using CFD simulations with varying number of wind directions, Frontiers in Built Environment, Frontiers in Built Environment 8, 110, https://doi.org/10.3389/fbuil.2022.858067, 2022. a
Hågbo, T.-O., Giljarhus, K. E. T., and Hjertager, B. H.: Influence of geometry acquisition method on pedestrian wind simulations, J. Wind Eng. Ind. Aerod., 215, 104665, https://doi.org/10.1016/j.jweia.2021.104665, 2021. a
Hargreaves, D. and Wright, N. G.: On the use of the k–ε model in commercial CFD software to model the neutral atmospheric boundary layer, J. Wind Eng. Ind. Aerod., 95, 355–369, 2007. a
Jørgensen, H., Nielsen, M., Barthelmie, R., and Mortensen, N.: Modelling offshore wind resources and wind conditions, in: Proceedings (CD-ROM) Copenhagen Offshore Wind, 26–28 October 2005, https://backend.orbit.dtu.dk/ws/portalfiles/portal/107919793/Modelling_offshore_wind_resources_and_wind_conditions.pdf (last access: 10 January 2023), 2005. a
Kaimal, J. C. and Finnigan, J. J.: Atmospheric boundary layer flows: Their structure and measurement, Oxford University Press, https://doi.org/10.1093/oso/9780195062397.001.0001, ISBN 9780195062397, 1994. a
Palma, J. M. L. M., Silva, C. A. M., Gomes, V. C., Silva Lopes, A., Simões, T., Costa, P., and Batista, V. T. P.: The digital terrain model in the computational modelling of the flow over the Perdigão site: the appropriate grid size, Wind Energ. Sci., 5, 1469–1485, https://doi.org/10.5194/wes-5-1469-2020, 2020. a, b, c, d
Lalic, B. and Mihailovic, D.: An Empirical Relation Describing Leaf-Area Density inside the Forest for Environmental Modeling, J. Appl. Meteorol., 43, 641–645, https://doi.org/10.1175/1520-0450(2004)043<0641:AERDLD>2.0.CO;2, 2004. a
Menke, R., Vasiljević, N., Mann, J., and Lundquist, J. K.: Characterization of flow recirculation zones at the Perdigão site using multi-lidar measurements, Atmos. Chem. Phys., 19, 2713–2723, https://doi.org/10.5194/acp-19-2713-2019, 2019. a
Queck, R., Bienert, A., Maas, H.-G., Harmansa, S., Goldberg, V., and Bernhofer, C.: Wind fields in heterogeneous conifer canopies: Parameterisation of momentum absorption using high-resolution 3D vegetation scans, European J. Forest Res., 131, 165–176, https://doi.org/10.1007/s10342-011-0550-0, 2012. a
Quimbayo-Duarte, J., Wagner, J., Wildmann, N., Gerz, T., and Schmidli, J.: Evaluation of a forest parameterization to improve boundary layer flow simulations over complex terrain. A case study using WRF-LES V4.0.1, Geosci. Model Dev., 15, 5195–5209, https://doi.org/10.5194/gmd-15-5195-2022, 2022. a
Richards, P. and Hoaxey, R.: Appropriate boundary conditions for computational wind engineering models using the k−ε turbulence model, in: Computational Wind Engineering 1, edited by: Murakami, S., Elsevier, Oxford, 145–153, https://doi.org/10.1016/B978-0-444-81688-7.50018-8, 1993. a
Rodrigo, J., Santos, P., Chavez, R., Avila, M., Cavar, D., Lehmkuhl, O., Owen, H., Li, R., and Tromeur, E.: The ALEX17 diurnal cycles in complex terrain benchmark, J. Phys.-Conf. Ser., 1934, 012002, https://doi.org/10.1088/1742-6596/1934/1/012002, 2021. a
Salmon, J., Bowen, A., Hoff, A., Johnson, R., Mickle, R., Taylor, P., Tetzlaff, G., and Walmsley, J.: The Askervein Hill Project: Mean wind variations at fixed heights above ground, Bound.-Lay. Meteorol., 43, 247–271, https://doi.org/10.1007/BF00128406, 1988. a
Schmidt, J., Peralta, C., and Stoevesandt, B.: Automated generation of structured meshes for wind energy applications, Open Source CFD International Conference, London,United Kingdom, https://www.researchgate.net/publication/271754239_Automated_generation_of_structured_meshes_for_wind_energy_applications (last access: 10 January 2023), 2012. a, b
Sørensen, N. N., Bechmann, A., Réthoré, P.-E., Cavar, D., Kelly, M. C., and Troen, I.: How fine is fine enough when doing CFD terrain simulations, in: EWEA 2012-European Wind Energy Conference & Exhibition, European Wind Energy Association (EWEA), 1167–1172, https://backend.orbit.dtu.dk/ws/portalfiles/portal/7951156/How_fine_is_fine.pdf (last access: 10 January 2023), 2012. a, b
van der Laan, M. P., Kelly, M., Floors, R., and Peña, A.: Rossby number similarity of an atmospheric RANS model using limited-length-scale turbulence closures extended to unstable stratification, Wind Energ. Sci., 5, 355–374, https://doi.org/10.5194/wes-5-355-2020, 2020. a
van der Laan, M. P., Kelly, M., and Baungaard, M.: A pressure-driven atmospheric boundary layer model satisfying Rossby and Reynolds number similarity, Wind Energ. Sci., 6, 777–790, https://doi.org/10.5194/wes-6-777-2021, 2021b. a
Vasiljević, N., L. M. Palma, J. M., Angelou, N., Carlos Matos, J., Menke, R., Lea, G., Mann, J., Courtney, M., Frölen Ribeiro, L., and M. G. C. Gomes, V. M.: Perdigão 2015: methodology for atmospheric multi-Doppler lidar experiments, Atmos. Meas. Tech., 10, 3463–3483, https://doi.org/10.5194/amt-10-3463-2017, 2017. a
Vassallo, D., Krishnamurthy, R., Menke, R., and Fernando, H.: Observations of Stably Stratified Flow through a Microscale Gap, J. Atmos. Sci., 78, 189–208, https://doi.org/10.1175/JAS-D-20-0087.1, 2020. a
Yang, Y., Gu, M., Chen, S., and Jin, X.: New inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary layer in computational wind engineering, J. Wind Eng. Ind. Aerod., 97, 88–95, https://doi.org/10.1016/j.jweia.2008.12.001, 2009. a