Multi-element ducts for ducted wind turbines: a numerical study

Multi-element ducts are used to improve the aerodynamic performance of ducted wind turbines (DWTs). Steady-state, two-dimensional computational fluid dynamics (CFD) simulations are performed for a multi-element duct geometry consisting of a duct and a flap; the goal is to evaluate the effects on the aerodynamic performance of the radial gap length and the deflection angle of the flap. Solutions from inviscid and viscous flow calculations are compared. It is found that increasing the radial gap length results in an augmentation of the total thrust generated by the DWT, whereas a larger deflection angle has an opposite effect. Reasonable to good agreement is seen between the inviscid and viscous flow calculations, except for multi-element duct configurations characterized by large flap deflection angles. The viscous effects become stronger at large flap deflection angles, and the inviscid calculations are incapable of taking this phenomenon into account.

1 Introduction

Ducted wind turbines (DWTs) represent an interesting technological solution for increasing energy extraction with respect to conventional horizontal-axis wind turbines (HAWTs) for a given rotor radius and free-stream velocity (de Vries, 1979). DWTs consist of a rotor and a duct; the role of the latter is to increase the mass flow rate through the rotor relative to a similar rotor operating in the open atmosphere, thereby increasing the generated power. There is more than one explanation for how this occurs. One explanation, as stated by van Bussel (2007), is that the duct forces an expansion of flow downstream of the turbine beyond what is attainable for a bare wind turbine. This provides reduced pressure on the downstream of the turbine, thereby increasing the total mass flow through the turbine. A second explanation, as argued by de Vries (1979), is that if the sectional lift force of the duct is directed towards the turbine plane, then the associated circulation of the duct induces an increased mass flow through the turbine. This solution is particularly suited for urban areas where the radius of the rotor is a constraint and the free-stream velocity is low due to the presence of buildings. van Bussel (2007), using a one-dimensional momentum theory approach, found that the maximum power coefficient obtained by a DWT can exceed the Betz limit by a factor of 2.5. The best aerodynamic performance for a DWT can be achieved by increasing the duct expansion ratio (Liley and Rainbird, 1956; Foreman et al., 1978; Loeffler and Vanderbilt, 1978; Foreman and Gilbert, 1984; Samson and Katebi, 2014). As a drawback, for a duct with a large ratio of duct outlet to rotor area, flow separation along the duct inner walls might be present (Aranake et al., 2015; Dighe et al., 2018a). Tang et al. (2018) experimentally investigated the effects of variable duct expansion ratios on the aerodynamic performance of DWTs. They found that increasing the expansion ratio over a certain limit results in a power reduction. This was linked to the appearance of flow separation within the inner walls of the duct. An alternative approach to improve the aerodynamic performance is to increase the camber of the airfoil, used as a cross section of the duct, until separation occurs behind the rotor plane (Dighe et al., 2019). Since, as expected, flow separation has an undesired effect, solutions to prevent flow separation via active boundary layer control techniques have been proposed (Igra, 1977; Foreman et al., 1978; Abe and Ohya, 2004). However, the performance benefits were limited by the cost of the active system and its installation.

Another possible solution to improve the aerodynamic performances of DWTs, which has only been explored to a very limited extent, consists of using a duct with a flap (i.e., a multi-element duct). The flap is realized as a secondary duct with a small chord airfoil cross section mimicking high-lift devices for aircraft wing (see Fig. 1). A first theoretical and experimental analysis of DWTs with a flap was carried out by Foreman et al. (1978) and Igra (1981). The latter found that the addition of a flap improves the DWT aerodynamic performance by 25 % with respect to a single duct. The flap inhibits flow separation along the inner duct wall and increases the camber of the equivalent airfoil, thus being beneficial for the aerodynamic performances of the DWT (Dighe et al., 2019). The literature, however, is missing a detailed parametric study that investigates the effect of flap installation setting, i.e., the radial location and its angle of attack, on the total power generated by a DWT.

Figure 1 Schematic of flow around a multi-element ducted wind turbine.

The goal of this paper is to conduct a parametric study to investigate the effect of the installation settings of the flap on the aerodynamic performance of a multi-element DWT. This is performed using computational fluid dynamics (CFD). To this aim, a reference multi-element duct is selected and the rotor is simulated by a uniformly loaded actuator disk (AD) model. The rest of the paper is organized as follows. Section 2 reports the most relevant nondimensional coefficients adopted for characterizing the multi-element duct–AD model. Section 3 details the verification and validation of the numerical methods against experimental findings. Section 4 describes the computational settings and parameters with a brief description of the numerical methodology. The dependence of the nondimensional coefficient on the flap installation settings is discussed in Sect. 5. Finally, the most relevant results are summarized in the conclusions.

2 Multi-element duct–AD flow model

The incompressible flow past a wind turbine is computed by substituting the rotor with an AD of infinitesimal width. The AD exerts a uniform thrust force T_AD per unit area. Then, the nondimensional thrust force coefficient is

$$\begin{array}{}\text{(1)}& C_{\mathrm{T},\mathrm{AD}}={\displaystyle \frac{T_{\mathrm{AD}}}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{U}_{\mathrm{\infty }}^{\mathrm{2}}{S}_{\mathrm{AD}}}},\end{array}$$

where ρ is the fluid density, U_∞ is the free-stream velocity and S_AD is the AD area.

T_AD is obtained by forcing a uniform pressure drop across the AD, $T_{\mathrm{AD}}=\mathrm{\Delta }p\times S_{\mathrm{AD}}$. The pressure drop Δp is taken from experiments and is given as an input parameter to the numerical simulations. The mean velocity across the AD, U_AD, is obtained by integrating the difference of the streamwise velocity component across the AD surface U_x:

$$\begin{array}{}\text{(2)}& U_{\mathrm{AD}}={\displaystyle \frac{\mathrm{1}}{S_{\mathrm{AD}}}}\underset{S_{\mathrm{AD}}}{\oint }{U}_x\mathrm{d}S.\end{array}$$

Then, the power coefficient is

$$\begin{array}{}\text{(3)}& C_{P_{\mathrm{o}}}={\displaystyle \frac{P_{\mathrm{o}}}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{U}_{\mathrm{\infty }}^{\mathrm{3}}{S}_{\mathrm{AD}}}}={\displaystyle \frac{U_{\mathrm{AD}}}{U_{\mathrm{\infty }}}}{C}_{\mathrm{T},\mathrm{AD}}.\end{array}$$

For a multi-element duct–AD configuration, additional thrust forces are exerted by the duct and the flap. Then, the total thrust force T is the vectorial sum of the AD thrust force T_AD and the axial thrust force exerted by the duct T_D and the flap T_F. It can be written as

$$\begin{array}{}\text{(4)}& T=T_{\mathrm{AD}}+T_{\mathrm{D}}+T_{\mathrm{F}}=T_{\mathrm{AD}}+T_{\mathrm{M}}.\end{array}$$

The total thrust coefficient is then defined as

$$\begin{array}{}\text{(5)}& C_{\mathrm{T}}=C_{\mathrm{T},\mathrm{AD}}+C_{\mathrm{T},\mathrm{M}},\end{array}$$

where C_T,M is the multi-element duct thrust coefficient.

To highlight the relative contribution of the multi-element duct thrust T_M and the AD thrust T_AD to the total thrust T, a dimensionless thrust factor τ is introduced (Bontempo and Manna, 2016):

$$\begin{array}{}\text{(6)}& \mathit{\tau }={\displaystyle \frac{T_{\mathrm{M}}}{T_{\mathrm{AD}}}}={\displaystyle \frac{C_{\mathrm{T},\mathrm{M}}}{C_{\mathrm{T},\mathrm{AD}}}},\end{array}$$

so that the total thrust coefficient can be written as

$$\begin{array}{}\text{(7)}& C_{\mathrm{T}}=(\mathrm{1}+\mathit{\tau })C_{\mathrm{T},\mathrm{AD}}.\end{array}$$

Following Bontempo and Manna (2016), the normalized axial velocity at the AD for a ducted configuration can be also expressed as a function of the thrust coefficient:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{U_{\mathrm{AD}}}{U_{\mathrm{\infty }}}}={\displaystyle \frac{\mathrm{1}+\mathit{\tau }}{\mathrm{2}}}\left(\mathrm{1}+\sqrt{\mathrm{1}-C_{\mathrm{T},\mathrm{AD}}}\right).\end{array}$$

Using Eqs. (3) and (8), the power coefficient of the multi-element duct–AD model considering S_AD as the reference area can be written as

$$\begin{array}{}\text{(9)}& C_{\mathrm{P}}={\displaystyle \frac{\mathrm{1}+\mathit{\tau }}{\mathrm{2}}}\left(\mathrm{1}+\sqrt{\mathrm{1}-C_{\mathrm{T},\mathrm{AD}}}\right)C_{\mathrm{T},\mathrm{AD}}.\end{array}$$

In Eq. (9), C_P indicates the power coefficient of the multi-element duct–AD model. The above relation is also valid for a simple AD model setting τ=0:

$$\begin{array}{}\text{(10)}& C_{P_{\mathrm{o}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}+\sqrt{\mathrm{1}-C_{\mathrm{T},\mathrm{AD}}}\right)C_{\mathrm{T},\mathrm{AD}}.\end{array}$$

Equations (9) and (10) can be used to evaluate the contribution of the multi-element duct through a power augmentation parameter r:

$$\begin{array}{}\text{(11)}& r={\displaystyle \frac{C_{\mathrm{P}}}{C_{P_{\mathrm{o}}}}}=\mathrm{1}+\mathit{\tau }=\mathrm{1}+{\displaystyle \frac{C_{\mathrm{T},\mathrm{M}}}{C_{\mathrm{T},\mathrm{AD}}}}.\end{array}$$

Equation (11) states that r for a multi-element duct–AD model is proportional to the ratio between the multi-element duct thrust coefficient C_T,M and the AD thrust coefficient C_T,AD. Thus, if τ>0, then a higher power coefficient can be obtained for a multi-element DWT in comparison to a HAWT with the same rotor. The performance coefficients described above were evaluated by means of axial momentum theory (AMT) for DWTs by van Bussel (2007). However, the AMT cannot be used to estimate the performance of the duct–AD model for a prescribed C_T,AD and a given duct geometry. This problem can be solved using numerical solutions based on panel and RANS (Reynolds-averaged Navier–Stokes) methods (Bontempo and Manna, 2016; Dighe et al., 2019).

3 Numerical validation

For validating the numerical methods, experimental data reported by Igra (1977) were simulated. Igra's experiments were conducted in the subsonic wind tunnel of Israel Aerospace Industries (formerly Israel Aircraft Industry); this tunnel has a large test section that measures 3.6 m × 2.6 m.

Eight different geometries were investigated experimentally, but only two geometries are used for the validation study. The two geometries are a duct–AD model with $C_{\mathrm{T},\mathrm{AD}}=\mathrm{0.434}$ (Model B) and a multi-element duct–AD model with $C_{\mathrm{T},\mathrm{AD}}=\mathrm{0.550}$ (Model C (ii) + flap). A schematic of the cross sections of the two geometries is shown in Fig. 2. The longitudinal cross section of the duct and of the flap is a NACA 4412 airfoil. The leading edges of both duct geometries are identical. For Model C (ii) + flap, the trailing edge of the duct is radially stretched, resulting in a duct expansion ratio $\frac{A_{\mathrm{e}}}{A_{\mathrm{AD}}}=\mathrm{1.84}$; this ratio is 1.71 for Model B. The flap chord measures 35 % of the duct chord length, c, and the deflection angle θ=30^∘ with respect to the free-stream direction. The experimental dataset consists of static pressure distribution at different axial and radial positions, as well as forces generated by the duct and the flap surfaces. During the experiments, the inflow velocity was set at U_∞=32 m s⁻¹, corresponding to $Re\approx \mathrm{4.5}\times {\mathrm{10}}^{\mathrm{5}}$. Following Igra (1977), the wall interference and blockage correction can be ignored.

Figure 2 A schematic cross-sectional layout of the three-dimensional experimental model used for the numerical validation study.

For numerical validation, three numerical methods are considered: a two-dimensional panel method based on Euler equations, a two-dimensional steady RANS method based on time-averaged RANS equations and a two-dimensional unsteady RANS method based on time-filtered URANS equations. A detailed description of the differences in the governing equations is beyond the scope of the current discussion. The reader may refer to Versteeg and Malalasekera (2007). For both steady and unsteady RANS methods, the k−ω SST (shear stress transport) model is used as a turbulence model. Preliminary investigations showed good agreement with the experiments (Dighe et al., 2018a, 2019). A numerical stability analysis is performed in the context of time-filtered URANS simulations. The physical time step corresponding to a Courant–Friedrichs–Lewy (CFL) number of 1 in the finest mesh resolution level is 1.67×10⁻⁵ s. More details on the grid settings and the boundary conditions will be presented in Sect. 4

In Fig. 3, the power augmentation factor r is plotted as a function of the inflow yaw angle α. CFD results, obtained using panel, steady RANS and unsteady RANS methods, are compared with the experimental data. A very good agreement between the CFD simulations and the experimental findings is found for Model B. On the other hand, the deviation between the CFD and the experimental findings is larger for Model C (ii) + flap, in particular for α≠0^∘. The discrepancies might be due to three-dimensional effects not accounted for in the two-dimensional simulations. At non-yawed conditions (α=0^∘), however, the CFD results agree well with the experimental findings, wherein the maximum deviation of the CFD results from the experimental data is 1.2 % and 4.3 % for Model B and Model C (ii) + flap, respectively.

Figure 3 Effect of yawed inflow on the power augmentation factor: comparison between the experiments (Igra, 1977), panel method and RANS method for different duct configurations.

The differences in the CFD results can be explained by looking at the flow field. Figure 4 shows the contours of nondimensional axial velocity $\frac{U_x}{U_{\mathrm{\infty }}}$ computed with the three numerical methods for Model C (ii) + flap. The velocity contour from the panel method is plotted on the left; steady RANS is in the center, and unsteady RANS is on the right. A clear difference in the flow field between the panel and the RANS (steady and unsteady) methods could be identified from the contour plots. Neglecting viscosity, as in the panel method solution, the flow remains attached over the suction side of the duct. As a result, the magnitude of velocity on the suction side, and ultimately the value of r (Fig. 3), is larger for the panel method solution in comparison to the RANS (steady and unsteady) method solutions. On the other hand, the flow fields obtained using steady and unsteady RANS methods are almost identical. Both the steady and unsteady RANS solutions show flow separation along the inner walls of the flap. Subtle differences appear in the flow separation region, where the velocity contour patterns differ slightly in their spatial organization. In the URANS solution, the turbulent flow structures, which evolve in time and space, are explicitly computed. These flow structures are temporally averaged in the RANS solution. The net result of such different formulations explains the difference in the value of r in Fig. 3 calculated using steady and unsteady RANS methods.

Figure 4 Velocity contours colored with normalized free-stream velocity obtained using the (a) panel method, (b) RANS method and (c) URANS method for Model C (ii) + flap.

Although URANS simulations increase the level of description of the unsteady flow due to the multi-element duct–AD interaction, the computing cost incurred by going from RANS to URANS does not justify the scope of the current study, in which the effects of distributed AD loading, wake rotation, divergence and inflow yaw angle are totally ignored.

4 Numerical approach

In this section, the numerical methods employed will be briefly described. For an in-depth description, the reader can refer to Dighe et al. (2019).

4.1 Panel method

A two-dimensional potential flow panel method has been used to compute the steady iso-entropic incompressible flow field around the multi-element duct–AD model following the work of de Oliveira et al. (2016). The governing flow equations are a simple form of the Euler equations. The AD is represented by a pair of symmetric counter-rotating vortices. The duct and the flap geometries are defined using a distribution of vortices located on the panels to reproduce the desired cross-sectional shape. A uniform distribution of vorticity on the panels is assigned by assuming the Kutta condition. The assumption of uniform vorticity distribution over the panels represents a simplification of real physics and prevents flow separation on the multi-element duct surface, even for larger pressure gradients. The duct and the flap surface discretization is based on the constant spacing approach. The streamwise discretization in the near and far wake is nonuniform, with initial panel length equal to 1.0 % of duct chord length c, just behind the AD, and increasing gradually in length as the wake expansion settles further downstream (see Fig. 5).

The panel method is particularly appealing for routine design analysis due to its short execution time. A typical converged panel method solution is obtained in roughly 0.05 h on a multicore workstation desktop computer.

Figure 5 Panel distribution along the duct surface and the wake region used for the inviscid panel method calculations.

4.2 RANS method

A commercial CFD solver, ANSYS Workbench^®, has been used for a complete viscous solution of steady incompressible flow around the multi-element duct–AD model. The governing flow equations are the Reynolds-averaged Navier–Stokes (RANS) equations. The 2-D computational domain is shown in Fig. 6, where the distances from the AD location to the domain inlet and outlet are 12c and 24c, respectively. The computational grid consists of quadrilateral cells with a minimum y⁺ value of 1 on the duct walls. Boundary conditions are a uniform velocity at the inlet, zero gauge static pressure at the outlet, and no-slip walls for duct and flap surfaces. A symmetric boundary condition is applied along the center-line axis, while a fan boundary condition is used for the AD. Experimental data from the wind tunnel (Dighe et al., 2018b) in the form of pressure drop against velocity normal to the AD plane were extrapolated to determine the input parameter for the fan boundary condition. RANS solutions are considerably more reliable and accurate than the panel method solution but at the expense of computational cost. A typical converged RANS solution with approximately 0.1 million mesh elements is obtained in roughly 0.5 h on a multicore workstation desktop computer.

Figure 6 Computational domain showing the boundary conditions employed. The length is indicated in terms of duct chord length c (representative, not to scale). The computational grid along the duct and flap leading edge is zoomed in to show the C-grid and boundary layer refinement.

RANS solutions are sensitive to the discretization of the computational domain. For the present computations, a C-grid structured zonal approach is chosen (Fig. 6), which proved advantageous in the case of a curved boundary (duct and flap leading edge). The C-shaped loop terminates in the wake region. Grid independence analysis has been carried out using three grid sizes; the refinement factor in each direction is approximately 1.5. The refinement factor is defined as the rate at which the grid size increases far from the object.

The multi-element duct thrust force coefficient, C_T,M, is taken as a reference for the convergence analysis. The results of the grid independence study are shown in Table 1. Convergence is reached for the medium grid, which is then used in the rest of the paper.

Table 1 Grid statistics for the grid independence study of the reference case.

5 Results and discussion

5.1 Multi-element duct geometry

In the following sections, the effects of flap installation settings on the aerodynamic performance of the multi-element duct–AD model are described. The multi-element duct–AD configuration investigated in the present work is shown in Fig. 7. The longitudinal cross section of the duct is a DonQi D5 airfoil; the profile is chosen based on the duct shape parametrization study conducted by the authors (Dighe et al., 2019). For the DonQi D5 duct an optimal $C_{\mathrm{T},\mathrm{AD}}=\mathrm{0.7}$ was obtained. This value is employed throughout the present discussion. A NACA 4412 longitudinal cross section, measuring 0.35c, is chosen for the flap following Igra (1977). The flap installation settings are the radial gap ζ and the deflection angle θ. The radial gap ζ, indicated as a percentage of duct chord length c, is defined as the distance from the trailing edge of the duct to the leading edge of the flap. A positive value of radial gap (ζ>0) indicates that the leading edge of the flap is positioned below the trailing edge of the duct. A positive deflection angle (θ>0) corresponds to a downward flap deflection, whereby the angle is defined relative to the free-stream direction. The axial gap between the trailing edge of the duct to the leading edge of the flap is zero based on the findings of Igra (1981). The numerical study is performed at a fixed Re of 4.5×10⁵ as in the experiments. In Sect. 5.2 and 5.3, the changes in the aerodynamic performance coefficients with respect to the flap's geometric orientation are quantified.

Figure 7 A schematic cross section of the multi-element duct–AD model with the variable flap parameters used for the flap installation study.

5.2 Duct force coefficient

Contours of the multi-element duct force coefficient C_T,M, obtained from panel and RANS methods, are shown in Figs. 8 and 9, where C_T,M as a function of the radial gap ζ and the deflection angle θ are reported. The figures show that C_T,M increases for larger ζ. Conversely, C_T,M decreases with increasing θ. The maximum C_T,M obtained from both the numerical methods lies in the same region, i.e., ζ≈5 % and θ≈10^∘. The differences between results obtained using the panel and RANS methods are highly contrasting for θ≥60^∘.

Figure 8 Effect of the variable radial gap and deflection angle of the flap on the duct thrust force coefficient using the panel method.

The differences can be explained by looking at the flow field. Contours of nondimensional axial velocity $\frac{U_x}{U_{\mathrm{\infty }}}$ from both methods are reported in Fig. 10a–f. Results from the panel method are plotted on the left, while the ones from RANS are on the right. Contours for the no-flap configuration are shown in Fig. 10a and b. Two flap settings, in order to explain the aerodynamics behind Figs. 8 and 9, are shown: ζ=5 % and θ=10^∘ in Fig. 10c and d as well as ζ=5 % and θ=70^∘ in Fig. 10e and f.

Figure 9 Effect of the variable radial gap and deflection angle of the flap on the duct thrust force coefficient using the RANS method.

Contour plots show a higher velocity at the rotor plane for the configuration with flap in comparison with the no-flap configuration. This is due to the additional aerodynamic thrust force generated by the flap. The presence of a radial gap between the duct and the flap accelerates the flow over the flap. This reduces the pressure recovery demands on the multi-element duct, thereby reducing flow separation. Obviously, flow separation is seen for RANS contours only. The overall integral contribution of the viscous forces increases the C_T,M magnitude in the RANS solutions relative to the panel solutions, a trend that can be clearly observed by comparing Figs. 8 and 9. For the flap configuration with ζ=5 % and θ=70^∘ (Fig. 10f), the flow over the flap's inner walls separates completely. The separation along the inner walls of the multi-element duct reduces the C_T,M, which rapidly becomes large and negative at higher flap deflection angles as seen in Fig. 9. For panel solutions, however, the drop in the C_T,M magnitude for higher flap deflection angles is gradual (see Fig. 8) because viscous effects are neglected (see Fig. 10e).

Figure 10 Velocity contours colored with normalized free-stream velocity obtained using the (a) panel method with no flap, (b) RANS method with no flap, (c) panel method with ζ=5 % and θ=10^∘, (d) RANS method with ζ=5 % and θ=10^∘, (e) panel method with ζ=5 % and θ=70^∘, and (f) RANS method with ζ=5 % and θ=70^∘.

5.3 Power augmentation

Figures 11 and 12 represent contours of power augmentation factor r using the panel and RANS solutions, respectively, as a function of radial gap ζ and deflection angle θ. Recall from Eq. (11) that an r gain for a multi-element duct–AD model can be attained by increasing the C_T,M magnitude for a constant C_T,AD. Evidence of this is provided in Figs. 11 and 12, which exhibit the r maximum in the same region of C_T,M maximum, as in Figs. 8 and 9. Then, the maximum power augmentation factor r_max=1.25 and 1.38, obtained for panel and RANS calculations, respectively, corresponds to ζ≈5 % and θ≈10^∘.

Figure 11 Effect of variable radial gap and the flap deflection angle on the power augmentation factor using the panel method.

Figure 12 Effect of variable radial gap and the flap deflection angle on the power augmentation factor using the RANS method.

A comparison of r_max obtained from CFD simulations with the one-dimensional axial momentum theory (van Bussel, 2007) is carried out. The analysis of AMT is not reported herein for the sake of brevity; a description of the AMT theory applied is reported in Appendix . The velocity and pressure values within a DWT, using the AMT approach, are determined by applying the relations of duct expansion ratio β and back-pressure velocity ratio γ. With the AMT approach, the maximum power augmentation factor of r_max=1.50 for ζ≈5 % and θ≈10^∘ is obtained. In comparison to the CFD methods, the AMT approach does not take into account the effects of multi-element duct geometry and the nonlinear mutual interaction between the multi-element duct and the AD. As a result, the values of C_P and the power augmentation factor r, determined using the AMT approach, are overestimated in comparison to the CFD methods shown in the above discussion.

6 Conclusions

In this work, the aerodynamic performance of a multi-element DWT is studied using a numerical approach. To this aim, two-dimensional numerical calculations using the panel method and the RANS method are employed. A simplified AD model is used for simulating the rotor. Based on the existing studies conducted by the authors, the multi-element duct geometry consists of a DonQi D5 airfoil and a NACA 4415 airfoil for the duct and the flap cross sections, respectively. To validate the numerical methods, the present simulations are compared with similar experimental data. In order to deepen the design principles of multi-element ducts, the effects of radial gap ζ and the flap deflection angle θ on the global performance of DWTs are investigated. Clear trends of the multi-element duct thrust force coefficients C_T,M and the power augmentation factor r are observed across a range of multi-element duct configurations. An increase in the flap deflection angle θ results in a decrease in C_T,M, whereas an increase in the radial gap ζ shows an increase in C_T,M. The analysis of flow field shows that flow separation in the multi-element duct inner walls increases for higher values of θ. This phenomenon determines the reduction in C_T,M and ultimately the augmentation factor r. As expected, the RANS method is more suitable for representing solutions for highly deflected flap configurations. The viscous effects become stronger at higher flap deflection angles, and the panel method is inherently incapable of taking this into account. Regarding prediction of the near-optimal multi-element duct configuration, both the numerical methods show good agreement. Moreover, in comparison to the AMT approach, the CFD method fully takes into account the mutual interactions between the multi-element duct and the AD.

Appendix A: r_max calculation using AMT

The velocity and pressure relations within a DWT are obtained using the AMT approach (see Fig. A1). Based on the geometry under consideration, the duct expansion ratio β=2.25 is established and the AD inside the duct operates at $a=\frac{\mathrm{1}}{\mathrm{3}}$. V₀=6 m s⁻¹, corresponding to the numerical study reported above. Having said that, the velocity at the duct exit equals

$$\begin{array}{}\text{(A1)}& V_{\mathrm{3}}=(\mathrm{1}-a)V_{\mathrm{0}}=\mathrm{4}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}.\end{array}$$

Using the diffuser area ratio β, the velocity at the AD location equals

$$\begin{array}{}\text{(A2)}& V_{\mathrm{1}}=\mathit{\beta }(\mathrm{1}-a)V_{\mathrm{0}}=\mathrm{9}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}},\end{array}$$

and the value of the back-pressure velocity ratio γ yields

$$\begin{array}{}\text{(A3)}& \mathit{\gamma }={\displaystyle \frac{V_{\mathrm{3}}}{V_{\mathrm{0}}}}=\mathrm{0.67}.\end{array}$$

Finally, the power coefficient considering AD surface area returns

$$\begin{array}{}\text{(A4)}& C_{\mathrm{P}}=\mathit{\beta }\mathit{\gamma }{\displaystyle \frac{\mathrm{16}}{\mathrm{27}}}=\mathrm{0.89},\end{array}$$

where $\frac{\mathrm{16}}{\mathrm{27}}$ is the theoretical limit for maximum achievable performance for a bare wind turbine, most commonly addressed as the Betz limit. Then, the maximum power augmentation factor r_max equals

$$\begin{array}{}\text{(A5)}& r_{\max }={\displaystyle \frac{C_{\mathrm{P}}}{C_{P_{\mathrm{Betz}}}}}=\mathrm{1.50}.\end{array}$$

Figure A1 Pressure and velocity relations in a DWT using AMT.