A CFD Model for Spatial Extrapolation of Wind Field over Complex Terrain—Wi.Sp.Ex

Abstract

High-resolution wind datasets are crucial for ultra-short-term wind forecasting. Penetration of WT installations near urban areas that are constantly changing will motivate researchers to understand how to adapt their models to terrain changes to reduce forecasting errors. Although CFD modelling is not widely used for ultra-short-term forecasting purposes, it can overcome such difficulties. In this research, we will spatially extrapolate vertical profile LIDAR wind measurements into a 3D wind velocity field over a large and relatively complex terrain with the use of stationary CFD simulations. The extrapolated field is validated with measurements at a hub height of three WTs located in the area. The accuracy of the model increases with height because of the terrain anomalies and turbulence effects. The maximum MAE of wind velocity at WT hub height is 0.81 m/s, and MAPE is 7.98%. Our model remains accurate even with great simplifications and scarce measurements for the complex terrain conditions of our case study. The models’ performance under such circumstances establishes it as a promising tool for the evolution of ultra-short-term forecasting as well as for the evaluation of new WT installations by providing valuable data for all models.

1. Introduction

Wind energy forecasting enables energy management techniques to evolve, optimise the usage and storage of all available forms of energy [1], and act fast in extreme conditions. Forecasting the energy production in short time horizons in various energy sectors combined with knowledge of demand and storage capabilities is fundamental for the success of a dynamic pricing system based on an SG similar to larger time frames [2,3,4]. Thus, the need for ultra-short-term (≤15 min) wind energy forecasting has emerged. The task of predicting the E-yield of a WT is a difficult task, as wind movement in short time scales and/or over complex terrain is chaotic. Two main types of models are used for wind forecasting: statistical models and physics-based models.

Statistical models (like ANN) are capable of accurate and fast predictions that can be used for operational purposes. Some models directly predict the E-yield [5,6] instead of wind speed, such as the models presented by [7]. These models need to be specifically adjusted for different locations while training, so such a model requires large datasets. WFs that are already in operation for a long time without significant infrastructural changes can provide such data, but new wind farms need to start measurement campaigns long before the installment of the first WT. This kind of model needs to be retrained each time a change is made at the location of interest, such as new installments of WTs or buildings, which is time-consuming and may be crucial to the operation of a WF.

Physics-based models, such as CFD models [8,9,10,11] or the model presented by [12], can provide accurate estimations of complex wind fields, and their accuracy does not depend on changes in the environment, in contrast to statistical models. Wind flows are sensitive to terrain changes and obstacles [13], especially near buildings and urban areas. For these kinds of ultra-short-term predictions, models with errors below 10% are needed. Unfortunately, CFD models require great computational power. The time required to compute time-dependent simulations is usually so long that they are not widely used for ultra-short-term forecasting.

The effect of obstacles on a wind flow has been thoroughly studied by researchers [13,14]. Even relatively small terrain changes (a small shed that is constructed within the WT vicinity, Figure 1) will have a non-negligible impact on the overall wind field. Given that any 3D wind flow profile has a complex interaction with the blades, the E-yield is expected to be impacted by terrain changes, especially when elevation fluctuations are present. The importance of studying and modelling these effects and the physics that governs them is essential in the renewable wind energy sector. Our motivation is to explore the capability of CFD modelling to overcome the consequences of terrain changes. We believe more research will be motivated in that direction in the future.

This paper focuses on empowering small, low-cost wind farms to operate with existing systems while enabling research and development using physics-based synthetic datasets. We present a fast and accurate CFD model based on stationary equations that spatially extrapolates scarce wind measurements to provide high-resolution input datasets for ultra-short-term forecasting [15], even over complex terrain. The forecasting method chosen by any wind farm can be coupled with our proposed extrapolation method. This approach remains cost-effective and can leverage new measurement techniques.

Acknowledging the impact of complex terrain on LIDAR accuracy, as noted by [16,17,18], we discuss three simulations over CRES to understand model behaviour under varying conditions. The first examines model behaviour when extrapolating a 3D wind field from a 1D field with simplified inlet conditions to prioritize computational efficiency. The second focuses on extrapolating 1D LIDAR measurements to a 2D wind field, which serves as input for the third simulation, which aims to reduce errors caused by initial simplifications.

2. Data and Experimental Setup

CRES is located at a complex terrain at Lavrio, Greece (37.76, 24.06). Measurements used for the simulations came from a campaign that took place between 22 October and 24 October 2010. Wind measurements at nine different heights (40 m, 54 m, 78 m, 100 m, 120 m, 140 m, 160 m, 180 m, and 200 m) from a wind LIDAR are provided. The LIDAR wind speed measurement errors at CRES are 4% to 6%, according to [16].

Horizontal wind speed measurements at hub height from three different wind turbines located at the WF were available at the time that the LIDAR campaign took place and were provided by CRES. Wind speed on WT hub height is measured by anemometers mounted directly on the hub (cup anemometers or ultrasonic anemometers). Due to the influence of the WT itself on the airflow, the directly measured wind speed data are corrected online by the WT manufacturers based on pre-existing calibration data. The calibration data are obtained from a separate mast of equal height positioned at a distance of 2–3 times the WT diameter. The characteristics of the wind turbines can be seen in Figure 2. The location of the LIDAR, the wind turbines, and the terrain characteristics can be seen in Figure 3.

A CAD file provided by CRES was used to create the 3D representation, as seen in Figure 3, via COMSOL 5.6. The internal area of the 3D representation is the area of interest, which contains the WTs, and the LIDAR is located at the inlet plane. The external area is used solely to expand the simulation domain and extrapolate LIDAR measurements along the inlet plane. The simulation domain represents an area with a size of 2 km × 2.5 km × 0.6 km. The average elevation is 66 m, and the standard deviation is 33 m. The surface has a maximum elevation of 136 m and a minimum elevation close to sea level.

3. Model Description and Methodology

Our aim is to to develop and validate a base CFD model that can spatially extrapolate a wind field in under 5 min with an accuracy below 10%, even with extremely limited wind and temperature data. This is important for generating high-resolution wind datasets when key information is missing, as in our case study. Due to the unavailability of geometry files for the installed WTs, they were not included in the simulations, and the domain is assumed to be obstacle-free.

We conducted simulations to understand the impact of different obstacles and simplifications of our study. A non-rotating wind turbine of similar size to the ones of our case study can leave a trace up to 270 m. In Figure 4, we see stationary simulation differences in wind speed magnitude before and after importing a WT (non-rotating) as well as the iso-surfaces of 1, 2, 3, and 5 m/s. The 1 m/s iso-surface expands approximately 250 m behind the WT (six times the WT’s rotor diameter) and is 50–80 m wide even if the WT is relatively thin and not rotating. Figure 1 illustrates the effect of a small shed on a wind field over flat terrain. The influence of a similar structure located 200 m away from a wind turbine could extend in all directions and potentially reach the rotor area. These simulations expose the impact of different types of obstacles on the behaviour of wind flow over large areas.

COMSOL 5.6 was used for modelling and running the simulations with respect to the experiment’s geometrical and physics-based representation. It was selected due to its ability to create multi-physics simulations and customise the base model and the geometry with ease, without limitation, according to the needs of the user. The description of the model created for this study is separated into three main parts: the stationary hypothesis, the geometrical representation, and the stationary–turbulence CFD model.

3.1. Multi-Physics Framework

The hypothesis that stationary solutions of CFD equations can represent the flow over non-controlled and complex environments under specific conditions is made. Based on the principles of calculus used in CFD and the finite-element theory [21,22], we want to asses the assumption that:

The temporal evolution of airflow within a volume can be represented by a series of stationary states. As the volume increases, a stationary flow should effectively represent the average flow over a longer duration. Conversely, as the mean velocity of the airflow within the volume increases, a stationary flow should represent the averaged flow over a shorter duration.

The inlet is selected be perpendicular to the mean wind direction over the selected averaging period when a simulation is run, especially when other boundaries are considered as open, such as in our case study.

The k-$ϵ$ turbulent flow model [23] is used for the simulations because of the lower time consumption and negligible differences with other available turbulent models. The equations solved by this interface are the RANS equations for conservation of momentum and the continuity equation for the conservation of mass. Turbulence effects are modelled using the standard two-equation k-$ϵ$ model with realizability constraints. Wall functions are used to model the flow near walls. All parameters are set to automatic values apart from the inlet velocity field and temperature. Pressure is set to 1 atm.

Three hundred and seventy-five simulations (representing 3750 min) were run to validate the performance of the model. Because of the number of simulations run for this study and time constraints, we selected a fixed position for the inlet. Due to the setup of the WTs and the presence of open boundaries, the maximum absolute deviation of measured wind direction with respect to the normal vector to the inlet plane was approximately 30°. This way, an air parcel that starts from the side boundaries with an average absolute direction of 30° hardly reaches the WTs or the LIDAR position. The presence of WTs was excluded from the presented simulations, acknowledging the influence of this simplification on the results [11]. Small buildings in the area were also not included in the simulations. The simulations were simplified based on the hypothesis that the flow is isothermal (uniform distribution, $dT/dt=0$).

There are different types of mesh structures used in CFD models that use different geometrical shapes. COMSOL 5.6 gives the option to use tetrahedrons, hexahedrons, triangular prisms, and pyramid shapes. For the geometrical discretisation of the simulation domain, a tetrahedral non-uniform grid structure was selected, and the domain was separated into two different regions, as seen in Figure 5. The selection of unstructured tetrahedral mesh helps to represent the terrain anomalies more accurately without raising the numbers of elements needed. Automatic normal-size mesh was applied to the inner region that contains the instruments and the WTs. Thus, the mesh of the inner region had a maximum element size of 73.5 m and a minimum element size of 21.9 m. An automatic extra coarse mesh was applied to the outer region used in the EAS, resulting in a mesh with a maximum element size of 219 m and a minimum element size of 54.8 m.

The meshing of the simulation domain was primarily determined by the low resolution of the inlet data, computational time constraints, and assumptions made at the inlet of the EAS. The simulation domain was separated into two regions. The inner one, which is the location of the WF, is denser than the outer one. The outer region was created to capture the flow characteristics via the EAS. A dense mesh (even a uniform one) applied over the outer region could not produce more accurate results over the inner area than a non-uniform coarse mesh over the expanded area that gradually becomes denser as we reached the inner region. The height of the domain was set at 0.6 km according to results from [24]. The assumption that the wind speed does not change along the inlet at a given height combined with a dense grid and the complexity of the terrain slows down the convergence of the simulation and propagates the errors created by the assumption, leading to unsafe results (simulated flow experiences a strong shock at the first interaction with the terrain due to the homogeneity of initial conditions over a complex terrain).

A smooth transition from a coarser mesh dictates the simulated flow to absorb the shock produced by the surface anomalies due to the smoother terrain curvature constructed through the meshing procedure. The distance from the inlet area is more than 1 km away from the inner region, which is more than five times the height from the lowest to the highest altitude of the terrain. The top boundary of the whole domain is the coarsest one, as it is more than 450 m above the highest altitude of the surface, and the flow there tends to be laminar.

3.2. Simulations and Methodology

3.2.1. IAS and EAS

Our goal is to assist wind farms that have a limited amount of information needed for proper operation and research by providing them with a solution that can produce reasonably accurate results. In our case study, we divided the profile vertically in three sections ($z<40$ m, 40 m $\le z\le$ 200 m and $z>$ 200 m). Missing data for $z<40$ m and for $z>200$ m were the sole reason for such a division. Both IAS and EAS were based on the same assumptions and used the same data as inlet conditions. The inlet of the expanded area was located 1.3 km away from the LIDAR, which was positioned at the IAS inlet. The average distance of the WTs from the IAS inlet was 566 m, as seen in Figure 6.

A ramp function ($U_{ramp}\left(z\right)$, z$ϵ$$\mathbb{R}$) from the COMSOL 5.6 interface was used to extrapolate LIDAR wind measurements $(U_{LIDAR}(x,z,t_n)$ vertically below 40 m (parameters: location = −10, slope = 1/30, cutoff = 1; smoothing parameters: size of transition zone at start = 0, size of transition zone at cutoff = 50, number of continuous derivatives = 2). The ramp function graph is demonstrated in Figure 7. In the absence of data below or above a certain height, researchers must adapt the ramp function to their specific case. We chose the simplest approach due to this uncertainty, in order to demonstrate the performance of the base model we present in this article. Constant extrapolation was used above 200 m. The EAS goal was to physically extrapolate the wind profile vertically and horizontally over the E-IAS(/IAS) inlet plane. To generate the 2D inlet plane wind field ($U_{inlet}(x,z,t_n)$), the velocity along the X-axis of the inlet was assumed to be uniform, and the vertical velocity component was set to zero.

If $U_{LIDAR}(x,z,t_n)=(u_{x-LIDAR},u_{y-LIDAR},u_{z-LIDAR})(x,z,t_n),nϵ[0,375)=C\subset \mathbb{N}$ and $U_{inlet}(x,z,t_n)=(u_{x-inlet},u_{y-inlet},u_{z-inlet})(x,z,t_n),nϵC$, then $\forall nϵC$:

$$\begin{array}{c}\begin{array}{c}\partial U_{inlet}(x,z,t_n)/\partial x=(0,0,0)\\ u_{z-inlet}(x,z,t_n)=0\end{array}\end{array}$$

(1)

and

$$\begin{array}{c}\begin{array}{c}{U}_{inlet}(x,z,t_n)=U_{ramp}\left(z\right)\ast U_{LIDAR}(40,t_n),z<40\\ U_{inlet}(x,z,t_n)=U_{LIDAR}(z,t_n),40\le z\le 200\\ U_{inlet}(x,z,t_n)=U_{ramp}\left(z\right)\ast U_{LIDAR}(200,t_n),200<z\end{array}\end{array}$$

(2)

The assumption that the averaged wind flow would have similar behaviour if the terrain was homogeneous matches wind LIDAR measurement assumptions above complex terrain [18]. The aim of EAS simulations was to capture the wind flow characteristics at the inner area’s inlet (wind velocity and direction distribution along the 2D inlet).

3.2.2. E-IAS

After the EAS results were created, a 3D matrix that contains the simulated wind velocity vectors $U_{EAS}(x,z,t_n)=(U_{x-EAS},U_{y-EAS},U_{z-EAS})(x,z,t_n)$ at the inlet of the examined area of our study (where the LIDAR is located) was created. A deviation matrix was created for 0 m $\le z\le$ 600 m and −500 m $\le x\le$ 500 m:

$$U_{dev}(x,z,t_n)=U_{EAS}(x,z,t_n)-U_{LIDAR-sim}(z,t_n)$$

(3)

where $U_{LIDAR-sim}$ represents EAS simulated LIDAR velocity at measurement heights, extrapolated with constant extrapolation below 40 m and above 200 m. To extrapolate the LIDAR measurements along the inlet plane, we created a new matrix that contains the extrapolated inlet wind field:

$$U_{LIDAR-ex}(x,z,t_n)=U_{dev}(x,z,t_n)+U_{LIDAR\ast }(z,t_n)$$

(4)

Here, $U_{LIDAR\ast }$ represents LIDAR measurements extrapolated with constant extrapolation below 40 m and above 200 m. The vertical velocity component is set to 0 for the E-IAS inlet field as for IAS and EAS. The E-IAS simulation is run over the area of interest (internal domain in Figure 3), which is 1/4 the volume of the expanded area.

4. Results

In this section, we will discuss the results from IAS, EAS, and E-IAS and compare their performance. Simulation results are compared with WT and LIDAR data after being corrected with linear regression. Performance of state-of-the-art models in ultra-short-term forecasting are presented and indirectly compared with our extrapolation model’s performance to assess its evolution and usage capabilities.

4.1. IAS

Figure 8, Figure 9 and Figure 10 show velocity magnitude plots and their linear regression scatter plots. The MAE of wind velocity estimation is 0.88 m/s at Enercon WT’s hub height, 0.76 m/s for NEG WT’s hub height, and 0.56 m/s for Vestas WT’s hub height. Enercon’s MAPE is 8.26%, NEG’s is 7.94%, and Vestas’s is 7.09%, as seen in

Table 1. IAS-WT corrected velocity magnitude statistics.

Velocity Magnitude [m/s]
	Enercon	NEG	Vestas
MAE	0.88	0.76	0.56
MSE	1.29	0.98	0.6
RMSE	1.14	0.99	0.77
MAPE [%]	8.26	7.94	7.09

. This is expected due to their locations with respect to LIDAR positioning, as seen in Figure 3.

4.2. EAS

Table 2. EAS-LIDAR and WT corrected velocity magnitude statistics.

	LIDAR Levels									WTs
	40 m	54 m	78 m	100 m	120 m	140 m	160 m	180 m	200 m	Enercon	NEG	Vestas
MAE	0.41	0.37	0.30	0.22	0.17	0.13	0.10	0.09	0.11	0.85	0.85	0.67
MSE	0.30	0.23	0.15	0.08	0.05	0.03	0.02	0.01	0.02	1.20	1.13	0.77
RMSE	0.55	0.48	0.39	0.29	0.22	0.16	0.12	0.11	0.13	1.10	1.06	0.88
MAPE [%]	4.28	3.70	2.81	2.05	1.56	1.17	0.86	0.76	1.00	7.78	8.78	8.41

shows that MAE values for LIDAR wind velocity magnitude estimations across all available heights span from 0.09 m/s to 0.41 m/s, with corresponding MAPE values ranging from 0.76% to 4.28%. For 90% of all LIDAR level estimations, absolute errors remain below 0.7 m/s, and absolute percentage errors remain under 10%. Higher errors were detected when abrupt shifts in wind speed or direction were observed (ramp events). LIDAR estimation errors decrease with height due to the diminishing influence of terrain curvature and friction. At the Enercon WT’s hub height, the MAE of wind velocity estimation is 0.85 m/s, aligning with the MAE for NEG, whereas Vestas exhibits a MAE of 0.67 m/s. MAPE values for Enercon, NEG, and Vestas are 7.78%, 8.78%, and 8.41%, respectively. Results are presented in Figure 11, Figure 12 and Figure 13.

4.3. E-IAS

Table 3. E-IAS-WT corrected velocity magnitude statistics.

Velocity Magnitude [m/s]
	Enercon	NEG	Vestas
MAE	0.81	0.76	0.56
MSE	1.08	0.93	0.58
RMSE	1.04	0.96	0.76
MAPE [%]	7.62	7.98	7.08

reveals that the MAE of wind velocity estimation is 0.81 m/s for Enercon WT, 0.76 m/s for NEG WT, and 0.56 m/s for Vestas WT. Corresponding MAPE values are 7.62%, 7.98%, and 7.08%, respectively. In Figure 14, Figure 15 and Figure 16, linear regression scatter plots suggest that extrapolating 1D measurements to a 2D inlet wind field reduced the error dependency on wind speed (w/s) for both Enercon and NEG simulations. This finding is promising, considering the E-IAS objective to improve estimations for locations not directly downstream of the LIDAR. Figure 14a–Figure 16a depict the X-axis velocity component (green), demonstrating that a wind parcel passing the LIDAR location typically travels southwest, often bypassing the Enercon WT.

4.4. Discussion

EAS LIDAR estimations are sufficiently accurate despite the simplifications of the study. Lack of measurements close to the ground level, where turbulence increase is vital, leading to higher errors. The extrapolation method used to estimate wind field below 40 m greatly affects the accuracy from ground level up to 54 m. Errors decrease as we estimate wind speed at higher heights, but, at 200 m, which is the highest measurement point, errors increase even if the flow becomes less turbulent. This increase in estimation errors highlights the importance of sufficient measurement points against the extended use of extrapolation. Because measurements below 40 m are much more important than measurements above 200 m for our case, we estimate that their use could significantly increase the accuracy of EAS.

As shown in

Table 4. Velocity magnitude estimation statistics comparison table.

		MAE	RMSE	MAPE [%]
	Enercon	0.88	1.11	8.26
IAS	NEG	0.76	0.99	7.91
	Vestas	0.56	0.6	7.09
	LIDAR 40 m	0.41	0.55	4.28
	LIDAR 54 m	0.37	0.48	3.70
EAS	Enercon	0.85	1.10	7.78
	NEG	0.85	1.06	8.78
	Vestas	0.67	0.88	8.41
	Enercon	0.81	1.01	7.62
E-IAS	NEG	0.76	0.96	7.98
	Vestas	0.56	0.76	7.08

, EAS results for ENERCON WT are better than IAS because the inlet is positioned further away from the wind turbines. The offset distance in the EAS combined with the unstructured mesh absorbs the shock of the simulated flow before it reaches the inlet of the IAS. This not only results in more minor errors but also removes bias dependency on wind speed.

E-IAS initial conditions depend on EAS. This translates to error dependency from EAS simplifications and slightly more than two times the computational time of IAS. E-IAS outperforms EAS and IAS due to the absence of measurement offsets and the elimination of flow shock in the simulated flow’s initial conditions at the inlet plane. This is achieved by implementing the curvature effect on the flow through physics-based extrapolation via EAS, rather than using a ramp function as in IAS and EAS. The E-IAS results strongly suggest that incorporating more measurement locations at the inlet will further enhance the model’s accuracy.

Although direct comparison between different model types is not typically recommended, we assess our model’s performance against state-of-the-art forecasting models in the literature to gain a nuanced understanding of its capabilities. This non-dogmatic evaluation aims to determine if this base model can be adapted and expanded for ultra-short-term wind power generation forecasting. As the core physics-based concept of the model will remain unchanged for this purpose, we do not anticipate any significant alterations in its behaviour.

In Figure 17, we can see the basic statistics of state-of-the-art models as presented by [25]. We see that the MAPE of our model is in the middle of the MAPE range of the other models at the first forecasting step. This is promising due to the simplifications adapted in our simulations. The MAPE of all the forecasting models at the second forecasting step is greatly increased, ranging between 8.09% and 14.89%, with the [25] proposed model almost doubling its MAPE to 8.09%. At the third forecasting step, their model sets the lower limit for MAPE at 11.73%, further raising it to more than 40%.

The behaviour of all models shows the difficulty in estimating such a volatile atmospheric variable in smaller time frames, as well as that our model has the capacity to create accurate estimations. Generating 3D high resolution datasets from 1D measurements with low errors is useful for any forecasting model like the ones in the comparison table. The creation of reliable extrapolated wind datasets is essential to any stage of a WF, from the instalment location of the WTs to the operation of the WF and the optimisation of a smart grid.

5. Conclusions

The presented model efficiently reconstructs a 3D wind field at 29,928 estimated locations using only 9 vertical measurement points per time step, achieving this in under 2 min. This represents a remarkable 3325-fold increase in data resolution, which is particularly valuable for enhancing ultra-short-term forecasting applications. Providing researchers with a physics-based model capable of rapidly overcoming terrain changes without disrupting energy production is crucial. We consider physics the foundation of research and development, and the accuracy of physics-based models is directly linked to the quality of input data. Consequently, high-resolution data inputs are expected to inherently yield more accurate results.

The average computational time needed for a single IAS simulation is 43 s; for EAS it is 61 s. The E-IAS study indicates that higher-resolution measurements could increase the accuracy of the model and reduce extreme errors, but the amount of time needed for the E-IAS is more than 100 s, which is much longer than the IAS due to extrapolation calculations and dataset loading. Results from the E-IAS model could be used to create extrapolation look-up tables that can be used in ultra-short-term forecasting models, increasing their accuracy without consuming much time. Statistical models can be retrained after terrain changes by reconstructing any given wind field without the need for new datasets.

Despite the simplifications and the lack of high-resolution measurements [14,26], the E-IAS model estimated wind speed at WTs’ hub height with relatively small errors, most of the time. Seventy percent of the absolute errors were below 1 m/s, whereas the absolute percentage errors were below 10%. Fifty percent of the absolute errors were below 0.71 m/s, and the absolute percentage errors were below 7.6%. High absolute errors occurred mostly at high wind speeds, whereas high error percentages were detected at low wind speeds. The highest errors were produced when extreme changes in wind speed and direction occurred.

To properly validate the model, more wind measurements are needed along the inlet, especially close to the ground level. Different vertical profile extrapolation methods of LIDAR measurements are expected to increase the quality of the results. Other meteorological variables like temperature, pressure, and humidity are also necessary for the calculation of atmospheric stability. This is essential for the creation of realistic simulations to understand the model’s behaviour under different atmospheric conditions. The accuracy of the estimations is already sufficient, considering that the model extrapolated a scarce set of vertical measurements (1D) to a 3D wind field over a large complex terrain area.