Introduction

The growing need for sustainable energy generation has led to the development of huge wind turbine rotors. However, this is not the only way to generate more electrical power from wind. The design and optimization of mini and micro wind turbines is also becoming an essential element to this field of research due to the increasing interest in smart grid generation and widespread generation of electrical power. It means that it is possible to generate great electrical power through large rotors or through the use of numerous mini and micro rotors, which are designed and optimized for local electrical power generation. Basically, the widespread generation of electrical power is a very smart way to produce energy from wind because it can reduce the visual and environmental impact, maximizing annual energy production in the same way as the large rotors do.

However, design and optimization of mini and micro wind turbines require an in-depth study of fluid dynamics due to low Re effects like transition, laminar separation and laminar bubbles.

These phenomena are not of the utmost importance for the large rotors due to operation in high turbulent flow conditions; however, they are quite important in low Re flows, because they greatly influence performance and thus have to be taken into account in the aerodynamic modeling of micro rotors.

In this paper the authors describe the development strategy for high accurate design of micro wind turbines, operating in low Re flows, using 1D BEM modeling, 3D CFD computations and wind tunnel experiments. Following previous works [1–6], the 1D BEM model was used to generate a first geometrical draft, defined as the twisting and tapering action of the blade to maximize energy production, based on wind distributions. The widely known problems with 1D BEM model are the necessity for accurate experimental data of airfoil aerodynamic coefficients and the limitations of the 1D modeling in reproducing 3D effects.

To overcome the limitations of 1D modeling, the authors developed a 3D CFD model, using the ANSYS Fluent solver and a four equations RANS transition turbulence model, therefore optimizing it for wind turbine applications [7]. The CFD model was validated, showing excellent compatibility with experimental data. Subsequently, this model may be used to analyze the effects of the aerodynamic coefficients on the 1D model. Moreover, the wind tunnel experiment results were compared to the numerical results in terms of torque, power and power coefficient. Therefore, in designing the rotors with the use of both models, it will then be possible to greatly reduce the prototyping phase and wind tunnel experimentation, thus saving time and costs.

The entire procedure was validated through the design of two micro rotors to scale, using the 1D BEM model and analyzing them with the 3D CFD model. The larger prototype was made with wood whereas the smaller one was made with ABS in a 3D printer. Both prototypes were tested in the (0.5 × 0.5 × 1.2 m) wind tunnel using an AC brushless servo system to measure torque and rotational speed with a pc—Labview PCB interface. A hot-wire anemometer and a five hole Pitot probe were used to measure the flow speed in the open test section wind tunnel. In this way it was possible to verify the blockage effects on rotor performance as well.

Mathematical models

Review of the in-house 1D BEM model

The 1D BEM based models are powerful tools for obtaining design drafts in terms of diameter, twist, taper and rotational speed thus quickly optimizing rotor performance of HAWTs.

The in-house 1D code developed by the authors (see [1–6]), is based on BEM theory and Glauert propeller theory with suitable modifications for application to wind turbines.

In particular, 3D effects like centrifugal pumping were taken into account using a suitable mathematical representation of the lift and drag coefficients [6] whereas, for the tip losses, a modified Prandtl function was implemented [1, 2, 4, 5].

The correct evaluation of axial and tangential induction factors was obtained comparing experimental measurements and numerical simulations as well [1].

The code was validated using experimental data from literature [9], comparing generated and simulated power and torque as a function of wind speed [1]. In Fig. 1 a power comparison is presented for the NREL Phase VI wind turbine. Comparing numerical and experimental data [9] an excellent compatibility can be seen.

Fig. 1
figure 1

NREL phase VI. Comparison between experimental and simulated power

Full size image

As the 1D model was validated, it was used in this paper in a reverse way, defining a power target, based on a hypothetical PDF wind distribution. First of all the design was optimized maximizing the C l/C d ratio, it was then refined to optimize the power curve and the performance for maximizing annual energy production [5]. The final design was obtained in terms of twist and taper, choosing the optimal airfoil as well.

The most significant problem was the lack of airfoil experimental data in scientific literature, above all at very low Reynolds numbers. As BEM models are strictly dependent on the experimental airfoil coefficients, it was not a simple task to adequately scale-down the prototype fluid dynamic design for the wind tunnel size requirements. This issue was successfully resolved thorough the development of both a CFD and experimental approach. The 3D CFD model was optimized for transitional flows at very low Re and the results were tested in the wind tunnel.

Review of the CFD 3D model

CFD modeling is widely used in industrial and scientific research, especially with the superior performance of calculators in modern technology. This has led to the possibility of simulating numerous important aerodynamic phenomena related to wind turbines, thus aiding in the improvement of researchers’ designs and optimizations.

However, it is well known that the best results are obtained when large rotors are simulated. This is because RANS “full turbulent” modeling works very well at higher Reynolds numbers, reproducing fully developed turbulence [10–16].

As the mini and micro rotors operate at low Reynolds numbers, the fully turbulent modeling lacks the capability to adequately simulate this kind of flow-field and rotor performance. Laminar bubble, laminar separation and boundary layer transition are quite important phenomena that must be captured by the CFD model to reproduce the actual physical aerodynamic behavior. In general, the traditional full turbulent modeling leads to an overestimation of fluid-dynamic torque and power. This is due to the earlier separation of the laminar boundary layer in incipient stall region compared to a fully developed turbulent flow [17–22]. Moreover, the need to scale-down the prototypes for wind tunnel testing, led to the necessity to further examine the study of laminar and transitional flows.

For these reasons the authors developed a strategy to generate a HAWT high accurate 3D CFD model [7], using the correlation based transition model of Menter [20, 23, 24], calibrated for airfoils which operate at low Re number flows.

To compare the results and to validate the procedure, a model of the NREL PHASE VI wind turbine was first implemented [7], calibrating the model using experimental data found in scientific literature [9]. In this way it was also possible to qualitatively compare the results with the 1D BEM calculations (Fig. 1).

A 3D CAD was used to accurately reproduce the rotor geometry, the file was then imported into the ANSYS Workbench multiphysics platform. The computational domain was generated using the ANSYS CAD interface ‘Design Modeler’.

Spatial discretization was optimized doing a grid independence study, focusing on turbulence model limitations [23, 24] and geometric complexity. Specifically a y + < 1 and a maximum skewness value of 0.85 were reached. The mesh was later converted by Fluent solver from a tetrahedral to a polyhedral geometry. This led to a considerable reduction in the number of cells. As this type of cell considerably reduces mesh skewness, better alignment of the flow with the cell faces were obtained. Interpolation errors and false numerical diffusion decreased as well.

The innovative concept was the use of the transition SST turbulence model. To calibrate the transition model for wind turbine applications, a long process of optimizing the local correlation variables was carried out. Several simulations on typical wind turbine airfoils (S809, NASA ls421, NACA 4415) were performed using Fluent 2D. Inlet and outlet turbulence boundary conditions were also optimized (turbulent intensity and turbulent viscosity ratio).

In detail, the local correlation variables F length, Re θc , Re θt were modified through the use of a UDF, written in C language and interpreted by the Fluent solver [7].

A moving reference frame model was used for rotation in a steady state, pressure based, coupled solver. This was found to be the best balance between reasonable calculation time and accuracy of results. This was an important achievement in the obtaining of a quick response in design process and scientific research.

As the CFD 3D calculation results were in close proximity to experimental data, the methodology of generating and optimizing the CFD 3D model was considered valid. In Fig. 1 a comparison between Transition CFD calculations, 1D BEM results and experimental data [7] of generated power for the NREL PHASE VI turbine is also reported.

The NREL PHASE VI operates at Reynolds numbers of approximately 1 million. To demonstrate the validity of the methodology even at lower Reynolds numbers and validate the process of designing micro rotor for wind tunnel testing, two micro rotors with different scales were designed. They were subsequently analyzed using the proposed CFD modeling. The rotors were finally constructed and tested in the wind tunnel, owned by the University of Catania. In this way, it was possible to identify both the advantages and the limitations of the different numerical approaches, as related to the absence of experimental data for low Re flows, as well.

Design and numerical analysis of micro HAWTs

Geometrical design of prototypes

The objective of this paper was to design and analyze the performance of two scaled micro rotors to fully examine both the capabilities and the limitations of 1D BEM and CFD 3D modelling for low Re applications. The rotors were then tested in the wind tunnel to evaluate the blockage effects as well.

Through previous works by the authors [3, 4], it was possible to choose the geometric characteristics of the turbines (rotor diameter, number of blades, airfoils, chord, pitch, twist and rotational speed). It was also possible to evaluate the forces acting on the blades, torque and power at the rotor shaft.

Choosing a NACA 4415 airfoil, the rotors were first optimized maximizing the C l/C d ratio, and hence the power coefficient. To make the machines suitable for a wide range of wind conditions, the twist was modified to flatten the power curve, thus increasing annual energy production [3, 4].

One of the objectives was to examine the effects of the use of aerodynamic coefficients related to higher Reynolds numbers. As the only complete available experimental data were those reported in [25], these aerodynamic coefficients were used in the design process. These data refer to a Reynolds number of approximately 1 million while the prototypes operate at Re < 80,000.

Final features of the rotors are presented in Table 1. The rotors were named as rotor A, for the largest one, and rotor B, for the smallest one.

Table 1 Experimental rotors features
Full size table

Rotor A was made with wood, while rotor B was built in a 3D printer with ABS material. The construction clearly followed the twist and taper laws, defined by 1D BEM calculations, as reported in Fig. 2. Although rotor B was half the size of rotor A, the differences in twist and taper laws were related to the different operative conditions, defined in 1D BEM calculations. In this way it was possible to further examine different fluid dynamic behaviors.

Fig. 2
figure 2

Twist and taper design: rotor A (left) and rotor B (right)

Full size image

The prototypes are presented in Fig. 3. Both rotors were built with three separated blades, assembled with a central pentahedron as a shaft connection system and two plates for longitudinal structural connections (exploded view Fig. 3, right). An aerodynamic nose cone completed the assembly.

Fig. 3
figure 3

Fully assembled rotor A (left) and exploded view of rotor B (right)

Full size image

Numerical analysis of rotor performance

Once the features of the rotors were defined, they were analyzed both with 1D BEM and CFD 3D models, comparing the results in terms of generated fluid dynamic power.

1D BEM performance simulation analysis was quick and simple to implement. Once the geometric characteristics and the rotational speed were known, wind speed was simulated from 5 to 25 m/s, thus obtaining power curves and power coefficient curves in very few seconds. The results are presented in Fig. 4 for both types of rotors along with the CFD 3D calculation results for immediate comparison.

Fig. 4
figure 4

Calculated fluid dynamic power comparison and Reynolds number trends for rotor A (left) and rotor B (right)

Full size image

The 3D CFD models instead, require all the steps presented in paragraph 2.2 as reported in a previous work by the authors [7]. Specifically, an accurate grid independence study was carried out in both cases, reaching a y + < 1 in all simulations, as required by the transition turbulence model.

The optimal local correlation variables for very low Re flows were found after a lengthy calibration process. Several simulations on NACA 4415 airfoil were performed using Fluent 2D and the local variables were modified until the numerical aerodynamic coefficients matched the experimental one. Inlet and outlet turbulence boundary conditions were optimized as well (turbulent intensity and turbulent viscosity ratio) [7]. As complete experimental data for Re < 105 were not available in scientific literature, the calibration process results were compared to the data reported in [26].

All of the simulations were performed on a workstation HP z820, with 2 Intel Xeon E5-2695 and 6 cores for each, 128 GB of RAM memory and a NVIDIA TESLA K20c for GPU computing. To speed up calculation time, the grid was partitioned using a METIS technique on the 24 threads available. The use of the coupled solver allowed Fluent to enable the GPU computing, speeding up the calculations as well. The main features of the CFD models are summarized in Table 2.

Table 2 CFD 3D setup features
Full size table

The simulation results are presented in Fig. 4, as a generated fluid dynamic power comparison between both models. The substantial differences are quite obviously due to the use of inadequate aerodynamic coefficients and limitations of 1D modelling. Specifically, a general overestimation of the BEM calculated power at lower wind speed and an underestimation at higher wind speed are evident.

It was possible to evaluate both of these effects on performance overcoming the 1D modelling limitations using a more powerful 3D approach. This will be fully assessed in the following paragraphs where experimental and numerical results will be compared to validate procedure and hypotheses.

Experimental setup: wind tunnel features and measuring chain

Experimental tests on both types of rotors were carried out in the subsonic wind tunnel of the Department of Industrial Engineering (University of Catania, Italy). The wind tunnel was designed and calibrated by the authors [27], the main geometrical features are presented in Fig. 5.

Fig. 5
figure 5

Wind tunnel geometrical features (units in mm)

Full size image

It is a closed circuit wind tunnel with a test section of 0.5 × 0.5 × 1.2 m and it is possible to perform tests with an open or closed test section. There are four corner vanes to correctly guide the flow. An axial fan of 5.5 kW power peak at a maximum rotational speed of 1400 r/min is used to generate the flow. The fan speed is regulated using an inverter from 0 to 50 Hz. A settling chamber and a honeycomb with three screen levels are used to reduce turbulence in the test section. The tunnel was tested and calibrated by the authors following what is reported in [28–33]. The maximum flow speed possible is 30 m/s with closed test section and 27 m/s in the open test section configuration. Turbulence Intensity, measured at the contraction nozzle section was lower than 0.5 % at maximum speed [27]. The misalignment of the flow-field was also studied to adequately align the flow-field with the rotational axis of the rotors.

The measuring chain included: the use of a five hole Pitot probe to determine the real flow speed and direction; a hot wire anemometer to confirm the measured wind speed; a thermometer and a barometer to define the environmental conditions before each test. A Sanyo Denki AC brushless servo system was used for torque, power and rotational speed measurement with a National Instruments PCB–LabVIEW interface for acquisition and signal conditioning [8].

Specifically, the Sanyo Denki servo motor R2AA04010FXH1CM was controlled by the servo amplifier RS1A01AA through the National Instruments DAQ NI—6008 USB PCB to connect the servo amplifier with a pc where the LabVIEW interface allowed generation control signal and data acquisitions. A regenerative resistor was connected to the servo amplifier to dissipate braking energy. In Fig. 6 an outline of the control system is reported.

Fig. 6
figure 6

Torque, power and rotational speed control system

Full size image

Essentially, in Fig. 6, the control strategy was generated in LabVIEW. The signal was passed through the PCB DAQ to CN1 connection of the Servo Amplifier, that directly controlled the brushless generator in terms of rotational speed and braking torque. In this way the motor worked like a dynamic brake (generator) and the encoder dissipated the braking energy in the regenerative resistor. As the rotors were directly locked into the shaft, knowing the conversion efficiency of the brushless generator, the fluid dynamic torque was obtained in the PC LabVIEW Controller [8].

The control strategy was to fix the rotational speed and to vary the wind speed. In this way the system reacted by braking continuously to maintain the dynamic balance, dissipating surplus energy in the resistor, indicating the relative torque and hence generated power.

An example of the connection of rotor B to the generator shaft is presented in Fig. 7.

Fig. 7
figure 7

Connection between micro rotor B and generator shaft

Full size image

Analysis of the results: experimental validation

Experiments were first carried out for rotor A in open test section configuration. Five tests were carried out to verify the dispersion of the measurements. For each test, the rotational speed was fixed to the on-design value (n = 1000 r/min), through the servo system control while the flow speed was gradually increased from 0 to 26 m/s and then decreased. Resistive torque was obtained directly on the LabVIEW interface and therefore generated power was calculated.

The experimental and numerical power data were compared and reported in Fig. 8.

Fig. 8
figure 8

Experimental, CFD 3D and 1D BEM power comparison for rotor A

Full size image

It was possible for the rotor A measurements to be affected by blockage effects. This was due to an overly large blockage ratio, \( B_{\text{r}} = \frac{{A_{\text{wt}} }}{{A_{\text{ts}} }} = 0.318 \). To further examine this issue and validate numerical models, the experiments were replicated for rotor B (\( B_{\text{r}} = \frac{{A_{\text{wt}} }}{{A_{\text{ts}} }} = 0.159 \)), using the same methodology. The results are presented in Fig. 9 where the same generated and calculated power comparison is proposed.

Fig. 9
figure 9

Experimental, CFD 3D and 1D BEM power comparison for rotor B

Full size image

Due to blockage effects, the experimental results for rotor A underestimated the actual power, as evidenced in Fig. 8. In an open test section, in fact, the divergence of the flow field, occurring at a short distance from the contraction nozzle and the wake expansion effects caused by the size of the larger rotor (A), lead to a slowing of the fluid flow and to a misalignment of the flow field [28]. Hence, the experiment for rotor A was useful in the verification of the limitations of rotor dimensions to avoid blockage effects.

The above is supported by the results presented in Fig. 9. In this case, in fact, CFD results were in close proximity to experimental data and blockage effects were definitely of minor importance as rotor B had a diameter which was half the diameter of rotor A. The optimized transition CFD 3D model was therefore considered valid even at very low Reynolds numbers and it can be used to analyze, and finally optimize, rotor geometry in a wide range of transitional flow conditions. It was useful in assisting wind tunnel testing in the proposed design process as well.

To further confirm the validation, a power coefficient comparison between CFD calculations and experimental data is presented in Fig. 10. The chart refers to rotor B at n = 2450 r/min.

Fig. 10
figure 10

Power coefficient vs λ comparison between CFD and experimental results (rotor B)

Full size image

The range of validation was limited by the maximum possible rotational speed of the servo system (n = 4500 r/min). For this reason, in Fig. 10, the experimental data are related to low λ. However, all the following considerations were made in the validation range of the CFD model.

Specifically, through the use of CFD post-processing, referring to the validated rotor B results, it was possible to analyze the discrepancies between both numerical models.

A comparison of calculated flow-field at different radial stations and different wind speeds was obtained using ANSYS CFD post-processor and is presented in Fig. 11. The sections were obtained using cylindrical surfaces at a distance of 0.026 m (near the hub), 0.06 m (medium radius) and 0.1 m (near the tip). Three different wind speeds were examined as well: low (10 m/s), medium (20 m/s) and high (30 m/s). By doing so, it was possible to analyze a wide range of fluid dynamic rotor behaviors.

Fig. 11
figure 11

Calculated flow-field at different radial stations r = 0.026 m (a), r = 0.06 m (b), r = 0.1 m (c) for different wind speeds

Full size image

At a lower wind speed (Fig. 11, 10 m/s) the flow was mostly attached along the blades. Only laminar separation bubbles were detectable and this transitional behavior was the reason why the 1D BEM calculations overestimated power at lower wind speeds, as seen in Fig. 9. The excellent capabilities in both transitional modeling and behavior of the proposed CFD 3D model were clearly evident.

More interesting was the fluid dynamic behavior at medium and higher wind speeds. The flow in fact is fully separated along the blade, already at 15 m/s. The phenomenon is shown in Fig. 11 (20; 30 m/s). The transition between attached and separated flow is abrupt as the leading edge stall characterizes low Reynolds flows. This is well captured by the CFD model. Specifically, the influence of the hub at r = 0.026 m, where the flow is accelerated and the stall vortex is flattened, is also evident. Moreover, extrapolating the sectional AoAs, post-processing CFD calculations [34, 35], it was possible to highlight a particular fluid dynamic behavior that influences rotor performance, as evidenced in Fig. 9.

In Fig. 12, a comparison between numerical predictions of AoA trends along the blade, for three different wind speeds, is presented. The way in which the 1D BEM model underestimates angles of attack along the entire blade has been highlighted. It is therefore possible to conclude that lift and drag coefficients will also be different than the actual ones.

Fig. 12
figure 12

Comparison of CFD 3D and BEM calculated AoAs at different wind speeds (rotor B)

Full size image

In fact, through the calculations of the 3D aerodynamic coefficients, from CFD simulations, and comparing them with 2D data, the effects of the centrifugal pumping can be evaluated. As is known from literature [35], the 3D coefficients are much higher than the 2D values even in deep stall conditions.

For example, at 20 m/s, r = 0.06 m (medium radius) and AoA = 46°, the calculated 3D aerodynamic coefficients were C l = 1.49 and C d = 1.42 (medium Reynolds number along the blade was 33,000) while corresponding experimental 2D values at Re = 106 are C l = 1.39 and C d = 1.32. This indicates that centrifugal pumping effects are clearly evident in deep stall conditions as well, increasing forces on the blades compared to 2D conditions.

In this way, two issues were highlighted. This is of the utmost importance because 1D BEM code is a fundamental tool in wind turbine design and optimization. Therefore, the joint use of a validated CFD 3D model and wind tunnel experiments will allow the BEM code to be refined.

Certainly, the errors introduced through the use of inadequate aerodynamic coefficients lead to errors in the optimal design and modeling. The authors will try to resolve this issue in future works, generating a database of 2D aerodynamic coefficients at very low Re numbers through the use of experimental and CFD calculations. However, the increase of aerodynamic coefficients, due to centrifugal pumping effects, even in deep stall conditions, suggests that it is necessary to further study this phenomenon. It is also important to take this phenomenon into account in the BEM code, modifying the proposed empirical equations [6] for centrifugal pumping in low Re conditions. This will also be examined in future works.

An example of radial flow along a blade was obtained using the ANSYS CFD-Post processing and is presented in Fig. 13 for a wind speed of 20 m/s, in deep stall condition. A centrifugal vortex is highlighted along the entire blade and this supports the hypothesis that centrifugal pumping heavily influences fluid dynamic behavior, even in fully stalled conditions.

Fig. 13
figure 13

Streamlines of velocity highlighting radial flows along a blade (rotor B; V w = 20 m/s)

Full size image

Conclusions

In this paper a design strategy for low Re HAWTs is proposed and verified. Through the joint use of 1D BEM code, 3D CFD modeling and wind tunnel experiments, the authors examined the advantages and limitations of numerical modeling, thus resolving issues related to very low Reynolds number flows, leading the way for future research.

Using results from previous works, an in-house 1D BEM model was applied to obtain a first geometrical draft of two micro rotors for wind tunnel applications, thus optimizing the twist and taper of the blades and fixing the on-design rotational speed.

To fully examine the effects of the Reynolds number on performance, the authors developed a 3D CFD model, using and optimizing a SST transition turbulence model. The CFD model was validated in a previous work and was applied to verify its capabilities in reproducing very low Reynolds number flows as well.

The rotors were simulated with both of the numerical codes and the results were compared to wind tunnel experiments, carried out in the subsonic wind tunnel, belonging to the University of Catania.

In addition, rotor dimensions were chosen to verify the widely known blockage effects in the open test section of the wind tunnel. The largest rotor (rotor A), had a diameter which was twice the diameter of the other smaller rotor (rotor B). From the experimental results, an underestimation of power for rotor A, due to high blockage ratio, was verified. For rotor B, however, the power comparison showed a good correlation with 3D CFD computations, demonstrating good predictive capabilities of the modified transition turbulence model, even at very low Re conditions.

As highlighted by the CFD computation results, errors were induced in the 1D model through the use of unsuitable aerodynamic coefficients, related to higher Reynolds numbers. This led to errors in the optimal design and modeling. This is an important issue with the 1D BEM modeling that must be resolved to improve the good reliability of this code at low Reynolds number ranges as well. With the support of the proposed CFD model and wind tunnel experiments, a database of aerodynamic coefficients at low Re will be generated for implementation into the BEM model.

CFD results also showed that centrifugal pumping effects heavily influence fluid dynamic behavior along blades, increasing the 3D aerodynamic coefficients, even in deep stall conditions. This proved to be true not only in incipient stall regions (stall delay). In fact, comparing calculated angles of attack and the 3D aerodynamic coefficient values, it was possible to deduce that the combination of the widely known sensitivity to flow separation at low Re with the centrifugal effects increase the aerodynamic forces more than at higher Reynolds numbers.

Both of the aforementioned effects will be further examined in future works to define new equations. This definition is necessary to accurately take into account the 3D radial flow inside the 1D numerical code, based on BEM theory.