Maximal power per device area of a ducted turbine

Abstract. The aerodynamic design of a ducted wind turbine for maximum total power coefficient was studied numerically using the axisymmetric Reynolds-averaged Navier–Stokes equations and an actuator disc model. The total power coefficient characterizes the rotor power per total device area rather than the rotor area. This is a useful metric to compare the performance of a ducted wind turbine with an open rotor and can be an important design objective in certain applications. The design variables included the duct length, the rotor thrust coefficient, the angle of attack of the duct cross section, the rotor gap, and the axial location of the rotor. The results indicated that there exists an upper limit for the total power coefficient of ducted wind turbines. Using an Eppler E423 airfoil as the duct cross section, an optimal total power coefficient of 0.70 was achieved at a duct length of about 15 % of the rotor diameter. The optimal thrust coefficient was approximately 0.9, independent of the duct length and in agreement with the axial momentum analysis. Similarly independent of duct length, the optimal normal rotor gap was found to be approximately the duct boundary layer thickness at the rotor. The optimal axial position of the rotor was near the rear of the duct but moved upstream with increasing duct length, while the optimal angle of attack of the duct cross section decreased.



Introduction
The power output of a wind turbine can be augmented by surrounding it with a duct, typically referred to as a ducted wind turbine (DWT), a diffuser augmented wind turbine, or a shrouded wind turbine. The effect of the duct is to increase the mass flow rate through the rotor. For a given rotor area, significantly more power can be obtained for a DWT compared to an open wind turbine. However, by adding a duct, the total area of the device facing the wind direction is increased. If the power produced per total projected frontal area of the device is calculated for DWTs, often values closer to that of open wind turbines are found (van Bussel, 2007). When nondimensionalized by the kinetic power available in a unit area of freestream, the power per rotor area and power per total area of the DWT are referred to as rotor and total power coefficients and are denoted by C P and C P,total respectively. For an open rotor turbine, these two power coefficients are equal. Achieving values of C P,total greater than the Betz-Joukowsky limit (Okulov and van Kuik, 2012) of 0.593 for a DWT is significant as it means a DWT can capture more power per unit area of the device than an open rotor turbine. Optimization studies of DWTs have shown that DWTs can achieve values of C P,total beyond the Betz-Joukowsky limit (Aranake and Duraisamy, 2017;Bagheri-Sadeghi et al., 2018). C P,total not only is a useful metric to compare DWTs with open rotor wind turbines but could be an important design objective in certain problems like fully integrated DWTs for sustainable buildings (Ishugah et al., 2014;Agha and Chaudhry, 2017) or other applications where a designer seeks to maximize the power output from the limited space allocated to a wind turbine.
Since the experimental demonstration of the power augmentation provided by shrouding wind turbines (Sanuki, 1950), numerous studies on the design and optimization of the DWTs have been carried out (Lilley and Rainbird, 1956;Igra, 1981;Loeffler, 1981;Gilbert and Foreman, 1983;Georgalas et al., 1991;Politis and Koras, 1995;Hansen et al., 2000;Phillips et al., 2002;Ohya and Karasudani, 2010;Aranake and Duraisamy, 2017;Venters et al., 2018; Bagheri-  Sadeghi et al., 2018). Only a few (Aranake and Duraisamy, 2017;Bagheri-Sadeghi et al., 2018;Venters et al., 2018) have used C P,total as a design metric, while most other studies have focused on maximizing C P . The design variables of a DWT include the rotor blade design, its axial location (z rotor ), tip clearance or rotor gap ( n), the angle of attack (α), length (l), and shape of the duct cross section. These design variables can be seen in the schematic shown in Fig. 1 with the rotor replaced by an actuator disc.
Many of the numerical studies use axisymmetric computational fluid dynamics (CFD) models (Loeffler, 1981;Georgalas et al., 1991;Politis and Koras, 1995;Phillips et al., 2002;Aranake and Duraisamy, 2017;Venters et al., 2018;Bagheri-Sadeghi et al., 2018). If the rotor blades are modeled as an actuator disc, the thrust coefficient can be considered a design variable, and different rotor loadings can be represented by different values of the thrust coefficient (Loeffler, 1981;Hansen et al., 2000;Phillips et al., 2002;Venters et al., 2018;Bagheri-Sadeghi et al., 2018). Similarly, some experimental studies replace turbines with screens of different porosity to study DWTs at various rotor loadings (Igra, 1981;Gilbert and Foreman, 1983). Most simplified theoretical models of DWTs indicate that the optimal power output is achieved at a thrust coefficient of 8/9 similar to open wind turbines ( van Bussel, 2007;Jamieson, 2008;Bontempo and Manna, 2020). Bontempo and Manna (2020) reviewed various theoretical models of ducted wind turbines and concluded that they are all equivalent and that their apparent differences are due to different choices of flow parameters used to characterize the effect of the duct (e.g., the exit pressure coefficient, Foreman et al., 1978; or extra back pressure velocity ratio, van Bussel, 2007). They also show that these models are insufficient in predicting the optimal design of ducted turbines as they neglect the dependence of the flow parameters on the thrust coefficient. However, they demonstrated that C P,total greater than the Betz-Joukowsky limit can be achieved at C T = 8/9. Bagheri-Sadeghi et al. (2018) reported an optimal value of thrust coefficient close to this value when optimizing for C P,total using Reynoldsaveraged Navier-Stokes (RANS) simulations with an actuator disc model. However, McLaren-Gow (2020) performed axisymmetric inviscid simulations of DWTs with various duct shapes with an actuator disc and concluded that the value of C T to maximize C P is lower than 8/9.
As most studies focus on maximizing C P , with a few exceptions (Georgalas et al., 1991;Politis and Koras, 1995;Venters et al., 2018;Bagheri-Sadeghi et al., 2018), the axial location of the rotor is usually fixed at the smallest cross section of the duct where the maximum velocity is assumed. In our previous paper (Bagheri-Sadeghi et al., 2018), we included the axial position of the rotor as a design variable and compared the optimal designs for maximum C P and C P,total . We observed that the optimal axial location of a rotor to maximize C P or power is close to the duct throat. However, when optimizing for maximum C P,total , the optimal axial position moves further downstream of the duct throat, which for a given rotor size results in a significantly smaller total area of the device.
The effect of the angle of attack of the duct cross section has been included in most studies (Georgalas et al., 1991;Politis and Koras, 1995;Aranake and Duraisamy, 2017;Venters et al., 2018;Bagheri-Sadeghi et al., 2018). The results of these studies show that power output is considerably sensitive to the angle of attack of the duct cross section and that more power is obtained by increasing the angle of attack of the duct cross section up to the point where the flow separates inside the duct. As noted by Bagheri-Sadeghi et al. (2018), the optimal design of a DWT is on the verge of flow separation, which is often accompanied by a sharp decrease in the power output. Therefore, the accuracy of CFD simulations significantly depends on the accurate prediction of flow separation. The k − ω shear-stress transport (SST) turbulence model (Menter, 1994) is more accurate than the k− models in prediction of the flow separation for flows with adverse pressure gradients and thus is preferred when RANS simulations are used to study DWTs (Hansen et al., 2000;Bagheri-Sadeghi et al., 2018). This sharp drop in the power, which results in a discontinuity in the objective function, has implications in the choice of optimization method too, as it renders methods assuming a smooth objective function ineffective (Bagheri-Sadeghi et al., 2018).
There are a few studies on the optimization of the shape of the duct cross section for optimal C P,total such as Aranake and Duraisamy (2017), but the design space is limited by fix-ing some other design variables. For instance, Aranake and Duraisamy (2017) use a penalty function to avoid large values of the thrust coefficient, and the axial location of the rotor, the chord length of the duct cross section, and the rotor gap were not introduced as design variables. In most studies, a high-lift airfoil with the suction side inside the duct is used as the cross section of the duct (de Vries, 1979). A high-lift airfoil shape creates circulation and thereby increases the mass flow rate through the duct. Further increases in C P,total could be possible by delaying boundary layer separation, e.g., by using multi-element or slotted ducts (Igra, 1981;Gilbert and Foreman, 1983;Phillips et al., 2002;Hjort and Larsen, 2014), as the optimal design tends to be on the verge of flow separation (Bagheri-Sadeghi et al., 2018). Large flanges are often used at the exit of ducted turbines to further lower the pressure at the exit plane of the duct and increase the swallowing effect. Limacher et al. (2020) showed that large flanges lead to reduced values of C P,total .
The effect of the rotor gap as a design variable is considered in a few studies (Politis and Koras, 1995;Venters et al., 2018;Bagheri-Sadeghi et al., 2018), and although optimum rotor gaps have been examined (Bagheri-Sadeghi et al., 2018), no conclusions about the optimal rotor gap for optimal C P,total have been obtained to our knowledge. Lastly, the effect of the duct length for a given rotor diameter on the optimal design of a DWT is only considered in a few works (Georgalas et al., 1991;Politis and Koras, 1995;Ohya and Karasudani, 2010;Venters et al., 2018). Within the range of duct lengths that seem practical (lengths smaller than the rotor diameter), studies suggest that C P can be increased monotonically by increasing the duct length. However, we are not aware of any studies on the effect of the duct length on the optimal design for C P,total .
This study investigates the effect of the duct length on the optimal design for maximizing the total power coefficient (C P,total ) of a DWT having the Eppler E423 airfoil as the cross section. This entails identifying whether there is an optimal duct length and how the optimal design variables of a DWT change as the duct length varies. The results show that there is an upper limit to C P,total for a DWT, which is similar to the Betz-Joukowsky limit for open rotors. The result established here is specific to the Eppler E423 used for the duct cross section, but a similar result should hold for other duct cross sections as well. Additionally, the results indicate that the optimal rotor gap is close to the boundary layer thickness at the rotor independent of the duct length. Furthermore, this paper illustrates how the optimal axial position and the angle of attack of the duct change with increasing duct length.
The paper is organized as follows: Sect. 2 discusses the CFD model used for RANS simulations including a validation study of the actuator disc model used followed by the details of the DWT design parameterization and the optimization method used. The optimization results and how the optimal design changes with the duct length are discussed in the third section. This section also involves a sensitivity analysis of C P,total of the optimal design, as well as a comparison of the power per unit device area vs. rotor thrust of the optimal DWT and an open rotor.

Methods
Ansys Fluent 17.1 was used to solve the incompressible RANS equations with the k − ω SST turbulence model (Menter, 1994). The computational domain used is shown in Fig. 2, which extends to max(4D, 15c) upstream of the rotor and max(8D, 25c) downstream of it, where max(x, y) is the greater of x and y, D is the rotor diameter, and c is the chord length of the duct cross section. The flow was considered axisymmetric, and the rotor was modeled as an actuator disc. The pressure drop across each cell of the actuator disc was where ρ is the air density and C T,rotor is the thrust coefficient based on the axial velocity (V z ) at the rotor. The value of C T,rotor was defined in the fan boundary condition of Ansys Fluent. Using C T,rotor as a design variable means that the pressure drop across the actuator disc was not constant. However, as it is easier to interpret, all the results are given in terms of the thrust coefficient defined as where T is the thrust force on the rotor, V ∞ is the freestream velocity, and A rotor is the swept area of the rotor. At the inlet, the freestream values of the turbulence variables were set as ω ∞ = 5V ∞ D and k ∞ = νω ∞ × 10 −3 , where ν is the kinematic viscosity of air. The flow field of RANS actuator disc simulations are sensitive to freestream values of turbulence variables (Bagheri-Sadeghi et al., 2020), and the values selected correspond to recommended values by Menter (1994). Ansys Fluent uses a cell-centered finite-volume method. The pressure-based solver with the coupled algorithm and Fluent's second-order accurate schemes were used for all flow variables. The values of rotor power output and thrust were monitored to ensure iterative convergence.
The grid near the duct, which uses both structured and unstructured elements, is depicted in Fig. 3. The average nondimensional wall distance of the first grid point in the boundary layer mesh was y + 1 ≈ 1, and the aspect ratio of the first element on the airfoil was set to about 20. The thickness of the boundary layer mesh was set equal to min(δ, 0.95 n/c), where δ is the thickness of the boundary layer over the airfoil estimated from the flat-plate correlation δ = 1.1 0.16c where the 1.1 factor accounts for the longer curved surface of the airfoil compared to its chord (White, 2006). This prevents the actuator disc from penetrating the boundary layer mesh, which can result in the failure of boundary layer mesh generation. The fan boundary condition in Fluent requires identifying the direction of the positive pressure jump. The use of triangular elements as shown in Fig. 3 between the boundary layer mesh and the outside structured mesh made the grid generation more efficient. However, with the triangular elements used on the fan boundary condition, the specified direction of the fan boundary condition randomly changed from case to case, and sometimes the actuator disc became undesirably distorted. To fix this issue, a thin structured quadrilateral grid was created just downstream of the fan boundary, which is seen in Fig. 3. The width and cell sizing of this thin grid are scaled with c/D to prevent large cells near the boundary layer mesh for smaller ducts.
Two metrics were used to characterize the performance of a DWT: first, the power coefficient based on the swept area of the rotor, and, second, the power coefficient based on the exit area of the duct, (3) C P,total is a measure of the performance for a given total cross-sectional area of the device, whereas C P is the performance for a given rotor cross-sectional area.
In order to validate the actuator disc model, the axisymmetric actuator disc without a duct was simulated in the domain shown in Fig. 2. Figure 4 compares the axisymmetric RANS actuator disc simulation results on three different grids with the 1-D actuator disc momentum theory. The coarse, medium, and fine grids had about 1.6×10 4 , 6.5×10 4 , and 2.6 × 10 5 cells. The only noticeable difference between the three grids is observed at values of C T close to 1. At lower values of C T the RANS actuator disc model and the momentum theory results are visually indistinguishable. For heavily loaded actuator discs with C T > 1, the simple 1-D momentum theory is not valid and empirical correlations are often used (Glauert, 1935;Sørensen et al., 1998).

Optimization
The performance of the DWT was considered to be a function of a number of design variables, mentioned in the introduction, including the chord length of the duct cross section (c), the thrust coefficient based on the axial velocity at the rotor (C T,rotor ), the angle of attack of the duct cross section (α), the normal rotor gap ( n), and the axial position of the rotor (z rotor ). These are shown in Fig. 1. The optimization problem was to maximize C P,total as a function of normalized design variables c D , C T,rotor , α, n c , and z rotor l , where l is the duct length as shown in Fig. 1. Although the design variable was C T,rotor , all the results are presented in terms of the easier-to-interpret thrust coefficient based on the freestream velocity C T = The nondimensionalization of the rotor gap by chord length instead of rotor diameter was done to help the optimization process as the optimal normal rotor gap seemed to scale with the chord length of the duct cross section in general. The results given in the next section support this scaling. When c/D is included as a design variable, one can fix the value of Re c or Re D . If Re D is fixed, the smaller values of c/D result in values of Re c lower than the operating range of the airfoil for which Reynolds number dependency can be expected. Additionally, the larger values of c/D result in high values of Re c , which require more computationally expensive boundary layer meshes. For this reason, Re c was fixed and Re D was allowed to vary. Re c was set to a value of 3 × 10 5 . For the Eppler E423 airfoil used here, the operating range extends to Re c = V ∞ c ν as low as 1.4 × 10 5 . For lower Reynolds numbers, large flow separation is observed before the airfoil can attain its design maximum lift coefficient (Selig et al., 1996;Jones et al., 2008). For the range of chord lengths studied, Re D varied from 6.0 × 10 5 to 6.0 × 10 6 . The axisymmetric actuator disc model used here is insensitive to the Reynolds number once Re D is greater than 2000 (Sørensen et al., 1998;Mikkelsen, 2004). Within this range of Re D , the value of C P only changed by about 0.07 % when using the fine grid of the actuator disc validation study at C T = 8/9. Therefore, fixing Re c should minimize any Reynolds number effects and keep the computational expense of cases with larger c/D manageable. To determine Reynolds number sensitivity, the optimization was repeated at Re c = 1.2×10 6 using the optimal design at Re c = 3×10 5 as the starting point.
The Hooke and Jeeves direct search optimization technique (Hooke and Jeeves, 1961;Kelley, 1999;Kolda et al.,   2003) was employed in this study. Our optimization study (Bagheri-Sadeghi et al., 2018) concluded that the flow inside the duct of an optimal design for either C P or C P,total as the objective function is on the verge of separation. The flow separation inside the duct can result in a significant drop in power output, and therefore the objective function can be considered almost discontinuous close to such optimal design points. The performance of optimization methods that rely on objective function gradients is affected by the presence of such a discontinuity. The Hooke and Jeeves method, however, is less affected by such discontinuities as it does not assume a smooth objective function and only uses the objective function values to identify if a better design point is found or not.
The optimization technique starts by modifying the initial design, one design variable at a time. These exploratory moves in the design space are called steps in coordinate directions. All the design variables were scaled by their initial values, and the initial step size in each coordinate direction was set to 5 %. Based on the success or failure of these steps in the coordinate directions of the five-dimensional design space, the algorithm creates pattern directions, moves the base point, and increases or decreases the step sizes. The stopping criterion was N i=1 1 N | x i x i,0 | < 0.1, where x i is the step size in each coordinate direction, x i,0 is the initial step size, and N is the number of design variables. This means that the initially small step sizes should on average reduce by an order of magnitude for the optimization to stop. Figure 5 shows the convergence of the optimization technique towards the optimal design. The optimization was repeated with a different initial design, and the optimization approached the same optimal design point.

Results and discussion
The optimization was performed at a constant Re c = 3 × 10 5 as described in the Methods section. The optimal design achieved C P,total = 0.69 at c/D = 0.18, which corresponds to l/D = 0.15. The values of design variables and duct length at this optimal design are given in the second row of Table 1. The optimal design was also simulated on a finer grid with about 1.11 × 10 6 cells compared to 2.8 × 10 5 cells of the original grid. There was only a 0.03 % difference between C P,total values. Also, simulating on a larger domain extending 1.4 times further in each direction, the value of C P,total changed by about −1 %, which was considered sufficiently accurate for the optimization study carried out here.
To verify that there is indeed a maximum in C P,total with duct length, optimizations were carried out at several other fixed c/D values at Re c = 3 × 10 5 . The optimal designs at these fixed values of c/D are shown in Table 1 as well. Figure 6 shows the power coefficients of the designs for optimal C P,total . These additional optimization studies confirm that there is an optimal duct length for a DWT that maximizes C P,total and that this maximum is greater than the Betz-Joukowsky limit. Thus, for applications desiring the greatest power per unit device area, the duct length should be around 15 % of the rotor diameter. As the results here are specific to Eppler E423, further studies are needed to verify this conclusion. The values of power coefficient based on the rotor swept area (C P ) are also shown in Fig. 6 for the designs of Table 1. The values of C P seem to increase almost linearly with c/D. This suggests that significantly larger values of C P can be obtained by increasing the duct length, but this will result in a lower C P,total .
Both here and in Bagheri-Sadeghi et al. (2018) the values of optimal C T were close to 0.9. For all optimal designs of Table 1, the value of the thrust coefficient is close to C T = 8/9, which is the optimal thrust coefficient of open wind turbines at the Betz-Joukowsky limit. This is in agreement with the momentum analysis and CFD studies of DWTs, which conclude that the optimal C T values of open and ducted wind turbines are similar ( van Bussel, 2007;Jamieson, 2008;Bagheri-Sadeghi et al., 2018).
The geometry and streamlines of the first three configurations in Table 1 are shown in Fig. 7b-d along with those for an open rotor (in Fig. 7a). The optimal angle of attack of the duct cross section decreased with increasing duct length. For the open rotor, the streamlines close to the tip are at an angle of about 30 • . Furthermore, at c/D = 0.05 almost all of the duct cross section can be considered to be in the vicinity of the strong divergence of streamlines close to the rotor tip. This explains why for the small c/D = 0.05 the flow is still attached at α = 46 • . The presence of the small duct with z rotor / l = 0.98 adds extra suction inside the duct without significantly increasing the total area of the device and achieves about 10 % more power per unit device area than an open rotor (i.e., 10 % more than the optimal power output of an actuator disc RANS simulation, which was about 0.6 as shown in Fig. 4). For the optimal design with a variable c/D (i.e., c/D = 0.18), the increase in C P,total is 15 % over an open rotor, and the angle of attack is reduced to 31 • . Note that the angle of streamlines of the actuator disc with respect to the horizontal decreases further downstream as the rotor wake recovers in Fig. 7. Similarly, further upstream of the rotor for the actuator disc case, the angle of streamlines decreases. At greater duct lengths, the airfoil of the duct cross section becomes less influenced by the presence of the actuator disc and hence the maximum α without a large flow separation decreases. This trend continues for c/D = 0.5 (see Table 1).
The optimal normal rotor gap ( n/c) decreased with increasing duct length, but the ratio of the rotor gap to the estimated boundary layer thickness at the rotor ( n/δ rotor ) stayed nearly constant. A smaller rotor gap means a smaller exit area of the duct, and therefore it helps to increase C P,total . On the other hand, the presence of the rotor gap results in an annular jet of high-velocity air which imparts momentum The optimization performed at fixed c/D at Re c = 3 × 10 5 . b The optimization performed with c/D as a design variable at Re c = 3 × 10 5 . Figure 6. The power coefficients of designs for optimal C P,total .
to the boundary layer and helps the flow stay attached. A n that is too small or too large weakens this annular jet. The jet is easier to see in the contour plot of c/D = 0.35 in Fig. 7d. For all of the optimal cases shown in Table 1, n/δ rotor was close to 1. When the chord is small relative to D, the dominant length scale is c. This determines the flow scales as well as the boundary layer thickness. Also, note that because Re c is held constant, the boundary layer thickness at the trailing edge δ scales linearly with c. Therefore, as optimal n/δ rotor is nearly constant, optimal n/c only slightly decreases as the optimal rotor position is moved further upstream, resulting in a smaller boundary layer thickness at the rotor (δ rotor ) compared to the trailing edge. This justifies the scaling of n by the chord length of the duct.
Similar to Bagheri-Sadeghi et al. (2018), the design for optimal C P,total resulted in downstream placement of the rotor. The optimal rotor position moved further upstream in the duct as c/D increased. Note that at greater duct lengths the annular jet is stronger as can be seen in Fig. 7. This further upstream placement of the rotor means that the annular jet formed between the rotor tip and duct can exchange momentum with a larger portion of the boundary layer, which may help the DWT attain more power per total unit device area without flow separation.
The sensitivity of the total power coefficient of the optimal design to different design variables x i is shown in Fig. 8. The greatest sensitivity is to the thrust coefficient of the rotor which matches the results of previous studies (Venters et al., 2018;Bagheri-Sadeghi et al., 2018) and illustrates the impor-tance of rotor design in achieving optimal power output from a DWT (a rotor design approach based on the blade element momentum method using axisymmetric RANS actuator disc simulations as input is discussed by Kanya and Visser, 2018). For α and z rotor / l, an increase in the design variable causes flow separation inside the duct and a significant drop in the power output. This is demonstrated by an increased sensitivity of C P,total to increase in α and z rotor / l compared to their reduction. The sensitivity to the reduction in z rotor / l is mainly driven by changes in the total area of the device rather than changes in the power output. Therefore as concluded in Bagheri-Sadeghi et al. (2018) a smaller DWT with similar power output (i.e., a greater C P,total ) can be designed by placing the rotor further downstream of the duct. The observation that the optimal design is on the verge of flow separation with respect to α agrees with the results of Bagheri-Sadeghi et al. (2018) obtained for the design for optimal C P . The flow of the design for optimal C P in Bagheri-Sadeghi et al. (2018) is on the verge of separation with respect to rotor gap and shows little sensitivity to z rotor / l, whereas the optimal design for C P,total shows the least sensitivity in Fig. 8 to n/c, and a slight increase in z rotor / l leads to flow separation. Figure 9 compares the C P,total vs. C T curves of the optimal DWT with the RANS actuator disc model. Note that for an open rotor turbine C P,total = C P . Close to the optimal design, at a given C T the DWT produces about 15 % more power per unit device area. However, at lower values of C T , this increase in C P,total becomes smaller (e.g., at C T ≈ 0.55 the DWT produces only 6 % more power per total device area).  In other words, at lower rotor loadings, the power output of the optimal DWT per total device area approaches that of an open rotor turbine. Generally, when using actuator disc simulations, the increase in velocity of a DWT compared to an open wind turbine should be fixed so the ratio of C P values (and therefore C P,total values in Fig. 9 as the A total is fixed) should be independent of C T (Hansen et al., 2000). Indeed, if the curves of Fig. 9 were plotted for another DWT design not on the verge of flow separation (e.g., the optimal design but at a reduced α) the ratio of the power coefficients would stay constant. The constant ratio of the power coefficients is not seen here at the optimal design because the optimal design is on the brink of flow separation, and at lower values of C T the rotor is less effective in keeping the flow attached and the induced velocity decreases. As the thrust coefficient reduces, the high-speed annular jet becomes weaker and therefore cannot keep the flow attached, and gradually the region of flow separation extends further upstream. The flow separation with thrust coefficient reduction is gradual rather than Figure 9. Comparison of C P,total vs. C T for RANS actuator disc model and the DWT for optimal C P,total . the abrupt separation which can be observed by slightly increasing the angle of attack or z rotor / l in Fig. 8.

Reynolds number sensitivity
To examine Reynolds number sensitivity, the simulation of the optimal design was repeated at Re c = 1.2 × 10 6 . The Reynolds numbers used in this study correspond to typical values for DWTs, which typically target residential applications. For example, at a wind speed of V ∞ = 10 m s −1 at standard sea-level temperature and pressure, the Re c = 3 × 10 5 design corresponds to c = 0.44 m, D = 2.4 m, and P = 2.3 kW, and Re c = 1.2×10 6 corresponds to c = 1.76 m, D = 9.8 m, and P = 38 kW.
The result of optimization at this higher Reynolds number is shown in the second row of Table 2. The maximum C P,total slightly increased to 0.70. The optimal value of c/D did not change to the precision of the optimization. The optimal values of n/c and z rotor / l also did not change, and the optimal C T and α only varied slightly. This confirms that at Re c = 3 × 10 5 the Reynolds number dependency was small.
The contour plot of eddy viscosity ratio µ t µ for the optimal design is shown in Fig. 10a. On the suction side of the duct at Re c = 3 × 10 5 the turbulence model is activated near the leading edge of the duct. This suggests that the optimal design should be insensitive to Reynolds number. The eddy viscosity ratio at Re c = 1.2 × 10 6 is shown in Fig. 10b. At this Reynolds number the turbulence model is activated slightly further upstream and the power output increased to C P = 0.82. The small increase can be explained by the reduced flow separation near the trailing edge at this greater Reynolds number, which increased the effective exit area of the duct. Also, note that the turbulence model on the pressure side becomes activated further upstream as well. However, this should not matter in identifying the optimal design as the flow stays attached on the outside of the duct. The k − ω SST turbulence model, with the apparent transition near the leading edge, indicates that the flow is entirely turbulent over the airfoil, and therefore the results should not change significantly even at higher Reynolds numbers.

Conclusions
The optimal design of a ducted wind turbine with the Eppler E423 airfoil as its cross section was investigated using CFD simulations of axisymmetric RANS equations with the k − ω SST turbulence model and an actuator disc. The total power coefficient characterizing the power output per total area of the device (C P,total ) was used as the design objective. The design variables included the chord length and the angle of attack of the duct cross section, the thrust coefficient, the rotor axial position, and the tip clearance of the rotor.
The results demonstrate the existence of a maximum C P,total for ducted wind turbines, which sets an upper limit for ducted wind turbines similar to the Betz-Joukowsky limit for open rotors. For the Eppler E423 airfoil, this maximum was obtained for a duct length of about 15 % of the rotor diameter. With this duct length a C P,total of 0.70 was obtained, which exceeds what can be obtained with an open rotor by 16 %. This is the first time that an optimal duct length has been identified, although the optimization was for C P,total and not C P . The value of C P increased almost linearly with duct length over the range investigated.
Additionally, the results of optimization at fixed c/D suggested that the optimal value of the design variables can sig- nificantly change with the duct length. In agreement with previous theoretical and numerical studies, for all duct lengths the optimal thrust coefficient (C T ) was almost 0.9, which is similar to open rotor turbines. The results also showed that the optimal design for C P,total was most sensitive to the thrust coefficient of the rotor, which indicates the importance of proper rotor design. At lower-than-optimal thrust coefficients, the power per unit device area (C P,total ) of the optimal DWT design gradually approached that of an open wind turbine. This can be explained by considering that the optimal design was on the verge of flow separation and the rotor became less effective in keeping the flow attached as the thrust coefficient decreased. The optimal angle of attack of the duct cross section decreased significantly with increasing the duct length. Additionally, the optimal design was on the verge of flow separation with respect to the angle of attack of the duct cross section and the axial position of the rotor.
The optimal normal rotor gap was close to the boundary layer thickness at the rotor. Therefore, the optimal normal rotor gap scaled proportional to the chord length of the duct cross section as the turbulent boundary layer thickness almost linearly increases with the chord length of the duct cross section. This gap is needed to create the high-velocity annular jet, which helps keep the boundary layer attached.
The optimal rotor position was at the rear of the duct, but at greater values of duct length it moved further upstream in the duct. This further upstream position was more effective at eliminating flow separation and hence allowed greater values of C P,total . Data availability. Data are available upon request from the corresponding author.
Author contributions. NBS contributed to the methodology, ran the simulations, post-processed the data, and wrote the first draft of the paper. BTH supervised the study and contributed to the conceptualization, methodology, writing, and revision of the paper. KDV contributed to the conceptualization and revision of the paper.
Competing interests. The authors declare that there is no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Review statement. This paper was edited by Jens Nørkaer Sørensen and reviewed by Peter Jamieson, David Wood, and one anonymous referee.