Validation of the actuator disc approach using small scale model wind turbines

The aim of the present study is the validation of the implementation of an actuator disc model in the CFD code PHOENICS. The flow behaviour for three wind turbine cases is investigated numerically and compared to wind tunnel measurements: (A) the flow around a single model wind turbine, (B) the wake interaction between two in–line model wind turbines for a uniform inflow of low turbulence intensity and (C) the wake interaction between two in–line model wind turbines at different separation distances in a uniform or sheared inflow of high turbulence intensity. This is done by using Reynolds av5 eraged Navier–Stokes (RANS) and an actuator disc (ACD) technique in the computational fluid dynamics code PHOENICS. The computations are conducted for the design condition of the rotors using four different turbulence closure models and five different thrust distributions. The computed axial velocity field as well as the turbulent kinetic energy are compared with hot wire anemometry (HWA) measurements. For the in–line two model wind turbine cases, the thrust coefficient is also computed and compared with measurements. The results show that for different inflow conditions and wind turbine spacings the proposed 10 method is able to predict the overall behaviour of the flow with low computational effort.


Introduction
The study of wake properties is important for assessing the optimal layout of modern wind farms.Wind turbine wake development may be studied using field experiments, smallscale wind tunnel measurements or numerical simulations with computational fluid dynamics (CFD).There are several advantages of CFD over field experiments and small-scale wind tunnel measurements, e.g.no violation of similarity requirements, control over inflow conditions and information about the relevant parameters, e.g.wind velocity, over the entire flow.However, as CFD results are sensitive to the experience and knowledge of the user of the CFD code and to the numerous computational parameters and assumptions involved, it is imperative to perform validation studies.Previous work on validating CFD wake models using a wind turbine tested in wind tunnels has been presented by Simms et al. (2001) and by Schepers et al. (2012).These studies demonstrated that there was a significant deviation between the various prediction tools and the wind tunnel measurements.Similar results for a small-scale model wind turbine are reported by Krogstad and Eriksen (2013) and by Pierella et al. (2014), indicating the importance of validating existing wind turbine modelling tools and methodologies.
Advanced methods of wake modelling with CFD may be implemented by using large-eddy simulation (LES) techniques in which the wind turbine forces may either be prescribed with an actuator line (ACL) method or an actuator disc (ACD) method.Work along these lines has been performed by numerous researchers such as Breton et al. (2014), Nilsson et al. (2015), Churchfield et al. (2012), Andersen et al. (2015), Calaf et al. (2010) and Ivanell et al. (2007).
Although LES provides high-fidelity results comparable to field measurements, the computational requirements of the method (Churchfield et al., 2012;Laan et al., 2015b) are still too expensive and therefore not yet suitable for engineering Published by Copernicus Publications on behalf of the European Academy of Wind Energy e.V.
N. Simisiroglou et al.: Validation of the actuator disc approach practices of wake computations for whole wind farms.A less computationally expensive alternative to LES are Reynoldsaveraged Navier-Stokes (RANS) simulations.RANS simulations have been used with the ACD method to simulate wind turbine wakes by numerous researchers, e.g.Laan et al. (2015a), Prospathopoulos et al. (2011), andEl Kasmi andMasson (2008).
The aim of the present study is the validation of the implementation of an ACD model in the CFD code PHOENICS (Spalding, 1981).In order to do so, computational and experimental results are compared for three cases.Case A consists of a single wind turbine in a low-turbulence-intensity environment with a uniform wind inflow.Case B is composed of two wind turbines positioned in-line in the same low-turbulence-intensity environment with a uniform wind inflow.Case C again uses two wind turbines positioned inline, but in this case multiple inflow conditions are studied for different spacings of the wind turbines.This is done to investigate the influence of the inlet conditions on the wind turbines' thrust.
As this method is intended to be used for industrial purposes, it therefore needs to provide accurate and reliable results with low computational effort.The simulations are performed according to the "Blind test 1", "Blind test 2" and "Blind test 4" invitation workshops organized by NOWITECH and NORCOWE (Krogstad et al., 2011;Pierella et al., 2012;Saetran and Bartl, 2015).The goal of these three workshops is to serve as an ideal test case for CFD tools by providing detailed measurements of the thrust coefficient and the wake properties behind the rotor both in terms of mean flow and turbulence kinetic energy within a controlled wind tunnel environment.Note that this work is an extension of the preceding work of Simisiroglou et al. (2016).
The paper unfolds as follows: Sect. 2 presents the experimental set-up of the workshops, in which the three test cases are outlined.This is followed by a description of the numerical method and of the computational settings used to perform the simulations.The results from the numerical simulations are introduced and discussed in Sect. 3. Lastly,in Sect. 4 the main conclusions of this study are presented.

Experimental set-up
The experiments are performed in the large closed-return wind tunnel facility at the Norwegian University of Science and Technology (NTNU).The test section for all three cases has the width (W ) and length (L) dimensions of W × L = 2.710 m × 11.150 m; Fig. 1.To maintain a zero pressure gradient and maintain a constant velocity along the streamwise direction, the height H of the wind tunnel increases from 1.801 m at the inlet to 1.851 m at the outlet.Velocity measurements are performed using both hot-wire anemometry (HWA) and laser Doppler anemometry for verification purposes.The tip Reynolds number for these three cases is approximately Re c,tip ≈ 10 5 for the upstream wind turbine.This tip Reynolds number is based on the velocity of the tip and the chord length at tip.For full-scale experiments a typical tip Reynolds number is on the order of 10 6 .The air density ρ is equal to 1.2 kg m −3 .In all cases the results are only considered where the wind turbines operate at their design condition, i.e. tip speed ratio of 6 (TSR = 6).
In case A, the three-bladed wind turbine is positioned at a distance of 3.660 m from the inlet.The model wind turbine has a tower that consists of four cylinders of different radii.The hub height is H hub = 0.817 m and the rotor radius is R = 0.447 m.The rotor blades are designed to produce a constant pressure drop across the rotor, which resembles a uniformly distributed thrust, when operating at their design condition (Krogstad and Eriksen, 2013).The airfoil used is the NREL S826 and to increase the Reynolds number of the blades, a chord length of approximately 3 times longer than normal was used.The blades have a circular shape close to the nacelle primarily to allow them to be attached to the hub.The transition from the airfoil section of the blade to the circular section is abrupt.An asynchronous generator of 0.37 kW is located under the tunnel floor and is connected to the wind turbine rotor by a belt located behind the tower.The total blockage effect, defined herein as the fraction of the total tower-and rotor-swept area to the wind tunnel cross section, is approximately 13 %.As a result the flow will be impacted by the walls and this interference will lead to artificial speed-up effects.The stream-wise inlet velocity is U ref,A = 10 m s −1 and the stream-wise turbulence intensity at the turbine position is I u,A = 0.3 %.A thin boundary exists near the wall of the wind tunnel.This boundary layer has been measured in an empty tunnel using pitot tubes for four distances downstream of the wind tunnel inlet, i.e. 1. 80, 4.50, 6.30 and 8.10 m.Further information on the details of the experimental investigations are reported by Krogstad and Adaramola (2012) and Krogstad and Lund (2012).
For case B the stream-wise inlet velocity and turbulence intensity are similar to case A. Here, two in-line wind turbines horizontally centred in the wind tunnel are investigated, where the downwind wind turbine is the same one as used in case A. The upstream wind turbine hereafter is always referred to as T 1 while the downstream wind turbine is referred to as T 2 .Both wind turbines rotate in a anticlockwise direction as seen from the inlet and are three bladed with the same blade geometry and airfoil, i.e. the NREL S826 airfoil.As the nacelle diameter of T 1 is somewhat larger than T 2 , the turbines have slightly different rotor diameters.The rotor radius of the upstream wind turbine, T 1 , is R T 1 = 0.472 m and the radius of the downstream wind turbine, T 2 , is R T 2 = R = 0.447 m.For T 1 as opposed to T 2 , the belt connected to the 0.37 kW asynchronous generator is located within the tower.To calculate the thrust force of the turbines, they are mounted on a six-component force bal-  ance.Further information on the experimental investigations is reported by Pierella et al. (2012Pierella et al. ( , 2014)).Case C is divided into two sub-cases, C1 and C2.For both sub-cases the same wind turbines as in case B are used with a hub height of 0.827 m instead of 0.817 m.The distinction between sub-cases is made because the wind turbines are exposed to different inflow conditions in terms of the wind velocity profile and turbulence intensity as seen in Table 1.For both sub-cases the upstream wind turbine is positioned at 4R from the inlet.Looking at each case individually, case C1 has a uniform inflow velocity of U ref,C = 11.5 m s −1 measured at the inlet and a turbulence intensity of 10 % measured at the first wind turbine position.The turbulence in the wind tunnel is created by a bi-planar grid built from wooden bars installed at the inlet.To estimate the effect of unintended stream-wise velocity gradients, horizontal stream-wise velocity values at four positions downstream of the inlet are measured in an empty domain.For case C1 the thrust values along with the velocity and turbulence kinetic energy at three distances downstream of T 1 are measured.Regarding case C2, a sheared inflow is considered with a turbulence intensity of 10.1 % at the hub height position of T 1 .The inlet velocity as a function of height z is expressed by the power law used for atmospheric flows, which is given as where a is the shear exponent equal to 0.11, the reference height is z ref = 0.827 m and the reference velocity is For C2, similar to case C1, empty-domain measurements are conducted at the same four positions from the inlet for the horizontal stream-wise velocities and turbulence intensity.The position of the second wind turbine is fixed to 10.36R and wake measurements for stream-wise velocity and turbulence kinetic energy are taken at a downstream distance of T 1 equal to 5.54R.

Numerical method
The simulations are performed with the commercial CFD code PHOENICS in which the RANS equations are solved using four different turbulence models.The turbulence models are (1) the standard kε (Launder and Spalding, 1974), (2) the re-normalization www.wind-energ-sci.net/2/587/2017/Wind Energ.Sci., 2, 587-601, 2017 group (RNG) k-ε (Yakhot and Smith, 1992), (3) the KL kε (Kato, 1993) and (4) Wilcox's k-ω turbulence model (Wilcox, 1988).The flow variables are stored in a uniform fully structured staggered grid and the Cartesian coordinate system is used.The SIMPLEST algorithm (Spalding, 1980) is used to solve the RANS equations, and the hybrid differencing scheme (Spalding, 1972) is used to discretize the convective terms.The diffusion terms are discretized using the central differencing scheme.In the computations, the wind tunnel conditions are replicated accordingly for each case and a zero static pressure is applied at the outlet plane.
The lateral, top and bottom faces of the domain are set to be impermeable and a wall function method according to Launder and Spalding (1974) is employed to introduce the effects of the wind tunnel walls into the numerical simulation.This particular method is preferred for its advantages in terms of low computational requirements and storage needs.The wall function method for a flow in local equilibrium obeys the relations where U p is the absolute value of the velocity parallel to the wall at the first grid node, U * is the friction velocity calculated as U * = √ τ w /ρ, κ is the von Kármán constant equal to 0.41, E is a roughness parameter dependent on the wall roughness taken equal to 8.6 for smooth walls, Y + is a dimensionless near-wall quantity for length determined as ν , ν is the turbulence viscosity, Y is the distance of the first grid node to the wall and C µ is a dimensionless constant equal to 0.09 in the standard k-ε turbulence model.If the wall is considered to be rough (not smooth) then the roughness parameter E is a function of the Reynolds roughness number defined as Re r = U * h r ν , where h r is the sand grain roughness height.The relation between the roughness parameter E and Reynolds roughness number Re r follows the empirical laws proposed by Jayatilleke (1966): where b = 29.7, a = (1+2x 3 −3x 2 ) and x = 0.02248•(100− Re r )/Re 0.564 r .
For the simulations no tower or hub effects are considered.The presence of the rotor is modelled using an ACD method based on the 1-D momentum theory.The thrust force F i of each individual cell of the disc is calculated according to where U 1,i is the velocity of the flow at the individual cell numbered i of the disc, α i is the axial induction factor calculated for each individual cell of the disc, A i the surface area of the cell facing the undisturbed wind flow direction and C T (U 1,i ) is a modified thrust coefficient curve dependent on the velocity at the disc.The modified thrust coefficient curve is created in a preprocessing step by replacing the undisturbed wind velocity values of the thrust coefficient curve with the wind velocity values at the disc U 1 .To do this Eq.( 9) is used, where C T is the thrust coefficient for the respective undisturbed wind velocity U ∞ .
The total thrust force applied to the flow is calculated by summing the individual thrust forces according to F tot = i F i over the disc area.This total thrust force may then be distributed in different ways over the disc.In this work, apart from using Eq. ( 8) as it is to prescribe the forces in each individual cell, referred to as the undistributed thrust, four different thrust distributions are tested: a uniform, a polynomial, a triangular and a trapezoidal distribution.Their equations are presented in Table 2. Figure 2 presents a normalized plot of all four distributions along a diameter of the disc.The uniform distribution is chosen to match the thrust distribution of the actual rotor of this case.Full-scale wind turbines however have a zero thrust value at the hub and at the tip of the blades.The polynomial distribution, which is a fourth order polynomial, is intended to respect this by having a zero thrust at the hub and at the tip of the disc.The triangular distribution is designed to have a zero thrust at the hub and to linearly increase the thrust force along the radius, up to the tip of the disc.Lastly, the trapezoidal distribution is set up to resemble the thrust distribution produced using the ACL method presented in Sarmast et al. (2012).While it is possible to determine a thrust distribution given the rotor geometry and airfoil data through a blade element momentum theory, it is somewhat impractical for industrial applications.Airfoil data of commercial wind turbines are generally not available to the typical industrial user.The purpose of testing different thrust distributions with this ACD method is that these will probably produce different wake properties, e.g. with respect to the velocity deficit and turbulence kinetic energy of the wake.Two questions thus arise, i.e. which thrust distribution within this ACD method better captures the wake produced by a wind turbine and up to which distance does the thrust distribution have an effect on the wake?Here it should be noted that the primary goal is not to isolate the influence of the thrust distribution on the wake flow, as the total thrust For the trapezoidal distribution there is no force applied in the region of 0 ≤ r < 0.2R.over the disc will intrinsically vary depending on the thrust distribution used within the method.Here the goal is to investigate the effect the ACD method with different thrust distributions has on the wake flow.
For the first part of the simulations (case A) the numerical domain was defined according to the wind tunnel geometry as reported in Krogstad et al. (2011).Initially, emptydomain simulations were conducted to assess the extent of unintended stream-wise gradients for the mean velocity and turbulence parameters.For this purpose horizontal profiles of U , k and ε are extracted at the inlet, turbine location and x/R = 10 downstream of the turbine position.As the roughness height of the wind tunnel walls is not known a priori, a comparison between the experimental boundary layer profile and the simulated boundary layer profile for different roughness height values is conducted in a trial and error fashion until the appropriate value for the roughness height is found.When considering the wind turbine in the simulation, the computed results are compared against the HWA measurements for the normalized axial velocity U/U ref,A and normalized turbulence kinetic energy k/U 2 ref,A at the three downstream positions mentioned in Table 1 along the horizontal line through the centre of the wake in the crosswise direction.
For case B the domain geometry and the positioning of the wind turbines are in accordance with the invitation sent out by Pierella et al. (2012).The equilibrium wall function method for smooth walls, that is E = 8.6, is used to introduce the effects of the wind tunnel walls into the numerical simulation; this applies as well for case C. The computed results when the ACDs are considered are compared to HWA measurements for the normalized axial velocity U/U ref,B and normalized variance of the axial velocity component u 2 /U 2 ref,B at the three downstream positions shown in Table 1 along the horizontal line through the centre of the wake in the crosswise direction.Further, the thrust coefficients of the two wind turbines are compared with the experimental results, where A is the rotor cross section of each individual wind turbine.Here it should be noted that even though the thrust coefficient curve is an input to the simulation, the thrust coefficient value applied against the flow depends upon the velocity at the disc, which changes as the simulation progresses.Lastly, for case C the domain geometry and the positioning of the wind turbines are in accordance with the invitation sent out by Saetran and Bartl (2015).Prior to the simulations with the ACD, empty-domain simulations are performed for the sub-cases (C1, C2).This is to match the inlet wind profile and turbulence intensity with the experimental measurements at four downstream positions from the inlet, that is x/R = 4.00, 9.54, 14.36 and 22.00.When considering the wind turbines in the numerical simulation via the ACD method, the computed results are compared with the HWA measurements for the normalized axial velocity U/U .

Grid convergence analysis
A grid independence study is carried out according to the recommended procedure of Roy (2003) for mixed-order schemes.For these simulations a uniform grid is used based on the cells per rotor diameter.dence study.Even though the grid independence study is performed solely for case A it is considered to apply for the other cases as well.
According to Roy (2003) the series that represents the discrete solution for each grid is given by where f k is the discrete value solution of grid k, g i is the ith order error term coefficient and h k is a measure of the grid spacing.The three unknowns (f exact , g 1 and g 2 ) may be found by expanding Eq. ( 10) for three consecutive grids and by solving the resulting three-set equation.
The spatial discretization error is calculated according to the following formula: where f k is the discrete value solution of grid k and f exact is an approximation to the exact solution f exact , which is found by disregarding the higher-order terms of Eqs. ( 11) to (13).
The normalized magnitudes of the first-and second-order error terms and the magnitude of their sum (mixed order) is given by where g 1 and g 2 are approximations to g 1 and g 2 , which are found as mentioned previously by solving and disregarding the higher-order terms of Eqs. ( 11) to (13).

Results and discussion
A summary of the cases for which results will be presented is shown in Table 4.For cases A and B, empty-domain re-sults for the axial velocity extracted at cross-sectional horizontal profiles at the inlet, turbine location and at a position x/R = 10 downstream of the turbine location show an approximately 2.7 % increase in the axial velocity.The turbulence parameters k and ε, however, decrease steadily from the inlet to x/R = 10, which is due to the lack of a turbulence generating mechanism along the domain, e.g.shear.The decrease in k and ε along the empty domain is 5 orders of magnitude lower than their average value when an ACD model is present in the computations.Figure 3 presents the spatial discretization error results obtained for the normalized axial velocity profiles at three distances downstream of the wind turbine position for case A. The error is estimated to be less than 2.4 % for finest grid (grid 1).Therefore, for the purpose of this investigation, a uniform grid resolution of 40 cells per rotor diameter is found suitable for all cases.Also shown in Fig. 3 are the normalized magnitudes of the first-and second-order error terms and of their sum, which is given by Eq. ( 15).
Figure 4 illustrates the computed stream-wise velocity results for two different wall function values, against velocity measurements conducted in the wind tunnel with pitot tubes.A quite good match between experimental measurements and the computed results exists when considering a smooth wall, i.e.E = 8.6.Therefore, the equilibrium wall function for smooth walls is used in the simulations, as it is found to have the best agreement with the measurements.By setting a sand grain roughness height other than that for a smooth wall causes the discrepancy between the simulated boundary layer development and the measurements to increase for the axial velocity profile.
Figure 5 illustrates the normalized axial velocity contours for cases A and B. The polynomial thrust distribution is used along with the k-ε turbulence closure model.It is clearly seen that the method reproduces what is expected, that is, by positioning a second turbine in the wake of the first, the axial velocity of the flow is further reduced.This reduction is due to the further energy extraction of the second wind turbine from the mean flow.The dashed lines in Fig. 5 indicate the positions at which flow values are extracted and compared with the HWA measurements.
For case A, the computed results are validated against HWA measurements for the normalized axial velocity and normalized turbulence kinetic energy; see Fig. 6.These results are computed using different turbulence models and the uniform thrust distribution.To investigate the influence of the thrust distribution on the wake development, simulations using the k-ε turbulence model with different thrust distributions were conducted; results are shown in Fig. 7.
In Fig. 6a it is observed that the k-ε and the KL k-ε turbulence models produce results similar to the measurements, with the KL k-ε model being less diffusive than the k-ε model in the crosswise direction.Apart from the undistributed thrust, all thrust distributions used in this study assume axisymmetry; therefore, the simulated profiles are sym-Table 4. Overview of the cases for which results will be presented.

Case Grid and wall
Thrust distribution used for Turbulence model used for Thrust coefficient Empty-domain function study different turbulence models different thrust distributions comparison study metrical to the rotor centre.However, this is not the case with the measurements, which exhibit asymmetric profiles as seen in Fig. 6a.According to Krogstad and Eriksen (2013) this asymmetry may be produced by the slowly rotating tower wake as seen at the downstream position of x/R = 10, for example.As the ACD is non-rotating in the simulations, it is expected that the predictions will not capture this asymmetry in the wake.Further, as the effects of the nacelle and tower are not considered, it is also anticipated to find small deviations of the predictions from the measurements in the near wake.On average the blockage effect is captured by the simulations as seen in Fig. 6a.This effect is apparent outside the wake region (|y/R| > 1.5) where the simulated and measured normalized axial velocity values are higher than one.Considering the normalized turbulence kinetic profiles in Fig. 6b, the shape of the profiles is not successfully predicted by any of the turbulence models.The k-ω turbulence model tends to over-predict the turbulence kinetic energy production in this environment with low background turbulence.As a result, the simulated wake recovery in comparison to the ve-locity measurements is too high.The discrepancy between the measured profile shape of the turbulence kinetic energy from that predicted at the downstream position of x/R = 2 for the wake region of |y/R| < 0.5 is mainly due to the presence of the nacelle and abrupt change of the blade shape from the airfoil profile to a cylinder near the nacelle.
When keeping the turbulence model constant and changing the thrust distribution, it is observed in Fig. 7 that the effect of the thrust distribution is pronounced in the nearwake region and diminishes further downstream.Porté-Agel et al. (2011) also observed that the effect of representing the forces of a wind turbine differently, such as by a rotating or non-rotating ACD or ACL was more pronounced in the near-wake region, rather than in the far-wake region.Regarding the turbulence kinetic energy, all thrust distributions seem to capture the position of the tip vortex apart from the polynomial.The increased turbulence production due to the breakdown of the tip vortex at the x/R = 2 position is not captured by any combination of thrust distribution and turbulence model.This is also observed by other researchers  such as Réthoré et al. (2014), who concluded that the ACD method lacks the ability to simulate the turbulence structures present in the near-wake region.Sumner et al. (2013) results also show that in an environment with low background turbulence intensity there seemed to be a perceptible dependency of the wake development on the turbulence closure used, in terms of velocity deficit and turbulence kinetic energy.
For case B when using the undistributed thrust, the normalized axial velocity and stream-wise variance of the velocity at three positions downstream of the second turbine are shown in Fig. 8 for different turbulence models.In Fig. 9, the effect of the thrust distribution on the wake downstream of the second wind turbine is investigated by varying the thrust distribution while keeping the same turbulence model.Results of the thrust coefficient values for the upstream and   5 for when the k-ε turbulence model is used.
For case B in which the wakes of both wind turbines interact, similar to the results of case A, the k-ε and the KL kε turbulence models produce axial velocity results in agreement with the measurements.The KL k-ε seems however to underestimate the normalized stream-wise variance of the velocity.The k-ω turbulence model here again over-predicts the wake recovery and overestimates the normalized stream-wise variance of the velocity.Conversely, the RNG k-ε underpredicts the wake recovery and underestimates the normalized stream-wise variance of the velocity.The blockage effect is on average captured for all turbulence models, apart from when the k-ω turbulence model is used, and thrust distributions are as seen in Figs. 8 and 9.The wake expansion is accurately predicted when using the k-ε turbulence model.Further, when keeping the turbulence model constant and changing the thrust distribution (Fig. 9), it is observed that the effect of the thrust distribution is less pronounced in the case with two in-line wind turbines than for the case with a single wind turbine due to the higher turbulence diffusion.The wake predicted when using the undistributed or uniform thrust distribution seems to be in closer agreement with the measurements.This is possibly due to the fact that these distributions produce a fairly constant pressure drop over the disc, as do the blades when operating at their design condition.The thrust coefficient values summarized in Table 5 for the different thrust distributions and the k-ε turbulence model agree quite well with the measured data.There is, on average, a 5 % difference between the measured thrust and the results for the upstream wind turbine and less than a 10 % difference for all cases concerning the downstream wind tur-     5 for cases C1 and C2; the k-ε turbulence model is used here.
The empty-domain simulations presented in Fig. 10 give a reasonably good agreement between simulations and measurements for both the uniform and the sheared inflow condition.Concerning the sub-case C1, the simulated streamwise velocity is in quite good agreement with the measurements for all turbulence models and thrust distributions; see Figs. 11 and 12.In this case of an environment with high background turbulence, the k-ω turbulence model shows remarkably better agreement with the measurements when compared to the previous cases.The shape and level of the normalized turbulence kinetic energy is now captured by all turbulence models, though small differences are ob-served between the simulated wakes when different turbulence models are used.This finding is in agreement with results from the studies performed by Laan et al. (2015b) and Sumner et al. (2013), in which small differences between the simulated wakes are also found when different RANS turbulence models are used in an environment with high background turbulence intensity.It should be recalled that emptydomain simulations were performed to match the background turbulence intensity with the experimental measurements when using different turbulence models for all cases prior to the simulations with the ACD.It appears that for the cases with high turbulence intensity this procedure has a significant effect on the computational results compared to the cases with low turbulence intensity.As the purpose of using different turbulence models in this study is to investigate the effect of the turbulence model with its defined constants on the wake development, it is crucial to set the background turbulence intensity throughout the domain in accordance with the experimental set-up when using different turbulence models.This is achieved here by varying the inlet turbulence parameters (k, ε or ω).In this way the background turbulence intensity is similar when using different turbulence models and the effect of the turbulence model on the wake development may be clearly accounted for.When the cases with higher background turbulence intensity are considered, it is observed that the effect of the thrust distribution is less apparent further downstream compared to the lower turbulence intensity of cases A and B. Similarly, for sub-case C2 (see Figs. 13 and 14) the axial velocity in the wake is predicted quite well for all turbulence models and thrust distributions.Because the measurement position is in the near-wake region, the effect of the thrust distribution is apparent on the turbulence kinetic energy and axial velocity profile.From Table 5 it is found that thrust coefficient values are estimated on average to have less than a 10 % difference from the measured values for both sub-cases when using the k-ε turbulence model.
Lastly, in terms of computational time or CPU hours, herein defined as the number of CPUs × wall clock time needed to perform the simulation, results are shown in Table 6.These results present the CPU hours needed to perform the simulations using this method and a LES with the ACL method described in Sørensen et al. (2015) for case A. It is found that the RANS-ACD method is significantly faster in simulating this one wind turbine case compared to the LES-ACL method.Although the LES-ACL method provides high-fidelity results comparable to the measurements, the computational requirements of this method, as for now, are still too demanding to make it usable for wake modelling in industrial applications.

Conclusion
The main conclusions of this study are summarized as follows: (i) the present results, considering the simplicity and low computational needs of the method, generally show satisfactory agreement between the simulations and the measurements used for both the set-up with one wind turbine (case A) and the set-up with two in-line wind turbines (cases B and C).(ii) The effect of using different thrust distributions on the profiles is generally present in the near wake and fairly absent in the far wake.Moreover, the impact on the near wake is more pronounced for the set-up with a single wind turbine than in the set-up with two wind turbines.
(iii) The uniform and undistributed thrust distributions generally outperformed the other distributions in terms of the estimated wake.Please note however that the uniformly distributed thrust might not be the best suited when considering near-wake effects if a full-size wind turbine that typically has a non-uniform thrust distribution is modelled.(iv) Changing the turbulence model has a noticeable impact on the wake development in the cases with low background turbulence intensity.When using the k-ε and KL k-ε turbulence models, the velocity results are in agreement with the measurements, but this is generally not the case with the k-ω turbulence closure models.(v) When considering the cases with high background turbulence intensity, small differences in www.wind-energ-sci.net/2/587/2017/Wind Energ.Sci., 2, 587-601, 2017 the wake development are found by changing the turbulence model.The results, however, are quite sensitive to the inlet conditions of the turbulence parameters used for the simulations.Therefore, depending on the turbulence model, the turbulence parameters at the inlet should be carefully considered as to represent the background turbulence experimental conditions.(vi) The wake in terms of the velocity and turbulence profiles was captured more accurately in the cases with high background turbulence than in the cases with low background turbulence.This method has shown to give reliable results for a number of different wind flow conditions and separation distances with respect to the case with a single in-line wind turbine and the case with two in-line wind turbines.However, it has not been validated yet for wind turbines operating in a situation in which only a part of the rotor is in the wake of the upstream wind turbine (partial wake situation).Moreover, it has also not been validated against operational data measured within existing wind farms operating in full-scale atmospheric conditions.Therefore, future research will focus on validating the method against data retrieved from operating wind farms.Cases in which wind turbines operate partially in the wake of the upstream turbine will be of special interest as well.

Figure 1 .
Figure 1.Illustration of wind tunnel layout for (a) the set-up with one model wind turbine (case A) and (b) the set-up with two in-line model wind turbines (case B).The three downstream positions x/R = 2, 6, 10 and x/R = 2, 5, 8 are where the measurements are extracted; radius R = 0.447 m.

Figure 2 .
Figure 2. Normalized plot of all four distributions along a diameter of the disc.
ref,C and normalized turbulence kinetic energy k = k/U 2 ref,C , where U ref,C = 11.5 m s −1 .The thrust coefficients of the two turbines are calculated as C T = 2F tot ρU 2 ref,C A

Figure 3 .
Figure 3. Normalized axial velocity and spatial discretization error computed behind a single model wind turbine (case A) for three grids at (a, d) x/R = 2, (b, e) x/R = 6 and (c, f) x/R = 10 for the k-ε turbulence model and the undistributed distribution.

Figure 4 .
Figure 4. Normalized axial velocity for different wall functions plotted against the experimental data for four downstream distances from the inlet.

Figure 5 .
Figure 5. Normalized axial velocity contours for the k-ε turbulence model using the polynomial thrust distribution.(a) One model wind turbine, case A, and (b) two in-line model wind turbines, case B.

Figure 6 .
Figure 6.(a) Normalized axial velocity and (b) normalized turbulence kinetic energy computed behind a single model wind turbine (case A) at x/R = 2, 6 and 10 for different turbulence models using the uniform distribution.

Figure 7 .
Figure 7. (a) Normalized axial velocity and (b) normalized turbulence kinetic energy computed behind a single model wind turbine (case A) at x/R = 2, 6 and 10 for different thrust distributions and the k-ε turbulence model.

Figure 10 .
Figure 10.Empty-domain stream-wise velocity and turbulence intensity results for case C1 (a, b) and case C2 (c, d) for four vertical profiles downstream of the inlet when using the k-ε turbulence model.

Figure 11 .
Figure 11.(a) Normalized axial velocity and (b) normalized turbulence kinetic energy computed for the two in-line model wind turbines (case C1) at x/R = 5.54, 10.36 and 17.00 downstream of the first wind turbine for a separation distance of 18.00R.Different turbulence models are used along with the undistributed thrust.

Figure 12 .
Figure 12.(a) Normalized axial velocity and (b) normalized turbulence kinetic energy computed for the two in-line model wind turbines (case C1) at x/R = 5.54, 10.36 and 17.00 downstream of the first wind turbine.Different thrust distributions and the k-ε turbulence model are used.

Figure 13 .
Figure 13.(a) Normalized axial velocity and (b) normalized turbulence kinetic energy computed for the two in-line model wind turbines (case C2) at x/R = 5.54 downstream of the first wind turbine.Different turbulence models are used along with the undistributed thrust.

Figure 14 .
Figure 14.(a) Normalized axial velocity and (b) normalized turbulence kinetic energy computed for the two in-line model wind turbines (case C2) at x/R = 5.54 downstream of the first wind turbine.Different thrust distributions and the k-ε turbulence model are used.

Table 1 .
Overview of cases.

Table 2 .
Thrust distributions over the disc, where r is the distance from the centre of the disc and R is the radius of the disc.Note that f on the left-hand side of the equations has dimensions of force per unit area.

Table 3 .
Grid levels and size.

Table 6 .
Averaged computational effort in CPU hours to perform the simulations of case A.