A semi-parabolic wake model for large offshore wind farms based on the open source CFD solver

Wake effect represents one of the main sources of energy loss and uncertainty when designing offshore wind farms. Traditionally analytical models have been used to optimize and estimate power deficits. However these models have shown to underestimate wake effect and consequently overestimate output power [1, 2]. This means that analytical models can be very helpful at optimizing preliminary layouts but not as accurate as needed for an ultimate fine design. Different techniques can be found in the literature to study wind turbine wakes that include simplified kinematic models and more advanced field models, that solve flow equations with different turbulence closure schemes. See the review papers of Crespo et al. [3], Vermeer et al. [4], and Sanderse et al. [5]. Purely elliptic Computational Fluid Dynamics (CFD) models based on the actuator disk technique have been developed during the last years [6–8]. They consider wind turbine rotor as a disk where a distribution of axial forces act over the incoming air. It is a fair approach but it can still be computationally expensive for big wind farms in an operative mode. With this technique still active, an alternative approach inspired on the parabolic wake models [9, 10] is proposed. Wind turbine rotors continue to be represented as actuator disks but now the domain is split into subdomains containing one or more wind turbines. The output of each subdomain is mapped onto the input boundary of the next one until the end of the domain is reached, getting a considerable decrease on computational time, by a factor of order 10. As the model is based on the open source CFD solver OpenFOAM, it can be parallelized to speed-up convergence. The near wake is calculated so no initial wind speed deficit profiles have to be supposed as in totally parabolic models and alternative turbulence models, such as the anisotropic Reynolds Stress Model (RSM) can be used. Traditional problems of elliptic models related to the estimation of the reference wind speed at each rotor position are mitigated due to the semi-parabolic algorithm. The model has been validated at the ECN test farm and at the offshore Horns Rev wind farm with significant results and also have been compared to other wake models. 1. Description of the semi-parabolic model The semi parabolic wake model is an in-house code based on the open source CFD solver OpenFOAM [26] and specifically designed for offshore wind farms. This is an Open Access article distributed under the terms of the Creative Commons Attribution License 2.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Article available at http://www.itm-conferences.org or http://dx.doi.org/10.1051/itmconf/20140206002 ITM Web of Conferences The model is inspired on the parabolic approach of other models such as UPMPARK [9] and WindFarmer [1, 2] but using the actuator disk technique to represent the wind turbine, instead of wind speed deficits as used in fully parabolic models [9, 10]. This way the model is solved elliptically but based on several subdomains instead of a single one, reducing thus significantly computational time. 1.1 Surface boundary layer modeling The free stream flow upstream of the wind turbines at the inlet corresponds to the fully developed vertical profiles in the surface boundary layer (SBL), characterized by a non-uniform shear boundary layer flow. OpenFOAM was adapted in order to solve the mean wind speed components and turbulence intensity according to the Monin-Obukhov theory, based on [11] and on the expressions that describe the SBL flow as presented in [12]. Firstly, turbulent viscosity is assumed to vary linearly with height:


ITM Web of Conferences
The model is inspired on the parabolic approach of other models such as UPMPARK [9] and WindFarmer [1,2] but using the actuator disk technique to represent the wind turbine, instead of wind speed deficits as used in fully parabolic models [9,10].This way the model is solved elliptically but based on several subdomains instead of a single one, reducing thus significantly computational time.

Surface boundary layer modeling
The free stream flow upstream of the wind turbines at the inlet corresponds to the fully developed vertical profiles in the surface boundary layer (SBL), characterized by a non-uniform shear boundary layer flow.OpenFOAM was adapted in order to solve the mean wind speed components and turbulence intensity according to the Monin-Obukhov theory, based on [11] and on the expressions that describe the SBL flow as presented in [12].
Firstly, turbulent viscosity is assumed to vary linearly with height: where t = turbulent viscosity (m 2 /s 3 ) u * = friction velocity (m/s) = von karman constant (= 0.4) z = height above ground (m).Neutral atmosphere is supposed, so that thermal effects are neglected, together with the Coriolis effects.Assuming shear stress to be constant over the surface boundary layer, a logarithmic velocity profile is defined at the entrance of the domain as: where z 0 = equivalent roughness length (m) and is Karman's constant (= 0.41).Assuming equilibrium between production and dissipation rates of the turbulent kinetic energy, the vertical profiles for the turbulent quantities can be derived as: for the turbulent kinetic energy, that turns out to be constant with height, where C is a constant of the standard k-turbulence model.On the other hand: for the dissipation rate of the turbulent kinetic energy.

Turbulence modeling
Turbulence closure is made for two turbulence models: the standard k turbulence model (isotropic) [13] and the Reynolds Stress Model (RSM, anisotropic) based on Launder and Gibson approach [14].
The first one is used on the conventional approximation of Richard and Hoxey [4] of getting the equilibrium of the set of constants (5) such that the inlet profiles (2) to (4) represents an exact solution of the Navier Stokes equations.
First Symposium on OpenFOAM in Wind Energy  This imposes the next expression to be fulfilled reaching a situation of equilibrium between the inlet profiles, the turbulence model and the wall: On the other hand, the next inlet profiles for uniform values of the turbulence components as defined in [7] are imposed when the RSM turbulence model is used instead: Since the constants of the model are the ones by default and a zeroGradient condition at the wall is imposed as a first approach for the Reynolds stress tensor components, the inlet profiles are not in full equilibrium with the wall.An equivalent set of constants similar to (5) should be adjusted for that purpose, although this task is out the scope of this work.Figures 1a and 1b shows respectively the difference between the inlet and the outlet profiles over flat terrain (length = 1000 m, width = 100 m, height = 100 m) using the k and RSM turbulence models.

Actuator disk model
From the linear momentum theory in its simplest form, it can be derived that the axial force from which the wind turbine extracts the kinetic energy from the incoming air flow, is just a function of the local induction factor or alternatively, of the thrust coefficient for the upstream wind speed [6][7][8].The wind turbine is thus considered as an actuator disk upon which a uniform distribution of axial forces, defined as momentum sink terms, is applied.These local axial forces can be prescribed over the disk area as: where: D = rotor diameter (m 2 ) V inf = upstream wind speed (m/s), obtained at this case from the upstream mast C t = thrust coefficient, obtained from the C t -velocity curve for the corresponding value of V inf = air density (kg/m 3 ).The local prescribed force is to be applied on a volume of cells defining the rotor disk.Total thrust is divided in the cells in the rotor from the ratio between the volume of each cell related to the total volume of the rotor disk.This is: where: dF i = thrust over cell i (N) dV i = volume of cell i (m 3 ) V = volume of rotor disk (m 3 ).

Numerical approach
The solution algorithm consists of a decomposition of the domain into a finite number of adjacent subdomains as shown in Figure 2.These sub-domains are solved elliptically in the axial direction, using the output variables of the previous sub-domain to be mapped onto the inlet of the subsequent sub-domain.This is done till the last wind turbine is detected at the end of the process.This way the computational time becomes significantly lower in comparison to the solution of an equivalent single domain by means of an elliptic approach.The model distinguishes thus between two types of subdomains: a free-stream-type if no wind turbine is present and a wind-turbine-type if one or more wind turbines co-exist in the same sub-domain.
The splitting algorithm of the domain into sub-domains preserves a minimum separation of 2D distance between the wind turbine and the inlet as well as between the last wind turbine positions and the outlet.This way, the flow solution ensures pressure jump at the position of the disks and no singularity is reached close to the boundaries.Nevertheless, the model is open to change the upstream distance by the user, although 2D is highly recommended.
A grid independence analysis, executed for an equivalent elliptical case in [7], reached to the conclusion that an horizontal (axial and transversal) resolution of 0.1D was found to be optimal.Using then this resolution, free-stream sub-domains of 4 cells and wind-turbine sub-domains of 40 cells (20 cells upstream and 20 cells downstream) were created.
Boundary conditions for each sub-domain can be summarized as follows: Inlet: Fixed values, mapped from the outlet of the previous sub-domain (except for the first sub-domain where free-stream SBL velocity and turbulence profiles -Eqs.( 2) to (4) -are assigned).

06002-p.4
First Symposium on OpenFOAM in Wind Energy Outlet: zeroGradient condition, which means no variation in the direction normal to the boundary.

FS WT FS FS FS FS FS FS FS FS FS T W T W
Lateral: symmetryPlane condition, which means zero gradient and no component of the velocity normal to the boundary.By default, the semi-parabolic model keeps the lateral boundaries at 8D apart from the last turbine.
Ground: Modified wall functions according to Blocken et.al. [15], ensuring thus continuity between the log-law velocity at the first cell next to the wall and the SBL log profile (2).
Top: fixed velocity and turbulence values, equivalent to fixed shear stress in the surface boundary layer.By default, the top is located at a height of 10 H above ground, where H is turbine hub height.
The current version of the model is still based on a first order upwind numerical scheme.It starts solving the first free-stream sub-domain and then flow variables are mapped sequentially from one sub-domain to the subsequent one.When the simulation finds a wind-turbine type sub-domain, the calculation is made in two phases.The first one solves the sub-domain without wind turbines till convergence is reached such that reference wind speeds at the positions of the wind turbines are estimated.These reference wind speeds are then used to calculate the momentum sinks at the position of the disks.A second calculation is then made with the sinks terms as described in Equations ( 5) and ( 6), activated at those cells coincident with the positions of turbines.
The current model is able to solve the domain for one wind speed and one wind direction.An extension of the model to solve the whole range of wind speed and wind direction bins is under development.

Set up of the model
As in similar models, the semi-parabolic model has to be first set-up by filling an input-data file containing the main parameters related to the free-stream flow as well as wind farm features.The next table shows a description of the parameters which should be assigned before launching the model.
The input-data file is located in a directory called input/, that must also contain the wind turbine coordinates [UTM_X UTM_Y] and two text files, one containing the power curve [WS POWER] and the thrust coefficient [V Ct].
The execution of the model consists of the three typical blocks: pre-processing, run and postprocessing.The three of them are totally automatized and minimum user interaction is required.Number of iterations for wind-turbine-tupe sub-domain.
Pre-processing is internally organized at the same time in two parts: 1. Pre-process 1: A splitting algorithm first rotates the wind farm to the wind direction specified in the input file such that the wind farm is aligned to the horizontal axis.The inlet flow is thus always aligned with the horizontal direction so that it is not necessary to decompose the wind speed vector.Once aligned, the code generates a text file with 5 columns related to the main information about the sub-domains to be simulated: • Relative number of sub-domain • Label of wind turbine as numbered in the original file with the UTM coordinates.If the sub-domain is free-stream, label is equal to −1.
Apart of this text file, some other auxiliary OpenFOAM files are generated, such as files for the automatic generation of the grids, for the creation of the cell zones corresponding to the rotor disks, for the assignment of disks features and for the decomposition process in case that more than one processor is available and a parallel simulation wants to be launched.

Pre-process 2:
The second pre-process generates a shell script with all commands to run the sequence of subdomains.The script follows the logical order, in case of each wind-turbine type sub-domain: • Map inlet boundary conditions from previous sub-domain • Grid generation • Run simulation with wind turbines deactivated • Read reference wind speed • Run simulation with wind turbines activated.
And in case of each free-stream type sub-domain: • Map inlet boundary conditions from previous sub-domain The script is launched and the simulation starts and continues till the last wind turbine is found.
Each simulation corresponding to a different sub-domain is stored in a different folder such that any information can be retrieved afterwards.
The post-processing is carried out through another script that directly access and reads the reference wind speeds as well as turbulence kinetic energy at the disks position.Using the values of wind speed, output power is derived for each wind turbine from its power curve.

06002-p.6
First Symposium on OpenFOAM in Wind Energy

Decrease in computational time
Parabollization of wake models is justified for far wake calculations inside big wind farms located in flat or moderately complex terrain.This brings as a consequence that the time needed to inverse smaller matrixes of any system is considerably lower.This means that it is then more profitable to run and solve a sequence of several small matrixes rather than solving one single huge matrix, leading to the same level of accuracy and convergence.
In order to illustrate this, a simple simulation corresponding to the Sexbierum case [16] (case of one isolated wind turbine) is run both in elliptic (one single domain) and in semi-parabolic mode (11 free-stream-type sub-domains +1 wind-turbine-type sub-domain).In order to make them comparable, they have been solved in one single processor and using the same boundary conditions, grid resolution and numerical schemes.
Table 2 shows that for this test case, the semi-parabolic model is able to solve the wake flow more than 4 times faster than the elliptic model.This ratio is expected to be higher for bigger domains.

Results
The semi-parabolic wake model has been validated in two phases.The first one demonstrates that wake flow solved using this approach is comparable to the elliptic case in terms of axial wind speed and turbulence.
The second phase demonstrates the ability to reproduce the power deficit along rows of wind turbines on two existing wind farms: the ECN test wind farm (coastal, quasi offshore) and the offshore wind farm Horns Rev.

Wind speed and turbulence profiles
In order to ensure the equivalence between both approaches, the transversal evolution of axial wind speed and turbulence kinetic energy of the Sexbierum test case [7,16] at 3 different downstream positions [2.5D, 5.5D and 8D] have been compared.
This has been done for the elliptic and the semi-parabolic versions of the model, and using both turbulence models: the isotropic k and the anisotropic RSM.
Figures 3 to 6 show the transversal profiles of wind speed deficit against lateral position [m] at the three positions mentioned before as described by the k turbulence model.As it can be observed, the semi-parabolic model reproduces accurately the maximum wind speed deficit at the centre of the wake as well as the wake width at downstream positions.This implies that the model can be applied with high degree of realibility to wake merge environments at wind farm areas.
In the same way, transversal profiles of wind speed deficit were extracted and compared to the elliptic case when using the anisotropic RSM model.Although a slight discrepancy was observed at 5.5D and 8D downstream at the centre of the wake, the flow is equally maintained for both cases, ensuring thus the validity of the model for the anisotropic treatment of the turbulence.
This small discrepancy between the parabolic and elliptic calculations is mainly apparent for the RSM.It may be due to the modeling of the pressure-strain interaction, whose role is to redistribute the  turbulent kinetic energy over the three normal stresses.This term is only in the RSM and not in the k-.It has the pressure, and pressure variations near the wind turbines may affect the parabolic approximation.With a larger size of the WT sub-domain this discrepancy is reduced.

Validation at the ECN test farm
Once the model has proved to be valid for an isolated wake case, a first result is obtained for a small wind farm like the ECN test farm [17][18][19], located at flat coastal terrain and composed by 5 wind turbines Nordex N80 of 2.5 MW in a row.Separation between wind turbines is 3.8D.
The main variables of interest consist of the power ratio of each wind turbine normalized with the power of the first one.Inlet conditions are considered, that are neutral atmosphere with a free-stream turbulence intensity equal to 7%, wind speed of 8 m/s and wind direction of 278 • , which fits with the exact row case where the 5 wind turbines are aligned in the same row.
This gives as a result maximum values of power deficits at the positions of the wind turbines.Since no thrust measurements were registered below 5 m/s, a similar value was considered from the cut-in wind speed in order to estimate the power ratio produced by the RSM turbulence model.

06002-p.8
First Symposium on OpenFOAM in Wind Energy  The first result is an overestimation of the power ratio at the second turbine according to the k turbulence model.This is an expected result for 3.8D separation since for isolated wakes in the Sexbierum test case and for similar distances [in the range 2.5D-5.5D] in the surrounding of the transitional region between the near wake and the far wake, an actuator disk approach based on an isotropic turbulence model leads to an overestimation of the power deficit at these downstream positions.
Since the load at the second actuator disk is less than expected the thrust and consequently the axial momentum extracted is smaller so wind speed obtained at the subsequent rotor disks tend to be gradually uniform.This gives as a result accurate agreement at wind turbines located downstream of the second wind turbine, at least for 3.8D distances and for the k turbulence model used at this case.On the other hand, the RSM turbulence model produces a general underestimation of the power ratio, a deeper impact at the second turbine and a recovery at the subsequent turbines.

Validation at the Horns Rev offshore wind farm
In order to complete the validation process of the semi-parabollic model, a real operating offshore wind farm like Horns Rev [1,2] was used to analyze the wake flow behavior when using this approach.Horns Rev is composed by an array of 8 × 10 Vestas V90 m [2 MW] wind turbines, with a hub height and rotor diameter of 80 m.
The flow case named as 1.8.1 was considered, which corresponds to a free-stream wind speed of 8 ± 0.5 m/s and a wind direction of 270 • ± 2.5 • at hub height (flow and wind farm aligned).The 06002-p.9 free-stream turbulence intensity measured at M2 met mast (located upstream of the wind farm) at 62 m height is considered as 9.53%.Due to symmetry for the considered flow case, just three wind turbines rows were simulated and the central row was post processed (experimental values correspond to the average of central rows).
Figure 7 shows the axial evolution of the power ratio corresponding to the central row of turbines normalized to the first one.
As it is observed, the semi-parabolic model shows a decreasing tendency with a very good agreement for the first four turbines and an underestimation of the power ratio for the last group of turbines.
Results from a fully parabolic model approach [1] have been included in the same graph in order to demonstrate that the semi-parabollic model is in the state-of-the-art of wake models focused on the solution of big offshore wind farms.The results from the elliptic model [20] based on the 2D reference wind speed (i.e.wind speed read from a cell 2D upstream of each rotor disk) are also included.
A second run was carried out using the anisotropic RSM model and compared to the previous run based on the k turbulence model.
Results on power ratio in figure 9 show a deep impact of the wake produced by the first turbine in turbine 2, an important recovery between the 2nd and 3rd turbines and a uniform tendency for the remaining turbines downstream.
This recovery produced by the semi-parabolic RSM was also reproduced by the elliptic model as it was showed in [20] for the same case and also observed experimentally by Smith and Taylor in [21].A possible explanation could be based on the fact that the flow behind the first turbine and consequently impacting on the second turbine is more turbulent that the flow behind the second turbine, leading to a higher recovery of the flow in comparison to other positions.
Figure 10 partly explains the difference between both turbulence models.Higher turbulence levels are simulated by the isotropic model leading thus to a higher recovery of the flow and higher power ratio, as it was depicted in Figure 6.
For the subsequent turbines, the thrust at the rotors estimated from the wind speed at the centre of each wake reaches equilibrium with the anisotropy of the turbulence, leading to a uniform power value from the 3 rd turbine onwards.The tendency fully agrees with the experimental values of the power ratio.This result points out that the RSM turbulence model and the anisotropic treatment of the free-stream and wake turbulence can represent a promising alternative at the CFD simulation of wind turbine wakes.

Advantages, shortcomings and further improvements
The main advantages of using the semi-parabolic model for the solution of wake flow in big offshore wind farms can be summarized as follows: • The near wake is calculated instead of modeled because the rotor is specifically solved by imposing the momentum sink terms at the rotor disk.This way, it is not necessary to impose self-similar profiles of wind speed deficits and added turbulence adding at the position of the expanded wake at 2.5D downstream.• The existing problem on the elliptic models related to the estimation of the reference wind speed at the position of each rotor disk is mitigated since the disks are not present and they do not create any influence upstream when the reference wind speed has to be estimated.• As the model is based on an open source CFD solver, the solution can be parallelized on an unlimited number of processors, leading to an important decrease in computational time.This decrease is even higher for the semi-parabolic model which has demonstrated to show an additional decrease in computational time of the order of 10 in comparison to purely elliptic models.• Taking advantage of this computational speed-up, alternative anisotropic turbulence models such as the Reynolds Stress Model (RSM) can also be used.The use of this model has shown to be promising with a tendency in the experimental validation even better than the isotropic models.This represents an important step forward on the solution of the anisotropy of the wake and atmospheric flow with an open source CFD solver.

06002-p.12
First Symposium on OpenFOAM in Wind Energy On the other hand, there are particular identified shortcomings and still some improvements to be made: • The model, as proposed here, can deal with an unlimited number of parallel rows.The main limitation imposes a separation of no more than 4 diameters for two consecutive wind turbines in the axial direction for any inflow, which is mostly fulfilled in current offshore wind farms.• The outlet boundary condition imposes a zero gradient in the stream-wise direction for all the variables.Since this boundary can be valid for pressure, a more appropriate condition should be implemented and developed for velocity and turbulence quantities.A zero second-derivative could be a nice approach since a linear tendency has been mostly observed at the end of the near wake region.• More accurate numerical schemes such as second order upwind or even cubic approaches could be implemented and validated.• As it has been demonstrated, the anisotropic RSM model is still under development and the set-up of the free-stream wind flow model has to be improved.Appropriate modifications on the turbulence constants should be made in order to ensure equilibrium between the inlet profiles, the wall boundary condition and the turbulence model.The implementation of the RSM turbulence model to simulate SBL flows is then preliminary and it should be taken as a first start point for future development.

Figure 1a .
Figure 1a.Simulation of the surface boundary layer flow over flat terrain based on the k turbulence model.

Figure 1b .
Figure 1b.Simulation of the surface boundary layer flow over flat terrain based on the RSM turbulence model.

Figure 8 .
Figure 8. Power ratio normalized to the first turbine.Flow case 1.8.1 and k turbulence model.

Figure 9 .
Figure 9. Power ratio normalized to the first turbine.Flow case 1.8.1 and k /RSM turbulence model.

Figure 10 .
Figure 10.Turbulent kinetic energy values at wind turbine positions for the semi-parabollic model.Flow case 1.8.1 and k /RSM turbulence model.

Table 2 .
Computational times for the Sexbierum case.Elliptic vs semi-parabollic models.