Meso-Scale CFD Simulation for Wind Resources : A Case Study of Complex Mountainous Terrain

Land based site selection for a wind farm has some challenging criteria, namely, cost for electricity generation and distribution, acquiring ownership of the site, potential barriers from various laws and permits, security concerns, access issues, feasibility of accommodation, etc. However, wind resource assessment deems the first criterion to rule out a site before other criteria can play roles in the selection process. In this paper, a Computational Fluid Dynamics (CFD) study has been performed on a complex mountainous terrain near a shore in the west coast of the US to assess the wind resource in order to spot potential suitable sites for wind turbines. Average wind speed at a height of 10 m at the centre (44°22′12.0′ ′ N, 123°59′24.0′ ′ W) of the chosen region under study has been compared with the simulated data for validation. Results from the study, which yields a continuous map of flow field variables, have revealed much more detailed features than the available state-wise wind maps. For example, it has revealed as high as 147% variation in wind speeds and 438% in wind power, making it possible to choose suitable sites without the need for, or perhaps in advance of, expensive direct measurements. This type of analysis may help in preliminary assessments and expedite the site selection process.


Introduction
Energy demand is increasing around the world due to industrialization and development, and fossil fuels have been the major player to provide this energy.However, its negative impact on environment and thus climate change has been a topic of discussion for decades.Burning fossil fuel produces CO 2 , CO, and other toxic gases that affect health and environment such as causing respiratory diseases and global warming, respectively.
Renewable Energy Sources (RES) are not only an alternative to conventional fossil fuel but can also provide power for the electricity deprived regions (where electricity supply from national grid is limited for installation costs and maintenance) as in many rural and hilly areas in many countries and areas hit by natural disasters.
The clean energy market has been expanding gradually and has managed to reduce profits of fossil fuel companies.It has grown nearly six times compared to its business in 2004 by totalling about $329 billion in investment [1].Awareness of governments around the world and support from corporations and the public alike have made it possible to cut down clean energy cost over years.Government incentives have played a major role in this scenario.For example, large-scale solar energy and wind energy cost have reduced by 60 percent and 40 percent, respectively, since 2008 [2].People living in rural and mountainous areas not connected with the central grid have pursued solar energy as a quicker and less expensive option.Corporations have come forward by cutting down heat-trapping gas emissions and putting a price on carbon [1].
Wind is one of the RES that is always available free and wind power plants do not emit any air or water but noise pollution.Cumulative wind capacity around the world has been reported in 2018 to have seen a growth from 23,900 megawatts in 2001 to 539,581 megawatts in 2017 with an average annual increase of 25 percent over the past decade [3].
The crucial factor in developing any large-or small-scale wind farms is to assess wind resources accurately and characterize them properly.An assessment of wind resource in 2009 [4] has reported that the land-based wind energy has a potential of generating power of 10,500 GW (gigawatt) at 80 m and 12,000 GW at 100 m with a capacity factor of at least 30 percent in the USA.
The National Renewable Energy Laboratory (NREL) has been a dedicated US government organization to develop and manage wind farms through its research and data collection, where US based characterized wind map are available.Figure 1 from NREL is a presentation of a characterized wind map based on the annual average wind speed at 80 m height above sea-level.This map has been constructed with a low resolution (2.5 km ≈ 1.24274 mile) data recording and then the data have been interpolated to a finer resolution [5].This type of categorized wind map helps find target regions for consideration within seconds to build large-scale wind farms.State-wise wind maps are also featured in NREL website.Figure 2 presents a similar wind map of California, USA based on the annual average speed at 30 m.It has been developed similarly with recorded data of spatial resolution of 2.0 km and then filled up with finer resolution via interpolation algorithm.This map definitely gives a good colored characterized wind speed distribution over the state.However, it still lacks the necessary resolution needed to accurately spot potential regions with suitable wind speed if a small region (a city, a county, or a village) is picked for consideration of small scale installation.To find wind resource sites, GIS (Geographic Information System) based MCDM (Multi Criteria Decision Method) is a very useful tool if numerous data points (e.g., terrain, wind resource, etc.) are available [7,8].Categorized wind map may not be available for all locations especially for remote and complex mountainous terrains.Computational Fluid Dynamics (CFD) software packages have been employed to analyze micro-scale and meso-scale sites [9,10] to characterize winds and model wind farms [11,12].However, repetitive wind analysis may become a common practice in the near future as the world is leaning towards renewable energy and urbanization keeps changing the demography and city structures every few years.To adapt the situation, small and large scale power plants (solar, wind and other renewable energy sources) are likely to be constructed in a large number.A detailed high resolution wind map is undeniably an essential tool in this scenario.
This paper uses OpenFOAM (Open source Field Operation And Manipulation), an open-source CFD software package, to generate various fluid flow property maps (e.g., velocity, pressure, turbulence kinetic energy, etc.) over a meso-scale 17.5 km × 11.12 km complex mountainous terrain in Oregon as a case study in search of potential small-scale wind turbine (especially Vertical Axis Wind Turbine) candidate sites where high resolution wind data are not currently available.Use of meso-scale simulations such as this may not provide the granularity of micro-scale simulations but do improve upon, with respect to resolution, the currently available wind maps from NREL, for example.It may be possible to use such a tool to improve resolution on national-scale wind maps.

Computational Domain
The domain has been modeled as a rectangular prism, the six sides of which correspond to the phenomena in the real site, i.e., the top as sky, the bottom as ground, the inlet as the face the wind is entering through, the outlet as the face the wind is leaving through, and the sides as the faces which experience no variation in their normal direction (properties remain same on both sides) (Figure 3).The origin "O" and the Cartesian coordinate system shows the orientation of the domain.The scale of simulation investigated herein can be assumed of having no Coriolis force effects [13].The equations of motion are governed by steady state Reynolds-Average Navier-Stokes (RANS) equations in the domain.The fluid is wind at the standard atmospheric temperature which can be assumed to be incompressible.Other thermo-physical effects (e.g., gravity) are negligible.The chosen site is designated by latitude ψ and longitude θ and elevation h, which is the height above sea-level.This information needs to be transformed into points (or vertices) in three dimensional space (x, y, z) for blockMesh to use in formation of cell blocks.The coordinate of each point in the computational domain is bounded by (ψ min , ψ max ) and (θ min , θ max ).Following the origin in Figure 3, mapping of vertices can be described as follows: where R = 6.371 × 10 6 m is the radius of the earth.The domain is almost rectangular and the top boundary is chosen at a height, H, far enough away from the ground not to affect the solution.
Blocks are formed from the coordinates of points (or vertices) with height from z = h to z = H while the subdivided cell volumes along the height vary with an expansion ratio, ∆ r , as follows: where ∆ h is the end cell thickness and ∆ s is the starting cell thickness in the z-direction.Intermediate cells have linearly interpolated thickness.

Governing Equations
The equation of conservation of mass is written as and the equation of conservation of momentum can be written as where u i is the time-averaged fluid velocity and two tensors-S ji and τ ij are the mean strain rate tensor and the specific Reynolds stress tensor, respectively.τ ij is derived from the fluctuation of velocity.These tensors can be expressed as:

Turbulence Model
A turbulence model needs to be adopted to address this specific Reynolds stress tensor, τ ij .Several models have been developed which base on Boussinesq hypothesis that assumes τ ij varies linearly with the mean velocity gradients and thus introduces two other quantities in the expression of τ ij : (a) turbulence kinetic energy, κ; and (b) kinetic eddy viscosity, ν T .
In this paper, the κ − model has been used as used in other larger scale simulations such as wind turbine simulations [14,15].The standard κ − model is a two-equation model as it requires solution of two transport equations to achieve closure and deliver accurate enough predictions.It brings in another term, , the rate of dissipation of κ per unit mass due to viscous stresses to establish relationship between ν T and κ.
where C µ = 0.09 is a constant.The stead-state transport equations for κ and can be written as: where C 1 = 1.44,C 2 = 1.92, σ κ = 1.0, and σ = 1.3 are typical recommended values for this turbulence model.

Solution Algorithm
An algorithm is needed to solve the system of equations with proper boundary and initial conditions.SIMPLE (Semi-Implicit Method for Pressure Linked Equations), a widely used Navier-Stokes Equations' solver, has been chosen for the current steady state fluid flow case in question.OpenFOAM has named the solver "simpleFoam".

Boundary Conditions
The boundary conditions mimic the phenomena of the chosen site such as the air flows in the x-direction as it enters the region of interest but its speed varies with height.Therefore, a power-law relationship describing the speed of air with height is employed to model the wind speed profile.This power-law relationship is frequently used to measure power densities of wind and macroscopic wind velocity [16,17].The relationship can be shown as: where a is a constant equal to 1 7 or 0.143, u r is a reference wind speed measured at a height, h, and z is the elevation in meters.The chosen value as 1/7 is aimed at coinciding with flat near-ocean region [18] as the inlet.Obviously, choosing other values of "a", such as 0.1, will definitely change the inlet wind profile.However, the results, which rely more on flow field comparisons than on absolute value, should remain valid regardless of (reasonable) choice of "a".All other boundaries but the inlet are set to a zero velocity-gradient boundary condition that acts on the normal direction of boundary surfaces: where n j is the outward pointing surface normal.Pressure is forced to have a value of zero at the outlet.
On all other boundaries, a zero gradient boundary condition normal to the surface is invoked.
An important point is to be noted at this stage is that the unit of the pressure is m 2 /s 2 due to its division by density, ρ, that facilitates the solution of the governing equations in the incompressible solvers.To initialize, pressure p = 0 is chosen arbitrarily to be used everywhere since only the relative pressure difference is taken into account.Turbulence kinetic energy, k, can be expressed as follows: where I is the turbulence intensity.The other boundaries have similar zero-gradient boundary condition normal to their surfaces.
Rate of turbulence dissipation can be written as: where l is the turbulence length scale that approximates to the average of the length of cell volumes in ground discretization.Other boundaries are again enforced to have zero-gradient boundary condition normal to their surfaces.

Simulation Software
The governing equations along with sufficiently described boundary and initial conditions are solved using OpenFOAM [19], a finite volume based C++ toolbox capable of solving a wide range of multi-physics problems including CFD.The domain is descritized by blockMesh, a native meshing tool that comes with the OpenFOAM package.The vertices of grid elements are used to form blocks and then cell volumes are created from the subdivision of those blocks.

Case Study and Validation
The chosen site is of length 17.507 km (8.75 miles) and width 11.12 km (6.67 miles) approximately and is shown in Figure 4 and the elevation data have been collected from Google [20].The ground of domain in Figure 3 shows the topography of the chosen mountainous terrain.The redder locations have greater elevation while the bluer areas are closer to sea-level.
The air properties have been taken as the standard values at room temperature and at sea-level such as air density, ρ = 1.2 kg/m 3 , and dynamic viscosity, µ = 1.8 ×10 −5 kg/ms.The turbulence intensity, I, is taken as 1% or 0.01.The real time velocities at six points (plus four points around the validation point) on the area under investigation have been collected from OpenWeatherMap [22] for 24 h starting on 24 June 24 2017 at 2:35 a.m.(UTC).The reference velocity, u r , is calculated as 4.1 m/s.This value comes as an average of the recorded velocities at six observation points, as shown in Figure 5.The distribution of velocities in Figure 5b is of Point 1 only.The other points show similar trends.The inlet reference velocity is taken as the average of speeds at all six points over the period of 24 h.Despite the availability of standard registered data (i.e., wind maps), preference was given to record the wind speeds to observe if there was any considerable difference and the calculated velocity is (Figure 6) reasonably close to the standard.The computational domain is made up of 200 × 125 nodes.The average side length of cell volumes gives an estimated turbulence length scale of 44.12 m.The ceiling of the domain in Figure 3 is at 1000 m away from sea-level.A grid independence test (Figure 7) has been conducted to ensure that results do not vary with grid size in the domain.Several simulations have been run using various mesh sizes and velocity in the x-direction has been compared with the recorded data as a measurement of success.The simulation result with 1.92 million cells seems to match with the collected data perfectly, but u x does not remain consistent with increased number of cells (for example, 3.2 million).Therefore, it can be inferred that this particular result is very prone to errors because of comparatively low resolution of the computational domain and thus should be avoided.The mesh size that is chosen eventually to conduct the study is of 1 × 10 7 cells.
The wind is never unidirectional and uniform.However, the recorded velocities have orientations around 350 • ∼ 360 • .Therefore, the wind velocities have been just time-averaged to determine the inlet velocity and it has been considered to be unidirectional for simplification.In addition, the record does not include wind speed in z-direction as it is assumed to be very low and cannot impact this study to a considerable level.
The OpenWeatherMap [22] collects data from airport, public and private weather stations and they are not necessarily calibrated for the same height which may explain some of the variations in Figure 5b.However, most of these are calibrated to a height to 10 m.Data for each station have not been tracked and recorded separately.Therefore, chances of stations' data for a height other than 10 m remain, although these stations can be considered small in number.However, averaging data reduces this error to a negligible level.Each block in the domain has 100 cells in z-direction with a cell expansion ratio ∆ r = 750 and two cell volumes in x and y directions.
The simulation was then run using a linear 2nd order, cell-centered finite volume spatial discretization with a relative residual tolerance less than 10 −4 for u, v, w, p, and κ solvers.

Results and Discussion
The objective of this study is to generate a wind map of a region to spot locations as potential wind resource sites.It also varies to some extent, if not significantly, at different elevations depending on the topographical features of the area.Therefore, wind maps of Oregon, USA (where the area under investigation is located) at 30 m and 80 m have been the bases for comparison with the flow properties of air from the simulation.

Wind Flow Properties at 30 m
Figure 8 presents a characterized wind map of Oregon, USA from NREL [23].The red dot on the map (shown with the arrow for better visualization) is the region of interest for the current study.It can be easily concluded from the map that most of the state has an annual average wind speed of 4-5 m/s.In addition, the resolution has been interpolated from a 2.0 km grid to finer grids.Therefore, trying to characterize wind on any particular small area using this wind map is challenging and results will be prone to errors.The flow velocity in x-direction (u x ,) from left to right is positive.Figure 9 shows u x distribution of the region at 30 m above ground.The inflow comes from the left and being a flatter region; this portion of the region has mostly a similar velocity as the inlet.However, there are considerable number of locations where it drops too low (even negative) and, at some other locations, rises to high.Any blockage in the flow path in the form of hills or cliffs can make the flow stop and turn, which reduces the speed.On the other hand, air flows faster through a narrow passage (such as valley) or over water body as the friction of water is less than that in solid surfaces.
The flow velocity components are combined together to get its magnitude that has been produced in Figure 10a.NREL suggests that at 30 m height a velocity magnitude of 4 m/s or higher is suitable for wind farm [24].The pie chart of area distribution (Figure 10b) shows that almost 67% area of the region meets the minimum target.Other categorized wind speed zone can by pointed out and its area share can be known from Figure 10.
The inflow is completely unidirectional (u y = 0, u z = 0 at inlet).Thus, u y and u z components come from wind shear, turbulence, wake propagation, change in flow path due to blockages, etc. Figure 11a reveals spots with positive and negative u y where bluish and reddish regions are two extremes respectively.This gives the idea about where the turning of flow direction is likely to occur and thus topographical feature of those spots can be guessed.Likewise, Figure 11b reveals spots where vertical turnings of flow occur as it passes over and through the terrain.Slopes will likely produce positive u z while cliffs or blocks will generate negative u z .Turbulence causes fatigue loading on turbine blades, i.e., time varying torque and thrust that reduces turbine's longevity.Turbulent flow has random wind velocities and is quantified by turbulent kinetic energy (TKE), κ. Figure 12 shows distribution of κ at 30 m.As shown in the figure, the vicinity near the inlet has higher κ because of the viscous effect from solid surface friction and thick turbulent boundary layer.As the flow nears the mountainous terrain, the eddies decay as they approach solid surfaces of the terrain.Wind power density is to | u| 3 and thus depends on the magnitude of velocity and is considered as an important criterion for assessment of a region for potential wind resource.Here, wind power density has been calculated as 1  2 ρ| u| 3 per m 2 area.Figure 13a has been produced to observe wind power distribution over the region.The pie chart in Figure 13b shows percentage of area share of wind power densities.The maximum value that can be found as 294.1 W/m 2 which is 5.4 times higher than the minimum required average wind power density (54.675W/m 2 ).

Wind Flow Properties at 80 m
Oregon's wind map from NREL [25] is shown again in Figure 14 at 80 m and, at this elevation, 6.5 m/s or higher wind speed is suitable for wind farms [26].The blue dot shows the location that our study has focused on.Again, this map gives an overall idea about the average wind speed throughout the state but lacks the details to use it for meso-or micro-level location analysis for potential wind resource.
Since the inflow is in +x direction, this velocity component, i.e., u x , is larger than the other two and its distribution is shown in Figure 15.At this elevation, areas with higher velocities (blue regions) are larger than those at 30 m and areas with lower velocities (red regions) are lesser than those at 30 m.It is obvious that blockage in the flow path decrease with the increase of height.
Similarly, u y (Figure 16a) and u z (Figure 16b) maps show how these velocity components change with elevation.Both bluish (+u y ) and reddish (-u y ) regions reduce in size in comparison with the map for 30 m.It is true for u z as well.The reason definitely can be attributed to less deflection and thus less blockage in the flow path.
Velocity magnitude and its pie chart have been produced in Figure 17.A barchart (Figure 18) summarizes the comparison of wind speeds at two different elevations.The percentage of areas with higher speeds are larger at 80 m.The velocity found at 80 m is also higher than that at 30 m.
Turbulent kinetic energy (TKE) map in Figure 19 for the height of 80 m shows that regions with higher turbulence have increased in size in comparison with the map at 30 m.The region near the inlet experiences more turbulence as wind shear affects flow at different layers within thick turbulent boundary layer.The mountainous terrain has also higher κ since there are less solid surfaces for the eddies to decay, rather wind shear contributes to the eddies.
Figure 20 shows the wind power distribution on the region at 80 m.The pie chart (Figure 20b) shows percentage of area share of different wind power densities.To compare it with wind power density at 30 m, Figure 21 has been added.The magnitude of maximum wind power density is higher at 80 m.Regions with higher wind power have more percentage in comparison to those at 30 m. NREL recommends a site should have as low as 185 W/m 2 wind power density to qualify while we get a wind power density as high as 305 W/m 2 , which is 1.65 times higher.

Conclusions
Site selection for wind farms involves a few important criteria such as cost for power generation and distribution, ownership of the site, limitations due to laws and permits, security, accessibility, feasibility, and so on but the major one is the availability of wind source.This study shows a CFD approach to identify possible site selection.
The current study has adopted an approach to model a meso-scale simulation with an open-source CFD software, OpenFOAM, and has utilized the standard κ − turbulence model.By using terrain elevation data from Google and the wind profile (power law), steady-state solution has been simulated and then flow field variables have been analyzed.To better understand the effect of altitude on the variables better, results have been analyzed at two different altitudes (30 m and 80 m).
The results demonstrate that the state-wise wind map may be useful for large wind energy projects, but may not be as useful for small scale installations.The area investigated herein is rectangular and spans 17.507 km in length and 11.12 km in width.Due to the resolution of typical wind maps, it is difficult to gauge how small-scale topographical variations can contribute to wind energy availability.
By observing the wind maps of Oregon from NREL, the region under the current study is assumed to have an annual average wind speed of about 4.5∼5.0m/s at 30 m and 6.0∼6.5 m/s at 80 m , however, our study reveals much more variation in wind speed and wind energy density.Our results indicate several locations within this region that are unsuitable for wind energy conversion due to relatively low wind speeds.Conversely, several other locations are highly recommended as candidate sites where average wind speeds are up to 147% higher than the NREL average at 30 m and 75% higher at

Figure 1 .
Figure 1.Annual average wind speed of the US at 80 m [5].

Figure 3 .
Figure 3. Computational domain and the colored ground is showing elevation above sea-level.

Figure 7 .
Figure 7. Grid independence test and validation.

Figure 8 .
Figure 8.Average annual wind speed of Oregon at 30 m.

Figure 9 .
Figure 9. Wind speed in x-direction at 30 m.
(a) Wind power density.(b) Wind power density's relative distribution.

Figure 14 .
Figure 14.Average annual wind speed of Oregon at 80 m.

Figure 15 .
Figure 15.Wind speed in x-direction at 80 m.

Figure 16 .
Figure 16.Velocity components in yand z-direction at 80 m.

Figure
Figure Wind speed's relative distribution by areas at two elevations.
(a) Wind power density.(b) Wind power density's relative distribution.

Figure 21 .
Figure 21.Wind power density's relative distribution by areas at elevations.