Numerical Investigations of Wake Expansion in the Offshore Wind Farm Using a Large Eddy Simulation

Abstract

Due to abundant wind resources and land saving, offshore wind farms have been vigorously developed worldwide. The wake of wind turbines is an important topic of offshore wind farms, in which the wake expansion is a key issue for the wake model and the layout optimization of a wind farm. The large eddy simulation (LES) is utilized to investigate various offshore wind farms under different operating conditions. The numerical results indicate that it is more accurate to calculate the wake growth rate using the streamwise turbulence intensity or the total turbulence intensity in the environment. By fitting the results of the LES, two formulae are proposed to calculate the wake growth rate of the upstream wind turbine. The wake expansion of the downstream wind turbine is analyzed, and the method of calculating the wake growth rate is introduced. The simulation indicates that the wake expansion of the further downstream wind turbines is significantly reduced. The smaller lateral distance of wind turbines in the offshore wind farm has the less wake expansion of the wind turbines. The wake expansion under different inflow wind speeds is also analyzed, while the wake expansion of wind turbines under more complex conditions needs to be further studied.

1. Introduction

Wind energy as a kind of renewable energy has been increasingly developed in many countries in the recent years. Due to the limited resources on land, offshore wind power is becoming a new direction of wind energy research for the future. Compared with onshore wind power, offshore wind power has the merits of higher wind speed, smaller wind vertical shear and longer annual utilization hours. Therefore, offshore wind power generation would become the main field of wind energy utilization.

Generally, people pay more attention to the power of a wind farm, and the wind turbine wakes have significant impacts on the power of the wind farm. Because the upstream wind turbine has a loss in the wake velocity, while the downstream wind turbine is affected by such a wake, the inflow conditions become worse, resulting in a reduction in the output power. Studies have shown that this power reduction can reach 25% of the total output power in the wind farm [1]. Therefore, it is important to accurately predict the wake flow and to optimize the layout of a wind farm so as to improve the output power of the wind farm. In order to accurately predict the wake, a lot of research has been conducted and many analytical wake models have been developed. Among them, the Jensen model proposed by Jensen [2] and the Gauss model proposed by Bastankhah and Porté-Agel [3] are commonly adopted, and they will be briefly introduced in the next section. In each of the models, the wake growth rate is a key parameter.

In subsequent studies, many researchers proposed new models based on the two analytical models. Tian et al. [4] improved the Jensen model and proposed the 2D_k Jensen model. This model modifies the wake growth rate in the Jensen model and uses the ratio of the wake turbulence intensity to the free stream turbulence intensity as the coefficient. Sun and Yang [5] proposed a new 3D wake model, which can calculate the vertical wake velocity distribution using the improved wake growth rate of Tian et al. [4]. Ishihara and Qian [6] introduced a novel wake model based on the Gaussian model, which considers the influence of the free stream turbulence intensity and the thrust coefficient, and the influence of thrust coefficient was added into the wake growth rate of the model.

In addition to the analytical wake model, the tool of computational fluid dynamics (CFD) is also extensively used in studying wind turbine wakes, typically the large eddy simulation (LES), which can accurately describe wind turbine wakes, compared with the Reynolds Average Navier-Stokes (RANS) model. Troldborg et al. [7] used the LES to research the wake of wind turbines and obtained the effect of inflow turbulence on the wake. Troldborg et al. [8] utilized the LES to investigate the wake interaction between two wind turbines under multiple inflow conditions. Wu and Porté-Agel [9] performed an LES to study the structure and characteristics of the wind turbine wake under four different aerodynamic roughness lengths. Zhou et al. [10] applied exact geometry to perform an LES and researched the influence of diverse inflow conditions on aerodynamic load and wake characteristics.

At present, the wind turbine wakes can be simulated accurately by using the LES technique. However, in the analytical wake model, it is difficult to accurately express the wake growth rate in theory, which can be obtained by fitting high-fidelity LES data. Niayifar and Porte-Agel [11] fitted the LES data to obtain a linear relation between the streamwise turbulence intensity and the wake growth rate, but the formula was not applicable to the condition of a streamwise turbulence intensity lower than 0.065. Ishihara and Qian [6] used LES data to fit the relationship between the wake growth rate, the thrust coefficient and the streamwise turbulence intensity. Using the previous LES data, Cheng et al. [12] obtained a linear relationship between the wake growth rate and the lateral turbulence intensity.

Obviously, the wake growth rate is a significant parameter to the analytical models, including the Jensen model, the Gaussian model and other modified models. This parameter is necessary to quantify the wake expansion and the velocity deficit. Therefore, it is important to accurately determine the wake growth rate for analyzing the wind turbine wake and optimizing the wind farm layout. The earlier studies have fitted the formulas of wake growth rate based on a single wind turbine with a fixed inflow wind velocity. However, the investigation of wake expansion and the validation of fitted formulas in the wind farm, especially in the offshore wind farm that has low turbulence intensity, are insufficient. The wake growth rate in offshore wind farms needs to be further studied. In the present work, the LES coupled with the actuator line model is utilized to simulate the offshore wind farms. In the simulations, different surface roughness is implemented to generate inflow conditions with different turbulence intensities. Based on the simulation results, the effect of different turbulence intensities on the wake growth rate is analyzed in the offshore wind farms. Then, two formulae of the wake growth rate are proposed by fitting the LES data, which are subsequently examined in the offshore wind farm and compared with the formulae introduced in earlier works. The wake expansions with different distances of turbine arrays and inflow wind speeds are also analyzed.

In the present paper, Section 2 presents an overview of the Jensen model and the Gauss model, Section 3 introduces the numerical techniques and the definition of wake growth rate, Section 4 presents the LES results and the study of the wake growth rate and the final conclusions are given in Section 5.

2. An Overview of the Jensen Model and the Gauss Model

The Jensen model [2] is an analytical wake model commonly used in many studies and commercial programs, in which it was supposed that the wake velocity deficit is in the shape of a top hat, and the velocity deficit is expressed by

$$\frac{\Delta U}{U_{\infty }}=(1-\sqrt{1-C_T})/{(1+\frac{2k_{wake}x}{D})}^2$$

(1)

$$\Delta U=U_{\infty }-U_{wake}$$

(2)

where, $U_{\infty }$ is the wind velocity of incoming flow; $U_{wake}$ is the wake velocity; $C_T$ is the thrust coefficient; $D$ is the diameter of the wind turbine; $k_{wake}$ is the wake growth rate; and $x$ is the downwind distance. Jensen assumed that the wake width increases linearly with downstream distance, and the wake growth rate $k_{wake}$ is constant. Jensen proposed $k_{wake}$ = 0.1 [2]. The value of $k_{wake}$ was corrected in subsequent studies; it is 0.075 in the onshore environment [13] and 0.04 or 0.05 in the offshore environment [13,14]. In addition, for the logarithmic wind profile, Frandsen [15] proposed a formula to calculate $k_{wake}$:

$$k_{wake}=\frac{0.5}{\ln (z_h/z_0)}$$

(3)

where, $z_h$ is the hub height of the wind turbine; and $z_0$ is the surface roughness length.

The Jensen model assumes that the velocity loss is in the shape of a top hat, which underestimates the velocity loss at the wake center and overestimates the velocity loss at the wake edge. Bastankhah and Porté-Agel [3] proposed a self-similar Gaussian distribution to better represent the velocity deficit profile in the wake. The normalized velocity deficit is expressed by

$$\frac{\Delta U}{U_{\infty }}=(1-\sqrt{1-\frac{C_T}{8{(\sigma /D)}^2}})\exp \left(-\frac{r^2}{2{\sigma }^2}\right)$$

(4)

$$\Delta U=U_{\infty }-U_{wake}$$

(5)

where, $r$ is the radial distance to the wake center; and $\sigma$ is the standard deviation of the Gaussian velocity deficit profiles. The wake width is also assumed to increase linearly, as expressed by

$$\frac{\sigma }{D}=k^{*}\frac{x}{D}+\epsilon$$

(6)

where, $k^{*}$ is the wake growth rate, which is different from $k_{wake}$; and $\epsilon$ is the value of $\sigma /D$ as $x$ approaches zero. In the Gauss model, the wake growth rate $k^{*}$ is also a key parameter to determine the region and velocity deficit of wakes, which hinges on the turbulence intensity of the incoming flow.

3. Numerical Model

3.1. Large Eddy Simulation (LES)

LES is a method based on spatial filtering operation, which separates large eddies from small eddies through a filter. In this way, large eddies can be simulated directly, and small eddies can be replaced by the subgrid-scale model. The filter width is defined as $\Delta ={(\Delta x\Delta y\Delta z)}^{1/3}$, where $\Delta x$, $\Delta y$ and $\Delta z$ are the grid lengths in the $x$, $y$ and $z$ directions, respectively. The filtered equations are derived from the incompressible Navier-Stokes equations. The filtered continuity and momentum equations are [16,17,18]

$$\frac{\partial {\overline{u}}_i}{\partial x_i}=0$$

(7)

$$\frac{\partial {\overline{u}}_i}{\partial t}+\frac{\partial ({\overline{u}}_i{\overline{u}}_j)}{\partial x_j}=-\frac{\partial \widehat{p}}{\partial x_i}-\frac{\partial {\tau }_{ij}^D}{\partial x_j}-\frac{1}{{\rho }_0}\frac{\partial p_0(x,y)}{\partial x_i}+F_{ext}$$

(8)

where, the upper bar represents the spatially resolved components; $t$ is time; $u$ is the wind speed; ${\rho }_0$ is the reference air density; ${\tau }_{ij}^D$ is the deviatoric part of the wind stress tenor; $F_{ext}$ is the external forces in the wind field; and $\widehat{p}$ is the modified pressure, defined by

$$\widehat{p}=\frac{\overline{p}(x,y,z,t)}{{\rho }_0}-\frac{p_0(x,y)}{{\rho }_0}+\frac{{\tau }_{kk}}{3}+gz$$

(9)

where, $\overline{p}$ is the mean pressure; and $p_0$ is the static pressure. According to the Boussinesq eddy viscosity assumption, ${\tau }_{ij}^D$ can be defined by

$${\tau }_{ij}^D=-2{(c_s\Delta )}^2\left|\overline{S}\right|{\overline{S}}_{ij}$$

(10)

where, $c_s$ is the Smagorinsky constant; ${\overline{S}}_{ij}=(\partial {\overline{u}}_i/\partial x_j+\partial {\overline{u}}_j/\partial x_i)/2$ is the strain rate tensor after filtration; and $\left|\overline{S}\right|=\sqrt{2{\overline{S}}_{ij}{\overline{S}}_{ij}}$ is the norm of the strain rate tensor. The external force $F_{ext}$ includes the force exerted by the blades of wind turbines, the Coriolis force and the buoyancy force. The external force $F_{ext}$ is expressed by

$$F_{ext}=\frac{1}{{\rho }_0}{F}_i+g\left(\frac{\overline{\theta }-{\theta }_0}{{\theta }_0}\right){\delta }_{i3}-{\epsilon }_{i3k}{f}^{}{\overline{u}}_k$$

(11)

where, $F_i$ is the force exerted by the wind turbines, which can be calculated by the actuation line model; $g$ is the acceleration of gravity; ${\epsilon }_{ijk}$ is the alternating unit tensor; $\overline{\theta }$ is the resolved potential temperature; ${\theta }_0$ is the reference temperature, with the value of 300 k; ${\delta }_{i3}$ is the Kronecker delta; and $f$ is the Coriolis parameter.

The potential temperature can be obtained by combining the following potential temperature equation with Equations (7) and (8):

$$\frac{\partial \overline{\theta }}{\partial t}+\frac{\partial ({\overline{u}}_j\overline{\theta })}{\partial x_j}=\frac{\partial q_j}{\partial x_j}$$

(12)

where the temperature flux $q_j$ is defined by [19]

$$q_j=-\frac{{\nu }_t}{{\mathrm{Pr}}_t}\frac{\partial \overline{\theta }}{\partial x_j}$$

(13)

where ${\mathrm{Pr}}_t$ is the subgrid turbulent Prandtl number.

3.2. Wind Turbine Modeling

In the numerical simulation of wind turbines, in order to reduce the computational cost, the wind turbine is generally not directly modeled, but represented by a simplified model. Two universally-used models include the actuator disk model (ADM) and the actuator line model (ALM). The ADM [20] simplifies the wind turbine load to the actuator-disk surface. The actuator-disk surface refers to the area swept by the blade. The lift and drag are calculated by the blade-element momentum (BEM) theory [20] and are distributed to the actuator disc in a certain way. However, there is a defect in the ADM, which cannot simulate helicoidal tip vortices. The ALM proposed by Sørensen and Shen [21] does not have this defect and is widely used in the LES of wind turbines [22,23,24,25]. Therefore, the ALM is used in the present work.

The aerodynamic force of each blade element is expressed as $f_i^a(x_j,y_j,z_j,t)$, which represents the value of the blade element $j$ at the position $(x_j,y_j,z_j)$ and the time $t$. After the correction of regularization kernel, the sum of the aerodynamic forces of the N blade elements is expressed by

$$F_{\mathrm{i}}={\displaystyle \sum _{j=1}^N\frac{f_i^a(x_j,y_j,z_j,t)}{{\epsilon }_i{}^3{\pi }^{3/2}}\exp \left[-{\left(\frac{r_n}{{\epsilon }_i}\right)}^2\right]}$$

(14)

where, $r_n$ is the distance between the center of the grid and the point representing the actuator line; and ${\epsilon }_i$ is the parameter to adjust the regularized load distribution. In order to obtain more accurate results, the tower and nacelle aerodynamics are simulated. The wind turbine selected in the present work is NREL 5 MW [26], the rotor diameter is 126 m and the rotor shaft height of the wind turbine is 90 m.

3.3. Numerical Setup

The solver used in the present work is the Simulator for Wind Farm Applications (SOWFA), which is developed by the National Renewable Energy Laboratory (NREL) [16,27], Colorado, USA. The version is SOWFA-2.4.x. The continuity and momentum equations are discretized employing an unstructured, collocated, finite-volume formulation with second-order central differencing. Time discretization is the second-order backward. The solution algorithm is the Pressure Implicit Splitting Operation (PISO) [28]. Further details can be found from Churchfield et al. [16,29]. All the simulations are conducted on the National Supercomputer Center in Tianjin (TianHe-1A), using 280 processors. In this study, the computational domain size is 3040 m × 3000 m × 1000 m, and the initial grid size is 304 m × 300 m × 100 m. The simulations are carried out with two steps. Firstly, the precursor simulation without wind turbines is simulated to generate the turbulent flow field. In this calculation stage, the inlet air flows forward along the X-axis. The velocity is set as the desired mean wind speed $U_{\infty }$ at the height of the hub center, and small, divergence-free perturbations are superimposed on the flow for generating the turbulent atmospheric winds. All lateral boundaries are periodic. Then, the wind turbines are added to the generated flow field for the simulation. At this stage, the upstream boundary condition uses the results of the precursor simulation as input, and the downstream boundary condition has zero velocity gradient and temperature gradient. The grid near the wind turbines is refined. In order to reach a quasi-equilibrium condition, the precursor simulation is conducted for 24,000 s with a time step of 0.5 s and the wind farm simulation for 4000 s with a time step of 0.2 s. The flow is considered in a quasi-equilibrium state when the oscillation of the sampled velocity is about some mean value at half the boundary layer height [16]. The data of the last 4000 s in the precursor simulation from the most upwind lateral boundary are used as the input boundary conditions of the wind farm simulation. In order to ensure the accuracy of results, the last 1200 s of data are extracted for inspection.

Table 1. Numerical settings.

Parameters	Value
Computational domain size (m)	$3040\times 3000\times 1000$
$\Delta x$ (m)	2.5
$\Delta y$ (m)	2.5
$\Delta z$ (m)	2.5
Precursor simulation time (s)	24,000
Precursor simulation time step (s)	0.5
Wind farm simulation time (s)	4000
Wind farm simulation time step (s)	0.2
Number of cores	140
$f$ (s−1)	9.6 × 10−5
$c_s$	0.13
Molecular Prandtl number	0.7
Solver tolerance of convergence	1 × 10−6
Solver relative tolerance of pressure	0.01
Maximum number of iterations of convergence	1000

shows the settings of the numerical simulation, where $\Delta x$, $\Delta y$ and $\Delta z$ represent the grid sizes in the $x$, $y$ and $z$ directions of the refinement area, respectively.

In this study, at first, a single wind turbine is simulated for the grid inspection and model verification. See the next section for details. In the subsequent study, multiple wind turbines are simulated. The layout of the wind farm with 6D lateral spacing is shown in Figure 1. The origin of coordinates is located on the coordinate axes in the figure. The grid refinement areas are two equal boxes with a length of 22.5D, a width of 3.5D and a height of 1.8D, as shown in Figure 1. In the figure, D represents the wind turbine rotor diameter, the short thick line represents the wind turbine, WT0, WT1, WT2, WT3, WT4 and WT5 are the names of wind turbines, the two dotted boxes are grid refinement areas and the larger dotted line in the middle is the axis of symmetry. The size of the grids in $x$, $y$ and $z$ directions is all 2.5 m. In order to eliminate the mutual influence of two rows of wind turbines, the working conditions of three wind turbines are simulated with only WT0, WT1 and WT2. In these cases, the distance between the upstream and downstream wind turbine is 7D, and the grid refinement area is similar to the above. In the case of the wind farm with 3D lateral spacing, the arrangement is similar to the case with 6D lateral spacing, but the lateral distance between two arrays is modified from 6D to 3D in the $y$ direction. The grid refinement area is similar to the above. By changing the surface aerodynamic roughness height $z_0$, the inflow conditions under different turbulence intensities are obtained. The present work focuses on the offshore environment, in which $z_0=0.001$ represents the typical offshore conditions and $z_0=0.016$ represents the typical rough sea conditions [16]. The settings of the cases are listed in

Table 2. Case settings.

Case	${\mathit{z}}_{\mathbf{0}}$ (m)	N	${\mathit{U}}_{\mathit{\infty }}$ (m/s)	LD (D)
1	0.001	3	8	-
2	0.005	3	8	-
3	0.01	3	8	-
4	0.016	3	8	-
5	0.001	6	8	6
6	0.003	6	8	6
7	0.007	6	8	6
8	0.016	6	8	6
9	0.001	6	8	3
10	0.016	6	8	3
11	0.001	3	9	-
12	0.001	3	10	-
13	0.001	3	11	-
14	0.001	3	12	-

, where $N$ represents the number of wind turbines; $U_{\infty }$ represents the incoming wind speed at the height of hub center; and LD represents the lateral distance between wind turbines.

3.4. Grid Sensitivity Verification and Model Validation

In this study, the grid sensitivity verification and the model validation are carried out at first. A single wind turbine is simulated, and the computational domain, initial grid and simulation methods are the same as described previously. The wind turbine is located between WT0 and WT3 in Figure 1, and the grid refinement area is similar to the above. The surface aerodynamic roughness height is set to 0.001 m.

The results of Churchfield et al. [16] are selected for validation. Their research has been widely accepted and is used for the validation of the NREL 5 MW wind turbine. The present inflow conditions are basically the same as their inflow conditions, but an improved ALM is adopted, which takes into account the tower and nacelle aerodynamics.

In order to verify the grid independence, three different grid sizes in the refinement area are set up, as presented in

Table 3. Numerical setup with different grid resolutions (GR).

GR Case	$\mathbf{\Delta }\mathit{x}$ (m)	$\mathbf{\Delta }\mathit{y}$ (m)	$\mathbf{\Delta }\mathit{z}$ (m)	${\mathit{N}}_{\mathit{y}}$	${\mathit{N}}_{\mathit{z}}$	Grid Number	Time (Core Hours)	Mean Relative Error (%)
1	5	5	5	25	25	$4.2\times {10}^5$	1176	6.91
2	2.5	2.5	2.5	50	50	$3.2\times {10}^6$	1792	6.82
3	2.5	1.25	1.25	100	100	$1.3\times {10}^7$	2744	6.75

, where $\Delta x$, $\Delta y$ and $\Delta z$ represent the grid sizes in the streamwise, the lateral and the vertical directions of the refinement area, respectively; $N_y$ and $N_z$ represent the grid numbers in the lateral and the vertical directions in the diameter of the rotor, respectively. The time-mean wake velocity of the three cases on the horizontal profiles of the hub height are calculated. Three distances of 2D, 4D and 6D downstream of the wind turbine are selected and compared with the results of Churchfield et al. [16], by calculating the mean relative errors of GR1, GR2 and GR3, as shown in Table 3. The results show that the error decreases with the grid refinement, but the results are insensitive to the grid resolution within this resolution range. In order to obtain more accurate results and spend less computational cost, GR2 is chosen for the subsequent study. The turbulence resolution $M$ is calculated as $M=k_{SGS}/(K+k_{SGS})$ [30], where $k_{SGS}$ represents the sub-grid energy, and $K$ represents the resolved part of the energy. The turbulence resolution of all grids in GR2 was less than 0.1, that is, the percentage of turbulence kinetic energy resolved is greater than 90%. Figure 2 shows the comparison between the normalized time-mean streamwise velocity deficit of GR cases and the results of Churchfield et al. [16]. The error is probably due to the improved ALM adopted in the present work. The improved ALM can simulate the aerodynamic forces of the tower and the nacelle, which were not considered in the study of Churchfield et al. [16]. These forces affect the distribution of wake deficits, especially in the near field. The maximum Courant number (CFL) is less than 1. The maximum wall Y+ is between $1\times {10}^5$ and $2\times {10}^5$, which is a very large value. However, in many LES studies of large-scale wind farms, due to the stress and temperature flux models employed at the planetary surface, the magnitude of the Y+ value is common, and the simulation results have proven to be good [9,16,17,18,31]. Anyway, the present LES results exhibit acceptable accuracy.

3.5. Wake Growth Rate

The physical meaning of wake growth rate is not exactly the same in the different analytical wake models, which can be roughly divided into two categories: $k_{wake}$ in the Jensen model and $k^{*}$ in the Gauss model. Although the two types of wake growth rates have different physical meanings, they can be easily converted. In all cases, the value of $k_{wake}$ is about twice of $k^{*}$. In this study, the wake growth rate is defined the same as $k_{wake}$ [32]. The wake radius is expressed by

$$r_w=k_{wake}x+r_0$$

(15)

where, $r_w$ is the wake radius; $k_{wake}$ is the wake growth rate; $x$ is the downwind distance; and $r_0$ is the initial radius of the wake. In the analysis of numerical simulation results, the criterion proposed by Ge et al. [32] is adopted to obtain the physical boundary of wakes. The time-mean velocity deficit at the wake boundary is defined as 5% of the maximum velocity deficit, and the wake radius $r_w$ equals approximately $2\sigma$.

4. Results and Discussion

4.1. Upstream Wind Turbine

In order to study the wake boundary of wind turbines, the time-mean streamwise velocity and turbulence intensity of each case are calculated. The streamwise turbulence intensity is defined as $I_u={\sigma }_u/U_{hub}$, the lateral turbulence intensity is defined as $I_v={\sigma }_v/U_{hub}$ and the vertical turbulence intensity is defined as $I_w={\sigma }_w/U_{hub}$, where ${\sigma }_u$, ${\sigma }_v$ and ${\sigma }_w$ represent the standard deviations of the streamwise, the lateral and the vertical velocities, respectively, and $U_{hub}$ represents the incoming wind speed at the hub height. The standard ${\sigma }_u$ is calculated as ${\sigma }_u=\sqrt{{\displaystyle \sum _{t=1}^T{(U_t-\overline{U})}^2/T}}$, where $U_t$ represents the streamwise velocity at time $t$; $T$ stands for the length of time; $\overline{U}$ represents the average streamwise velocity over time. The data are taken every second. The standard deviations in other directions are calculated similarly. The total turbulence intensity is defined as $I=\sqrt{({\sigma }_u^2+{\sigma }_v^2+{\sigma }_w^2)/3}/U_{hub}$ [16,31,33], which can reflect the turbulence levels in the three directions. In this study, the cases of three wind turbines (Cases 1–4) are used for fitting, and the cases of six wind turbines (Cases 5–8) are used for verification. The data from 2D to 6D downstream of the wind turbine are used to calculate the wake growth rate. Figure 3 shows the time-mean streamwise velocity contours in the x-y plane of rotor center height for some cases. Figure 4 shows the contours of streamwise turbulence intensity in the x-y plane of the rotor center height for some cases. The short black lines represent the wind turbines.

As illustrated in Figure 3 and Figure 4, the turbulence intensity has an obvious influence on the wake structure of wind turbines. The higher turbulence intensity in the environment accompanies the faster wake recovery, which means a shorter distance is required for wake recovery. In the previous studies [6,11], the streamwise turbulence intensity was used to describe the wake growth rate. Some studies [12,33,34] have also shown that the streamwise turbulence intensity is insufficient to describe the wake development. In this study, the effects of streamwise turbulence intensity, lateral turbulence intensity and total turbulence intensity on wake growth rate are studied. The relationship between the wake growth rate and the streamwise, lateral and total turbulence intensities is fitted respectively. The results show that the expressions with the streamwise turbulence intensity or the total turbulence intensity can better reflect the wake growth rate. Based on the analysis and previous studies, two formulae for calculating the wake growth rate are proposed.

The fitting formula of wake growth rate of the upstream wind turbine with respect to turbulence intensity is derived as

$$k_{wake}=1.233I-0.024$$

(16)

where $I$ represents the total turbulence intensity in the free stream area at the height of the hub center. The applicable scope of the turbulence intensity is $0.04<I<0.07$.

The fitting formula of wake growth rate of the upstream wind turbine with respect to streamwise turbulence intensity is derived as

$$k_{wake}=0.861I_u-0.016$$

(17)

where $I_u$ represents the streamwise turbulence intensity in the free stream area at the height of the hub center. The applicable scope of the turbulence intensity is $0.05<I_u<0.08$. Figure 5 shows the wake growth rate of WT0 in the LES results and the fitted lines.

In the offshore wind farms, wind turbine wakes may be affected by lateral wind turbines. To prove the applicability of the formula in the offshore wind farms, six wind turbine cases are used for validation. Figure 6 and Figure 7 show the contours of time-mean streamwise velocity and streamwise turbulence intensity, respectively, in the x-y plane passing through the hub center of six wind turbines.

This study also examines the applicability of the wake growth rate formula proposed in the previous study for offshore wind farms. Frandsen [15] proposed a formula to calculate $k_{wake}$ in 1992, as presented in Equation (3). For the Gaussian model, Niayifar and Porté-Agel [11] obtained the linear relationship between $k^{*}$ and $I_u$ by fitting LES data, but this formula is only applicable to the range $0.065<I_u<0.15$. Fuertes et al. [35] fitted a different formula through measured data, expressed by

$$k^{*}=0.35I_u$$

(18)

Ishihara and Qian [6] put forward another form of empirical formula in the improved Gaussian model, expressed by

$$k^{*}=0.11C_T^{1.07}{I}_u^{0.2}$$

(19)

Cheng et al. [12] obtained the linear relationship between $k^{*}$ and $I_v$ by fitting previous LES results, and the expression is

$$k^{*}=0.223I_v+0.022$$

(20)

Table 4. Wind farm simulation results.

Case	${\mathit{C}}_{\mathit{T}}$	${\mathit{z}}_{\mathbf{0}}$	${\mathit{I}}_{\mathit{u}}$	${\mathit{I}}_{\mathit{v}}$	$\mathit{I}$	${\mathit{k}}_{\mathit{w}\mathit{a}\mathit{k}\mathit{e}}$
5	0.86	0.001	0.0559	0.0409	0.0447	0.03084
6	0.88	0.003	0.0632	0.0482	0.0513	0.03866
7	0.89	0.007	0.0689	0.0526	0.0544	0.04280
8	0.88	0.016	0.0746	0.0584	0.0607	0.04942

shows the wind farm simulation results used for validation, where $C_T$ and $k_{wake}$ are the calculation results of WT0. The validation results are listed in

Table 5. Validation of fitting formulae.

	Case 5		Case 6		Case 7		Case 8
Equation Number	${\mathit{k}}_{\mathit{w}\mathit{a}\mathit{k}\mathit{e}}$	Relative Error (%)	${\mathit{k}}_{\mathit{w}\mathit{a}\mathit{k}\mathit{e}}$	Relative Error (%)	${\mathit{k}}_{\mathit{w}\mathit{a}\mathit{k}\mathit{e}}$	Relative Error (%)	${\mathit{k}}_{\mathit{w}\mathit{a}\mathit{k}\mathit{e}}$	Relative Error (%)
3	0.044	42.67	0.049	26.75	0.053	23.83	0.058	17.36
18	0.039	26.46	0.044	13.81	0.048	12.15	0.052	5.22
19	0.105	240.47	0.110	184.53	0.114	166.36	0.114	130.68
20	0.062	101.04	0.065	68.13	0.067	56.54	0.070	41.64
16	0.031	0.52	0.039	0.88	0.043	0.47	0.051	3.20
17	0.032	3.76	0.038	1.7	0.043	0.47	0.048	2.87

, where $k^{*}$ is converted to $k_{wake}$ and the relative error is calculated based on the simulation results. It is observed that the error of the previous results is large, indicating their poor applicability in the environment with low turbulence intensity. This is probably because the fitting formula in previous studies was obtained under the condition of high turbulence intensity. The present fitting formulae have quite small errors, proving the applicability in the offshore wind farm with low turbulence intensity. From the fitting and validation results, the data of a total of eight cases (Cases 1–8) can well reflect the linear relationship between the wake growth rate and the turbulence intensity.

4.2. Downstream Wind Turbine

As mentioned in the previous section, after analyzing the simulation results, the wake growth rates of upstream and downstream wind turbines have different laws. It is found that the wake growth rate of the downstream wind turbine decreases distinctly compared with that of the upstream wind turbine, as shown in Figure 3. The downstream wind turbine is located in the wake of the upstream wind turbine, and its own wake is the superposition of the two wakes. Although the turbulence intensity of the downstream wind turbine is higher, the inlet wind speed of the downstream wind turbine is significantly lower than that of the upstream wind turbine. As shown in

Table 6. Time-mean streamwise velocity and streamwise turbulence intensity.

Case	${\mathit{U}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}$	${\mathit{I}}_{\mathit{u}}$
1	5.06	0.0907
2	5.15	0.0957
3	5.24	0.0967
4	5.38	0.1003

, $U_{mean}$ and $I_u$ are the time-mean streamwise velocity and streamwise turbulence intensity, respectively, upstream of the WT1 rotor center. As a result, the wake recovery speed of the downstream wind turbine is slow, and the wake growth is low. Through the study of the simulation results, it is found that the wake growth rate of the upstream wind turbine can be modified to obtain the wake growth rate of the downstream wind turbine. The wake growth rate of the first downstream wind turbine is about 0.65 times the size of the upstream wind turbine, while the wake growth rate of the second downstream wind turbine is close to zero.

4.3. Different Lateral Spacing

In offshore wind farms, wind turbines often interact with each other. To further study this interaction, the cases with small lateral spacing are simulated. Figure 8 and Figure 9 show the contours of time-mean streamwise velocity and streamwise turbulence intensity, respectively, in the x-y plane passing through the hub center of six wind turbines. It is observed that the wake growth rate of the wind turbine is significantly lower than that of the case with large lateral spacing. This is mainly because, when the lateral distance is too small, the interaction of adjacent wind turbines will significantly reduce the wake recovery speed and lead to a larger wake area. Therefore, when predicting the wake growth rate of offshore wind farms with small lateral spacing, the results need to be modified. The results of this study show that the modified result is about $0.7k_{wake}$.

4.4. Different Inflow Wind Speeds

In order to study the influence of different inflow wind speeds, the cases of different inflow wind speeds are simulated. The wind speed range in working conditions is 8–12 m/s, which is the mean range in most offshore environments of China, Europe and the United States [36]. Figure 10 and Figure 11 show the contours of time-mean streamwise velocity and streamwise turbulence intensity, respectively, in the x-y plane of rotor center height. It is observed that at the same surface roughness height, different inlet wind speeds have a similar environmental turbulence intensity. However, the wake growth rates of upstream and downstream wind turbines are not consistent with the proposed formulae. This is probably because wind turbines have different physical characteristics under different inflow wind speeds, such as the rotational speed and the thrust coefficient. Previous studies have found that there is a linear relationship between the wake growth rate and the turbulence intensity under different inlet wind speeds [11,35]. Therefore, it is considered that the wake growth rate and the turbulence intensity meet the linear relationship under different wind speeds. According to the LES results of this study, Equations (16) and (17) can be used to calculate the wake growth rate when the inflow wind speed is 8–9 m/s. However, when the inflow wind speed is 10–11 m/s, the wake growth rate calculated by Equations (16) and (17) needs to be modified to $1.6k_{wake}$.

5. Conclusions

The LES technique is used to simulate a series of offshore wind farms. The law of wake expansion of wind turbines is analyzed in offshore wind farms with different wind speeds, environmental turbulence intensities and lateral spacing. The equations of the wake growth rate of upstream wind turbines are proposed by fitting the LES data. Major findings are summarized as

(1)

The higher turbulence intensity in the environment accompanies the greater wake expansion of upstream wind turbines. There is a linear relationship between the wake grow rate $k_{wake}$ and the environmental turbulence intensities $I$ and $I_u$.

(2)

The downstream wind turbines have smaller wake growth rates compared with the upstream turbines, even though they are in the stronger turbulence produced by the upstream turbines. The wake expansion of the further downstream wind turbines is significantly reduced.

(3)

Although the nearby wind turbines are not in the wake expansion region of other turbines, their wake growth rate may still be affected. The smaller lateral spacing in the wind farm decreases the wake expansion of turbine arrays.

(4)

The linear relationship between the wake grow rate and the environmental turbulence intensities is different at different inflow wind speeds. The proposed formulae need to be modified for different inflow wind speeds.

(5)

The formulae of the wake growth rate and the calculation methods can obtain more accurate wake expansion and velocity deficit and, therefore, improve the precision of analytical wake models. However, the wake expansion of wind turbines needs to be further investigated under more complex conditions, such as different stability of atmospheric boundary layers, different wind turbine spacing and wind turbines in yaw conditions.