Characterization of aerodynamic performance of ducted wind turbines A numerical study

The complex aerodynamic interactions between the rotor and the duct has to be accounted for the design of ducted wind turbines (DWTs). A numerical study to investigate the characteristics of flow around the DWT using a simplified duct–actuator disc (AD) model is carried out. Inviscid and viscous flow calculations are performed to understand the effects of the duct shape and variable AD loadings on the aerodynamic performance coefficients. The analysis shows that the overall aerodynamic performance of the DWT can be increased by increasing the duct cross-sectional camber. Finally, flow fields using viscous calculations are examined to interpret the effects of inner duct wall flow separation on the overall DWT performance.

insights on the importance of the viscous effects and on the limitations of the inviscid flow model for the design of DWT.The rotor is modelled as an AD with radially uniform thrust coefficient.The following section discusses the governing performance coefficients for a duct-AD model.Section 3 describes the computational settings and parameters with a brief description of the numerical methods used.Section 4 explains the verification and validation studies.The duct shape modifications is carried out using the class-shape transformation (CST) method and is discussed in Section 5. Finally, some insights on the performance coefficients with respect to AD loading will be discussed in Section 6, together with flow analysis obtained using RANS simulations.

DUCT-AD FLOW MODEL
The incompressible flow past a turbine is modelled by a flat AD of infinitesimal width.The AD exerts a constant thrust force T AD per unit surface, which corresponds to a nondimensional thrust force coefficient: where  is the fluid density, U ∞ is the free stream velocity, and S AD is the AD surface area.
For a duct-AD configuration, an additional thrust force F x , exerted by the duct on the flow, appears.In order to quantify the relative contribution of the duct force F x on the AD, the duct force coefficient C x is defined: The presence of duct force F x results in an overall mass flow rate .m = S AD U AD , when calculated across the AD plane.Although the assumption of constant loading C T is employed, the free stream velocity computed at the AD plane cannot be regarded as radially uniform, especially for a duct-AD model where the accelerations are higher as compared with a bare AD.In order to account for the non-uniform flow in the radial direction, the mean AD velocity U AD is defined by integrating the differential terms for the local axial velocity U x across the AD surface: This leads to a power coefficient of the duct-actuator model using AD surface S AD as the reference area: It is worth mentioning that some studies adopt a different definition for power coefficient in which the reference area is taken at the duct exit section. 11As a result, the power coefficient obtained is smaller than the one obtained using Equation (4).So that, with the help of reference areas S AD and S exit , the difference in the power coefficients can be expressed as The velocity at the AD U AD and the duct force coefficient C x are obtained by postprocessing the numerical solutions.The normalized velocity component is then combined with the respective AD loading C T leading to the power coefficient C P of the duct-AD model as in Equation ( 4).

COMPUTATIONAL SETTINGS AND PARAMETERS
A commercial DWT model DonQi ® is selected as the reference case.The geometrical parameters of the reference DWT model are taken from the existing studies conducted by the authors. 12,13e computational set-up replicates the experiments conducted in the closed-loop open-jet (OJF) wind tunnel facility at the Delft University of Technology (see Figure 2).A 0.002-m-thick AD model with porosity  = 70% is used to mimic the rotor.It results in an equivalent thrust   The tip clearance between the AD tip and the duct nozzle surface is 2.0%c.The numerical study is performed at a fixed Re of 4.5 × 10 5 as in the experiments and various AD loadings 0 ≤ C T ≤ 0.89 are investigated.
In this section, the discretization details of the two computational codes used will be briefly described.For an in-depth description, the reader can refer to de Oliveira et al 14 and Dighe et al. 15

Panel method
A two-dimensional potential flow panel method is used to model the steady isentropic incompressible flow field around the duct-AD model following the work of de Oliviera et al. 14 The governing flow equations are a simple form of the Euler equations.The AD is represented as a nonconservative force term in the vorticity equation.The duct surface is defined using a distribution of vortex panels with the local slope tangential to the duct surface such as to reproduce the desired duct cross-sectional shape; see Figure 3.In essence, a uniform distribution of vorticity on the duct panels is assigned by assuming the Kutta condition.The Kutta condition states that the flow leaves the sharp trailing edge of the duct smoothly.The simplified assumption of uniform vorticity distribution over the panels represents a simplification of real physics and prevents separation on the duct surface even for larger pressure gradients.The duct surface discretization is based on the constant spacing approach.The streamwise discretization is non-uniform, with initial panel length equal to 1.0%c just behind the AD, and increasing gradually in length as the wake expansion settles further downstream.This approach allows to capture the solutions of interest and also to reduce the overall computational costs.Due to the non-linear nature of solutions obtained, the panel method, in any case, should be considered to be a successor of the analytical models used for the analysis of DWT.The panel method is particularly appealing for repeated analysis design arrangements due to its short execution time.

RANS method
A commercial CFD tool ANSYS Workbench ® has been used for a complete viscous solution of steady incompressible flow around the duct-AD model.The governing flow equations are the Reynolds-averaged Navier Stokes (RANS) equations in planar coordinates.The two-dimensional computational domain is shown in Figure 4, where the distance from the AD location to domain inlet and outlet are 12c and 24c, respectively.The computational grid consists of quadrilateral cells with maximum y + value of 1 on the duct walls (see Figure 5).Boundary conditions are the following: a uniform velocity inlet, zero gauge static pressure outlet, no-slip walls for duct surface, fan boundary condition for AD implementation, and symmetry for centre line axis (see Figure 4).The k- SST (shear stress transport) model is used as the turbulence model.Preliminary investigations show good agreement with the experiments. 13The RANS solution was obtained using the coupled algorithm; it offers robustness and faster convergence for steady state flows as compared to the segregated solution schemes.A least squares cell-based method is used to evaluate the pressure gradient, with pressure and momentum equations solved using a second-order upwind differential scheme.The convergence criteria is set to 10 −6 for all the residuals.A typical converged RANS solution (single duct-AD configuration) with approximately 0.1 million mesh elements is obtained in roughly 0.5 hour on a multicore work-station desktop computer.

GRID INDEPENDENCE STUDY AND SOLUTION VALIDATION
Solution verification and validation have been performed to ensure the accuracy of the present CFD approach.The results have been published previously. 13,15grid independence analysis for the RANS simulations has been carried out using three grid sizings where the refinement factor in each direction is approximately 1.5.The decay of the axial velocity distribution across the AD in the wake region, from the AD location to the outlet of the computational domain, is taken as reference.The results of the grid independence study are shown in Figure 6.The medium refined grid, consisting of approximately 137 600 elements, is selected for the rest of the study.
A validation study have been performed where the CFD results, obtained from the panel method and the RANS method, are compared against the experimental data by Dighe et al. 13 Firstly, the comparison of centre line axial velocity distribution inside of the duct is reported in Figure 7.
Error bars represent the fluctuations in readings from several steady measurements.The velocimetry measurements, especially in the region just upstream and just downstream of the AD, fluctuates because of the flow porous medium.The deviation from the experiment is 1.8% to 2.7% and 3.2 to 7.5% for the upstream and downstream locations, respectively.The reason is likely due to the blockage effect of the wind tunnel, which results in an increased velocity decay rate at the measurement locations. 16Nevertheless, a good agreement is observed between the CFD results and the experiment.Further explanations for the observed deviation and wind tunnel blockage correction methods are comprehensively presented in Dighe et al. 13 Secondly, Figure 8 shows the comparison of pressure coefficients c p measured along the surface of the duct.A detailed explanation of the experimental set-up is presented in Dighe et al. 13 The effect of AD is present at the suction side as a pressure jump at AD location.The offset in the locations of the AD between experiment and computations causes the differences in the location of the pressure jump.The AD position during the experimental set-up is at x/c = 0.28 due to the design limitations.Differently, in the computations, the AD is placed at the nozzle plane (x/c = 0.2) in agreement with DWT theory. 2 It is strongly believed that the AD location influences the solution validation; however, the overall computed c p shows trends in good agreement with the experiment.

DUCT SHAPE PARAMETRIZATION
The duct shapes, in a topological space, represent one-dimensional manifolds whose complete description requires large number of parameters.
The manifolds, however, share peculiar features: cross sections always consist of a closed curve forming a smooth leading edge and a cusped trailing edge.These features are common with airfoil sections, so the design space of all possible duct cross sections can be approximated with airfoil parametrization techniques.
8][19][20] Amongst these, the CST method 19,21 is known to provide nearly complete coverage of the design space. 22Furthermore, it allows progressive design refinement 23 and does not suffer from surface-waviness issues observed in spline based methods. 24This section represents the DonQi ® reference duct with CST parameters and explains the duct shape modifications in the resulting parametric space.Appendix A includes additional information about the method used for the duct shape parametrization.
The CST representation of the duct describes the relative thickness ( = y∕c) of its sides (suction and pressure) as functions of chord-wise distance ( = x∕c) that depend on i = 0 … M shape parameters (A suction i and A pressure i ) : Class functions (C) provide basic topological features, viz, round leading edge and cusped trailing edge, while shape functions (S) represent the specificities of each design.
The parametrization procedure for duct shapes in the current analysis should preserve the following features: leading edge position (which defines the inlet area ratio), trailing edge position (which defines the exit area ratio), and inner side thickness (defines actuator diameter and clearance).This makes it ideal to isolate the effect of camber on the overall performance. 15The DonQi ® duct airfoil (A ref i ) is chosen as the reference shape.In effect of duct shape parametrization, the camber of the DonQi ® duct profile is varied by changing its pressure side while leaving its suction side untouched.First, the two extreme variants of the reference shape corresponding to cambered plate profile (A plt i ) and symmetric profile (A sym i ) are defined, see Figure 9.
The second step consists in obtaining shape parameters for the modified duct profiles (A mod i ) by interpolating between design variants.This is achieved with a quadratic interpolant of a single variable (): ) .
The reference DonQi ® duct is recovered by setting  = 0, while values of  = 1 and  = −1 produce symmetric and thin plate variants, respectively.Equation ( 7) provides parameters for any design with  greater than −1, and duct profile coordinates can be reconstructed from parameters with Equation ( 6).This section provides an extensive analysis of the effects of the duct shape on duct force coefficient, see subsection 6.1, and power coefficient, see subsection 6.2.To this aim, seven duct geometries, shown in Figure 10, are adopted.Finally, some insights on the changes occurring to the performance coefficients are obtained through a detailed flow analysis in subsection 6.3.

Duct force coefficient
Figures 11 and 12 illustrate the correlation between the duct thrust force coefficient C x and the AD loading coefficient C T obtained using the panel and RANS methods, respectively.In general, for lower AD loadings (C T ≤ 0.6), the C x value increases with the increasing C T .Subsequently, a local C x maximum for the individual duct airfoils appear.The value of C x decreases for C T beyond the local maxima.Furthermore, the magnitude of C x increases with the duct cross-sectional camber.A good agreement is achieved between the panel and the RANS solutions, with two noticeable exceptions.Firstly, the magnitude of C x calculated using the panel method is lower than the one calculated with the RANS method across the entire C T range.Secondly, unlike the panel method solutions, the trend lines for duct airfoils D4, D5, and D6 obtained using the RANS method intersect at higher AD loadings (C T ≥ 0.6).In fact, the value of the local maxima attained for duct D6 using RANS calculations is lower than the one attained for duct D5.The main cause of the differences in the two numerical solutions is the viscous interactions between the duct surface and the AD accounted by RANS solutions only.To this aim, some physical insights on the local flow changes due to the duct camber will be discussed in subsection 6.3.Nevertheless, the comparison of the C x predictions by both the computational methods indicate that the non-linear mutual interaction between the duct and the AD is taken into account.

Power coefficient
Figures 13 and 14 represent the power coefficients C P , using the panel and RANS solutions, respectively, as a function of variable AD loading C T .
By analogy, the C P trend appears as a characteristic corollary of the duct force coefficient C x .The larger the C x , the higher the C P reached, and vice versa.Similar to C x solutions, the C P trend lines for duct airfoils D4, D5, and D6 obtained using the RANS method begin to intersect at higher AD loadings (C T ≥ 0.6).In the absence of viscous mutual interaction, as in the panel solutions, the magnitude of C x and ultimately the C P increase with the duct camber.Contrarily, for the RANS solution, the effect of camber on the C P displays an upper limit.Then, the maximum C P ≈ 0.82 is obtained for D5 for C T ≈ 0.75.Moreover, present results show that the maximum C P is obtained for a C T value, which is lower than 8/9, unlike the AMT. 1

Flow field analysis
In order to highlight the effects of the duct shape and the viscous interactions onto the performance coefficients of the duct-AD model, shown in subsections 6.1 and 6.2, a flow field analysis using RANS solutions is carried out.Velocity contours coloured with normalized free stream velocity  The onset of suction-side flow separation, which is earlier for duct D6 than D5, reduces the C x and ultimately the C P of the prescribed duct-AD configuration.

CONCLUSION
The investigation presented in this paper focused towards aerodynamic performance improvement of ducted wind turbines using a simplified duct-AD model.To this aim, numerical calculations using panel method and RANS method are performed.A commercial DWT model DonQi ® is used as the reference case.Unlike the AMT, the two numerical approaches fully takes into account the non-linear mutual interactions between the duct and the AD.To validate the numerical methods, emphasis has been given to the comparison of numerical results with the experiments.
In order to deepen the design principles of DWT, the effects of the duct shape onto the global performance are investigated.An improved duct shape parametrization using CST method is proposed.The study has pointed out that the reference DonQi ® duct is not suited for maximum achievable performance.The analysis has shown the possibility to significantly increase the DWT performance by increasing the duct profile camber and a correct choice of rotor loading.In more detail, the overall performance improvement directly corresponds to the dimensionless duct thrust force coefficient.To better highlight the differences between the RANS and the panel codes, solutions for C x and C P as a function of C T are shown.The RANS solutions, however, exhibit an upper limit of the camber extent for maximum achievable performance.This phenomenon is characterized by a rapid reduction of C x and ultimately the C P , for highly cambered profile.The analysis highlights the limitations of the panel method when applied to highly cambered duct profiles for DWT analysis.A detailed flow analysis using RANS method shows that highly cambered duct profiles are prone to boundary layer flow separation, which is considered nonoptimal for the overall performance.
The duct shape parametrization is suitable for ducts of general shape and is freely available on contacting the authors.

FIGURE 1
FIGURE 1 Schematic of flow around a ducted wind turbine (DWT)

FIGURE 2
FIGURE 2 Pictures showing the experimental set-up with hot-wire anemometer placed on the traversing system (left) and pressure taps along the duct surface (right) [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 3 FIGURE 4
FIGURE 3 Panel distribution along the duct surface and the wake region used for the inviscid panel method calculations

FIGURE 5
FIGURE 5 Computational grid along the leading and trailing edge of the duct

FIGURE 6
FIGURE 6 Grid independence study showing decay of centre line axial velocity measured across the computational domain [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 7
FIGURE 7 Comparison of the numerical results with the experiments.The measurements indicate the axial velocity in the streamwise direction inside of the duct.Error bars indicate the fluctuations in readings from several steady measurements.RANS, Reynold-averaged Navier Stokes [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 8 FIGURE 9
FIGURE 8 Comparison of the numerical surface pressure coefficient c p with the experiments (C T = 0.65).RANS, Reynold-averaged Navier Stokes [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 10
FIGURE 10 Duct profiles used in the numerical study [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 11 FIGURE 12
FIGURE 11 Effect of variable actuator disc (AD) loading on duct thrust force coefficient using panel solution [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 13 FIGURE 14
FIGURE 13 Effect of variable actuator disc (AD) loading on power coefficient using panel solution [Colour figure can be viewed at wileyonlinelibrary.com] are shown in Figure15.Contours are plotted on a plane close to the surface of the duct thus allowing a better interpretation of the flow field associated with duct-AD interactions.In this study, duct airfoils D4, D5, and D6 for C T = 0.6, 0.7, 0.8, and 0.9 are considered to explain the intersecting region in the performance trend lines, as in Figures12 and 14.It is convenient to start with duct profile D4 as a function of C T , more precisely observing from top to bottom.The velocity contours show that, with the increasing value of C T , the leading edge stagnation point traverses further up along the suction side of the duct.Consequently, the magnitude of velocity at the suction side of the duct decreases, and the magnitude of velocity on the pressure side of the duct starts to increase significantly.This phenomenon, occurring at the surface of the duct,

FIGURE 15
FIGURE 15 Velocity contours coloured with normalized streamwise velocity.The results are depicted for ducts D4, D5, and D6 (left to right) bearing C T = 0.6, 0.7, 0.8, and 0.9 (top to bottom) [Colour figure can be viewed at wileyonlinelibrary.com]