Ducted Wind Turbine Optimization and Sensitivity to Rotor Position

The design of a ducted wind turbine modeled using an actuator disc was studied using RANS CFD simulations. The design variables included the rotor thrust coefficient, the angle of attack of the duct cross-section, the radial gap between the rotor and the duct, and the axial location of the rotor in the duct. Two different power coefficients, the rotor power coefficient (based on the rotor swept area) and the total power coefficient (based on the exit area of the duct) were used as optimization objectives. The optimal value of thrust coefficients for all designs was nearly constant having a value between 0.9 and 1. 5 The rotor power coefficient was sensitive to rotor gap but was insensitive to the rotor’s axial location for positions ranging from upstream of the throat to nearly half the distance down the duct. Compared to the design that maximized rotor power coefficient, the design for maximal total power coefficient was characterized by a smaller angle of attack, a smaller rotor gap and a downstream placement of the rotor. The insensitivity of power output to the rotor position implies that a rotor placed further downstream in the duct could produce the same power with a considerably smaller duct exit area and thus a greater total 10 power coefficient. The design for that maximized total power coefficient exceeded Betz’s limit with a total power coefficient of 0.67. Copyright statement. TEXT


Introduction
A properly designed duct placed around a wind turbine can increase power output by increasing the mass flow rate through the rotor.Ducted wind turbines (DWTs) are also called diffuser-augmented wind turbines (DAWT) or shrouded wind turbines.(Lilley and Rainbird, 1956) performed a one-dimensional momentum analysis of DWTs and concluded that higher expansion ratios of the duct and more subatmospheric pressures at the exit plane of the duct result in higher power outputs.They also suggested wind tunnel tests with screens of different porosities to model the pressure drop across the rotor.Such experimental tests were performed by Igra (1976Igra ( , 1977Igra ( , 1981)), Foreman et al. (1978), Gilbert et al. (1978), and Gilbert and Foreman (1979).The negative effect of flow separation on the power output of DWTs was observed, and various methods of preventing separation were investigated.Also, experimental tests with real turbines were performed, and the power augmentation of DWTs was demonstrated (Igra, 1981;Gilbert andFore-man, 1979, 1983).As the duct can be considered an annular wing (de Vries, 1979) with higher lift, meaning more suction and circulation, high-lift airfoils were used from early experimental studies.
Using lifting line theory for the rotor and modeling the duct as a superposition of vortex and source rings, (Koras and Georgalas, 1988) and (Georgalas et al., 1991) studied the power output of DWTs with airfoil cross sections and large rotor gaps (the clearance between the tip of the rotor and the duct) as a function of several design variables including the angle of attack of the duct cross section, the chord length of the duct, the maximum camber of the duct cross section, and the relative position of the rotor with respect to the maximum camber point of the duct cross section.They found a linear increase of power with duct chord length and angle of attack of the duct cross section.They also concluded that the effect of rotor position on the power output was weak.(Politis and Koras, 1995) extended the previous work to DWTs with any rotor gap.
N. Bagheri-Sadeghi et al.: Ducted wind turbine optimization Axisymmetric computational fluid dynamics (CFD) models were used (Phillips et al., 1999, 2002and Phillips, 2003) to improve the design of the first full-scale DWT built (the Vortec 7).(Hansen et al., 2000) performed a CFD study of DWTs and used the k − ω shear stress transport (SST) turbulence model for the axisymmetric model as it is more sensitive to adverse pressure gradients (Menter, 1994) and can be more accurate in predicting flow separation.Another similar CFD study was performed by (Abe and Ohya, 2004), where effects of rotor loading and the incidence angle of the duct on power output of a flanged DWT were examined and compared with experimental data.(Ohya et al., 2012) and (Kardous et al., 2013) did further CFD simulations of the flanged DWT with the rotor modeled as an actuator disc and found good agreement with wind tunnel data.
Van Bussel (1999Bussel ( , 2007) ) analyzed DWTs using 1-D momentum theory and concluded that optimal coefficient of thrust in a DWT is similar to an open rotor equal to 8/9.He also concluded that experimental power coefficients based on the exit area of the duct (the total power coefficient) above 0.5 have not been achieved yet and that very significant back pressure reductions are needed to achieve values of total power coefficients significantly above Betz's limit.(Jamieson, 2009) also used a similar momentum analysis and derived the same value of 8/9 for optimal loading on the rotor and noted that it should be independent of duct design.(Werle and Presz, 2008) in another study based on 1-D momentum analysis found that the maximum attainable power from a DWT is determined by shroud force coefficient, C s = F s /T , where F s is the axial force on the duct (shroud) and T is the thrust of the rotor.(Hjort and Larsen, 2014) used an axisymmetric CFD model with an actuator disc modeling the wind turbine for a multi-element DWT.They characterized the performance of the DWT using power coefficients based on the exit area of the duct with values well above Betz's limit.(Aranake and Duraisamy, 2017) also utilized an axisymmetric Reynoldsaveraged Navier-Stokes (RANS) solver with an actuator disc model for the turbine to optimize the airfoils used for the duct cross section and blades and verified the result with 3-D simulations.(Venters et al., 2017) investigated the optimal design of a DWT using the same approach (i.e., using a RANS solver and actuator disc model).The design variables investigated were the rotor loading, the angle of attack of the duct cross section, the rotor gap, and the axial position of the rotor.They used a response surface fitted to a number of design point calculations, and then the surface was searched using the NLPQL (Nonlinear Programming by Quadratic Lagrangian) algorithm (Schittkowski, 1986).They concluded that rotor loading is the main factor defining the performance of the DWT with the coefficient of thrust almost constant (close to 1) for different duct sizes.The power output of DWT was sensitive to the angle of attack of the duct cross section.However, the results for effect of the rotor gap and axial position of rotor were not conclusive.This paper improves on the work of (Venters et al., 2017) with a more accurate CFD model, a direct optimization technique, and a wider range of design variables.One of the goals of this study is to continue the investigation of (Venters et al., 2017) into how the objective of the optimization changes the optimal design.Specifically, (Venters et al., 2017) examined two objective functions, the rotor power coefficient and the total power coefficient.Their results indicated that the optimal design changes significantly depending on the objective function, but the results for optimizing the total power coefficient did not converge to an optimal solution.The goal of this work is to identify an optimal configuration for this objective function.
The paper is organized as follows.The details of the CFD model along with an evaluation of two different pattern search optimization methods are given in Sect. 2. Optimization results with the objective of maximizing the rotor power coefficient are given in Sect.3, and the variation of the rotor power coefficient and flow field with different design variables is presented.In Sect.4, optimization results are presented for the objective of maximizing the total power coefficient.These results are compared with the goal of understanding how the optimal axial position of the rotor depends on the optimization objective.

Method
A two-dimensional axisymmetric numerical model was developed in ANSYS Fluent 17.1 to simulate the flow field of a DWT.The wind turbine rotor was modeled as an actuator disc with a pressure drop, p, given by where ρ is the air density and C T, rotor is the thrust coefficient based on the axial velocity, V z , at the rotor.The thrust force, T , is given by where D is the rotor diameter.The extracted power, P , is given by Clearly, with an actuator disc model, rotor blade efficiency losses are not considered.The design variables, shown in   (Selig et al., 1996).Two power coefficients, , were used as objective functions for the optimizations and to compare the performance of different DWT designs.The domain and mesh used for the simulations are shown in Fig. 2.These were defined to ensure mesh independence for power coefficients as the design variables were varied.The domain extended 15 duct chord lengths upstream of the rotor and 25 chord lengths downstream.Numerical tests showed that this domain size gave power coefficient values that were independent of the domain size to two significant digits.As all optimizations were done with the same domain, this was deemed large enough to accurately calculate changes in the solutions with the design variables.The shape of the domain at its top made a distinct transition between inflow and outflow boundaries, which eliminated convergence issues due to reverse flow through outlet boundaries.The reverse-flow issue occurred when the inlet and outlet were smoothly connected.The mesh of Fig. 2 consisted of about 500 000 elements.The duct boundary layer mesh had a growth rate of 1.1, and the first mesh point was set at y + ≈ 1.
The boundary layer thickness was calculated as a function of Re c for each case, and enough inflation layers were used to span the entire boundary layer.The quality-based smoothing option in Fluent was used to improve the mesh quality.ANSYS Fluent's k − ω SST turbulence model was used to solve the incompressible Navier-Stokes equations.The pressure-based solver was chosen with the coupled scheme used for the pressure-velocity coupling.Gradients were calculated using the Green-Gauss node-based method, and second-order discretization schemes were used for pressure, momentum, turbulent kinetic energy, and specific dissipation rate.The output power, thrust, and drag coefficient of the duct were calculated and monitored at each iteration to ensure convergence.

Optimization techniques
For most of the optimization results, a pattern search method (Powell, 1964) was used to find the optimal design of the DWT.Optimizations were first performed with C P as the objective function and then with C P, total as the objective function.The optimization for both objective functions started from the same set of design variables (C T, rotor = 0.816; α = 25 • ; r/D = 0.03; and z/c = 0.14).In our implementation of Powell's method a quadratic interpolation of the function values is used to identify the optimal step length to move the design point in the coordinate or pattern directions.The optimization was stopped when the improvements obtained from the optimization methods were within a specified tolerance.The termination criterion was C P, optimal −C P,0 C P,0 +C P, optimal < 0.005, where C P,0 is the initial value of C P at the beginning of a search cycle.The design variables were then varied to determine the sensitivity of the objective function to the design parameters in the vicinity of the optimal design point.
Powell's method is known to be slow converging for objective functions that are discontinuous.As shown in Sect. 3 and Sect.4, it was observed that the optimal design points were on the verge of flow separation along the duct airfoil and that separation was accompanied with a large drop in power output.Therefore, the objective functions were nearly discontinuous at the optimal design point.With such an objective function, the optimizer worked inefficiently in finding the optimal step length.Also, when the optimizer moved the design point close to a discontinuity, it moved away from that point very slowly.Figure 3a shows the history of an optimization using Powell's method with z/c = 0.05 fixed and the design variables being C T, rotor , α and r/D.The optimization method was stopped at about 100 iterations without meeting its termination tolerance.At that point the search algorithm was jumping around significantly.This is shown in Fig. 3b, which shows the search history of C T , α points.The optimal point is shown as a red triangle which occurs at C T ≈ 1 and α ≈ 27 separation where the optimizer had trouble finding the optimal step length and was stuck close to the function discontinuity.The maximum C p obtained by the search method was 1.031.
The same problem was subsequently approached with the Hooke and Jeeves method (Hooke and Jeeves, 1961) with the same termination criterion and starting point.The optimization history is shown in Fig. 4a.This time the optimization algorithm reached the optimal design in only 16 function evaluations and reached the termination criterion in 32 function evaluations.In addition, a better design with C P = 1.053 was found.The more efficient performance of the Hooke and Jeeves method can also be observed in Fig. 4b, which shows the (C T , α) search points.The optimum point is again shown as a red triangle and occurred at C T = 0.97 and α = 30 • .Because the Hooke and Jeeves method does not fit an analytic function to the function values, it did not face the same difficulty when it got close to sharp variations in the objective function.
Although the Hooke and Jeeves method was more efficient, unless stated otherwise, most of the results obtained below were found using Powell's method as this was the first method implemented.As this method did not always satisfy the optimization stopping criterion, we call the optimized de-signs "near-optimum" points.The search history shown in Fig. 3b is fairly typical for the cases shown below.For each search direction, function evaluations at design points in the search directions of roughly ±1 % were evaluated with no increases in the optimal value.This, however, does not preclude the possibility of a slow variation along a ridge in the optimization function, which is essentially the reason we see 5 % differences in the optimized design variables between the two different optimization methods for this example.In many of the figures below, individual design variables are varied around the optimal point to give a further sense of the sensitivity to the design variables.

Design for optimal C P
The middle column of Table 1 shows the near-optimal design found with Powell's method when optimizing for maximum C P .The design values are close to what was observed by (Venters et al., 2017).(Venters et al., 2017) used a smaller chord length for the duct (c/D = 22.5 %) and a different turbulence model (k − realizable) and obtained a maximal value for C P = 1.00 at C T = 1.08 and α = 37.5 • .Our results predicted a lower value of optimal C T and α, which could be because of the more accurate turbulence model as the k − ω  SST turbulence model is known to be more accurate in prediction of flows with significant adverse pressure gradient and flow separation (Menter, 1994).
The results for the variation of C P with C T are shown in Fig. 5.The highest C P in this plot is at C T ≈ 0.93.When using the Hooke and Jeeves optimization, optimal C T values very close to 1 were observed, which is closer to that observed by Venters.One-dimensional momentum analysis done by (van Bussel, 1999) and (Jamieson, 2009) predicted that the optimal C T for a ducted turbine would be independent of duct design and have a value of 8/9, which is the same as that of an open rotor.The plot of Fig. 5 also shows the curve for an open rotor as predicted by actuator disc theory.Similar to an open rotor, increasing the loading on the rotor beyond the near-optimal design point of the DWT reduced the mass flow rate through the rotor and thus its output power.Also similar to an open rotor, at loadings less than the near-optimal design point, the flow rate through the rotor was larger, but the pressure drop was too low to obtain optimal power.In the ducted case, however, the reduction in C T had an additional effect, which was to cause flow separation in the duct.As shown next, there is a strong coupling between the coefficient of thrust, the angle of attack of the duct, and separation.Increasing the angle of attack or decreasing the coefficient of thrust can lead to separation.The effect of changing α is shown in Fig. 6.When α was increased beyond the near-optimal design point, a large flow separation resulted, which was accompanied by a sharp decrease in the output power.The flow field of the near-optimal design is shown in Fig. 7.The effect of increasing α on the flow field is shown in Fig. 8. Comparing the two flow fields, it is apparent that the small increase in angle of attack leads to a large separated region at the trailing edge of the airfoil.The separated region effectively reduces the exit area area of the duct, resulting in the capture of a smaller upstream flow area and a smaller power extraction.Similarly, reducing α from the near-optimal design also resulted in a decrease of C P because of the decreased exit area.
Likewise, as shown in Fig. 9, if rotor gap, r/D, was increased beyond the near-optimal design point (while keeping other design variables constant), a large power drop was observed due to flow separation and the streamlines appeared similar to Fig. 8. Decreasing r/D also reduced the power output of the rotor.Reduction of the rotor gap results in a decrease in the exit area of the duct, which could be the reason for the reduction in power.The dependence of C P on the axial position of the rotor, z/c, can be seen from Fig. 10.As z/c was varied from the near-optimal design, the power output did not change significantly.To better understand the effect of axial location on the power output of the rotor, the design was optimized using the Hooke and Jeeves pattern search method at a number of fixed z/c values from 0.05 to 0.35.The results shown in Fig. 10 confirm that C P within the range of z/c values shown is not very sensitive to the axial position of the rotor.The higher values of C P shown are due to better performance of the Hooke and Jeeves search algorithm as discussed in Sect. 2. This result shows that one can place the rotor anywhere from upstream of the throat to halfway down the duct and obtain similar performance.

Design for optimal C P, total
The last column of Table 1 shows the near-optimal design parameters when C P, total was the objective function, and Fig. 11 shows the geometry and flow field of the near-optimal design.Compared to the design for optimal C P , when a DWT was designed for optimal C P, total , the values of α and r/D were decreased, whereas z/c was increased.All of these changes have a similar effect: to decrease the exit area of the duct, which is in the denominator of the objective function.
The value of C T of the near-optimal design (0.87) was close to the optimal C T when C P was optimized using Powell's method (0.93).It is also close to the optimal value for an open rotor which (van Bussel, 1999) and (Jamieson, 2009) predicted.However, there is some ambiguity in the preciseness of this value because both Venters and the Hooke and Jeeves optimization showed values near 1.00 when optimizing C P .
The variation of C P, total with z/c is presented in Fig. 12.All other design variables were fixed at the near-optimal design point for C P, total as given in Table 1 as z/c was varied.Since C P, total depends on both power output and the exit area of the duct, the values of C P at each design point are also shown so that variations due to changes in exit area or 2.17e+01 2.07e+01 1.96e+01 1.85e+01 1.74e+01 1.63e+01 1.52e+01 1.41e+01 1.30e+01 1.20e+01 1.09e+01 9.78e+00 8.70e+00 7.61e+00 6.52e+00 5.44e+00 4.35e+00 3.26e+00 2.17e+00 1.09e+00 0.00e+00 power can be better understood.Similar to Fig. 10 the power of the rotor, C p , is not very sensitive to z/c when z/c < 0.5.The exit area decreases as z/c is increased, which makes C P, total increase as the rotor is moved towards the exit of the duct.When z/c is increased past 0.5, the power extracted decreases, but C P, total continues to increase because of the decreasing exit area.Just past the optimal value of z/c the flow separates, leading to a sharp decrease in both C P and C P, total .
The best value obtained here for C P, total = 0.67 was above Betz's limit and was also higher than the previous result by (Venters et al., 2017) of 0.621.Thus, it is possible to extract more power per unit device area using a ducted turbine than when using an open rotor.This is in agreement with theoretical predictions by (van Bussel, 2007) at high back pressure reductions.To obtain this value of C P, total , the rotor must be at the rear of the duct.The optimal design of a ducted wind turbine characterized by the thrust coefficient of the rotor C T, rotor , the angle of attack of the duct cross section α, the rotor gap r/D, and the axial location of the rotor z/c was investigated.The optimal design was significantly different when different power coefficients C P (based on rotor area) and C P, total (based on the exit area of the duct) were used as design objectives.Compared to the design for optimal C P , the design for optimal C P, total resulted in a duct with smaller α and r/D and a rotor placed at the rear of the duct rather than towards the front.This type of design has been experimentally investigated in (Kanya and Visser, 2018).The design for optimal C P, total attained C P, total = 0.67, which was above Betz's limit.This optimal design was on the brink of flow separation; increases in α, decreases in C T , or increases in r/D all resulted in flow separation and a sharp decrease in power output.The Hooke and Jeeves optimization method was found to be more efficient in finding the optimal designs than Powell's method, which was attributed to this sharp variation in C P around the design point.
Although optimal design data obtained from Powell's method should represent characteristics of a good design based on C P or C P, total , due to convergence issues they should be considered approximate.Also at higher Reynolds numbers (i.e., with larger DWTs) the flow field and separation characteristics of the DWT may change, which in turn can change the characteristics of the optimal design.Additionally, including effect like swirl and center body in the CFD model should give more realistic results and can change the optimal design.
Fig. 1, were the thrust coefficient of the rotor C T, rotor = T 1 2 ρV 2 z A rotor , the angle of attack of the duct cross section α, the radial gap of the rotor r/D, and the axial location of the rotor z/c.Because the thrust coefficient based on the

Figure 2 .
Figure 2. The domain and mesh (a); the zoomed-in view of mesh showing the structured mesh upstream and downstream of the duct (b); and the more-zoomed-in view of mesh showing the boundary layer mesh, the structured mesh of the actuator disc, and the unstructured triangular mesh between duct and the actuator disc.

Figure 3 .Figure 4 .
Figure 3. Optimization of C P using Powell's method for z/c fixed at 0.05.The red triangle denotes the optimized solution.

Figure 11 .Figure 12 .
Figure 11.Flow field of a near-optimal design based on C P, total .

Table 1 .
Comparison of designs based on C P and C P, total Design based on C P Design based on C P, total Streamlines at α = 28.2• (the streamlines are colored based on velocity magnitude in m s −1 ).