*Article* **Time-Averaged Wind Turbine Wake Flow Field Prediction Using Autoencoder Convolutional Neural Networks**

**Zexia Zhang 1, Christian Santoni 1, Thomas Herges 2, Fotis Sotiropoulos <sup>3</sup> and Ali Khosronejad 1,\***


Richmond, VA 23284, USA; sotiropoulosf@vcu.edu **\*** Correspondence: ali.khosronejad@stonybrook.edu

**Abstract:** A convolutional neural network (CNN) autoencoder model has been developed to generate 3D realizations of time-averaged velocity in the wake of the wind turbines at the Sandia National Laboratories Scaled Wind Farm Technology (SWiFT) facility. Large-eddy simulations (LES) of the SWiFT site are conducted using an actuator surface model to simulate the turbine structures to produce training and validation datasets of the CNN. The simulations are validated using the SpinnerLidar measurements of turbine wakes at the SWiFT site and the instantaneous and timeaveraged velocity fields from the training LES are used to train the CNN. The trained CNN is then applied to predict 3D realizations of time-averaged velocity in the wake of the SWiFT turbines under flow conditions different than those for which the CNN was trained. LES results for the validation cases are used to evaluate the performance of the CNN predictions. Comparing the validation LES results and CNN predictions, we show that the developed CNN autoencoder model holds great potential for predicting time-averaged flow fields and the power production of wind turbines while being several orders of magnitude computationally more efficient than LES.

**Keywords:** convolutional neural network; wind turbine; wake flow predictions; large-eddy simulation

### **1. Introduction**

In a wind farm, turbine wake interactions cause power losses and may increase fatigue loads on downwind wind turbines [1,2]. Therefore, the accurate prediction of turbine wakes is an important consideration in wind farm layout optimization, which can improve the efficiency of power production and reduce the overall levelized cost of energy. As a result, extensive efforts have been made on analytical and numerical models for the estimation of turbines wake [3–6].

Due to the simplicity and low computational cost, engineering models are widely used to predict wake flows and optimize wind farm power production, especially in industrial applications. The very first and extensively studied model was proposed by Jensen [7]. This model was derived from mass conservation, assuming a top-hat shape distribution of velocity deficit in the wake. However, the top-hat wake shape assumption is an oversimplification of the actual wake flow, which can be represented more accurately by a Gaussian distribution [8–10]. Furthermore, more complex real-life characteristics of wake flows have also been considered to improve the accuracy and flexibility of the Gaussian models, including the double-Gaussian type velocity profile of the near wake [11,12], threedimensional effects [13,14], more accurate models for turbulence intensities [15], wind turbine yaw offset [16], atmospheric stability, and Coriolis force [17]. Although these models are efficient, the accuracy varies significantly from case to case [18,19], especially in the near wake region [12,20]. In addition, wake overlapping effects are not accurately described, as shown by Archer et al. [21].

**Citation:** Zhang, Z.; Santoni, C.; Herges, T.; Sotiropoulos, F.; Khosronejad, A. Time-Averaged Wind Turbine Wake Flow Field Prediction Using Autoencoder Convolutional Neural Networks. *Energies* **2022**, *15*, 41. https:// doi.org/10.3390/en15010041

Academic Editors: Luis Hernández-Callejo, Sergio Nesmachnow and Sara Gallardo Saavedra

Received: 15 November 2021 Accepted: 17 December 2021 Published: 22 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Compared to engineering models, computational fluid dynamic (CFD) models can provide a better physics-based description of the dynamics of the turbine wakes, such as wake meandering [22–24] and the effects of atmospheric stability [25]. Moreover, some CFD models can even take into account the effect of complex terrain topology [26–28] in addition to the turbine tower and the nacelle [29,30]. The prediction of the velocity deficit and turbulent kinetic energy using the CFD methods is more accurate than engineering models [31]. However, the CFD methods employing LES are computationally expensive, and their use in wind farm optimization is becoming prohibitively expensive.

The development of machine learning and artificial intelligence has encouraged researchers to explore data-driven models to predict the wake and power production of turbines in a wind farm. For example, Japar et al. [32] used five different machine-learning methods, i.e., linear regression, linear regression with feature engineering, nonlinear regression, artificial neural network (ANN), and support vector regression (SVR), to estimate wind turbine power production based on free stream wind speed, wind direction and the turbine position in the wind farm. Although the more elaborate models, i.e., ANN and SVR, have higher accuracy, they slightly deviate from the measured power production in the high wind speed case. Sun [33] developed an ANN to predict power production of wind farms that considers wake effects for varying wind direction, wind speed, and yaw angle. The trained model successfully predicted the power production and was used to optimize the yaw angle of each wind turbine. However, these power production models need a large amount of input parameter combinations for the training and may only be applied to a particular wind farm. When considering the effects of more parameters, for example, turbulence kinetic energy, the whole neural network has to be retrained.

Other data-driven machine learning models have focused on the prediction of the velocity deficit of the wake. Wilson et al. [34] used Random Forest and Multilayer Perceptron (MLP) models to interpolate and predict wind velocity in the turbine wake, but it cannot be applied in different wind fields. Ti et al. [35,36] developed an ANN to predict the velocity deficit and turbulent kinetic energy field in the turbine wake from the incoming wind velocity and turbulence intensity. However, the method requires a large amount of CFD simulation results for training. Yang [37] developed a neural network model to predict the instantaneous position of the meandering turbine wake, using the upwind velocity, turbine torque, and turbine thrust as input features. Zhang and Zhao [38] proposed a neural network combining different dimensionality reduction techniques to predict the velocity field of distributed fluid systems and applied it successfully to predict the flow field of both a single turbine and turbine arrays. King et al. [39] proposed a Gaussian Process (GP) model to correct wind turbine flow field predictions from low-fidelity models, e.g., RANS model, to high-fidelity models, e.g., LES model. Ali et al. [40] used a Long-Short Term Memory (LSTM) model to successfully predict wind velocity fluctuation at specific locations in the turbine wake for a long time period. Renganathan et al. [41] combined an MLP and GP with a Convolutional Neural Network (CNN) decoder to map the wind turbine operation parameters, such as inflow wind speed, turbulent intensity, turbine power generation, atmospheric-boundary layer (ABL) Richardson number, rotor angular speed, and pitch angle, to the wake flow field. In addition to wake reconstruction, data-driven methods have also been used to identify and characterize turbine wakes. Aird et al. [42] developed a mask Region based Convolutional Neural Network (R-CNN) model that identifies turbine wakes in Lidar scan images with high accuracy, even with some missing data points, and is also able to character wake shapes in its forming and dissipating.

Despite these contributions, the accuracy of the existing algorithms for velocity field predictions varies with flow conditions and wind farm layouts, limiting their application for wind farm optimization. In this study, we develop a CNN autoencoder model for generating 3D realizations of time-averaged turbulent wake flow of wind turbines at the Sandia National Laboratories Scaled Wind Farm Technology (SWiFT) site in Lubbock, Texas. The site includes three Vestas V27 wind turbines to investigate the performance of the downwind turbine versus the upwind turbine wake. SpinnerLidar measurements of the

upwind turbine wake at the SWiFT facility have been used to validate the LES results of our in-house virtual flow simulator code, Virtual Flow Simulator (VFS-Wind model). After these validation comparisons, a series of LESs of the SWiFT site were conducted to produce instantaneous and time-averaged flow field results to train and test the CNN. Subsequently, the so-trained CNN was employed to predict the time-averaged flow field of new wind conditions. A data augmentation technique is employed to handle the location sensitivity problems of the CNN. The CNN predictions for the validation test cases were compared against the simulation results of the separately done LES validation case not used in the CNN training. In addition, the predicted time-averaged flow field of wind turbines was used to predict the time-averaged power production of the wind turbines.

This paper is organized as follows. In Section 2, the governing equations of the numerical model and the computational details of the LES of the SWiFT site are presented. In Section 3, the CNN autoencoder algorithm is described, followed by the results and discussion in Section 4. Final remarks can be found in Section 5.

### **2. Numerical Methods**

### *2.1. Governing Equations*

Simulations of the SWiFT site are performed by the LES module of the in-house incompressible Navier-Stokes solver code–VFS-Wind model. In the VFS-Wind code, the incompressible turbulent flow is described by the filtered Navier-Stokes and continuity equation written in curvilinear coordinates given as [43]:

$$J\frac{\partial \mathcal{U}^j}{\partial \xi^j} = 0,\tag{1}$$

$$\frac{1}{J}\frac{\partial \mathcal{U}^i}{\partial t} = \frac{\xi^i\_l}{J}\left(-\frac{\partial}{\partial \xi^j}\left(\mathcal{U}^j u\_l\right) + \frac{1}{\rho}\frac{\partial}{\partial \xi^j}\left(\mu \frac{\mathcal{g}^{jk}}{J}\frac{\partial u\_l}{\partial \xi^k}\right) - \frac{1}{\rho}\frac{\partial}{\partial \xi^j}\left(\frac{\mathcal{\xi}^i\_l p}{J}\right) - \frac{1}{\rho}\frac{\partial \tau\_{lj}}{\partial \xi^j\_l} + f\_l\right), \tag{2}$$

where *J* = *∂* \* *ξ*1, *ξ*2, *ξ*<sup>3</sup> /*∂*(*x*1, *x*2, *x*3) is the determinant of the Jacobian of the geometric transformation of the Cartesian coordinates {*xi*} and the generalized curvilinear coordinates ' *ξi* ( *, ξ<sup>i</sup> <sup>l</sup>* = *∂ξ<sup>i</sup>* /*∂xl* are the transformation metrics, *gjk* = *ξ<sup>i</sup> l ξk <sup>l</sup>* are the components of the contravariant metric tensor, *ui* is the *<sup>i</sup>*th Cartesian velocity component, *U<sup>i</sup>* = \* *ξi <sup>m</sup>*/*J um* is the contravariant volume flux, *p* is the pressure, *ρ* is the fluid density, *μ* is the dynamic viscosity of the fluid, *τij* is the sub-grid stress tensor for LES, which is modeled using the dynamic Smagorinsky sub-grid scale (SGS) model [44], and *fl* is the body force.

### *2.2. Numerics*

In the VFS-wind code, the governing equations are discretized on a hybrid staggered/nonstaggered grid in space. The convective terms are discretized using second-order accurate central differencing. For the divergence, pressure gradient, and viscous-like terms, the discretization used is the second-order accurate, three-point central differencing method [45]. The time derivatives are discretized using a second-order backward differencing scheme [46]. The discrete flow equations are time-integrated using an efficient, second-order accurate fractional step methodology in conjunction with a Jacobian-free, Newton-Krylov solver for the momentum equations and a GMRES solver enhanced with the multigrid method as a preconditioner for the Poisson equation.

### *2.3. Actuator Surface Model*

The wind turbine blades and nacelle are modeled using the actuator surface method developed by Yang and Sotiropoulos [47]. This method describes the flow field on the background Cartesian mesh and the wind turbine on the Lagrangian mesh following the actuator surfaces. Velocities on the actuator surfaces are interpolated from the background mesh using the smoothed discrete delta function proposed by Yang et al. [48]. Actuator surfaces of wind turbine blades are represented by chord lines along the radial direction. Lift and drag forces on blades are calculated similarly to the actuator line model, using the blade element momentum theory. For the nacelle actuator surfaces, the normal component of force is calculated by reconstructing the wall-normal velocity near the actuator surface to satisfy the non-penetration constraint, and the tangential forces are computed as a function of the local incoming velocity and a friction coefficient that parameterizes the effects of near-wall turbulence and the effects of surface geometry. Then the counterforces to the flow field (*fl* in Equation (2)) are calculated by distributing the forces on the blades to the background mesh using the above-mentioned interpolation method. Details of the actuator surface method can be found in Yang and Sotiropoulos [47].

### **3. Computational Details of SWiFT Site Simulation**

The SWiFT facility, located in Lubbock, Texas, is an experimental site supported by the U.S Department of Energy to investigate turbine wakes and turbine-turbine interactions. It is comprised of three research-scale wind turbines and two meteorological towers (METs). The layout of wind turbines is as shown in Figure 1b. The wind turbines, Vestas V27s, have rotor diameters of *D* = 27 m and hub heights of 32.1 m. A nacelle mounted Technical University of Denmark (DTU) SpinnerLidar is installed on turbine T1 to measure the turbine wake. The two METs are located upwind against the predominant wind direction to measure the atmospheric inflow. Details of the SWiFT facility can be found in Berg et al. [49].

**Figure 1.** (**a**) is the diagram of the computational domain. (**b**) is the relative location of three turbines.

A series of LESs were performed to generate 3D flow fields of the SWiFT site for a variety of inflow conditions, which give rise to different turbine wake interaction configurations. The computational domain of the real-scale modeled SWiFT site is shown in Figure 1a. It has a length of *L* = 80*D*, a width of *W* = 24.9*D*, and a height of *H* = 37*D*. Free-slip boundary condition is applied to the top and periodic condition along the spanwise direction. A logarithmic law of wall boundary condition is applied to the ground, which is given by

$$
\mu = \mu^\*/\kappa \, \, \ln(z/z\_0),\tag{3}
$$

where *u*<sup>∗</sup> is the friction velocity, *κ* is the von Karman constant, and *zo* = 0.0037*D* is the surface roughness height. The outlet is given by the Neumann boundary condition, and the inlet is fed with a fully-developed turbulent flow generated by a precursor simulation as described below. Three wind turbines are located over 8*D* downwind from the inlet. Each turbine has a hub height of *h* = 1.19*D* and a rotor diameter of *D* = 27 m. Turbine 2 is 2.99*D* west and 0.2*D* south from turbine 1. Turbine 3 is 5*D* north from turbine 1. The arrangement of the turbines is shown in Figure 1b. To generate different wake conditions, we conducted simulations for four wind directions (150◦, 0◦, 330◦, and 274◦, taking south as 0◦) as shown in Figure 2. Although the wind directions do not perfectly align with the cardinal directions, for the sake of brevity, they would be referred to as North-East, South, South-West, and West, respectively, through the rest of the paper. The flow domain is rotated in different

wind directions to ensure the *x*-axis is always along the wind directions. In addition, five

**Figure 2.** Diagram of wind turbine configurations in different wind directions; (**a**) north-east (150◦), (**b**) south (0◦), (**c**) south-west (330◦), and (**d**) west (274◦). Solid lines marked using T1, T2 and T3 represent the location of three turbines. Dashed lines represent the wake of each turbine.

The computational domain is discretized with a grid resolution of Δ*x* = 0.177*D* and Δ*y* = 0.089*D,* along the windwise and spanwise directions, respectively. A constant resolution of Δ*z* = 0.089*D* was given along the wall-normal direction up to a height of 7.46*D* and a grid stretching up to the top of the domain with a final resolution of Δ*z* = 1.48*D*. Therefore, a uniform grid is obtained in the bottom area that the turbines are located. The details of the computational domain are shown in Table 1.

**Table 1.** Geometrical and wind characteristics of the SWiFT site model simulations used for the training and validation of the CNN. *U*∞ is the free-flow velocity. *Re* is Reynolds number. *D* is the diameter of the turbine rotor (=27 m). *Nx*, *Ny*, and *Nz* are the number of computational grid nodes in windwise, spanwise, and vertical directions, respectively. Δ*x*, Δ*y*, and Δ*z* are the special resolution in windwise, spanwise, and vertical directions, respectively.


Precursor simulation of a neutral atmospheric boundary layer has been performed to prescribe the velocity at the inlet of the SWiFT site simulations. The precursor simulation has the same grid size and boundary conditions as the SWiFT site numerical domain. The initial transient of the simulation was discarded. After the mean kinetic energy of the computational domain reached steady state, velocities at a plane located in the center of the channel were saved at a <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> 2.0 <sup>×</sup> <sup>10</sup>−4*H*/*U*. The time-averaged velocity profile and turbulent intensity of the precursor simulation are shown in Figure 3.

⎯ **Figure 3.** Ten-minute time-averaged (**a**) velocity profile and (**b**) turbulence intensity; (◦) met-tower measurements and (**–**) precursor simulation.

### **4. Validation of the Computational Model**

A numerical simulation of a single wind turbine has been performed to validate the actuator surface model of the SWiFT wind turbine (Vestas V27) [50]. Numerical results are compared against MET data, the turbine supervisory control and acquisition system (SCADA), and SpinnerLidar measurements performed at the SWiFT facility by Sandia National Laboratories and the National Renewable Energy Laboratory [51]. Details of the field experiment and data acquisition can be found in Herges et al. [51].

An initial comparison of the inlet velocity profile obtained from the precursor simulation and the measurements was performed. Measurements of the 10-min time-averaged velocity and turbulence intensity obtained from a MET 2.5*D* upwind from the wind turbine are compared against the precursor simulation in Figure 3. Although the velocity profile shows good agreement, the turbulence intensity of the precursor simulation shows to be slightly larger than that obtained from the met-tower measurements. Furthermore, the difference in the turbulence intensity seems to be amplified with height, which suggests a slight decrease in turbulent kinetic energy due to atmospheric stratification may have occurred. However, the difference in the turbulence intensity is negligible, around 1% at hub height and 2.7% at 2.15*D*.

The yaw misalignment of the wind turbine was recorded with the SCADA system. Due to the constantly changing wind direction and the turbine yaw controller not perfectly tracking the wind during the measuring campaign, a time-dependent yaw offset is prescribed to the wind turbine in the simulations to match the measured flow misalignment, shown in Figure 4. In addition, the high-frequency fluctuations of the measured yaw offset, due to the small-scale turbulent coherent structures, are filtered by applying a locally weighted scatterplot smoothing with a window size of 100s to the offset signal. Therefore, a smoother yawing is prescribed to the turbine in the numerical simulation.

**Figure 4.** Yaw offset in time between the wind and the turbine; ( ⎯ **–**) SCADA measurement and ( ⎯ **–**) prescribed to the turbine on the numerical simulation.

⎯ ⎯

The numerical results are compared against the 10-min line-of-sight velocity measurements in the wake of the T1 wind turbine obtained from the nacelle-mounted DTU SpinnerLidar device [51]. The SpinnerLidar performed 984 scans at a constant focus distance from the device every two seconds. In addition, measurements were performed at 1*D* to 5*D* behind the wind turbine by varying the focus distance. The SpinnerLidar cycled the focus distance from 1*D* to 5*D* every 30 s.

To mimic the SpinnerLidar measurements, the velocity field of the numerical results are decomposed into a line-of-sight velocity (*VLOS*) with vertex, or origin, at the location of the nacelle. For higher fidelity in the comparison to the measurements, the simulated *VLOS* field is sampled at the approximate spatial coordinates and time as the SpinnerLidar measurements. The scattered data from the numerical simulation and the SpinnerLidar are interpolated into a spherical surface mesh to compute the 10-min time-averaged *VLOS*. A comparison of the horizontal line-of-sight velocity profiles at hub height is shown in Figure 5. Computed velocities in the wake of the turbine are slightly overestimated compared to the SpinnerLidar measurements, specifically close to the center of the rotor. Although the nacelle is modeled and avoids the formation of an unphysical jet in the center of the rotor (commonly observed in standard actuator line models), the momentum deficit seems slightly under-estimated. However, the numerical results from VFS-Wind show good agreement with the SpinnerLidar measurements.

**Figure 5.** Time-averaged line-of-sight velocity: (◦) SpinnerLidar and (–) numerical LiDAR from the LES; (**a**) 1*D*, (**b**) 2*D*, (**c**) 3*D*, and (**d**) 4*D* behind the wind turbine.

### **5. CNN Autoencoder Model**

We employed a CNN autoencoder model to predict the time-averaged flow field by extracting the key flow field features from instantaneous LES results. The CNN algorithm was originally developed to handle image recognition and image classification tasks [52,53]. Therefore, the CNN has some inherent advantages in handling high-dimensional data: it generally consists of convolutional layers and down-sampling layers that can reduce dimensions of the input image and extract abstract features of the image; the weight sharing concept (it will be explained in the next paragraph) used in CNN allows it to handle high-dimensional data using less learnable parameters and avoid location sensitivity problem [52]. As a variant architecture, the CNN autoencoder (Figure 6) consists of an encoder, which extracts features from input data, and a decoder, which is an inverse of the encoder. As a result, such CNN would enable the image reconstruction from the extracted features. Because of its ability to reconstruct field data, the CNN autoencoder has become a popular tool in the field of fluid dynamics, as well [54–57]. For these reasons, we employed the CNN autoencoder model in this study.

**Figure 6.** Schematic of the encoder-decoder CNN. Feature maps are depicted as solid boxes. Convolutional layers are depicted as gray dash lines. *L* × *W* × *H* × *channels* represents the dimensions of each feature map. L, W, and H represent the resolution of the input image in windwise, spanwise and vertical directions, respectively. The layer type, filter size, and stride size of each layer are shown below it. Strides represent the movement step-size of the convolutional filter.

The architecture of the CNN autoencoder used in this work is illustrated in Figure 6. The encoder part includes three 3-dimensional convolution layers for extracting features and down-sample the input data. Each convolution layer includes multiple channels corresponding to different features to be learned. The convolutional layer embedded with a nonlinear activation function and bias operates as follows [55]:

$$q\_i^l = \sigma\left(k\_i^l \odot \boldsymbol{\sigma}^{l-1} + b\_i^l\right),\tag{4}$$

where *q<sup>l</sup> <sup>i</sup>* is the output of *i*th channel in the *l*th convolutional layer, *σ* is the rectified linear unit (ReLU) nonlinear activation function [58], *k<sup>l</sup> <sup>i</sup>* is the *i*th trainable convolutional filter, <sup>⊗</sup> is the convolution operator, *ql*−<sup>1</sup> is the input of the *<sup>l</sup>*th convolutional layer and *<sup>b</sup><sup>l</sup> <sup>i</sup>* is the *i*th bias.

In Figure 7, we demonstrate the concept of the convolution operation. In this figure, *x* is a 4 × 4 input image with paddings of zero, *k* is a 3 × 3 convolution kernel (or filter), while *y* is the output of the convolution operation. The convolution window is traversed through the padded input image in both horizontal and vertical directions, where the convolution operation with the filter is performed. The step size of each move is the stride—i.e., the step size from red square to orange square. In the convolution operation, the filter can extract features from the input image, and the learnable weights stay unchanged as the convolution window moves. As a result, the weights of the filter are shared through the whole input image, and because of the weight sharing only one filter with nine weights is required to extract a feature through the entire padded image—instead of 16 different filters corresponding to each output element. In this approach, the filter will be independent of its location leading to fewer learnable weights and thus a more efficient training process. The convolution operation consists of an inner product between the convolution window (i.e., the red dashed box in the input of Figure 7) and the filter to generate the corresponding output element (e.g., the red dashed box in the output of Figure 7). Since the convolution operation reduces the image size, paddings of zeros are used to control the output size. For instance, the output size of the original 4 × 4 input image in Figure 7 (solid boxes in *x*) is 2 × 2 (solid boxes in *y*). However, with paddings, the output has the same size as the original input. In practice, the convolution operation can work on both 2- and 3-dimensional inputs.

**Figure 7.** Schematic of the convolution operation. *x* is the input, *k* is the convolution kernel, *y* is the output. Zeros around the input are padding. The orange and red squares represent corresponding input and output cells of the first and second convolution operation, respectively.

Since the convolution operation only involves a linear transformation, an activation function is needed to provide the nonlinear transformation into the CNN. Compared to the sigmoid or hyperbolic tangent activation functions, a rectified linear unit (ReLU) in hidden layers can increase the computational efficiency of the machine learning algorithm [58]. The ReLU function is given by [58]

$$\sigma(\theta) = \max(0, \theta), \tag{5}$$

where *θ* is the result of convolution operation plus the linear bias.

The decoder contains two transpose convolution layers, which are the inverses of the convolution layers, and a convolution layer to up-sample the data and construct the flow field. The discrepancy between output and the target values during the training iterations is calculated using the mean square error (MSE) loss function. Then, the weights in convolutional filters and the biases are updated by the backpropagation algorithm to minimize the loss function [59].

To determine the parameters of the CNN architecture, i.e., the number of layers and channels, and the kernel and stride sizes, a series of parameter combinations are tested to ensure the highest accuracy of results and least number of learnable parameters. The padding sizes are elaborately determined to guarantee the correct output size.

### **6. Results and Discussion**

We carried out LES of the SWiFT site under four different wind directions and five different wind velocities to train and validate the predictions of the CNN. First, we carried out the LES for all cases and discarded the initial two flow-through times—i.e., the duration of time it takes for an air particle to travel through the wind farm. Subsequently, the numerical simulations were continued until the first and second-order turbulence statistics were fully converged. The convergence of the time-averaged flow field is determined using a time-history-analysis approach reported in Khosronejad et al. [60]. During the training process, the fully converged instantaneous flow fields were fed into the CNN at the input layer, while the time-averaged flow field was designated as the target of CNN at the output layer. We note that the samples used in the training and, consequently, validation processes are taken from smaller domains around each turbine. These subdomains are 6*D* long, 3*D* wide, and 1*D* high. The turbine is located 2*D* downwind from the inlet and centered along the spanwise direction (Figure 8).

**Figure 8.** Schematic of the local area around turbine of CNNs sub-domain. The dashed line shows the plane of the rotor.

A schematic of the training procedure is shown in Figure 9. The input of each sample is composed of five neighboring snapshots with 1000 time-step intervals to convey the information of flow field fluctuations—induced by the wake meandering and large turbulent structures from the incoming flow. To examine the effect of the number of input snapshots on the accuracy of the trained CNN's predictions (for the time-averaged results), we conducted a systematic analysis in which we used a different number of input snapshots and calculated the computational errors. The computational error refers to the difference between the CNN's and LES results of the training and validation cases. Our findings for the relation between the number of input snapshots and the corresponding computational errors are shown in Figure 10. As seen in this figure, after five snapshots, the computational errors seem to plateau, suggesting that five snapshots are enough to reconstruct the time-averaged flow field. The target value of the CNN at the output layer is a statistically converged time-averaged flow field. Therefore, each training sample consists of five instantaneous snapshots and the time-averaged flow field.

**Figure 9.** Schematic of the training procedure to develop the CNN. The instantaneous velocity fields are fed into the CNN as the input signal while the time-averaged windwise velocity field is enforced as the target of the output signal.

**Figure 10.** Mean square error as a function of the number of input snapshots; (**a**,**b**) present the computational errors of the CNN for the training and test dataset, respectively.

The case we considered for the training of the CNN includes the SWiFT site with the north-east wind direction and wind speed of *U*<sup>∞</sup> = 7 m/s. The learning rate of the CNN had an initial value of 0.001 with a decay rate of 0.7 in a step size of 500 training epochs. Overall, 99 samples were used in the training process which took 1500 epochs of iterations until the loss curve was fully converged. The prediction of the training case is compared with the LES result in Figure 11. The difference between CNN prediction and the LES result is presented in Figure 12. Velocity profiles from six spanwise cross-sections are compared in Figure 13. The high accuracy of the CNN prediction demonstrates that the CNN is well trained.

Then the CNN is validated using 19 cases. Similar to the training, five successive snapshots (i.e., five instantaneous flow field data taken from five successive time steps) are time-averaged to produce the limited time-averaged flow field and used as the input to the trained-CNN to generate the reconstructed time-averaged flow field. The CNN-predicted time-averaged flow fields are then compared against the time-averaged results of separately conducted LES. For brevity, we only present the comparison of the overlapped turbine wakes for the four cases, e.g., north-east with *U*<sup>∞</sup> = 15 m/s, south with *U*<sup>∞</sup> = 13 m/s, south-west with *U*<sup>∞</sup> = 11 m/s and west with *U*<sup>∞</sup> = 9 m/s. In Figure 14, we plot the windwise velocity contours from top views at the hub level and side views at the hub layer. In Figure 15, we plot the difference between the CNN predictions and LES results. Additionally, velocity profiles taken from six spanwise cross-sections are compared in Figure 16. In the north-east (Figures 14a, 15a and 16a) and south-west (Figures 14c, 15c and 16c) wind direction cases, the presented turbines are 6*D* downwind of the other turbines. As seen in these figures, the trained CNN has been able to accurately resemble the velocity deficit in both the upwind and downwind wakes. The largest discrepancy is observed for the cases with the south (Figures 14b, 15b and 16b) and west (Figures 14d, 15d and 16d) wind directions, with the latter showing the maximum discrepancy. The present turbines are 5*D* and 3*D* downwind of the other turbines in the south and west case, respectively. The difference in the LES versus the CNN computed velocity fields seems to be caused by the proximity of the wind turbines, i.e., the closer the two turbines, the greater the computational error of the CNN predictions. We note that the CNN results for the velocity fields upwind of the turbines (Figure 15b,d, and profiles I, II in Figure 16b,d) in both the south and west cases seem to be slightly over-estimated. However, the CNN could predict the wake of the centered turbine with great accuracy, suggesting that the relatively higher discrepancy upwind of the turbine is due to the location sensitivity of the trained CNN. Theoretically, because of the "weight sharing" feature of the CNN, a trained CNN should be able to reconstruct the turbine wake at any location within the domain. However, the turbines and their wakes have almost the same location in all the training samples and the upwind wake does not strongly affect the training subdomain. Therefore, the trained CNN deems sensitive to turbine location.

**Figure 11.** Contours of time-averaged velocity normalized with the free-flow velocity for wind turbines in the training case. Top view cross-sections are at hub-height and side-view cross-sections are at the rotor center. Contours are from the CNN, the CNN with data augmentation method (CNNDA), and the LES results.

**Figure 12.** Contours of velocity difference between CNN and LES results normalized with the freeflow velocity for wind turbines in the training case. Top view cross-sections are at hub-height and side-view cross-sections are at the rotor center.

**Figure 13.** Velocity profiles along the spanwise direction at I, II, III, IV, V, and VI in Figure 11. (–) LES, (Δ) CNN, and (◦) CNN with data augmentation.

**Figure 14.** Contours of time-averaged velocity normalized with the free-flow velocity for wind turbines. (**a**) wind turbine 2 in north-east case; (**b**) wind turbine 3 in south case; (**c**) wind turbine 3 in south-west case; (**d**) wind turbine 1 in west case. Top view cross-sections are at hub-height and side-view cross-sections are at the rotor center. Contours are from the CNN, the CNN with data augmentation method (CNNDA), and the LES results.

**Figure 15.** Contours of velocity difference between the CNN and LES results normalized with the free-flow velocity for wind turbines. (**a**) wind turbine 2 in the north-east case; (**b**) wind turbine 3 in the south case; (**c**) wind turbine 3 in the south-west case; (**d**) wind turbine 1 in the west case. Top view cross-sections are at hub-height, and side-view cross-sections are at the rotor center.

**Figure 16.** Velocity profiles along the spanwise direction at I, II, III, IV, V, and VI in Figure 14. (**a**) wind turbine 2 in the north-east case; (**b**) wind turbine 3 in the south case; (**c**) wind turbine 3 in the south-west case; (**d**) wind turbine 1 in the west case. (–) LES, (- -) Analytical model, (Δ) CNN, and (◦) CNN with data augmentation.

To overcome the location sensitivity problem, we used the data augmentation technique. In this approach, instead of having a fixed location, the subdomains around the turbines are randomly moved up to 5*D* upwind. We also increase the size of training samples by a factor of five. Using this approach, the training samples contain the turbines and their wakes which are randomly scattered in each subdomain. The training and validation results of the CNN trained using the data augmentation technique (CNNDA) are depicted in Figures 11–16. The velocity deficits in the upwind wakes of the south and west cases are predicted more accurately in comparison to the results without data augmentation (Figure 15b,d, and I in Figure 16b,d). However, the discrepancies of the velocities around the turbines are still large in Figure 15b,d, because in these cases, the interaction with upwind turbine's wake is stronger than the training case.

The analytical model developed by Bastankhah and Porté-Agel [9] is used to compare against the CNN model. In this analytical model, the velocity deficit in the wake of a turbine is described as [19]:

$$\frac{\Delta u}{\mu\_0} = \left(1 - \sqrt{1 - \frac{C\_T}{8\left(\frac{\sigma\_y \sigma\_z}{D^2}\right)}}\right) \exp\left(-0.5\left[\left(\frac{y}{\sigma\_y}\right)^2 + \left(\frac{z}{\sigma\_z}\right)^2\right]\right),\tag{6}$$

where the Δ*u* is the velocity deficit in the wake, *u*<sup>0</sup> is the mean wind velocity perceived by the wind turbine, *CT* is thrust coefficient, *σ<sup>y</sup>* and *σ<sup>z</sup>* are the wake widths in spanwise and vertical directions, respectively, which are given by [19]:

$$\frac{x\_y}{D} = k\_Y \frac{x}{D} + 0.2 \sqrt{\frac{1 + \sqrt{1 - C\_T}}{2\sqrt{1 - C\_T}}},\tag{7}$$

$$\frac{\sigma\_z}{D} = k\_z \frac{\mathbf{x}}{D} + 0.2 \sqrt{\frac{1 + \sqrt{1 - \mathbf{C}\_T}}{2\sqrt{1 - \mathbf{C}\_T}}},\tag{8}$$

where *ky* and *kz* are the wake growth rates in spanwise and vertical directions, respectively. In the training and validation cases, *CT* = 0.8. The wake growth rates *ky* = 0.075 and *kz* = 0.075 are obtained by fitting the training case. For the turbines in the wake, a wake superposition model is considered [61]:

$$
\mu = \mathcal{U}\_{\text{hub}} - \sqrt{\sum\_{i}^{n} \Delta u\_{i}^{2}},
\tag{9}
$$

where *Uhub* is the free stream velocity at the hub height, *u* is the velocity in the wake of current turbine, Δ*ui* is the wake velocity deficit of the *i*th turbine in stand-alone condition, *n* is the number of superposition wakes. The velocity profiles of the analytical model are presented in Figure 16. The analytical model has a good performance in far wakes. However, the velocity deficits are overpredicted in the near wake areas (III in Figure 16) and underpredicted in the near upwind areas (II in Figure 16). In comparison, the presented CNNDA model seems to have a better performance than the analytical model in predicting the velocity field.

The discrepancy between the CNNDA predictions and the LES time-averaged results were quantified using the coefficient of determination (*R*2), mean absolute error (MAE), root mean square error (RMSE), and the mean absolute relative error (MARE). These statistical error indices are defined as follows [62]:

$$R^2 = 1 - \frac{\sum\_{i=1}^{N} \left(\psi\_{i(CNN)} - \psi\_{i(LES)}\right)^2}{\sum\_{i=1}^{N} \left(\psi\_{i(CNN)} - \overline{\psi}\_{i(CNN)}\right)^2} \tag{10}$$

$$\text{MAE} = \frac{\sum\_{i=1}^{N} \left| \psi\_{i(CNN)} - \psi\_{i(LES)} \right|}{N},\tag{11}$$

$$\text{RMSE} = \left(\frac{\sum\_{i=1}^{N} \left(\psi\_{i(CNN)} - \psi\_{i(LES)}\right)^2}{N}\right)^{0.5},\tag{12}$$

$$\text{MARE} = \frac{1}{N} \sum\_{i=1}^{N} \frac{\left| \Psi\_{i(CNN)} - \Psi\_{i(LES)} \right|}{\Psi\_{i(LES)}},\tag{13}$$

where *ψi*(*CNN*) is the predicted value (of time-averaged windwise velocity component) using the CNNDA, *ψi*(*LES*) is the value obtained from the LES model, *ψi*(*CNN*) is the mean predicted value using the CNNDA, and *N* is the total number of samples, i.e., the total number of computational nodes in the subdomain surrounding the turbine.

The discrepancies between the CNNDA predictions and the LES time-averaged results for the four cases are presented in Table 2. The CNNDA predictions maintain high accuracy for the turbines located in the free flow for all wind conditions (i.e., various magnitudes and directions). The largest discrepancies are observed for the turbines located in the wake and, specifically, the west wind direction case, where the distance between the turbines is the smallest. Nevertheless, the *R*<sup>2</sup> for all cases are over 0.95 and the RMSE less than 3%, which is quite remarkable.

**Table 2.** Statistical error indices of the CNNDA relative to the LES results for different wind direction cases.


The computational efficiency of the proposed CNN model has a great advantage over the LES. For each case, the LES required over 8 <sup>×</sup> 104 CPU hours to generate the timeaveraged flow fields. In comparison, although the proposed CNN model requires about 10 CPU hours for training, the well-trained model only requires 50 s to reconstruct the timeaveraged flow field. Considering the LES requires 9.6 <sup>×</sup> 103 to generate the instantaneous flow fields for inputs, the total cost of the proposed CNN model is still 88% less than the LES.

Now we turn our attention to the possibility of using the proposed CNN to predict aerodynamic power production of the individual turbines using the predicted velocity fields as follows:

$$P\_{\text{aero}} = \frac{1}{2} \rho A u^3,\tag{14}$$

where *Paero* is the aerodynamic power production, *ρ* = 1.225 kg m−<sup>3</sup> is the air density, A is the frontal rotor area, and u is the mean time-averaged windwise velocity component over the rotor area. In Figure 17, we compare the aerodynamic power productions of the four turbines using the LES results and CNNDA predictions. Wind velocities of the four cases range from 7 m/s to 15 m/s. Even though the velocity field discrepancies between LES and CNNDA are satisfactory, all wind turbines across all wind velocities overestimated the predicted power productions. As expected, the largest discrepancy is at turbine 1 in the west wind direction case, where the velocity deficit around the rotor area is underestimated (for instance, II in Figure 16d).

**Figure 17.** Power production of turbines located in the wake.

### **7. Conclusions**

In this study, we examined the capability of the CNN autoencoder to reconstruct the time-averaged flow field around the wind turbines at the SWiFT facility and predict turbine power output. LES of the SWiFT facility with four different wind conditions (i.e., north-east, south, south-west, and west directions with wind speeds of 7, 9, 11, 13, and 15 m/s) were performed to generate training and validation data for the CNN. A six-layer CNN autoencoder was developed and trained using both instantaneous and time-averaged LES results around three individual turbines using the north-east wind direction case. The input of every sample is constructed using five instantaneous velocity fields to reflect the temporal variations of turbulent structures. Subsequently, the trained CNN is validated and compared with time-averaged results of the additional large-eddy simulations. Based on the findings of this study, the following conclusions can be drawn:


and the cost of LES to produce instantaneous flow field (i.e., 9.6 <sup>×</sup> 103 CPU hours) for the inputs of the CNN, the total cost of the proposed CNN is 88% less than that of the LES. Therefore, the proposed CNN algorithms could enable reliable predictions of the wake flow field at a fraction of the cost required by the LES.

(3) The CNN predictions for the aerodynamic power productions were in good agreement with the LES results, except for the turbine located in the near wake of the upwind turbine owing to an underestimation of the velocity deficit within the wake. Overall, the comparisons between the LES results and CNN predictions of the SWiFT wind turbines demonstrate the potential of the developed CNN autoencoder for predicting time-averaged flow fields and the power production of wind turbines while being several orders of magnitude less computationally expensive than high-fidelity numerical simulations.

**Author Contributions:** Conceptualization, Z.Z., A.K.; methodology, Z.Z., A.K., T.H.; software, A.K., F.S.; validation, C.S., T.H.; formal analysis, Z.Z., A.K., C.S., F.S.; investigation, Z.Z., A.K., C.S.; resources, A.K., F.S.; data curation, C.S.; writing—original draft preparation, Z.Z., C.S.; writing review and editing, A.K., F.S., T.H.; visualization, Z.Z., C.S.; supervision, A.K.; project administration, A.K., F.S.; funding acquisition, A.K., F.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Offshore Wind Research and Development Consortium (NOWRDC) under agreement number 147503.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the size of the dataset.

**Acknowledgments:** The computational resources were provided by the Civil Engineering Department, Stony Brook Research Computing and Cyberinfrastructure, and the Institute for Advanced Computational Science at Stony Brook University. Additional computational resources were provided by the Minnesota Supercompting Institute at the University of Minnesota. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

