Experimental and numerical investigations of actuator disks for wind turbines

The expenses of a wind turbine modeling in particular a wind farm modeling either numerically or experimentally cause to take the advantages of actuator disks (AD). Porous disks (PD) are used to simulate actuator disks in experiments especially for wake studies. A rotor of a wind turbine replaced by a PD can be modeled efficiently with less cost and process time. For a proper PD selection, some semi‐empirical equations are suggested for a range of solidity of 0.2‐0.6. These equations came from a number of tests. The PDs with different dimensions and solidity values were considered for the tests and models. The results of 2D and 3D numerical simulations and the experimental results show that at the solidity of 0.5, the coefficient of performance of a PD reaches approximately 16/27 which is equal to the Betz's limit; the relation between the coefficient of power and the solidity is fitted by a parabolic equation. The average error of the proposed equations for the power coefficient is reported 1.59%. Finally, based on the current results, some semi‐empirical equations are suggested to help the initial PD selection, and then, the cost of further studies may decrease.


| INTRODUCTION
Wind turbines are going to produce more than 800 GW power worldwide by the end of 2021 1 showing the importance of wind turbine analysis either numerically or experimentally. It is not easy to model the wake of a wind turbine using a real turbine especially for wind farms. For saving the process time and cost of the studies, the rotors of wind turbines can be replaced by porous disks.
The actuator disk theory is an essential tool for wake analysis. Wake studies can be variously found in the past researches such as Naderi and Torabi 2 and Nedjari et al. 3 They used AD theory and numerical simulations of the wake study. Castellani and Vignaroli 4 showed the possibility of using ADs for large offshore wind turbines. Experimental approaches for AD theory can be found in the study of Stevens et al. 5 They compared AD and actuator line theories with wind tunnel experimental methods. From their work, the actuator line theory could capture the near wake structures up to 3 diameters downstream of a single turbine more accurate than the others. Airflow measurements behind wind turbines were studied vastly in the literature. Magnusson and Smedman 6 analyzed a four-turbine wind farm and showed that the coefficient of thrust is a useful criterion to study wake characteristics in comparison with wind speeds.
One of the main parameters of wind turbines is solidity. It can be replaced by porosity in porous disks. Porosity is the ratio of void spaces to the total area, and solidity is the ratio of covered spaces to the total area. In 2003, Duquette et al 7 showed that at a constant solidity, increasing the blade number increases the power coefficient of SG6037 airfoil blade profile for a range of Reynold number of 10 5 to 5 × 10 5 . There are some researches that studied porous-disk performances. Aubrun et al 8 studied experimentally the coefficients of power and thrust for a range of the porosity from 55% to 70%. They used wire meshes as disks and tested them under a constant wind speed. Theunissen et al 9 modeled an offshore wind farm by circular porous disks experimentally. The shape of porous disks was different from the literature; they made some different disk layouts in order to attain the same downstream wakes that they expected from real turbines. They modeled an offshore wind farm by their nonrotating disks. Sforza et al 10 investigated the wind turbine wakes. They showed that a modeled wind turbine by a porous disk can show the same wake characteristics as a real wind turbine if the coefficient of the thrust is the same for both the real and the modeled turbines. Dighe et al 11 studied diffuser augmented wind turbines (DAWT) by the AD approach; they installed a porous disk in the middle of a diffuser. The porosity of 40% was selected since this amount of porosity gave the same thrust coefficients as the real rotor. They continued their studies and showed the effects of Gurney flaps on the diffuser augmented wind turbines by using the porosity of 70% and compared with the results of a case without any Gurney flap. 12 Ranjbar et al 13 modeled the AD in a DAWT numerically. They showed that the DAWT increased the coefficient of performance over the Betz limit. Auburn et al 14 compared a nonrotating turbine model with a rotating one with the same coefficient of thrust. The scale of the model was 1:300. They measured the velocity deficit of the wake for relatively low turbulence conditions. A numerical study was done by Crasto and Gravdahl 15 in order to analyze the pressure drop of a specific area of the wake. The pressure drop was highly related to the axial induction factor and thrust coefficient. Lignarolo et al 16 tested an AD model and porous disk using the particle image velocimetry (PIV) technique. Their results including the thrust, enthalpy, and power were very close to the real turbine. Aloui et al 17 used the PIV method with a proper orthogonal decomposition (POD) method for a porous disk with a porosity of 55%. They investigated the wake characteristics of the disk with the aid of porous disks; researchers have had progress for wind farm studies. Charmanski et al 18 studied the boundary layer of a wind turbine array. Each porous disk has a porosity of 51%. This porosity value was chosen based on the extracted drag coefficient from a real rotating wind turbine. Hamedi et al 19 introduced a semi-analytical model for velocity profiles at the wake of a wind turbine using the blade element momentum (BEM) method. The model of Hamedi et al was used by Naderi et al for a wind farm analysis. 20 Naderi et al 20 modeled the wake of a horizontal axis wind turbine located in Horns Rev offshore wind farm using an improved actuator disk model. They analyzed the wind farm in different wind directions. They claimed that their proposed model performs similar to the large eddy simulation (LES) method. Camp and Cal 21 studied numerically and experimentally by particle image velocimetry technique to investigate the effects of the mean kinetic energy transport. They recommended to use of an actuator disk to model the wind turbine rotor for their computational simulations. The aim of the paper was based on using actuator disk to discuss the Reynolds stresses and the turbulence intensity. Pirrung and Van Der Laan 22 proposed a modification for the tangential loads by applying a commonly used tip loss factor for AD simulations. The proposed normal and tangential forces lead to an agreement between the AD and the reference BEM simulations, and their proposed modification is easy to implement and comes at a negligible computational cost.
Porous disks are used instead of actuator disks in experiments especially for wake studies. When a flow passes through a porous disk, the kinetic energy of the flow is transferred to the turbulent kinetic energy; the axial velocity of the flow is reduced near the porous disk, and the pressure of the flow is dropped similar to the actuator disk theory. The velocity reduction and the pressure drop are the main concepts for modeling the wind turbine based on the actuator disk theory. The porous selection is based on the thrust coefficient of a wind turbine. For the most of previous studies, there is no discussion in details about the design of the porous disks. Some studies have used a trial and error method to find suitable porous disk to mimic the desire wind turbine. The focus of the current study is to propose some methods for selecting a 2D porous disk with a proper solidity (instead of porosity) and geometry. Solidity is a nondimensional term which will be introduced later and different from the conventional solidity. The study has been conducted both experimentally and numerically. A comparison between the experimental and numerical results shows the reliability of the results. Also, the results are confirmed with the AD theory.

RANJBAR et Al.
where h and G are dimensions of the PD, and N is the number of Blocks (Figure 1). This name (solidity) is selected because the closest known concept in the aerodynamics of wind turbines to this nondimensional parameter is solidity. The velocity deficit which describes a velocity reduction at the disk, defines as axial induction factor (a): Since for a PD, there is no more induction factor, then for the rest, the term "induction factor" will be used instead of the "axial induction factor" for simplicity. The term induction factor for PDs has been used by others. 16,17,19,20 In this equation, ΔV is the difference between the tunnel velocity and the local velocity in the proximity of a PD (ΔV = V tunnel − V AD ). Since t is the thickness of the PD (Figure 1), the Aspect Ratio (AR) can be introduced as: The last dimensionless parameter is Reynolds number (Re) which is expressed as: For the circular screens, the parameter G should be replaced by the diameter. (2) induction factor = a = ΔV V tunnel .
(4) Re = VG ; An actuator disk (AD) is a theoretical concept of an infinitesimal thickness that extracts energy from the flow (or injects energy into the flow), while a porous disk (PD) is a physical screen which transfers kinetic energy of the flow into the turbulent kinetic energy. It is important to discriminate the AD from the PD. In this study, a 2D square PD and a circular PD are used to simulate actuator disks. Both pressure drop and velocity deficit, which are two main phenomena according to the actuator disk theory, are visible when the flow passes the porous disk.

F I G U R E 1 Parameters of a PD
In the current study, the square and circle shape for the screens are used. Although the shape of the disk may affect the size of the wake, in Sections 6.5 and 6.6, it will be proven that the shape of the screens (circle or square) do not change the results extracted from this study. Actuator disks were made by Plexiglas screens. The Plexiglas screens were cut with a laser (±0.1mm) with Solidities (S) from 0.2 to 0.6. PDs were produced in long rectangular shapes (Blocks). The screens are classified in two types described in the following.
Actuator disk I: h = t = 3 mm that means the thickness and the width of the square Plexiglas screens are 3mm (Figure 3, left).
Actuator disk II: h = t = 5 mm with the thickness and the width of 5mm for each screen (Figure 3, right).
All tested PDs with their dimensionless parameters are listed in Table 1.

| NUMERICAL SETUP
The governing equations, continuity and incompressible Navier-Stokes equations, in tensor notations are where v, ρ, P, and µ denote velocity, fluid density, mean static pressure, and dynamic viscosity, respectively. In these equations, term G i (body forces) is neglected.

| Geometry and mesh generation
For 2D simulation of PDs, small squares are used to cause pressure drop and velocity deficit across PDs. Figure 4 shows the cross-sectional view (x-y plane) of the screen (2D).
For mesh independency, five different meshes from 15 × 10 3 to 10 5 elements were tried. Based on the mean velocity in the  proximity of the PD, the selected grid includes 6.8 × 10 4 elements. For the case of S = 0.3 (PD I), the size of the first layer is 0.01 mm with 1.2 growth rate and the maximum cell size of 35 mm. These settings have been used to generate unstructured grids in all cases with the solidity of 0.2 to 0.6 with y + < 8.

| Solver
For the boundary conditions, in Figure 5, the inlet and the outlet are assumed velocity inlet and pressure outlet, respectively. The upper and bottom boundaries are considered symmetric. The inlet velocity is 8 m/s. The governing Equations 5 and 6 have been solved based on the SIMPLE algorithm with the aid of k − turbulence model. Also, for the validation purpose, k − SST and Transition SST turbulence models were tested too. The second-order upwind spatial discretization is used. The flow is assumed 2D, steady, and incompressible. The convergence criteria were set to 10 -3 for all the residuals.

| VALIDATION
The results are validated by two approaches. As the first step, the CFD results should follow the actuator disk theory which will be discussed in Section 5.1. The numerical results are also compared with the available experimental data from the literature, Section 5.2.

| Actuator disk theory
The pressure drop, velocity deficit, and velocity profiles are compared with those of the Betz approach and momentum theory. In this part, all figures are plotted from numerical results.

| Pressure drop and velocity deficit
According to the Betz theory, the AD blocks the free stream and causes a pressure drop and a velocity deficit. Figure 6A confirms that pressure drops between upstream and downstream of the PD in the center line of the PD. Far from the PD in the downstream, the pressure reaches to the atmosphere  Figure 6B). The pressure coefficient is de- where p,p ∞ , , and U ∞ are the static pressure of the flow field, the static pressure of the free stream, the free-stream fluid density, and the free-stream velocity, respectively. Zero pressure coefficient (in Figure 6B) means the pressure of the evaluated point is equal to the freestream static pressure. Figure 7 shows the velocity deficit in the wake of the PD. The PD generates some vortices ( Figure 7A) which oscillates the velocity of the wake. Similar to the Betz theory, the velocity of the free stream decreases before the PD and after the PD, Figure 7B. The downstream velocity reaches to a value less than the free stream velocity showing the energy extraction (broken lines in Figure 7B). Very close to the boundaries of the PD, the pressure and velocity values are scattered because of the stagnation point. Since the scattered values do not provide any information, those values are deleted in Figures 6 and 7.
As the solidity increases, the pressure drops in the wake, Figure 8. For S greater than 0.4, the low pressure coefficient covers most parts of the wake significantly. These pressure differences result in a high thrust coefficient which should be considered in the selection of S values.
Increasing the solidity magnifies the velocity deficit in the wake, Figure 9. Similar to Figure 8, for S greater than 0.4, the significant velocity drop dominates most part of the wake According to the continuity equation, the downstream (A out ) area should be bigger than the upstream area (A in ). Figure 10 demonstrates the expansion of the streamlines clearly. Since the velocity deficit is significant in the wake (Figure 9), based on the continuity equation, the wake expansion was expected. Figure 11 and Figure 12 illustrate the velocity profiles for two different disks of the same solidity (S = 0.5) in the proximity (Δx = 5 mm before) of the PD. In both cases, the mean velocities are approximately equal with the difference of 0.046 m/s. The average normalized velocity in Figure 11 is 0.63 (5.023 m/s) and in Figure 12   magnitude of the velocity increases because of the edge effect of the two sides of the actuator disk; the zoomed area in Figure 10 shows the edges of the actuator disk. The mentioned figures agree with the actuator disk (momentum) theory approving the quality of the numerical and experimental setup. Also, in the next section, the agreement between experimental results and the numerical ones will show the reliability of both methods. A wake profile comparison in the distance of the 2D of the disk between the experimental data 8 and the present CFD results is shown in Figure 14. The experimental data were collected by a hot wire. 8 Also, k − SST and Transition SST turbulence models are added to Figure 14. The figure shows that the results of the k − turbulence model are closer to the experimental data compared to the curves from the k − SST and Transition SST turbulence models. The same results were reported by Ref. 9 for selecting the turbulent model. According to Figure 14, the results of two approaches agree tolerably.

| RESULTS AND DISCUSSIONS
The determination of the performance/power and the thrust coefficients is based on the theoretical relation using the measured velocities and the Betz relation. First, the induction factor (Equation 2) is calculated according to the measured velocities, and then, it is possible to obtain the performance/power C P = 4a (1 − a) 2 and the thrust coefficients C T = 4a (1 − a) . These equations are extracted based on ideal conditions such as frictionless flows and no rotational velocity component in the wake. It is assumed that the disk works like a drag device for reducing the velocity. It should be mentioned that the induction factor also was calculated by the pitot tube array based on the integration of the velocity along the disk described by Lignarolo et al. 16 The same results from integration method and equation 2 were obtained.

| Power coefficient, C P
The performance coefficients from the experimental results versus Reynolds number (Re) are shown in Figure 15. The

V/Vref X-Position
solidity of 0.5 with the induction factor of 1/3 for both PD I and PD II has the maximum C P close to the Betz limit, 16/27; the C P is almost constant with the variation of velocity (from 4 to 14 m/s or 5 × 10 4 < Re < 1.85 × 10 5 ). Any solidity value except 0.5 results in the C P reduction. After S = 0.5, the more velocity reduction or kinetic energy extraction in the wake does not have positive effects on power coefficient ( Figure 15). For low solidity values, the C P variation is significant in particular at low velocities (V < 7 or Re < 9 × 10 4 ). According to the Betz approach, the Cp is independent of Reynolds number, C P = 4a (1 − a) 2 . According to Figure 15, deviation from S = 0.5 results in Reynold dependency in particular for Re < 10 5 . It is concluded that a PD with the solidity of 0.5 can be used as an ideal wind turbine candidate that extracts 16 27 of the wind energy (Lanchester-Betz limit). This conclusion may be used for further study of ideal wind turbines. 23 6.2 | Thrust coefficient, C T C T is directly proportional to the solidity (up to 0.6 solidity), Figure 16. For all experimental results, an exception is seen in Figure 16 near S = 0.4 where the C T variation measured by the anemometer is significant. The same discrepancy for S = 0.4 using the anemometer will be seen in next section. For S = 0.4, the measurements were done by the pitot tube array; with this measurement, for S = 0.4, the results agree well with the other methods ( Figure 16). From the AD theory, for an ideal rotor, C T should be almost 0.89 which agrees well with the results of S = 0.5 in Figure 16. The same conclusion as Figure 15 shows that this solidity is qualified for ideal rotors. In the literature, C T has been used as an indicator for replacing the rotor of a wind turbine or impeller by a PD. With this assumption without calculating C T for each PD, the solidity may be used for PD selection to replace the rotor. More details about the results and Eq. 13 shown in Figure 16 will be provided in the coming sections.

| Induction factor, a
The induction factor is calculated based on Equation 2. The velocities are measured in the proximity of the disk (5 mm or x/d = 0.025) located in the center and downstream of the disk. The distance was selected 5 mm since the probe caused limitation for closer distances. Moreover, some tests were done with the pitot tube for x/d < 0.025 and the results were almost the same. The previous study of Aubrun et al 14 showed that x/d = 0.025 was very close to the disk and the velocity of x/d = 0.025 presented almost the velocity inside the disk. The works of Dighe et al 11,12 provided some discussion based on the results of the particle image velocimetry (PIV) technique. Figure 7 also provides that for x/d = 0.025 (∆x = 5 mm), the closer distances do not have any effects; because of satisfying the continuity equation, velocity variation is negligible in small lengths.
The velocity reduction (ΔV) presented in Eq. 2 should be a function of the solidity; then, the induction factor is a function of S. Figure 17 shows increasing the solidity increases the induction factor or velocity reduction. For S = 0.5, the results of the experimental and numerical approaches are very close. For PD I, the S = 0.5 from the experiment provides the induction factor (a) of 1/3, which agrees with the Betz limit for maximum power production, C P max = 16 27 = 0.592. PD I with small thickness or low aspect ratio can provide a condition which is close to the Betz assumption. It can be concluded that to reach the Betz limit both experimentally or numerically, the disk should be thin enough and the solidity should be 0.5. From the experimental results of PD I, the induction factor versus the solidity is almost linear except for S = 0.4 when the velocity was measured by the anemometer, but the results from the pitot tube array have a good agreements with the rest. The values of the induction factors are illustrated in Table 2. More discussion about Figure 17 and Equation 12 ( Figure 17) will be provided in the following section.

| Equation extraction
From the theory and the experimental results, extracting equations for relating the solidity to key parameters such as the coefficient of performance (C P ), the coefficient of thrust (C T ), and the induction factor (a) are considered. These equations may make more practical the PD selection for analyzing wind turbines, propellers, and wind farms numerically and experimentally as discussed previously in Sections 6.2 and 6.3. These equations can prevent using the trial and error approaches and make the study more costeffective. For this purpose, the first step is to fit a curve. In Figure 18, the best fitted curve is shown for the C p ; the priority is with the experimental data from PD I case. The trend is parabolic.
The second-order polynomial curve for the C P versus the solidity is No obstacle in the flow (S = 0) with zero extracted power (C P = 0) gives = 0. From the experimental results, the maximum C P of 16 27 is for S = 0.5; then. and (7) C P = S 2 + S + . Substituting the constant parameters in Equation (7) gives Eq. 10 which came from the theory and the experimental results shows the dependency of the C P to the solidity. This equation is validated for the solidity from 0.2 to 0.6. In Equation 10, the solidities greater than 0.69 are not acceptable. Equating C P = 4a (1 − a) 2 from the theory to Equation 10, provides a relation between the solidity (S) and the induction factor (a): From the theory of the AD, a < 1 2 gives Based on the numerical and experimental results, all of these equations are extracted from the 2D AD concept with all assumptions from the AD theory; the errors will be discussed later. According to the experimental results, case PD I has the best agreement with the Betz theory; see Sections 5.2, 6.1, 6.2, and 6.3. Eggleston and Stoddard 24 showed that the AD theory is not valid when 0.5 < a for thrust (C T ) calculation. 0.5 < a is equal to 0.69 < S (Equation 12), which means, for 0.69 < S, the equation of the C T (Equation 13) is not valid. A summary of PD I errors between the results of the experiments and extracted equations (Eq. 10 and Equation 13) is provided in Table 3. Average errors between the results of experiments and equations are 1.59%

T A B L E 3 Summary of relative errors
of PD I and 2.82% for C P and C T , respectively. Equation 10, 12, and 13 can help the PD selection for analyzing either wind turbines or wind farms numerically and experimentally. In this regard, Eq. 13 can be equal to the C T of the real wind turbine; it can prevent using the trial and error approaches for the PD selection and makes the study more cost-effective. The common methodology for PD selection was mainly based on trial and errors; suggested approach in the current study can help to reduce the number of trails and reduces the cost of the study. The suggested PD can be used for analyzing either a wind turbine or wind turbine clusters in a wind farm especially in yaw condition. 25 As a preliminary selection of a PD, the following steps may be used for either experimental or numerical studies to prevent trial and errors: -specify the thrust coefficient (C T ) of the desired rotor, -use Figure 16 or Equation (13) to get the solidity of the PD, -select the PD with the proposed geometry (Sections 2, 3, and 6) with minimum thickness, and -consider Sections 6.5 and 6.6 for 2D versus 3D studies.

| 2D versus 3D cases, CFD approach
For the 2D study based on the numerical approach, it was assumed that the PD is made of some squares shown in Figure 4. For the experimental approach, the squares were extended to generate a square screen, Figure 1. In this regard, the flow field in the middle of the screen is completely 2D. If the squares create a screen which generates a 3D flow field such as the screen shown in Figure 19 (the right bottom figure), the possibility of using the 2D results for 3D studies can be investigated. The screen shown in Figure 19 has the solidity of 0.5. The size of the squares in the mid-span of the screen is the same as the size of PD I. For the numerical study, a 3D structured mesh with the size of 5 million cells was generated, Figure 19. All the settings for the numerical process were the same as the 2D simulation settings. The pressure coefficient contours and the velocity magnitude contours are shown in Figures 20  and 21, respectively. The figures agree well with the AD theory similar to Figures 8 and 9 for S = 0.5. In these figures, just a quarter of the disk is shown for the sake of visibility.
A comparison between 2D and 3D analyses is summarized in Table 4. The deviations of the 2D and 3D power coefficients from the Betz limit are estimated 1.28% and 1.13%, respectively. Since the results of the 3D case agree well with the results of the 2D case, it can be concluded that the disk selection method which is based on the 2D approach can be applied for a 3D study of wind turbines. In Figure 21, the tip effects are visible. Newman 26 suggested that if the outer edge of the disk is changed to slightly conical shape, it is possible to delete the edge effects. Since the main approach of this study is based on the 2D flow field, just the mid-span of the screens is important regardless of the shape of the disk borders. The rectangle or circle shapes of the screens may change the size of the wake which is not considered in this study; for a circle screen, the edge G in Figure 4 which is an important parameter should be equal to the diameter of the circular disk. Figure 22 shows the velocity contours in different positions from the 3D numerical simulation. For the upstream of the PD, one diameter ahead of the disk, the flow field is not disturbed. One diameter after the disk, the wake is significantly affected and the velocity reduction is clearly visible. For five diameters behind the disk, the effects of the disk on the wake are almost disappeared.

| 2D versus 3D cases, experimental approach
For the PD I with the solidity of 0.5, a comparison between the disk with a square border and the disk with the circle border was done experimentally. The geometry of the circular disk is similar to the one modeled numerically in Section 6.5. The diameter of the circular disk is about 20 cm, Figure 23. For the sake of the structure, four radial lines are added in the disks (Figure 23). These radial lines do not change the value of the solidity significantly, and their effects are ignored.
The tests were continued by using a pitot tube array (Figure 23, left). The results from both measurement tools were almost the same. The experimental results are summarized in Table 5. The deviation of power coefficients from the Betz limit is less than 1.3%. It can be concluded that the disk selection method which was based on the 2D approach can be applied for a 3D study of wind turbines.
For some studies, the limitation of modeling either numerically or experimentally may need a disk with a circular boundary or a disk with a square boundary. The shape of the boundary which may cause 2D or 3D effects may cause some advantages for these studies. The comparisons in this section show that the results are almost independent of the boundary. This founding can help to reduce the cost of further studies.

| CONCLUSION
According to the AD theory, experimental and numerical methods were applied for the flow field modeling and simulating in the vicinity of PDs with the solidity of 0.2-0.6. Relations between solidity and some wind turbine parameters including induction factor, power performance, and thrust coefficient were investigated. The 2D and 3D numerical results agreed acceptable with the experimental ones especially for low ARs. From the comparison between numerical results and the experimental results, the k − turbulence models are more accurate than k − SST and Transition SST turbulence models. A low AR disk with narrow stripes could produce more homogeneous and smooth velocity profiles in the downstream of the wake. From the experimental results and the momentum theory, a second-order relation between the power coefficient, thrust coefficient, induction factor, and solidity was extracted: the equations were valid for the range of