OffWindSolver: Wind farm design tool based on actuator line/actuator disk concept in OpenFoam architecture

Wind energy is a good alternative to meet the energy requirements in some parts of the world; however the efficiency of wind farm depends on the optimized location of the wind turbines. Therefore a software tool that is capable of predicting the in-situ performance of multiple turbine installations in different operating conditions with reliable accuracy is needed. In present study wind farm layout design tool OffWindSolver is developed within the OpenFoam architecture. Unsteady PisoFoam solver is extended to account for wind turbines, where each turbine is modeled as a sink term in the momentum equation. Turbine modeling is based on actuator line concepts derived from SOWFA code, where each blade of the turbine is represented as a line. The loading on each line/blade of the turbine is estimated using the Blade Element Method (BEM). The inputs for the solver are tabulated airfoil aerodynamic data, dimension and height of the wind turbines, wind magnitude and direction. OffWindSolver is validated for a real wind farm – Lillgrund offshore facility in Sweden/Denmark operated by Vattenfall Vindkraft AB. Because of the scale of the computation, we only examine the effect of wind from one direction at one speed. In the absence of time dependent Marine Atmospheric Boundary Layer (MABL), a log wind profile with surface roughness of 0.04 is used at the inlet. The simulated power production of each turbine is compared to the field data and large-eddy simulation. The overall power of the wind farm is well predicted. The simulation shows the significant decreases of the power for those turbines that were in the wake.


Introduction
The flow around the turbine is three dimensional due to centrifugal and Coriolis forces. Turbines installed in the flow extract energy from the wind, which affects the flow behind the wind turbines and creates wakes and turbulence. The approaching wake from neighboring turbines has a lower velocity and a higher turbulence intensity, which reduces the power production and increases the unsteady loads. The wake is classified into two categories near wake and far wake. The near wake is the area just behind the turbine rotor which is approximately up to one to three rotor diameters downstream. The near wake region, dominated by intense turbulence, mainly depends on a type of turbine, number of blades, aerodynamic shape of the blades and dimensions of turbine blades. The far wake region is a region beyond the near wake. The near wake research is very turbine specific, while the far wake research is focused on the mutual interactions of the turbines, when they are placed in clusters, like wind ITM Web of Conferences farms [3]. The far wake influences the turbine performance and overall power output of a wind farm. Therefore the main attention should be focused on far wake models, wake-wake interference, turbulence models etc. [1,2]. The wakes studies have been a topic of intensive research, both experimentally and numerically during the last decades and have been reviewed by Vermeer et al. [3].
The wake model can be classified either explicit or implicit. The explicit or kinematic wake models are the simplest and use self-similar velocity profile [4][5][6][7][8][9]. These models are based on the closed form solution for the width of the wake and the mean velocity profile in the wakes. The power prediction is reasonably accurate with these models. These explicit engineering models are used very often in industries because of less computational time and approximately good power prediction. However, these simple models do not account for detailed physics associated with wake and wake-wake interaction. Furthermore, these models also require tuning of parameters and this tuning varies from case to case. The other wake models are called field model where the modified governing flow equations are solved. The field models solve the Parabolized Navier-Stokes (PNS) equations with a turbulence model as a closure. These models are slightly computationally expensive than the explicit model, but the predictions with these models are better compared to the simple explicit model. Drawback with the model is the tuning of parameters and also it requires good initialization of wake in a case with multiple wakes. Another advanced model is elliptic field /full CFD, which solves an averaged NS equation, where the turbines are modeled as a sink term in the momentum equation. For the sink term in the momentum equation, either actuator disk (AD) or actuator line (AL) approach is used. The AD approach assumes the turbine rotor as a porous medium, and the AL approach resolves each blade of the turbine as a line or surface. In the AL model the turbine blades are represented by lines upon which a distribution of forces acts as a function of local incoming flow and blade geometry [1,2]. A main advantage of actuator line model is the representing the blades by its airfoil data. That is why it makes the approach well suited for wake studies [10]. However the biggest challenge with the AL model is that it requires many input parameters mainly the aerodynamic coefficients of the blade airfoil. Therefore, in the absence of detailed information, AD approach is used.
In the present study, a methodology is developed for simulating the whole flow around the turbine, including the near and the far wake. A wind farm layout design tool within OpenFoam architecture based on the AD and AL model is developed. Unsteady PisoFoam solver is extended to account for the wind turbine wakes and its interaction. The developed tool OffWindSolver has two options for wake modeling (1) model based on the AD and (2) based on Actuator Line (AL). If a user has fewer parameters i.e. tip radius, hub radius, thrust/power coefficients, height of the tower etc. then the actuator disk (AD) model can be used. Otherwise if the user has detailed input parameter [11,12] then the AL model can be used. The biggest challenge with the AL model is the shape and profile of the blade which is normally proprietary of a turbine manufacturer and difficult to obtain.

Theory
Turbine efficiency remains a critical component of the overall economic justification [13] for a potential wind farm, however the interaction of multiple wind turbine is more critical. Therefore a requirement for prediction methodologies that is capable of addressing the in-situ performance of multiple turbines is justified. The tool should be reliable, accurate and computationally less expensive. Simulation tools based on the explicit models are sufficient for small to medium sized wind farm. However for large size wind farm modeling better prediction tools are needed. In the present study a tool based on CFD is proposed. Where the actual geometry of the blades and viscous flow around the blades are not resolved, and the turbines are modeled with momentum sink/source approach (see Eq. (1)). The sink term (S) as shown in Eq. (1) is modeled in two different ways: (1) Actuator Disk (AD) and (2) Actuator Line (AL).

04001-p.2
First Symposium on OpenFOAM in Wind Energy The generalized N-S equation is as follows.
Where is the density, u is the flow velocity, is the shear stress, x i, x j are the directions, and p is pressure.

Actuator disk model
Based on one dimensional momentum theory and assuming no losses Rankine showed that the air speed through a propeller disc is half the sum of the air speeds far upstream and far downstream. Based on the control volume approach the thrust T and power P can be expressed as [14]: Where, V up is the upstream velocity, a is the axial interference factor, A is the area of the disk and V is the volume of the cell, C p and C T are power and thrust coefficients respectively. The sink term in momentum equation is implemented with two approaches. First approach is based on thrust Eq. (2) and second approach is based on Erik Svenning actuator disk model [15].

Actuator line model
AL based on Blade Element method (BEM) was proposed by Sørensen and Shen's [1,2]. The main philosophies in the AL model are the representation of a blade as a line, and the projection of the point force as a volume force. Figure 1 shows the geometrical representation of the wind turbine and their blades surface. The approaching flow introduces local forces (Lift and Drag) depending on the relative velocity, twist angle of the blade, the blade chord and the local angle of attack. In the AL, the blade is represented as a line and the blade forces are projected on a line. The linear force on the each element of the blades/line is decomposed into normal and tangential components denoted by F n and F t [11,12,16].

ITM Web of Conferences
A blade profiles section A-A' is plotted in Figure 1. The fluid velocity relative to the blade, t is decomposed into a normal component U n and a tangential component U t . The local angle of attack ( ) is calculated from the U n , U t and twist angle ( ).
Based on the local angle of attack and relative velocity the lift (L) and drag (D) forces per unit span are calculated. Where Where, controls the width of Gaussian and |r| is the distance between blade point and CFD cell.

Implementation in OpenFoam
OffWindSolver is derived from PisoFoam. This is achieved by adding a sink term in the momentum equation. The momentum sink term is based on the SOWFA [11,12] philosophy. The existing SOWFA code is Large Eddy simulation (LES) tool and mainly based on AL. SOWFA is computationally demanding and requires a large number of input parameters, specially aerodynamic coefficient of the blades. SOWFA is a good research tool but applicability for industrial problem is limited due to computational constraint and model input requirement. However, OffWindSolver is based on PisoFoam and can be used for both LES and RANS. Furthermore, it has AD model which does not require many input parameters. In the AD approach, the major input to the model are thrust/power curve, tip and hub radius, wind direction and location of the wind turbine. In the actuator line model, the input parameters are turbines location, airfoil aerodynamic data, hub and tip radius, base location, turbine height, rotational speed and direction. The major outputs we get are the forces, torque, power etc. on the each turbine and total power output of the wind farms. The AL model is derived from the SOWFA code [11,12]. Main objective of OffWindSolver was to develop a tool for industrial applications.

Validation and model setup
The OffWindSolver has two options for modeling the wake, however in the present study only advanced actuator line model is validated. The actuator disk model will be validated in future. The simulated wind farm is Lillgrund offshore wind farm [17] Figure 2 shows the plan view of the distribution of the wind turbines in the farm. The position of wind turbines is given in the report from Jeppersson et al. [17]. In our modeling the case set up was taken from SOWFA [11,12]. The computational domain was rotated along the wind direction as shown in Figure 2 to simplify the inlet and outlet boundary conditions.

Computation domain and boundary condition
The computational domain used for the study is 4500 × 3500 × 800 m in axial, lateral and vertical direction respectively. The numbers of mesh points used in axial, lateral and vertical direction are 370, 200 and 80 respectively. The wind-wave interaction plays a major role while modeling the offshore wind turbines. A coupling of the wind-wave interaction module and the wake models will be carried out in future.
The sea surface waves are not modeled explicitly. The water surface is assumed as a flat surface without any roughness. A no-slip boundary condition is used on the bottom water surface. On the sides and top surface of computational domain free slip boundary condition is applied. At the exit of the computational domain outlet boundary condition is used. A logarithm velocity profile as shown in Figure 3 and given in Eq. (10) is used to approximate the marine atmospheric boundary layer (MABL).
Where is the von Kármán constant, z is height, z o is the aerodynamic roughness length and u * is the friction wind speed defined by √ / is the wind stress, and is air density. The turbulence model used for the study is the standard k-without modifying the dissipation rate equation [11,12,16]. The source term in turbulent equation due to rotation of blades and also due to meandering [12] is not accounted in the present study.
In a reality, turbulent kinetic energy and dissipation vary along the height. The variation in turbulent kinetic energy and dissipation depend on the atmospheric condition. However in present study a constant fixed value of turbulent kinetic energy and dissipation rate were used at the inlet of the domain. The inlet turbulent intensity was 0.15 and the turbulent viscosity ratio was 0.2. Furthermore standard wall functions were employed without correction to sea surface roughness [18]. 04001-p.5  In present study neither the rigorous grid sensitivity studies nor the effect of model parameter associated with AL were performed. These parameters will strongly influence the turbulence and the velocity deficit behind the wind turbines. Nevertheless, some simulations with different meshes were performed until the power outputs of turbines were independent of the mesh size.

ITM Web of Conferences
In CFD studies the spatial and temporal discretization are the main source of the discretization error. Generalized Richardson extrapolation [19] is normally used for estimation the discretization error. Performing such studies for a large computation domain such as wind farm is beyond the scope of present work. Main objective of the present work was to introduce the AL model within OpenFoam 04001-p. 6 First Symposium on OpenFOAM in Wind Energy  architecture and use this approach to simulate one real size wind park. The model parametric studies and Richardson extrapolation will be performed in future studies. Figure 4 shows the iso surface of vorticity magnitude for a turbine. It can be observed that the rotating turbine blades induce the rotational component in the downstream flow.

Results and discussions
This will not only reduce the power production of an upstream turbine but it will also induce the unsteady loads. The formation of vortex structures on the down stream wind turbine depends on the 04001-p.7 relative distance from upstream wind turbine and tip speed ratio [20]. The pressure distribution around the wind plant (see Fig. 5) is obtained with OpenFoam implemented solver, OffWindSolver. A pressure jump across the wind turbine is clearly visible. The pressure is higher ahead of the turbine and lower behind the turbine. Figure 5 shows the existence of the wakes behind the wind turbines. The wakes are clearly visible. Near the center of the rotor, it is noticed powerful pulsating vortex that is more powerful 04001-p.8

ITM Web of Conferences
First Symposium on OpenFOAM in Wind Energy for front turbines than all other members of the row. The pressure distribution looks much smeared because of the poor resolution of the grids; there was no attempt on local refinement of the mesh close to the wind turbine. This will be addressed in future studies. Figure 6 depicts the contours of streamwise velocity at the hub height. Near the center of the rotors there is a higher speed region of the flow. The reduced velocity from the upstream turbines reaches to the downstream turbine which subsequently reduces the powers of downstream turbines. The region of higher-speed fluid persist for a longer distance behind the front turbines than all other member of the row [12]. Figure 7 shows contour of resolved turbulent kinetic energy at the hub-height. Though not clearly, but still higher turbulent kinetic energy is visible on the edges of wake because of the higher velocity gradient in lateral direction.
Previous LES simulation [12] obsreved the greatest turbulent kinetic energy in the downstream turbines due to the meandering of the downstream wakes. Which is not seen in the simulation with the RANS based approach. Nevertheless, we still noticed a large value of turbulent kinetic energy in the rows (B and C), which has more wind turbines. Figure 8 shows velocity, turbulent kinetic energy and pressure contour plot in the vertical direction, the plan is in the center of the turbines in row D (refer to Fig. 2). Most of the same behavior seen in the horizontal contour planes is seen here. It can be observed a gradual increase of the velocity deficit with downstream distance. An increased turbulent kinetic can be seen on the edges of the wake as explained earlier and also noticed with LES [12], however noticeable increased in turbulent kinetic energy with downstream distance is not visible with RANS based approach, these results calls for accounting a source term in turbulence modeling because of meandering [12]. The pressure look more smeared therefore a grid refinement studies require to address this problem.
One of the targets of this study was to explore the ability of the different methods to evaluate the power production of the wind turbine. Our simulations are compared to the production data (presented 04001-p.9 ITM Web of Conferences by Dahlberg [21]) and LES simulation data (SOWFA solver [12]). Figure 9 shows the power produced by each turbine normalized by the power of the corresponding of the first turbine of the row. The power prediction with OffWindSolver compares well with SOWFA and also with field data.
At Lillgrund the turbines spacing is very low for today's standards: 4.3 rotor diameters along the rows and 3.3 diameters between rows. In consequence the power decreases dramatically because of the presence of the wake behind the turbines, which create velocity deficit.

Conclusions
An OffWindSolver within OpenFoam architecture is developed for wind farm layout design. AD and AL models are implemented in the OffWindSolver. The Model selection will be based on the available input information and the required accuracy. Discretization error because of mesh size and AL model parameter is not estimated here. These studies will be addressed in future work.
In present study advance AL wake model is validated for the Lillgrund wind farm. The overall wind farm power production is predicted well with the OffWindSolver. The effect of wakes on power production is dramatic: the wind turbines affected by upstream wake have only 30-40% of the power of the first turbines in the row. The turbulence kinetic energy is not predicted well due to the approximated turbulence intensity and the meandering.