Formulation, Validation, and Application of a Novel 3D BEM Tool for Vertical Axis Wind Turbines of General Shape and Size

: Low order models based on the Blade Element Momentum (BEM) theory exhibit modeling issues in the performance prediction of Vertical Axis Wind Turbines (VAWT) compared to Computational Fluid Dynamics, despite the widespread engineering practice of such methods. The present study shows that the capability of BEM codes applied to VAWTs can be greatly improved by implementing a novel three-dimensional set of high-order corrections and demonstrates this by comparing the BEM predictions against wind-tunnel experiments conducted on three small-scale VAWT models featuring different rotor design (H-shaped and Troposkein), blade proﬁle (NACA0021 and DU-06-W200), and Reynolds number (from 0.8 × 10 5 to 2.5 × 10 5 ). Though based on the conventional Double Multiple Stream Tube (DMST) model, the here-presented in-house BEM code incorporates several two-dimensional and three-dimensional corrections including: accurate extended polar data, ﬂow curvature, dynamic stall, a spanwise-distributed formulation of the tip losses, a fully 3D approach in the modeling of rotors featuring general shape (such as but not only, the Troposkein one), and accounting for the passive effects of supporting struts and pole. The detailed comparison with experimental data of the same models, tested in the large-scale wind tunnel of the Politecnico di Milano, suggests the very good predictive capability of the code in terms of power exchange, torque coefﬁcient, and loads, on both time-mean and time-resolved basis. The peculiar formulation of the code allows including in a straightforward way the usual spanwise non-uniformity of the incoming wind and the effects of skew, thus allowing predicting the turbine operation in a realistic open-ﬁeld in presence of the environmental boundary layer. A systematic study on the operation of VAWTs in multiple environments, such as in coastal regions or off-shore, and highlighting the sensitivity of VAWT performance to blade proﬁle selection, rotor shape and size, wind shear, and rotor tilt concludes the paper.


Introduction
The history of the Vertical Axis Wind Turbine technology is millennial. The very first machines, conceived to extract energy from the wind thousands of years ago, featured the axis of rotation normal to the wind direction and they were based on drag (see, e.g, the Panemone turbine). At the beginning of the 20th Century, two novel VAWT concepts emerged, i.e., the Savonious and the Darrieus turbines, which managed to exploit (marginally the former, mainly the latter) the novel idea of lift force. At the time of the patent presentation, the Darrieus concept promised to outperform any other wind energy conversion device available so far. The Darrieus turbine was later object of intense research, mainly carried out by the Sandia Lab group [1][2][3], in the second half of the 20th Century; however, when modern lift-driven wind turbines evolved from concepts to technology, VAWT eventually lost the competition against Horizontal Axis Wind Turbine (HAWT).
A comprehensive picture of the VAWT technology at the beginning of this Century is depicted in the book of Paraschivoiu [4]. This book clarifies the motivations for the diffi-cult development of VAWT technology, that lie in the flow complexity and aerodynamic forcing of a machine which features unsteady operation by design. However, the text also highlights several aspects that might make it advantageous with respect to HAWTs: insensitivity to yaw angle, low sensitivity to gust and skew, lower noise (due to lower optimal tip-speed ratios), the possibility to mount gearbox, electric generator and other systems on the ground. In the last decade, the development of VAWT technology has seen a significant acceleration, and recent works have enriched the experimental database on prototypes at laboratory scale [5][6][7][8][9][10][11][12][13], have improved the modeling techniques [14][15][16][17][18][19][20][21][22][23], and have investigated issues related to the machine operation [21,[24][25][26][27][28]. This impulse has improved the VAWT concept, making it progressively competitive against the conventional HAWT.
These considerations have led to reconsider the VAWT as a promising technology for the efficient energy harvesting in non-conventional environments, especially for smallscale urban/distributed installations. The same arguments, however, apply also in another, extremely relevant application branch, e.g., Floating Offshore Wind Turbines (FOWT), which experience a significant inflow complexity due to the movement of the platform. FOWTs operate under time-periodic skew and yaw conditions and, in the case of HAWTs, have to cope with large overturning moments produced by the heaviest components of the turbine, such as the drivetrain and the generator, which have to be placed high above the sea level. In such a context, the VAWT may be better suited to the emerging floating offshore market, as heavy components would all sit at the base of the turbine, and the rotor aerodynamics is less sensitive to wind skew and yaw with respect to HAWTs. VAWTs could thus potentially enhance the turbine performance and, at the same time, reduce the cost of FOWTs installation and maintenance, owing to the absence of high-speed shafts, yaw systems and nacelles which are all potentially subject to faults.
In light of the renewed interest in VAWT technology in such emerging contexts, the development of accurate and effective tools for the prediction of VAWT performance and load are now of interest. The needs of future applications will involve to consider general rotor shapes, with special focus on the Troposkein shape (and its modifications) and considering the wind shear and skew which might be relevant, respectively, for mini-scale distributed applications and large-scale floating applications.
To satisfy aforementioned requests, recent literature offers comprehensive BEM models, employing the Double Multiple Stream Tube approach, improved by analytical corrections as in [18,29], implementing different models to correct incidence due to flow curvature, dynamic stall effects and streamtube expansion for H-shaped rotors. The main differences between the study by Bianchini et al. and present 3D BEM code are that the former employs Prandtl-Lanchester correction to estimate performance losses at tips, the latter adopts Prandtl's tip losses correction adapted from HAWTs; regarding the dynamic stall model, results in [18] were obtained with Berg's model whereas this paper entails Strickland's one. Both BEM codes envisage an extended polar database for accurate post-and deep-stall aerodynamic coefficients with smoothing and model strut and pole passive torque. Flow curvature is accounted for by Migliore's virtual mapping method adding virtual incidence concept; only a simplistic correction by Goude is chosen to correct the angles of attack in this work. On this point, Bangga et al. [29] developed a novel simple method to consider the 'decambering effect' of airfoil in cycloidal motion and compared results with already existent methods; the dynamic stall used is the Beddoes-Leishman and streamtube expansion model; however, passive torque corrections were not considered.
The present paper documents the formulation, the validation and multiple applications of a novel BEM tool conceived to deal with 3D in-flow configurations coupled with 3D rotor architecture, and featuring 3D corrections to the base model considering the tip and strut loss contributions. The paper is structured as follows. At first the three exemplary VAWTs configurations used for validation and application purposes throughout the study are presented, then the novel BEM formulation model is described and its experimental validation is performed by resorting to wind-tunnel tests of the aforementioned VAWT prototypes. Then the BEM tool is applied to multiple case studies, at first comparing different rotor architectures, rotor scales and blade profiles, and then considering exemplary on-shore installations in coastal regions and off-shore installations of the floating DeepWind turbine.

Case Study: Turbine Models and Experiments
The present computational study was made possible by the availability of three laboratory-scales models of VAWTs of different shape and size, which provided the fundamental benchmark for experimental validation and were used as references to construct the application case-studies. In this section, these turbine configurations are described.

Micro-Generation VAWTs
The first class of VAWTs used in this work is composed by two prototypes representative of real-scale turbines for micro-generation (P max < 1 kW) featuring the same cross-section and the same blade profile but different rotor shape, namely H-shape rotor and Troposkein. The straight H-shape VAWT consists of three unstaggered symmetric NACA0021 blades, each one connected to the hub through two spokes. Representation of the turbine configuration along with the coordinate system and triangle of attack notation are depicted in Figure 1a and Figure 1b. Troposkein rotor exhibits hub-to-hub blades, without supporting struts and hence offering a lower aerodynamic drag with respect to the H-shape design (at least for small-scale VAWTs). The main geometrical characteristics of the considered turbines are reported in Tables 1 and 2. The rotor solidity reported in Table 2 is calculated according to Sandia Laboratories definition [1], resulting higher than that of the H-shape rotor (even though the solidity evaluated locally at the equatorial plane is comparable to that of the H-shape rotor).  The performance and wake shed by the turbines were experimentally characterized in the large-scale wind tunnel GVPM (Grande galleria del Vento del Politecnico di Milano [30]) in free-jet configuration, featured by a square cross sectional area of about 16 m 2 with the rotor placed in the center of the jet. The application of dedicated correlations [31] indicates a negligible free-blockage. The incoming flow features very low turbulence intensity (below 1%). Performance and load measurements were obtained by combining angular speed, torque and force measurements on the shaft of the turbine; velocity and turbulence measurements in the wake were also performed by traversing multiple hot wires, even though they are not used in this study due to the inherent limitations of BEM methods in reproducing the wake [18]. Full details on the geometrical features, on the installation in the GVPM, and on the performance/load data of the two turbines can be found in [12].

Deepwind Troposkein Prototype
The second class of VAWT considered in this work is based on the DeepWind concept, which is a floating VAWT concept for large scale applications (rated power from 5 MW to 20 MW). A laboratory scale demonstrator (1 kW) of the DeepWind turbine was tested in the GVPM for different wind speeds and rotor rotational velocities, also including skewed flow operation due to a 15 • -tilted rotor arrangement (representative of the expected inclination of the tower during floating operation). All the main geometric features are reported in Table 3; DeepWind turbine blades feature the DU06W200 profile, which was specifically designed for VAWT applications, relying on good aerodynamic performance as well as mechanical properties [32]. Full details on measuring campaign and rotor geometry are reported in [33].

3D BEM-DMST Model: Basic Formulation and High-Order Submodels
The in-house BEM tool presented in this work was conceived to assess the performance of VAWTs featuring any general shape of the rotor, employing the most recent improvements available in the literature, and with the aim of reproducing the effects of 3D flow around VAWTs. The code is based on the steady-state approach of the Blade Element Momentum Theory, reformulated according to the so-called Double Multiple Stream Tube (DMST) model inherited by Paraschivoiu [4]; the variable induction factor formulation improves the level of local load prediction accuracy by considering upwind and downwind induced velocities as a function of the rotation angle of each blade. The chosen reference frame is depicted in Figure 1. The basic formulation is enriched by sub-models and corrections making the DMST able to account for several 2D and 3D phenomena affecting VAWTs aerodynamics and performance. The following sub-sections present in detail all the corrections, which are now briefly listed.
First considering 2D flow features, the unsteady operation of the VAWT blades induces large variations in angle of attack on the airfoils with high peak values (both positive and negative); the operation at large angle of attack might make the blade to operate in the post stall region, for part of the revolution, requiring the extension of profile polars beyond the static stall limits. Beside this, the change of relative velocity along the chord produced by the blade cycloidal motion, requires the introduction of a so-called "flow curvature correction" [34]. The inherent unsteadiness of VAWTs makes the aforementioned incidence variation to occur in relatively short time scales and partly within the post stall region, resulting in a hysteresis loop of lift and drag forces, phenomenon known as "dynamic stall". Semi-empirical models have been treated in the literature and show a relative ease of implementation in DMST codes such as the one proposed by Gormont [35] for helicopter applications and its modifications by Berg [36], Strickland [2] and Paraschivoiu [37]. Further studies brought to indicial approaches, e.g., the most renowned proposed by Beddoes-Leishman [38]. These physical-based models consider distinct contributions for the unsteady aerodynamics, leading edge and trailing edge flow separation. The implementation of these methods still relies on the availability of empirical constants which are Reynolds-dependent: literature offers a wide range of data at high Reynolds regimes that infrequently find space in small scale VAWT applications featured by typical low Reynolds number levels, for that a lack of experimental tuning of the empirical constants can be shown.
Focusing on 3D flow features, we first note that the finite length of the blade leads to the spanwise change in the aerodynamic forces (and hence in the vorticity) and eventually results in the formation of the tip vortex, so a model has to be adopted to consider the effect, either in an integral or in a distributed way. Moreover, VAWTs comprise supporting structures such as the spokes and pole which contribute to power losses and are modelled as a substantial drag augmentation affecting the torque harvested by the rotor. The effects of struts are evaluated as function of the azimuthal angle. Losses caused by wake-blade interaction still remain at an higher order of resolution, unaffordable for BEM codes. However, the wake can be considered nearly vanished when the blade passes through the downstream region as inferred by Rezaeiha et al. [39] simulating the wake structure interaction using CFD. In the present code, the drag effect of the rotating pole is accounted for as a passive torque, reducing the overall performance.
Finally, the breakdown of momentum theory due to the turbulent wake state needs to introduce a correction of thrust at high induction factors. This is particularly experienced with high TSRs at the downstream actuator disc. The formulation obtained by Glauert's experimental fitting [40] is adopted implying a linear trend of thrust coefficient for values of axial induction factors a > 0.33 and assuming a value of C X0 = 1.816 at the maximum induction factor a = 1.

Extended Polars
Airfoil polars in the range [−180 • , +180 • ] consist of three regions: pre-stall, post-stall and deep-stall. The pre-stall is characterised by a constant lift coefficient slope, after which separation at trailing edge starts to take place. XFOIL panel code [41] is adopted to create the polar database at several Reynolds numbers, ranging 40-5000 ×10 3 , in the region of attached flow; the free transition model parameter is set to N crit = 9 under characteristic turbulence measured in the wind channel.
Viterna-Corrigan correction applies in the post-stall region, from static stall angle to the deep-stall angle where drag coefficient drastically increases. The deep-stall region is modelled by semi-empirical correlations reported in [42] both for symmetrical and cambered profiles which expressions hail from data fitting of several airfoils according to Reynolds number.
Finally, polars are smoothed before being applied as input for the BEM code, as it was proved to ensure better accuracy. Despite Reynolds regimes covered by large scale DeepWind5MW turbine is in the order of 10 7 , aerodynamic coefficients are retained approximately constant beyond 5 × 10 6 where constant self-similarity condition is achieved. It should be noted that, in constructing polars, uncertainty is mostly related to low-Re regimes where transition may occur, rather than high-Re where the flow is assumed completely turbulent. In Figure 2, lift and drag coefficients are reported in the low Re range 40-360 × 10 3 for both the airfoils analysed in the paper (NACA0021 and DU06W200).

Flow Curvature
The effect of flow curvature is introduced according to the Goude's formulation [43], which suggests to retain the original airfoil polar only changing the angle of attack. The modification in the attack angle magnitude is treated assuming two-dimensional potential flow and considering the airfoil as a flat plate. This approach produces a modification of the attack angle, resulting in a general increase on the upstream disc and in a reduction on the downstream one. Such a correction plays a key role for high-solidity rotors involving a shift of the optimal TSR. In Equation (1), x 0r is the non-dimensional hook point respect to the mid-chord (e.g., x r0 = 0 for blade attached at c 2 ).

Dynamic Stall
To account for dynamic stall, the code implements multiple semi-empirical models, namely the ones proposed by Gormont, Berg, Strickland and Paraschivoiu; however, the analyses reported in this paper are only conducted with the Strickland dynamic stall model, which was specifically formulated for VAWT applications, mounting symmetrical airfoils with 10% chord thickness. The dynamic stall calculation is dependent on the rate of change of the attack angleα which is computed for each time step (any blade position) as a first order backward difference. This may lead to a slight instability when angles of attack vary abruptly in correspondence of θ = 180 • (blade in the leeward phase).
Strickland's formulation is summarised in Equations (2)-(5), where γ is a semiempirical constant, usually ranging between 0.5 and 1.5, assumed to be equal to standard values suggested by Strickland [2] for herein BEM analyses; Sα is the sign of the derivative of the effective angle of attack; α 0 is the zero-lift angle. Besides, Equation (5) extends the original formulation to cambered airfoils, generally featuring a negative zero-lift angle.

Tip Losses
Tip losses are commonly referred to as the difference between two-dimensional and three-dimensional configuration of a lift-driven device due to pressure equalization across the suction and pressure side at the tip. The effect was first considered through tip loss factor F defined by Prandtl as the ratio between actual circulation (and hence lift) on the blade at the local spanwise section and the optimal circulation of Betz, assuming a helical rigid wake accounting for the rotational component without wake expansion. Prandtl's tip loss factor found in the literature is in the form of Equation (6) where F is a factor which tends to zero as the propeller radius r approaches the tip of the blade. Several expressions of F derived from Prandtl's theory for screw propellers (and excellent review is reported in [44]). By virtue of the different character of wake dynamics in VAWTs respect to HAWTs and propellers, Willmer's modified version proposed in [45] is adopted in this work, where g is defined in Equation (7) as a function of the vertical coordinate z and the pitch of the wake s.
The linkage between Blade Element Theory and Momentum Theory inherited by traditional DMST formulation is employed to reach convergence on induction factors a to both upstream and downstream actuator discs.
Once two-dimensional solution is converged, the tip loss factor F is applied to correct force coefficients including three-dimensional tip losses, following the approach proposed in [46].
Flow parameters are calculated according to the original DMST approach so that relative velocity and angle of attack are obtained for each blade position θ.
where subscript j ∈ [0, 2π], H is the blade semi-height, z is the vertical coordinate (z = 0 at equatorial plane) and V e corresponds to the inter-disc velocity, whereas V D assumes the values of V or V when respectively the upstream or downstream disc is of concern. The stagger and sweep angles, for unstaggered and straight blades, likewise for those examined in this research, β c and γ c are both null and δ is the angle between the blade normal and the equatorial plane.
The proposed distributed formulation of the tip loss correction is particularly relevant in this work, since it perfectly fits for the 3D character of present BEM code. Since the discretization along the turbine height is anyway necessary to consider general rotor shapes, it can be exploited to calculate F at each spanwise location. Conversely, BEM analyses found in the literature for VAWT frequently apply the Prandtl-Lanchester [47] as single overall correction to introduce the effects of the finite blade length and the induced drag losses due to tip vortices. The correction is applied to the midsection where only two-dimensional effects are assumed and predicts a loss of performance enforcing a depletion of angle of attack over the blade span. The hypothesis underlying Prandtl-Lanchester theory refers to elliptical loading under lifting line and linear aerodynamic behaviour of the airfoils, assumptions which do not hold in common VAWT operation for that post-stall conditions and dynamic stall hysteresis occurs. For these reasons, the proposed 3D formulation of the tip loss correction appears more physically consistent and of more general validity with respect to the lumped-parameter approach of the Prandtl-Lanchester one.

Strut Correction
Supporting struts play a paramount role for the stiffness and mechanical integrity of the structure in VAWTs. On the other hand, they have a detrimental effect on turbine performance, particularly at high TSRs when the flow is severely distorted by the interaction with non-aerodynamic spokes. With emphasis on the H-shaped VAWT investigated in this article, spokes are modelled as two radial flat plates with variable width as function of turbine radius. Geometric parameters of the model are available in [12] while the estimation of passive torque induced by strut drag force is evaluated according to Equations (13)- (16), as also suggested in [17]. The drag force contribution is calculated for discretised elements of the strut and the drag of the rectangular section is assumed according to [48] and equal to C Dst = 1.3 for Reynolds number in the range 10 4 -10 5 . The velocity approaching the actuator disc is interpolated at the mid point position of each discretised strut element to calculate the relative velocity W st . The overall average passive torque estimation over a rotation is estimated in Equation (16) by integrating over the strut length and over the azimuthal angle.
The strut geometry implementation benefits of the flexibility to deal with both straight horizontal arms and inclined ones.

Pole Correction
The rotating pole is modelled as a smooth stationary cylinder subjected to the wake velocity of the upstream disc at the centre of the rotor. The only effect modelled is the resistive torque exerted on the moving pole surface. The value of C D cyl = 1.15 is suggested in [4]. To simulate the overall impact of the resistive torque Equation (17) is used.
The integral value is subtracted to the aerodynamic torque generated by the turbine blades. In addition, axial thrust exerted on the pole is contributing to the overall axial force on the turbine.

Wind Shear
Realistic atmospheric wind shear profile exhibits wind direction variations throughout the atmospheric boundary layer (ABL) with significant change in turbulence, playing a paramount role in the characterisation of inflow conditions. Furthermore wind shear depends on numerous factors such as stratification, temperature, Coriolis force and many others. Hence, high and medium fidelity (Actuator Line Method) codes are particularly suited to address the problem of the effects on wind turbine performance and wake propagation or breakup [49,50], particularly when modelling wind shear as in [51]. The engineering codes can deal with simplified wind shear models consisting in power law profiles: in such a way only mean wind speed change is retained. The so-called Hellman's formulation for wind shear profile reported in Paraschivoiu [4] employs the power law defined in Equation (18) where exponent α w describes the typology of atmospheric boundary layer, V ∞ is deemed the mean flow velocity at equatorial plane which blows at z EQ level from the ground, z i is the local vertical coordinate and V ∞ i is the wind speed as function of z-coordinate (see Figure 1a).

Limitations of the Model
The in-house BEM code contains most of the sub-models in the literature, formulated to be applied to general VAWT architectures (hence, to the H-shaped and Troposkein rotors considered in this study).
As a first note, polar curves have a dramatic impact on the performance estimate and they have been found to be highly sensitive to procedure followed in their construction: this remains a critical step for all those algorithms which rely on aerodynamic look-up tables; experimental and numerical simulations, CFD or panel codes are very sensitive to set-ups (e.g., sensitivity to turbulence and surface roughness [52]), especially at low Reynolds number regimes.
Moreover, the dynamic stall model here employed is based on semi-empirical constants (e.g., γ L , γ D and K 1 ), which are obtained by experimental data fitting on a wide range of airfoil series, hence a margin of uncertainty has to be accepted since no dedicated experiments are available on airfoils such as NACA0021 and DU06W200. A study on the uncertainty propagation ascribed to dynamic stall constants can be found in [53]. Alternative to the present approach, several authors employ the Leishman-Beddoes model, since this formulation includes the effects of unsteady airfoil aerodynamics both within and beyond the stall limit, including the different effects of stall occurring at leading and trailing edge. Anyhow, this model is still based on empirical constants, which vary with airfoil shape and that are only available for high Reynolds number regimes: thus, so far the model does not find a proper application in the field of small-scale VAWTs.
Finally, the code does not consider stream tube expansion from the upstream to the downstream discs, thus the mass conservation through the actuator discs is not rigorously fulfilled. Despite this, small VAWTs are not particularly sensitive to flow expansion [18], hence the present formulation is deemed suitable for the validation of test cases and application to distributed scale machines (though simplified in the application to the large scale DeepWind turbine).
A great influence on micro and small VAWTs is represented by the role of turbulence. Turbulence intensity and scales dramatically affect VAWT performance, although it is very challenging to take into account turbulence effects in engineering codes. The way they can be addressed is to consider the effect on airfoil polars: (1) a parametric analysis of the airfoil performance based on several CFD computations varying the inlet turbulence level (2) a parametric analysis with panel codes such as XFOIL (used in this study) varying the number N crit which is the parameter linked to the turbulent intensity used in the e N method. Since BEM analyses, for the here reported validation, have been carried out with models tested in the PoliMi wind tunnel where turbulence level was very low (<1%), 3D BEM code primarily simulates study cases where turbulence impact is small as first instance. However, turbulence modelling in 3D BEM code is particularly suited for openfield applications such as the one treated herein. A future work is conceived to assess turbulence effects trough a sensitivity analyses on N crit value with XFOIL tool.

3D BEM Algorithm and Simulation SET-UP
This section is devoted to simulation set-up of the in-house 3D BEM code. An explanatory flow chart of the adopted algorithm is reported in Figure 3. By virtue of the general 3D shape of the VAWTs (and in the presence of wind shear), the numerical domain has to be discretised along three dimensions: along z-coordinate (vertical) and the azimuthal coordinate θ laying on the x − y plane.
A dedicated preliminary study was performed to identify the resolution required for obtaining predictions independent from the discretization. The analysis revealed that 80 streamtubes for each double actuator disc (∆θ ≈ 2.25 • ) are required to solve the DMST algorithm at each blade section ensuring good accuracy, especially if dynamic stall is activated. Moreover, the VAWT semi-height H has to be split in 10 slices with a grading of 1.2, the minimum blade-element thickness being 0.164c at the blade tip.
The convergence loop is based on the thrust coefficient difference evaluated by the blade element model and by the momentum theory. The residual tolerance is 0.01% of the thrust coefficient. When 100 iterations are exceeded, the induction factor is updated anyway. The algorithm is speeded-up by increasing or decreasing the induction factor according to thrust coefficient trend.
By virtue of the code structure, the VAWT performance outputs are obtained for each blade slice, allowing in-depth investigation of the turbine behaviour as a result of tip effect, blade curvature and wind shear. Tip, strut, and pole effects are introduced once the induction factors reach convergence. Therefore, 2D converged velocity field is evaluated at each strut element to calculate its effect on the overall performance. The 3D BEM code calculates averaged and phase-resolved quantities with the aforementioned settings in 30 s for each wind velocity investigated. Runtime is doubled when neither symmetry of wind profile nor the rotor specular shape can be exploited, e.g., DeepWind5MW study case, taking around 20 minutes to range the entire operating TSRs. The code has been written in MATLAB® and runs on one core Intel(R) Core(TM) CPU @ 1.30 GHz.
Tip Losses:

3D BEM-DMST Code Experimental Validation
This section discusses the experimental validation of the 3D BEM-DMST code for the three VAWT cases summarised in Table 4. VAWTs herein analysed were all tested in the GVPM at constant rotational speed, hence the undisturbed wind velocity changes covering all operating TSRs; no wind shear is obviously considered. The DeepWind turbine includes tests at different rotational speed in order to study the effect of Reynolds number. Besides, the DeepWind demonstrator envisages two configurations, namely standard upright and 15 • -tilted rotor. Airfoil shapes employed for the construction of the VAWTs are compared in Figure 4.

Microgeneration H-Shaped 200 W
The performance of the H-shaped VAWT (rated power 200 W) is first considered. The contribution of the 3D BEM sub-models are discussed with the support of Figure 5a,b to achieve a detailed comprehension of the sub-model effects in terms of power coefficient and power. A standard BEM curve ('BEM STD') gives, as expected, an overestimation of the performance; combined effects of dynamic stall and flow curvature ('+DS+FC' curve) shift C P peak at lower TSRs, which is a peculiar effect induced by flow curvature, and reduce power beyond 10 m/s. In details, flow curvature severely decreases turbine performance however dynamic stall mitigates its detrimental effect resulting in C P curve mostly affected at high TSRs, even though correspondent power depletion is limited. The effect of passive torque is remarkable ('+struts+pole') lowering the entire curve both in terms of dimensional and non-dimensional power, especially at low wind speeds; this is reasonable since the drag forces at low wind speeds 'weight' more than at high ones due to non-linear character of passive torque contributions. Same conclusion was found by CFD analysis carried out in [28] with the same H-shaped VAWT: C P difference at medium and high TSRs exceeds 0.1 percentage points between 2D and 3D simulation results. The activation of all sub-models gives 'BEM full 3D' curve which encompasses 3D tip losses: the poor C P reduction at low TSRs actually becomes as a non-negligible decrease in power from 10 m/s, resulting crucial for enhancing the accuracy of the prediction at high wind speed.   Figure 5a. However, in dimensional terms this difference appears practically negligible, being even within the experimental uncertainty as shown by Figure 5b; the good match of experimental data at high TSRs/low wind speed is somehow expected, due to the limited angles of attack and the inherent accuracy of the polars in the fully attached region. At high wind speeds, where detached flow, post-stall and deep-stall take place, the agreement is very good, thus showing that the treatment of the blade aerodynamics beyond the static stall is appropriate for this kind of machine. This might be due to the operation of the VAWT well beyond static stall angles of attack, for most part of the rotation, somehow limiting the hysteresis loops of lift and drag forces, or alternatively due to the reliability of the dynamic stall model.
Coherently with the performance prediction, in the same wind speed range the forces on the rotor tend to increase as depicted by Figure 5b, where the 3D BEM code shows a very good agreement with the experimental data, both in terms of dimensional and non dimensional values. A slight deviation of thrust force from experiments is observed at high speeds, which could be ascribed to empirical Glauert expression applied to model high loaded conditions with induction factors a > 1/3.
A more detailed investigation on the BEM code prediction capabilities on a phaseresolved basis was also possible for an operating condition close to peak C P , namely for λ = 2.4. For such condition the experimental evolution of the rotor torque coefficient was measured [54] and it is compared to the 3D BEM prediction in Figure 5d. The qualitative agreement is remarkable: the two peaks at 60 • and 90 • are well captured as well as the local minima, even though the numerical simulation tends to be slightly delayed. Conversely, minimum values show higher discrepancies which originate the slight difference in averaged C P (still anyway lower than the uncertainty band). Besides, the dashed line represents the torque produced by a single blade, being the overall torque obtained superimposing the effect of each blade with a phase shift of 120 • . This trend shows that the blade faces high drag force increase when static stall limit is exceeded, leading blade torque coefficient to drop down to negative values (θ = 120 • ). However, the corresponding active contribution of the other blades make the overall torque to remain always positive, as also found in the experiment.

Microgeneration Troposkein 450 W
The predicted performance of the Troposkein VAWT with a rated power of 450 W is compared to experiments in Figure 6, where λ is referred as equatorial TSR. The matching between the experiments and the prediction is not excellent as the one found for the H-shaped turbine, though remaining satisfactory especially in the peak C P region. The predicted peak C P value is in close agreement with the experimental datum and occurs, both in the experiment and in the prediction, in a relatively wide 'flat' region of the C P curve comprised between TSR vales of 3.1 and 3.5. As already observed for H-shaped VAWT, the non-dimensional analysis shows the highest gap of C P at high TSRs (wind speeds lower than 10 m/s); conversely the agreement once again looks better when considering dimensional quantities. For wind velocities exceeding the peak one, both shows a constant power production (Figure 6b), while experiments evidence a power decrease up to 15 m/s where a plateau is finally reached.
The thrust coefficient (Figure 6c) shows a very good agreement up to the optimal TSR value; for higher TSRs, experiments show a very slight C X reduction, while 3D BEM conversely predicts a slight increase followed by a plateau at λ = 4.  Overall, the 3D BEM code appears capable of quantitatively predicting Troposkein performance, despite the quality of the prediction reduces for high wind speed regimes. At first glance, one might have expected a level of agreement similar to the one obtained for H-shaped VAWT, since the two VAWTs share the size and airfoil shape, but the radical difference between the rotor architectures must be taken into account. Due to complex shape of Troposkein rotor, featuring bent blades towards the basement and the top of the pole, huge angles of attack are expected at very low peripheral speeds at blade tips also for relatively high vales of TSR at the equatorial section. This implies to include, for any point along the performance curve, local predictions at the blade tip suffering severe limitations of the models implemented and possibly available in the literature. Blade sections close to the pole behave more like struts and fluid structure interaction is poorly captured by 3D BEM model. This is a crucial feature of any VAWT not featuring a 3D rotor shape, and deserves a proper analysis. Assuming to cut the rotor with planes normal to the turbine rotational axis, the blade sections coincide exactly with the design airfoil shape (NACA0021) only at the equatorial plane. Otherwise, the intersections with the blade originate thicker 'virtual' airfoils approaching the tips. As a consequence, the assumption underlying the 3D BEM approach is to consider, at each blade section, the design airfoil shape and correspondent streamtube wind velocity equal to the normal component to the airfoil stacking line. The neglected effect of the span-wise velocity component in Troposkein rotors is unavoidable unless new aerodynamic coefficients are obtained for each 'virtual' airfoil shape. For the sake of simplicity and model conciseness the 'standard' approximation is considered in order to have affordable run times of the 3D BEM code.
A detailed spanwise/axial analysis of the rotor aerodynamics is reported in Figure 7a, where the incremental contribution of rotor slices to the C P and power are shown to highlight the impact of each section of the Troposkein rotor shape. Troposkein semi-span is ideally split according to vertical discretization, resulting in several plain slices for each rotor hemisphere, with a decreasing slice width towards the blade tip. The incremental C P − λ trend (Figure 7a) is determined by adding the contribution of each pair of symmetric slices incrementally, so that at z/H = 1 the power coefficient curve results the one of the entire rotor. The power coefficient plot shows, as expected, that the most 'efficient' section lays at the equatorial plane: the flow field is reasonably two dimensional and the Troposkein features the minimum solidity. It has also to be noted that the optimal C P of equatorial sections is found at TSR comparable with the one of the experiments. Progressively adding the sections towards the blade tips, the optimal C P values shift at higher TSRs. This effect stems from the fact that blade sections towards the tips operate at low local speed ratio (due to smaller radii) hence the optimal operating point inevitably lays at higher TSRs evaluated at the equatorial diameter.
The power curve in Figure 7b shows that the power contribution of the Troposkein slices is beneficial, i.e., it provides a net positive power output, up to z/H = 0.86 (with H the semi-span of the rotor). The regions closer the tip constitute, instead, a detrimental effect on the overall performance. The reason is found in high drag forces originated by low tangential velocities (at small radii) and hence high angles of attack. By the physical point of view the presence of highly distorted flow and resistive forces make the blade extremities inefficient.

DeepWind Demonstrator
The 3D BEM code is finally validated against the experiments performed on the DeepWind Demonstrator, which features a Troposkein rotor shape and the DU06W200 airfoil designed on purpose for VAWT applications. The demonstrator turbine is a downscaled model of a VAWT conceived for future floating offshore applications; its installation in the actual floating foundation foresees a backward tower tilt β equal to 15 • in standard operation. The DeepWind demonstrator was tested in the GVPM, both in upright and tilted configurations, and it is here compared to 3D BEM predictions. In the calculation, the rotor tilt is modelled considering the rotor in upright configuration subjected to a wind skew equal to the rotor axis tilt angle. Moreover, the assumption that the velocity component is normal to the rotor axis is made. Effective wind speed approaching the blade is therefore obtained by Equation (19).
The spanwise component of the free stream velocity is neglected, which appears reasonable considering the limited skew angles as imposed for 3D BEM simulations conducted in this study.
Comparisons between 3D BEM and experiments are presented for DeepWind rotational speeds of 200 rpm and 300 rpm in Figure 8a,b respectively, both in upright and tilted operating conditions. First considering the power coefficient at 200 rpm, the 3D BEM captures the trend of experiments and the peak C P values, with a light underestimate of the TSR of peak C P , for both upright and tilted rotor. It has to be noted that power coefficient estimates remain within the uncertainty bands for TSRs higher than those of the peak power. Results lose in accordance for low TSRs, where C P predictions appear very close to the experimental data but larger discrepancies appear in dimensional terms (with maximum differences between 3D BEM of about 10% at V 0 = 11m/s for both upright and tilted rotor). The reason for that may be ascribed to the different aerodynamic behaviour of blade sections. Experimental data seem to be affected by a local severe stall, which leads to a local drop in the power plot for a wind speed of 11 m/s, after that the effect disappears and the power rises again as wind speed increases. Conversely the 3D BEM power curves decrease monotonically after the peak at V 0 = 9 m/s. Polars of the 3D BEM have been smoothed and abrupt drops of lift and drag coefficients are avoided to ensure good convergence ability of the algorithm. As already discussed, dynamic stall predictions still remain affected by uncertainty at very low TSRs because of the inherent approximation of the implemented models. Both these features might explain the discrepancy between experiments and predictions at high wind speed.
Comparing Troposkein 450 W and DeepWind Demonstrator results, one may notice that 3D BEM predictions are in better agreement for the latter case although the architecture of the two machines and geometry are quite similar and the models applied are the same. Above all, the attention must be drawn on the fact that largest discrepancies are observed at high wind speeds (low TSRs) in the case of Troposkein 450 W and DeepWind Demonstrator 200 rpm in Figures 6b and 8b. Good accordance is found at low wind speeds up to the power peak, whereas discrepancies become higher above 10 m/s. When Reynolds number increases, a better agreement is found for DeepWind Demonstrator at 300 rpm. Therefore it is assumed that uncertainty related to low Reynolds numbers is evident in the comparison with experimental data.
The effect of rotor tilt is also of interest: at high wind speeds, upright and tilted rotor show minimal differences for both experimental and 3D BEM data. Since the rotor tilt involves a reduction of the normal wind velocity component, the angle of attack slightly decreases, but the effect is more noticeable at low wind speeds where a limited change of angle of attack leads to a remarkable change in lift, while drag is almost constant. At high wind speeds, where the turbine is operating in deep stall for the largest part of the rotation, small changes of angle of attack provide a drag contribution which nullify the increase of lift. Besides, 3D BEM tilted rotor results overestimate the power at high wind speed respect to experiments because the flow might be severely disturbed in the presence of combined tilted configuration and turbulence associated to higher wind speeds. The physical phenomenon can be imputed to high disturbed flow and vortex detachment, in fact differences in configurations are expected at low speeds where the flow is attached: in tilted layout every streamline (supposed to be all parallel to flow direction) will approach a thicker airfoil section due to the inclination of the rotor; this implies a completely different aerodynamic behaviour which should lead to modified polars. The effect of the combination of rotor tilt and turbulence could be accurately modelled only with a high fidelity tool such as CFD.
Better experimental accordance is found at the higher rotational speed of 300 rpm (Figure 8c). The agreement at higher Re numbers (see Reynolds regimes in Table 4 at optimal TSRs, sampled at equatorial plane), is justified by considering the challenge in accurately modelling the aerodynamic phenomena where laminar-to-turbulent transition plays a key role. The polar construction uncertainty reduces as Reynolds number increases because transition is less likely to occur. The power curve does not feature a saddle within the wind speed tested; however the gap between 3D BEM and empirical sampling increases with wind speed magnitude, although it remains within 5% of the experimental datum. The motivations are similar to the ones reported for the case at 200 rpm, that are all relaxed by the 50% increase in the Reynolds number.
In light of the three comparisons reported above, the in-house 3D BEM code is deemed to be validated to performance predictions at medium and high TSRs with low tilt angles. In the successive sections, the application of medium and large scale VAWTs will be examined.
However, as a final step of the validation, the tip loss model is considered in detail and discussed in the following subsection.

Analysis of the Tip Loss Model
One of the novel aspects of the present tool is the spanwise distribution of aerodynamic loss resulting from the finite length of the blades. This section aims at presenting a detailed analysis of the distribution of the tangential forces along the blade semi-span in order to evaluate the effect of the tip-loss model implemented in the 3D BEM code. The discussion covers the case of H-shaped 200 W and Troposkein 450 W VAWTs, as exemplary configurations of two different rotor shapes.
First considering the micro H-shaped rotor, tangential force coefficient C T distribution along normalized semi-span is shown in Figure 9 for the upstream and downstream region for TSR = 2.4 (close to C P peak). The x axis corresponds to the transversal coordinate, made dimensionless using the rotor radius (y/R = 1 corresponds to azimuthal angle θ = 0 • and y/R = −1 to θ = 180 • ). First of all, the effect of the tip loss factor F is observable as force coefficients tend to zero approaching the tip. Tangential force is predominantly generated by the first half of the rotation in the upstream disc and in the range −1 < y/R < −0.5 in the downstream one (180 o < θ < 235 • ). Contour plot in Figure 9 conveys the information about the sector where the turbine blade is producing useful power: the tangential force originates torque and power, nonetheless the link is 'weighted' by the relative velocity on the discs. Even though high values of C T are found in the downstream disc, the output torque is halved compared to torque coefficient of the blade in Figure 5d since the blade in the leeward phase is subjected to low relative velocities.
An example of C T distribution is provided also for Troposkein shape. The C T decrease due to finite blade length is barely visible since aerodynamics at tip is challenging and forces strongly depend on the vertical coordinate because velocity triangles change continuously along the span. A similar trend discussed for H-shaped rotor is valid at equatorial sections, while towards the blade tips angles of attack increase abruptly resulting in high drag forces. Moreover, dynamic stall activation originates non-physical oscillations in C T values around θ = 180 • . Therefore, local minima and maxima at the blade tips are not trustworthy. High fidelity models are only capable of capturing flow physics in such complex regions. Particularly for Troposkein rotors, the tangential coefficient is mostly provided in the upstream windward phase. As discussed, conclusions drawn for Troposkein 450 W rotor are also applicable to the DeepWind Demonstrator case, which features a Troposkein shape as well.

Validation of the Tip Loss Model
The previous analysis suggests that the tip loss model has a significant impact on the performance prediction, especially for H-shaped rotors. To validate quantitatively the action of the tip loss model, the CFD simulations presented in [16] were considered. As matter of fact, in [16] the authors proposed comparative 2D and 3D CFD simulations of a turbine very similar to the present H-shaped in Table 1, but featuring only one blade and an increased blade height to 1.5 m and operating at TSR = 3.3. The model did not include any strut, so the difference between 2D and 3D performance is only due to the tip loss effect. In particular, the overall C M obtained with the 3D CFD resulted about 10% lower than the one evaluated at midspan or equivalently with a 2D model. The same case set-up with a single blade model with longer span operating at TSR = 3.3 has been simulated with the present 3D BEM code to have a fully consistent comparison with CFD data shown in [16]. The following discussion aims at highlighting differences between low and high fidelity approaches, hence both fully 3D CFD and 3D BEM analyses will be normalised by their respective values at midspan. CFD simulation gives the trend of the torque coefficient along the blade semi-span shown in Figure 10. At the blade tip, due to pressure equalization, forces tend to zero and torque coefficient consequently. The average CFD torque coefficient over the blade span is 92% of the value at midspan. The study also compares the finite aspect ratio model by Prandtl-Lanchester carried out by Balduzzi et al. to account for tip losses when simple 2D BEM theory is applied at the equatorial plane. 3D BEM average torque coefficient over the blade is 89% respect to the value found at the equatorial plane. The trends are slightly different since CFD predicts that the blade works as the midspan section up to z/H = 0.6, conversely, 3D BEM solution deviates from 2D torque coefficient already at z/H = 0.4. The value of average torque coefficient for CFD and 3D BEM is reached respectively at z/H = 0.8 and z/H = 0.7 suggesting that the tip influence is limited to the tip region when CFD tool is adopted. MoreoverC M /C M 2D also assumes negative values in the close vicinity of the blade tip. Whereas the tip loss model entails a slightly more 'distributed' loss over the span and the sections operating in a pure 2D way are concentrated in the first half inboard blade (z/H < 0.4).
Overall, despite the inevitable local differences, and also considering the inherent complexity of the phenomenon, the tip loss model appears capable of estimating the overall effect of the finite blade length on the turbine performance, providing quantitatively comparable effects on the total torque of a VAWT blade preserving also the qualitative decrease towards the tip.

Applications
The successful double validation of the 3D BEM tool, against both experiment and high-fidelity CFD, leads to consider the code suitable for systematic applications to multiple test cases representing potential areas of VAWT use/installations.
The 3D BEM code is employed to assess VAWT performance assuming different machine scale, rotor shape, and blade profiles (mounting both NACA0021 and DU06W200). Comparative analyses have been carried out for the following three categories:
Small-scale VAWTs: H-shaped and Troposkein, referred to rated power range of 10-50 kW, for application in coastal regions; 3.
The proposed scenarios are summarised in Table 5.

Micro Scale VAWTs
VAWTs for microgeneration have been already parsed at validation stage. The aim of the present analysis is to evaluate the performance sensitivity to airfoil shape of two VAWTs for microgeneration in urban environment (micro rotors in Table 5). The selected airfoil has nearly the same thickness evaluated at 30% of the chord length; however DU06W200 is slightly cambered (0.5%) at 84% of the chord. Power coefficient and power comparison are shown in Figure 11. NACA0021 airfoil outperforms DU06W200 for both VAWT architectures, except for Troposkein rotor operating at low TSRs as confirmed by power graph in which DU06W200 slightly overtakes NACA0021 at wind speeds higher than 12 m/s. Figure 11 shows that Cp differences resulting from the application of the two profiles are lower in the case of the Troposkein rotor: this can be explained by the fact that DU06W200 outperforms NACA0021 at higher Reynolds numbers, as shown by Figure 2, and the Troposkein rotor benefits of a larger Reynolds number in the equatorial section (which is the section that harvests the largest amount of power). The differences smear out as Reynolds number increases, and hence more the aerodynamics of H-shaped rather than Troposkein showing the maximum performance gap at optimal TSR.

Small Scale VAWTs
Small scale VAWTs, referred to rated power below 50 kW, may be considered suitable for installations in coastal regions. The present analysis have been inspired by the work of Burlando et al. [55], with the aim of computing the Annual Energy Output (AEO) of two turbines installed in an harbour area under real atmospheric wind conditions. H-shaped and Troposkein architectures with comparable swept area were upscaled maintaining the same aspect ratio (height to radius) of the VAWTs tested at GVPM. For these analyses, wind shear modeling was introduced by assuming a proper wind shear coefficient α w (see Equation (18)). The wind shear coefficient was selected according to the study of Secci et al. [56]; they conducted a measurement campaign based on 1h-sampled wind data recordings over 6 years, in three coastal locations in Italy (Brindisi, Portoscuso and Termini Imerese). Measurements were made at 10 m and 50 m sensor heights above ground level in order to reconstruct the velocity profile and, from this, to derive an averaged wind shear coefficient.
The present 3D BEM analysis was performed with a wind shear coefficient equal to 0.271 representative of the 6-year average recorded in Portoscuso (Sardinia, Italy). Small scale VAWTs characteristic parameters are reported concisely in Table 6. For the present analysis, the same systematic approach used for the micro-scale VAWTs was adopted, considering rated power ranging between 10 and 20 kW ('Small rotors' in Table 6). H-shaped and Troposkein rotors were upscaled from their original dimensions in Tables 1 and 2 preserving unaltered solidity and H/D ratio of the original laboratory models. Moreover, rotational speeds and chord lengths were set identical for both upscaled turbine layouts. Figure 12 shows a sensitivity analysis to VAWT hub height and airfoil shape with wind shear in the case of H-shaped and Troposkein VAWT types. The analysis was made in terms of relative efficiency reduction, defined in Equation (20), which quantifies the difference in C P with respect to the case with the same airfoil shape and without wind shear.
Overall, wind shear effect appears limited when an imposed parabolic wind speed profile is considered, resulting in a reduction of the power coefficient, less than 1% for both H-shaped and Troposkein architectures. Whereas, when z EQ ranges from 10 to 20 m, keeping the same wind shear coefficient, the effect on the turbine C P is more visible: 2.3% power coefficient reduction at optimal TSR is observed for DU06W200, whereas a lower reduction equal to 1.2% for NACA0021. However this trend is not respected at offpeak efficiency points where NACA0021 shows higher performance reduction larger than 2%. The impact of wind shear on VAWTs is lower for the Troposkein configuration (0.8% C P loss at optimal TSR, λ = 3.5) as the extraction of power is focused on the equatorial section; H-shaped VAWT therefore turns out to be more influenced by the hub location within an atmospheric boundary layer even though the overall effect is mitigated by the tip losses. Figure 13 highlights the enhanced efficiency of DU06W200 for either H-shaped and Troposkein VAWTs, achieving a higher optimal tip speed ratio over NACA0021. Just after the peak efficiency, DU06W200 experiences a steeper decrease in C P than NACA0021, thus the latter still shows better performance at low wind speeds. Observing power curves, a plateau is gained at high wind speeds according to the attainment of high angles of attack, due to the constant rotational speed.

Large Scale DeepWind 5 MW
The DeepWind 5 MW turbine has been simulated with the 3D BEM code considering the rotor axis tilt of 15 • , as representative of FOWT actual operation. DeepWind 5 MW concept has been developed in the framework of a European project aimed to optimize aerodynamic and structural components in a coupled approach for floating applications [57,58]. The optimised DeepWind shape led to a new concept which modified the classical Troposkein shape in order to minimize the mechanical stresses (Table 7 and Figure 14).
Despite previous research studies focused on operational NACA0018, NACA0025 and DU12W264 blade shapes, herein the comparison is carried out with NACA0021 and DU06W200 profiles. For multi-MW VAWTs a control logic is required to exploit optimal TSRs, therefore a variable rotational speed is imposed. Merz et al. [59] proposed several operating schedules, among which constant rated speed (5.95 rpm) was selected: before reaching the wind speed of 8 m/s (in Table 7 referred as 'Wind speed @ rated rpm'), the rotational speed varies linearly with the wind velocity. In Table 7 the relevant geometric characteristics of the rotor are listed according to [60]. A wind shear coefficient of α w = 0.14 was assigned to model the atmospheric boundary layer, which represents a typical value for offshore applications. The wind speed distribution V ∞ is referred to z EQ = 65.55 m above sea level, where the DeepWind 5 MW exhibits the largest radius of 60.49 m, while blade attachment with the tower is placed at 15 m above sea level.
The analysis aims at establishing proven design choice at the preliminary stage and also at estimating the effects of main parameters on VAWT aerodynamic performance. In Figure 15, results of sensitivity analyses on airfoil shape, rotor tilt and wind shear are shown. Rotor tilt evidently affects turbine efficiency below the rated wind speed of 15 m/s with a reduction of power accordingly.
Wind shear, with reference to baseline condition, exceeds rated power because of higher wind speeds over the most part of the turbine height, provided that reference wind velocity is assumed to be at the equatorial plane.
The overall superposition of the effects of rotor tilt and wind shear yields to a C P reduction of almost 15% points from the baseline case (no wind shear, no rotor tilt). First, the predominant effect of the rotor axis tilt on the peak power coefficient, Figure 15a, and peak power, Figure 15d, is observed: the free stream reduction due to the inclination of the turbine tower is 10%. The introduction of a freestream wind velocity profile with α w = 0.14 originates a sensitive increase in turbine performance only because the rotor portion above the equatorial plane experiences a higher approaching velocity since the reference freestream speed is considered at the equatorial section. It is therefore noted that the effects of tilt and wind shear do not linearly combine and the overall performance is lower with respect to the one with only rotor axis tilt. Such result represents an important conclusion which might benefit from a systematic comparison with high-fidelity CFD simulations, as both effects need a complex treatment and involve a high level of resolution, which demand proper assessment. As a matter of fact, it should be noted that the Reynolds number is in the order of 10 7 , where no experimental data are available and polar construction is only based on XFOIL calculations. Correlated to that, DU06W200 profile overcomes NACA0021 because of enhanced aerodynamic performance of the former at high Reynolds number.
The outcomes at high wind speeds shrink lowering the TSR, and the differences among the different configurations above rated speed become less relevant. However, the impact of the modelling features below rated speed is significant and demonstrates the relevance of using proper tools for VAWT modelling also in the offshore context.

Conclusions
The paper presented a new formulation of well-known DMST algorithm employing already fine-tuned sub-models implemented into an up-to-date 3D BEM code, such as: The spanwise distribution of the corrections implemented led to develop a fully 3D BEM code, especially featuring the tip losses. The treatment of the correction factor F is therefore addressed and a novel implementation into the BEM algorithm is proposed. The overall effect of the tip loss model herein implemented shows good agreement with high-fidelity CFD.
The in-house code was shown to be a flexible and fast tool to investigate VAWT performance with no limitations in the rotor shape: H-shaped and Troposkein layouts have been specifically investigated in the paper. The research has benefited from a wide and comprehensive experimental assessment with data available in the literature, especially for micro-generation VAWTs at different level of detail. Excellent or very good agreement has been found in terms of averaged quantities as C p , power and thrust forces. Phase resolved analyses for H-shaped VAWT have shown a trend very similar to the experimental one, with nearly identical peak values. Rotor tilt implementation has also been validated against measurements performed for the DeepWind Demonstrator for different rotational speeds.
By virtue of the extended validation, the code has been applied to perform sensitivity analyses on rotor and airfoil shape according to diversified scenarios for VAWTs installation, from micro-scale or small-scale machines for distributed generation up to large scale offshore wind turbines. Analyses have also been performed varying the wind shear profile, considering a constant velocity at the equatorial plane, showing low effects on small Troposkein rotors and slightly higher impact on small H-shaped VAWTs. Conversely, the impact of wind shear has been found quantitatively significant for large-scale VAWT rotors, even in the less disturbed off-shore environment.
The analysis was also aimed at obtaining a fast-response in evaluating the sensitivity to geometric parameters in the preliminary design stages of VAWTs; for this reason, two airfoil shapes have been considered for selected test cases (symmetrical NACA0021 and DU06W200). It has been found that NACA0021 exhibits higher performance as long as the Reynolds number remains below 200'000, while for higher Reynolds number values the DU06W200 starts outperforming NACA0021. This result suggests a careful approach when designing a VAWT according to Reynolds characteristic regimes.
In light of such results, the here-proposed 3D BEM tool appears to be inherently suitable for estimating the energy harvested by a turbine during long-term operation: introducing a wind velocity distribution (e.g., Weibull distribution), or alternatively measured wind data at site such as those reported in [55], the 3D BEM code may allow predicting the Annual Energy Output with very limited computational cost. As a further application area, the tool developed herein may be exploited to perform a design optimization of the machine according to the specific site conditions. Many aspects still play a paramount role in trustful performance predictions: for instance, tests at low Reynolds number are still subjected to uncertainty due to laminarturbulent transition on which research efforts are underway. Moreover, polars are sensible to the detailed surface quality and roughness; this sheds light on the importance in having accurate preliminary databases to enhance prediction accuracy, otherwise widely precluded to low-fidelity models. Future works will consider the implementation of a physical-based dynamic stall model such as the one of Beddoes-Leishman, which is able to capture unsteady aerodynamics and separation physics even at small angles of attack. Eventually, the implementation of the streamtube expansion model will be also considered, in view of improving the accuracy of performance predictions for large-scale VAWT applications.
Author Contributions: Conceptualization and methodology, A.G.S., V.D. and G.P.; software, A.G.S.; formal analysis A.G.S. and G.P.; writing-review and editing, A.G.S., V.D. and G.P.; supervision, V.D. and G.P. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
a induction factor C D drag coefficient