Cfd Wake Modelling with a Bem Wind Turbine Sub-model

Modelling of wind farms using computational fluid dynamics (CFD) resolving the flow field around each wind turbine's blades on a moving computational grid is still too costly and time consuming in terms of computational capacity and effort. One strategy is to use sub-models for the wind turbines, and sub-grid models for turbulence production and dissipation to model the turbulent viscosity accurately enough to handle interaction of wakes in wind farms. A wind turbine sub-model, based on the Blade Momentum Theory, see Hansen (2008), has been implemented in an in-house CFD code, see Hallanger et al. (2002). The tangential and normal reaction forces from the wind turbine blades are distributed on the control volumes (CVs) at the wind turbine rotor location as sources in the conservation equations of momentum. The classical k − ε turbulence model of Launder and Spalding (1972) is implemented with sub-grid turbulence (SGT) model, see Sha and Launder (1979) and Sand and Salvesen (1994). Steady state CFD simulations were compared with flow and turbulence measurements in the wake of a model scale wind turbine, see Krogstad and Eriksen (2011). The simulated results compared best with experiments when stalling (boundary layer separation on the wind turbine blades) did not occur. The SGT model did improve turbulence level in the wake but seems to smear the wake flow structure. It should be noted that the simulations are carried out steady state not including flow oscillations caused by vortex shedding from tower and blades as they were in the experiments. Further improvement of the simulated velocity defect and turbulence level seems to rely on better parameter estimation to the SGT model, improvements to the SGT model, and possibly transient-instead of steady state simulations.


Introduction
Computational fluid dynamics (CFD) codes with advanced wind turbine models included are important for the prediction of the wind energy production when optimizing wind farm layout.Such codes can calculate the velocity defect caused by the front row of wind turbines, which reduces the wind energy potential for the wind turbines located downwind, and the entrainment of momentum into the wake, which increases the energy production potential for downwind turbines.

Statement of the problem
Modelling of wind farms using CFD, resolving the flow field around each wind turbine's blade and the interaction of wakes is still too costly and time consuming to be used in engineering calculations.One strategy is to use sub-models for the wind turbines representing them as reaction forces acting on the flow field within the conservation equations for momentum.Depending on the computational grid resolution and the distribution functions for these forces it may be necessary to include sub-grid models for turbulence production and dissipation within the governing equations for turbu-lence.This can be done in order to model the turbulent viscosity accurately enough to handle turbulent entrainment of momentum in wakes and their interactions in wind farms.

Previous work
There are several sub-models which can be used to approximate the reaction forces from the wind turbines acting on the flow field.A group of these sub-models are the actuator disc models with-and without rotational forces which may be uniform or distributed over the disc, the actuator line method, and the actuator surface model.The latter two may be most suited for simulating unsteady flow, however the increased accuracy of these methods requires more computational resources, see Sandersee et al. (2011).
A wind turbine model for energy production, based on the Boundary Element Method (BEM) theory, has been implemented in Fortran, see Sand (2011).It has recently been shown by Lu and Portè-Agel (2011), that the BEM theory can also be used as a basis for the actuator line method.This approach makes it possible to model some of the unsteady flow behaviour of the near wake flow field experienced by the downstream wind turbines.

Scope of work
The numerical implementation of the BEM theory, see Sand (2011), will be included as a wind turbine sub-model in an in-house CFD code, Hallanger et al. (2002).
The tangential and normal forces to the rotor plane, computed by the BEM theory, will be distributed over the control volumes (CVs) located at or close to the rotor plane of the modelled wind turbine.
The classical k − ε turbulence model of Launder and Spalding (1972) will be used, and sub-grid models of turbulent kinetic energy (TKE) and its rate of dissipation are implemented, see Sha and Launder (1979); Sand and Salvesen (1994).
The simulations will be carried out steady state and compared with experiments available in the literature.In particular, the simulations will be compared with the experimental measurement of flow velocity and TKE downstream (in the near wake) of a model scale wind turbine located in the Norwegian University of Science and Technology (NTNU) wind tunnel, see Krogstad and Eriksen (2011).

Governing Equations
Transient and stationary turbulent flows may be described by conservation equations for the variation of mass-weighted time-mean quantities, see Favre (1965).Essentially these conservation equations are assumed to have the same form as the time-averaged conservation equations describing incompressible turbulent flow.The problem is also governed by an equation of state and by boundary and initial conditions.

Conservation equations and turbulence model
Conservation of mass: where ρ is density, U is velocity and sub-index i refers to Cartesian space direction.
Conservation of momentum: where p is the pressure and the stress tensor σ ij is defined by The effective viscosity is defined by: where µ and µ t are the laminar and the turbulent viscosity, respectively.
The turbulence model used is the k − ε model, see Launder and Spalding (1972).In this model the turbulent viscosity is related to the turbulent kinetic energy k and its rate of dissipation ε by: where k and ε are governed by conservation equations.
Conservation of turbulent kinetic energy: where S k is additional source terms and G is production rate of TKE defined by: Conservation rate of dissipation rate of turbulent kinetic energy: Conservation of specific enthalpy: σ ϕ for ϕ = k, ε, h is the turbulent Schmidt number see Table 1.

Solution procedure for governing equations
The conservation and transport equations for the flow and the turbulence parameters can all be written in the form: Here ϕ is the variable considered and S ϕ is the source term of the corresponding conservation equation.
To solve these equations we use a finite volume method, see Patankar (1980) and Ferziger and Peric (1996).The calculation domain is divided into a discrete number of CVs.All variables are stored at the CV center (collocated variables).Equations of the general form (10) are integrated over the CVs using interpolation formulas for the variation of ϕ between the grid points.A second order scheme Van Leer (1974) is used for the convective fluxes, a second order central scheme is used for the diffusive fluxes.A fully implicit formulation of the discretized equations is used.The result is a set of algebraic equations, where each equation connects the value of a scalar variable in a point to its neighbour values.
The algebraic equations for each CV have the form: where sub-script P refers to the CV center, subscript L refers to the neighbour CV centers and N is the number of neighbours.N = 6 for a three-dimensional problem on a Cartesian grid.These algebraic equations are solved together with proper boundary and initial conditions by a matrix solver for all CVs within the computational domain.
The velocity components computed within each time step or iteration by solving the conservation equations of momentum are only used as a first guess at the velocity field.The final values of the velocity components, the density and the pressure fields are calculated using a variational procedure known as the SIMPLE Method, see Patankar and Spalding (1972).The conceptual idea behind this method is given by Chorin (1968).For details concerning an extension to compressibility, see Hjertager (1985).

Turbulence sub-grid model
The sub-grid models for generation of TKE and its rate of dissipation are basically taken from the flow model in large rod bundles by Sha and Launder (1979).It is assumed that the production rate of TKE, G, and its rate of dissipation ε are integrated over the wake of the wind turbine and distributed over the near field.The additional source terms for TKE and its rate of dissipation are in the distribution area set to respectively and The coefficient C T was estimated to be C T = 0.24 • C D from detailed simulations of wakes downstream cylinders, circular and quadratic, by Sand and Salvesen (1994).Here it is assumed that the drag coefficient C D can be replaced by the thrust coefficient (C T h ) of the wind turbine.The coefficient C ε was chosen to be 0.05 after comparisons between test simulations and the NTNU experiments by Krogstad and Eriksen (2011).
Equations ( 12) and ( 13) may need further improvements and adjustments of the parameters based on detailed simulations of wind turbines and experimental data from turbines on different scales.

Wind turbine sub-model
The wind turbine model is based on the BEM theory with the Prandtl tip speed correction and Glauerts correction of wind turbine thrust for high loading with the modification of Buhl (2005) to avoid oscillations in the computational procedure.In the following we outline the equations with reference to Figure 1.The trust force on each blade element is given by the normal force component to the rotor plane where sub-index i refers to blade element i, V rel to the relative velocity (seen from the from the moving blade element), ϕ the angle between the relative wind (wind seen from the moving blade element) and the rotor plane, c the cord length and dr radial blade element length.C l and C d are respectively lift and drag coefficients for the airfoil profile used for the blade element.
The torque about the rotor axis is given by the tangential force component multiplied by the elements radial position ) where r is the radial posititon.
The relative wind is given by and where ω is the angular velocity in rad/sec and V the undisturbed axial velocity.
The induced non-dimensional element velocity in axial direction is given by The induced non-dimensional element rotational velocity is given by The solidity is given by where N B is the number of blades.
The normal force coefficient to the rotor plane is given by The tangential force coefficient to the rotor plane is given by The Prandtl tip loss factor is given by The trust coefficient based on axial momentum theory and modified by the tip loss factor is given by For high loading, based on Glauberts empirical approach, Wilson (in Spera ( 2009)) recommends where a c ≥ 1/3.From the balance between blade element force from N B blades and momentum force normal the rotor plane substituting for C T hi from Eq. ( 25), defining and solving the second order equation gives for high loading the axial induction factor The minus sign is chosen when a i [a c , 1].In the simulations a c was set to 0.4.
The axial and tangential induction factors for each blade element are determined by iteration satisfying a tolerance of 0.001.The algorithm uses drag and lift coefficients from airfoil profiles as function of angle of attack (α i ) of the relative wind, and Reynolds number based on cord length and relative wind, see Sand (2011).The wind turbine blades cord length and twist angle (θ i ) distributions are approximated using the trapezoidal rule and the geometry of the wind turbine blades given in Krogstad and Eriksen (2011).
The wind turbine model is structured so that individual treatment of each blade element is possible.By doing this extension, as suggested in Moriarty and Hansen (2005), it will be possible to capture energy from non-uniform incoming wind fields.

Modeling of a wind turbine wake in a wind tunnel
The NTNU wind tunnel is simulated with the wind turbine model of Krogstad and Eriksen (2011) included.
The inlet velocity is uniform with U ref = 10m/s, an isentropic relative turbulent intensity of 0.3 % and turbulent integral length scale of 0.005 m.The boundary conditions at the walls were set to no slip.The boundary condition at the outlet was found by extrapolating using the gradient and the nearest cell value for all variables.
The maximum tip speed of the wind turbine model is 100 m/s, which is low enough to neglect Mach number (compressibility) effects.The temperature is assumed to be constant during the experimental recordings.In consequence the flow field is approximately incompressible, and the conservation equation for specific enthalpy is not solved.
The reaction forces from the wind turbine sub-grid model and the hub are included in the CFD model.The physical blockage of the wind turbine and the reaction force from the tower are neglected, assuming that the modifications to the flow field from these geometries are small compared with the effect of the wind turbine reacting forces on the flow field.
Since the simulations are carried out steady state, axial and rotational forces from the three blades are averaged over circumferential rings.The hub was represented with the drag force.The averaged forces from each circumferential ring are decomposed to give the forces acting on each CV representing the wind turbine rotor location and added as source terms in the momentum equations.
As input velocity to the momentum and turbulence sub-routines the upstream uniform velocity is used.This must be changed for simulation of wind parks when turbines may be laying in the wake of upstream turbines and the inlet velocities will vary over the rotor plane.

Parameters used in the simulations
The wind tunnel test section is 2.71 m wide, approximately 1.81 m high and 11.5 m long.
The wind turbine rotor center is located in the middle of the wind tunnel test section 0.817 m above the floor and 3.66 m downstream of the test section inlet.
The rotor diameter D=0.894 m.The hub has a diameter d ≈ 0.09m over an axial length of 0.37 m.The total length, including the hemispherical front and end, is 0.46 m.The hub was represented with a drag coefficient = 0.6, taken from Hoerner (1965), and with the projected area in the axial direction given byA h = π • d 2 /4.The 3 blades have the NREL S826 airfoil profile along the entire span.The wind turbine blade's twist and cord length distributions are given in Krogstad and Eriksen (2011).
In the CFD simulations the blade's twist angle and cord-lengths distributions were represented in the BEM routines with two different element resolutions.Both data sets represented the hub by one element, and the blade's twist and cord-lengths by 9 elements and 18 elements in the low-and high-resolution case respectively.Each blade element has equal radial length and the trapezoidal model was used to give the corresponding cord length and twist angle.
The lift and drag coefficients as function of angle of attack and Reynolds are computed using XFOIL, Drela and Youngren (2001).The four Reynolds numbers used in the CFD simulations by the BEM routines, based on cord length and velocity seen from the blade element, are: 40,000, 80,000, 100,000 and 200,000.The transition amplification factor N crit = 3.0 is used as recommended by Karlsen (2009).Three of these data sets are shown in Figures 18 and 19.
The numerical grid used in the simulations is Cartesian.In the length direction the grid is non-uniform.
The highest resolution is around the turbine where the grid is uniform.Towards the inlet and outlet the grid is stretched.In the two other directions the grid is uniform.
To investigate the grid dependence, initial computations with a sub-grid turbulence model are carried out using three different grid resolutions.The CV number in all three space direction are increased with 50 % between each resolution.Between each resolution the maximum and minimum axial velocities change less than 0.7 % along profiles 1, 3 and 5 rotor diameters downstream of the turbine.Since the changes in velocity are small, it is decided to use the medium grid with a resolution of 71 x 91 x 61 CVs corresponding to respectively the length, width and height of the wind tunnel in the simulations.
Plots of the simulations with tip speed ratio (TSR) equal to 6 are shown in Figure 2 and in Figure 3.
Figure 2: Simulated horizontal cut through the wake of the wind turbine model of Krogstad and Eriksen (2011).Inlet velocity is 10 m/s, and the tip speed ratio of the wind turbine is 6.The position of the turbine blades is given by a black line.

Numerical simulations
The tip speed ratio is defined as ω Three TSR; 3, 6 and 10, used in the experiments given in Krogstad and Eriksen (2011), are simulated.The velocity defect is defined as 1−U/U ref , where U is the time averaged velocity in the axial (x) direction in the wind tunnel.
The turbulent kinetic energy k is related to the fluctuating velocity components measured in the experiments by k = 1 2 • ((u x ) 2 + (u y ) 2 + (u z ) 2 ).The overbar denotes time average of the squared fluctuating velocity component u .
In the figures comparing the simulations with experiments x/D is the non-dimensional distance down- The power and thrust coefficients generated by the BEM theory are given in Figure 16 and Figure 17 where they are compared with experimental values given in Krogstad and Lund (2012).
Airfoil data used in the simulations are given in Figure 18 and Figure 19.

Discussion of simulation results
The simulations cover three tip speed ratios, 3, 6 and 10.The data shown are the velocity defect and the TKE visualized cross-wise of the wake 1D, 3D and 5D downstream of the wind turbine rotor plane.The simulations are compared with experiments carried out by Krogstad and Eriksen (2011).

Velocity defect for TSR=3
For TSR=3, the simulations do not agree well with the experiments regarding the velocity defect, see Figure 4.For this low TSR which is half of the optimal design TSR of the model wind turbine, see Krogstad and Lund (2012), the boundary layer on the blades separates and stalling occurs, decreasing the ratio between lift and        drag.Taking out the angle of attack from the BEM procedure used for the given inlet velocity (10 m/s) and relative turbulence intensity 0.3% the angles of attack for the 18 elements used are all, with exception of the tip element, in the stall range, see α(ii) in Table 2 and airfoil data in Figure 18 and Figure 19.
A stall model is not included in the simulations.It is therefore expected that the drag coefficients from the blade are not representative and should have been higher.The deviation of the velocity defect is largest in the center region, indicating that the drag coefficients for about half of the center elements should have been higher, see Table 2, where C d (ii) is limited by the value 0.251 given at attack angle 20 o .
The lower drag coefficient does also affect the power coefficient (C p ) distribution in the lower TSR range from 3 and down to 1 (see Figure 16, and the thrust coefficient (C T h ) distribution (Figure 17).An improvement may be possible by expanding the data sets for drag and lift as function of Reynolds number and attack angle to cover a range up to possibly 30 o .
In Figure 5 the sub-grid contributions S k and S ε are added respectively to the generation of TKE and its rate of dissipation.The maximum center values of the velocity defect decrease further, however at the shoulders of the distribution.i.e. closer to the end of the turbine blades the agreement with the experiments is improved.

Turbulent kinetic energy for TSR=3
The TKE in the wake for TSR=3 is given in Figure 6.The simulated k is much lower than the experimental values of Krogstad and Eriksen (2011).When the subgrid contributions S k and S ε are added respectively to the generation of TKE and its rate of dissipation, the results improve, see Figure 7. However the center values of the distribution in the wake are still low.
The numerical simulations are carried out steady state.As a consequence the vortex shedding from the tower and the helical vortex wake shed by the three blades will in the simulations not give an oscillating contribution to the flow field as the vortex structures would in the experiments as they were convected downstream past the measurement positions.In the experiments these vortex induced velocity oscillations were included in the instantaneous velocity measurements and will therefore contribute to the TKE when averaged over time.It is therefore reasonable to expect large differences between measured TKE and simulated turbulent kinetic energy by steady state methods.

Velocity defect for TSR=6
For TSR=6, the simulations agree well with the experiments regarding the velocity defect, see Figure 8.The best result is obtained 3D downwind of the rotor plane were the results fit very well and have the same maximum value as the experiments.At 1D the velocity defect deviates much in the center of the wake were it is too low.This may be a consequence of not including the wind turbine tower in the simulations since the tower generates a wake on its own.The TSR=6 is the design TSR of the wind turbine blades, see Krogstad and Lund (2012), and for the wind speed 10 m/s the wind turbine model is operating below or at rated speed.In this case the boundary layer on each blade stays attached to the blade.
In the CFD simulations the angles of attack of the flow field, seen from the blade elements at convergence of the induced velocities, are given in Table 3.The corresponding lift and drag coefficients are all well within the Reynolds number range chosen for the airfoil data.All the attack angles for the blade elements are much smaller than the approximately 14 degrees where stalling starts to develop with decreasing lift coefficient and increasing drag coefficient as function of further increase in the angle of attack.
From Figure 16 and Figure 17 we see that both the power coefficient C p and the thrust coefficient C T h agree very well with the value measured by Krogstad and Lund (2012) for TSR=6.It should be noted that blockage effect of the flow field caused by the turbine standing in the wind tunnel test section is not included in the BEM method used here to generate the results given in Figure 16 and Figure 17 and in Table 3, however the thrust axially will affect the flow field in the CFD simulations, so blockage is partly taken into consideration in the CFD simulations.
In Figure 9 the sub-grid contributions S k and S ε are added respectively to the generation of TKE and its rate of dissipation.The maximum center values of the velocity defect increase in this case.However, all the velocity defects seem to become approximately of the form of the velocity defect at 5D.This indicates that the sub-grid contribution to the turbulent viscosity is too large.

Turbulent kinetic energy TSR=6
The TKE in the wake for TSR=6 is given in Figure 10.The simulated k is much lower than the experimental values of Krogstad and Eriksen (2011) at 1D downwind of the rotor plane.At 3D and 5D the results are better but still far off.When the sub-grid contributions for S k and S ε are added the results improve, see Figure 11.TKE at 1D downwind is now equal to the

Velocity defect for TSR=10
For TSR=10, the simulations agree well with the experiments regarding the velocity defect, see Figure 12.The best result is obtained 3D downwind the rotor plane were the simulated defect fit very well.At 1D the velocity defect agrees surprisingly well with the experiments at the center values, and even with the shape of the experimental velocity defect.However the peak value of the velocity defect at the position 0.5D crosswise is approximately 10% lower than the experiments.This seems to be in agreement with the thrust coefficient C T h , see Figure 17 for TSR=10.At 5D the maximum size of the velocity defect is similar to the experiments, but the crosswise mixing of flow is less than in the experiments, since the negative local torque dM (ii) in the center for ii=3 and 4, see Table 4 which generates the low value of the velocity defect for 1D and 3D has not disappeared at 5D in contrast to the experiments.
In the CFD simulations the angles of attack of the flow field, seen from the blade elements at convergence of the induced velocities, are given in Table 4.The corresponding lift and drag coefficients are all well within the Reynolds number range chosen for the airfoil data.All the attack angles for the blade elements are much smaller than the approximately 14 o were stalling starts to develop with decreasing lift coefficient and increasing drag coefficient with further increase in the angle of attack.
The TSR=10 is not the design TSR of the wind turbine blades, see Krogstad and Eriksen (2011), and for the wind speed 10 m/s the wind turbine model is operating below or at rated speed.However the high TSR gives a very high loading on the turbine blade and the induced axial velocities are close to 1 giving local thrust coefficients close to 1.44.How well the airfoil data fits in this case, where parts of the blades lose angular momentum to the flow field while other parts of the blades catch angular momentum from the flow field, is uncertain.One might expect that the opposite circulation sign of neighbour elements will generate additional turbulence.
In Figure 13 the sub-grid contributions S k and S ε are added respectively to the generation of TKE and its rate of dissipation.The maximum center values of the velocity defect increases in this case at all wake positions (1D, 3D and 5D).However all the defects seem to become approximately of the form of the defect at 5D.Further sideways all the velocity defects decrease.
We conclude that the sub-grid contribution to the turbulent viscosity is too large.

Turbulent kinetic energy TSR=10
The TKE in the wake for TSR=10 is given in Figure 14.The simulated TKE is much lower than the experimental values of Krogstad and Eriksen (2011).When the sub-grid contributions S k and S ε are added respectively to the generation of turbulent kinetic energy and its rate of dissipation the results improve, see Figure 15.However the center values of the distribution in the wake are still low.See also the corresponding section for TSR=3 to the effect of vortex shedding.

Conclusions and Recommendations
A numerical implementation of the generalized Blade Element Momentum (BEM) theory in Fortran, see Sand (2011), has been included as a wind turbine sub-model in an in-house CFD code, Hallanger et al. (2002).
The tangential and normal forces to the rotor plane acting on the wind turbine blades are computed by the BEM theory.These forces are distributed as reaction forces acting on the CVs located at the rotor plane of the wind turbine modelled, i.e. added as sources in the conservation equations of momentum at this location.
Steady state CFD simulations, with the wind turbine sub-model included, were compared with flow and turbulence measurements carried out in the wake of a model scale wind turbine located in the NTNU wind tunnel, see Krogstad and Eriksen (2011).
The velocity defect simulated without including subgrid models of the source terms in the k and ε equations compares well with the experiments of Krogstad and Eriksen (2011), the exception being for TSR=3 where stalling occurs because of high angle of attack.The airfoil data for drag and lift used in the BEM procedure did not cover this case properly.The simulations for this case (TSR=3) should be repeated with better suited airfoil data (expended to cover attack angle up to possibly 30 o ).The velocity defect, simulated when sub-grid models are included in the source terms of the k and ε equations, compares less well with the experiments.The increased turbulence viscosity caused by the sub-grid model seems to smear the velocity defect distribution indicating that the parameters in the subgrid model need adjustment through comparison with results from detailed simulations of the wind turbine and its wake.
The TKE simulated without including sub-grid models of the source terms in the k and ε equations does  Krogstad and Eriksen (2011) and is far too low.The TKE, simulated when sub-grid models are included in the source terms of the k and ε equations, compares much better with the experimental data than without sub-grid models.However the simulated TKE is still low.
It should be noted that the simulations were carried out steady state so flow oscillations caused by vortex shedding from tower and blades are not included in the simulated TKE as they were in the experiments.
Further improvement of the simulated velocity defect and the TKE seem to rely on: detailed simulations of the wind turbine in order to estimate better parameters in the sub-grid models, an improved subgrid model, and possibly transient-instead of steady state simulations.
It should be noted that blockage effect of the flow field caused by the turbine standing in the wind tunnel test section is not included in the BEM method used here to generate the results given in Figure 16 and Figure 17.In general for an object in a channel there is a contribution to the blockage from the object, the wake of the object and from the boundary layer along the channel walls.In consequence the blockage will act as if the flow velocity were higher at the inlet of the wind tunnel test section than measured, in particular when the tip speed ratio is large enough to generate high thrust coefficient Rae and Pope (1984).However the thrust axially will affect the flow field in the CFD simulations, so blockage is partly taken into consideration in the CFD simulations.
The TSR=10 is not the design TSR of the wind turbine blades Krogstad and Lund (2012), and for the wind speed 10 m/s the wind turbine model is operating below or at rated speed.However the high TSR give a very high loading on the turbine blades and the induced axial velocities are close to 1 giving local thrust coefficients close to 1.44.How well the airfoil data fits in this case, where parts of the blades loses angular momentum to the flow field while other parts of the blades catches angular momentum from the flow field, is uncertain.One might expect that the opposite circulation sign of neighbour elements will generate additional turbulence, see Table 4, Figure 16 and Figure 17.It may be possible that XFOIL, which is used to generate the airfoil lift and drag data, is not including this rotational effect.
CFD codes with advanced wind turbine models included are important for the prediction of the wind energy potential when optimizing wind farm layout.The latest version of the wind turbine model includes a yaw model.This functionality may be an advantage for reducing wind stress on wind turbines downstream the first wind turbine line.way (RCN).NORCOWE is a consortium with partners from industry and science, hosted by Christian Michelsen Research.The information about wind turbine blade profiles and coefficients from John Amund Lund is gratefully appreciated.The help from Thomas Hansen in generating drag and lift coefficients using XFOIL is also acknowledged.Discussions and comments from collages within the CMR organisation and participants in NORCOWE are gratefully appreciated.

A. Power and Thrust of Wind
Turbine from Airfoil Data A.1.Power and thrust coefficients of wind turbine The power coefficient of the wind turbine model is defined as where the power is P = ω • M and M is the torque summed over all blade elements on all blades, about the rotor axis.The torque from one blade element is given by Eq. ( 15).The power coefficient as a function of TSR is given in Figure 16.
The thrust coefficient of the wind turbine model is defined as where F n is the normal component of the wind force summed over all blade elements on all blades.The normal force from one blade element is given by Eq. ( 14).
The trust coefficient as a function of TSR is given in Figure 17.

A.2. Lift and drag coefficients used in the simulations
The lift and drag coefficients are given as functions of angle of attack and Reynolds number in respectively Figure 18 and Figure 19.The drag force D is the force from the wind acting on the foil in the direction of V rel .The lift force is normal to V rel , see Figure 1.

Figure 1 :
Figure 1: Loading forces from the wind acting on blade element i, αu i is angle of attack, θ i is blade element twist angle.

Figure 3 :
Figure 3: Simulated rotational velocities 2D downstream the wind turbine model of Krogstadand Eriksen (2011).The cut is crosswise the wake.Inlet velocity is 10 m/s, and the tip speed ratio of the wind turbine is 6.Maximum rotational velocity is 1.27 m/s.The velocities are proportional to the vector length.

Figure 16 :
Figure 16: Simulated Power Coefficient as Function of TSR using the BEM Procedure Compared with Experiments, Krogstad and Eriksen (2011).NREL-S826 Airfoil Profile, U=10 m/s.

Figure 17 :
Figure 17: A Simulated Thrust Coefficient as Function of TSR using the BEM Procedure Compared with Experiments, Krogstad and Eriksen (2011).NREL-S826 Airfoil Profile, U=10 m/s.

Figure 18 :
Figure 18: Lift coefficients of the NREL-S826 airfoil as function of attack angle and Reynolds number.N crit = 3.0.

Figure 19 :
Figure 19: Drag coefficient of the NREL-S826 airfoil as function of attack angle and Reynolds number.N crit = 3.0.

Table 1 :
Constants appearing in the governing equations.

Table 2 :
Parameters used within the CFD code by the blade elements, U ref = 10m/s and TSR=3, generated blade element force dFn(ii) and blade element torque dM (ii) both from one blade.

Table 3 :
Parameters used within the CFD code by the blade elements, U ref = 10m/s and TSR=6, generated blade element force dFn(ii) and blade element torque dM (ii) both from one blade.ii

Table 4 :
Parameters used within the CFD code by the blade elements, U ref = 10m/s and TSR=10, generated blade element force dFn(ii) and blade element torque dM (ii) both from one blade.