DOI: 10.3303/CET2188149 Paper Received: 28 May 2021; Revised: 9 August 2021; Accepted: 9 October 2021 Please cite this article as: Nikou T., Papadopoulos C., Vlahostergios Z., Misirlis D., Yakinthos K., 2021, Investigation of the Aerodynamic Characteristics of a Small Horizontal Wind Turbine Blade with Open-Source Software Using Advanced Turbulence Modelling, Chemical Engineering Transactions, 88, 895-900 DOI:10.3303/CET2188149 CHEMICAL ENGINEERING TRANSACTIONS VOL. 88, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Petar S. Varbanov, Yee Van Fan, Jiří J. Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-86-0; ISSN 2283-9216 Investigation of the Aerodynamic Characteristics of a Small Horizontal Wind Turbine Blade with Open-Source Software Using Advanced Turbulence Modelling Thomas Nikoua, Charalampos Papadopoulosb, Zinon Vlahostergiosa,*, Dimitrios Misirlisc, Kyros Yakinthosa a Laboratory of Fluid Mechanics and Hydrodynamic machines, Department of Production and Management Engineering, Democritus University of Thrace, Xanthi, 67100, Greece b Aristotle University of Thessaloniki, Department of Mechanical Engineering, Laboratory of Fluid Mechanics and Turbomachinery, Thessaloniki, 54124, Greece c Department of Mechanical Engineering, International Hellenic University, Serres, 62124, Greece zvlachos@pme.duth.gr Harvesting wind energy with wind turbine utilization is very important for the renewable sources industry. The maximization of a wind turbine efficiency is closely linked to the air flow characteristics around the wind turbine blades. As a result, the optimization of the flow development around the blades in respect to aerodynamic lift increase and aerodynamic drag reduction is the key factor for achieving this goal. In the current work, the characteristic airfoil of the blade of a Small Horizontal Wind Turbine (SHWT) is examined regarding its aerodynamic performance with the use of Computational Fluid Dynamics (CFD). Due to the high air flow volatility the modelling of the flow filed is very challenging while the accurate description and modelling of the effect of turbulence on the aerodynamic performance is critical for the SHWT overall design and efficient operation. For this purpose, more advanced turbulence models such as low-Reynolds number cubic eddy viscosity models are adopted in a free commercial CFD software and the results of the computations are compared with available experimental data of the lift and drag aerodynamic coefficients for various angles of attack of a selected airfoil profile of the SHWT blade. The study provides valuable information for a more accurate flow development prediction around the wind turbine blade and it shows that the adoption of more advanced turbulence models has the potential of providing SHWT designs with increased aerodynamic efficiency. 1. Introduction The usage of clean energy and especially the utilization of wind energy is very important for a sustainable environmental future. Especially during the last decades, in the renewable energy industry, wind turbines have taken an important role in the electric energy production. However, the construction of wind turbines with large pylons and increased rotor diameters, although provide increased energy generation, have an important drawback regarding their operational costs and maintenance. Additionally, the need for easily distributed electric energy to consumers and the appearance of more and more industrial renewable energy prosumers, UNIDO (2015), make the need of smaller, efficient, easily transferred and deployed wind turbines profound. Small Horizontal Wind Turbines (SHWT) are characterized by a mean rotor diameter up to 10 m, being able to produce up to 10 kW of electric power and operate in relatively small Reynolds numbers ranging from 200,000 to 350,000. The improvement of the SHWT overall efficiency is mainly focused on the aerodynamic shape optimization of the wind turbine blades for improving their aerodynamic efficiency. The operational environment of a SHWT is very challenging for blade shape optimization since the SHWT is installed inside the atmospheric buffer boundary layer region where strong shear phenomena and flow unsteadiness are present, leading to turbulence generation and boundary layer separation. If these flow conditions are not taken into consideration during the blade design phase, an overall poor wind turbine performance will finally occur. These flow conditions are also very challenging for Computational Fluid Dynamics (CFD) and the turbulence models that should be 895 used for modelling the aerodynamic flow during the design phase of the blades. Currently, the SHWT blade design methods still are based on airfoil shapes initially developed for higher Reynolds number, Van Treuren (2016). The laminar boundary layer characteristics are not sufficiently considered in the blade aerodynamics design phase resulting in not optimized SHWT designs. Significant effort is needed for improving the reliability of these methods including lower fidelity Blade Element Momentum (BEM) methods and CFD computations. For the latter, the use of more advanced turbulence models that treat better the complex transitional flow characteristics is a challenging topic. These models are more complicated but are able to better capture the complex flow development and compute more accurately the aerodynamic characteristics of the SHWT blades. Various turbulence modelling methodologies are presented in the literature such as the works of Balduzzi et al. (2018), where the γ-Reθ transitional turbulence model was used, and Lin and Sarlak (2016), where the k-kl-ω transitional model was assessed together with the γ-Reθ transitional model. Blade geometry optimization with accurate turbulence modelling is also reported, as presented in Jeong and Ha (2021). Furthermore, hybrid methodologies that combine turbulence modelling and Large Eddy Simulation (LES), Di et al. (2019), and pure LES methodologies, Li and Yang (2021), are also used for detailed design of wind turbine blades. A performance assessment regarding the accurate flow modelling of four widely used linear eddy viscosity turbulence models for an upscale wind turbine blade, can be found in Muiruri et al. (2019). A similar effort towards this direction is presented in the current work, but with the use of advanced turbulence models in order to assess their performance for use in SHWT design for low Reynolds number conditions. Two low-Reynolds number eddy viscosity turbulence models are used, the linear Launder and Sharma model as presented in Mathur and He (2013) and the non-linear model of Lien et. al (1996), with the open source CFD software OpenFOAM version 8.0 (2020). These models are not so widely used for SHWT design due to their modelling complexity which leads to numerical oscillations. To the authors knowledge their use in SHWT blade designs is not yet properly evaluated. The results of the computations are compared with available experimental data regarding the aerodynamic coefficients of a selected wind turbine blade airfoil. The study in the long run will help to the overall design blade optimization of SHWTs leading to the overall efficiency increase of wind turbines. 2. Airfoil geometry and turbulence modelling procedure Regarding the airfoil geometry that was selected for the computations, the NRL S834 arfoil, Somers (2005), low Reynolds number airfoil was chosen for the study. The specific airfoil belongs to a typical group of airfoils specifically used for SHWT blades. A characteristic view of the airfoil geometry and the computational domain with the boundary conditions is given in Figure 1. Figure 1: NRL S834 airfoil geometry and domain boundaries The numerical study is focused on the flow modelling of NRL S834 airfoil for two Reynolds numbers based on the airfoil chord with values 200,000 and 300,000. These values correspond to typical operational Reynolds numbers of a SHWT. Additionally, two turbulent intensities were considered for both Reynolds numbers with values 0.7 % and 7 % which correspond to the minimum and maximum turbulent intensities presented during a SHWT operation. The experimental measurements for validating the turbulence models were performed in a closed circuit low speed wind tunnel by the same research group and the aerodynamic lift and drag coefficient values were measured with an aerodynamic force balance as described in Papadopoulos et. al (2019). 2.1 Turbulence modelling For the modelling procedure, more advanced modelling approaches were selected in order to capture the flow development around the airfoil, especially near regions where the boundary layer has transitional characteristics. Two low Reynolds number turbulence models were adopted based on the same modelling principles regarding the modelling of the turbulence dissipation 𝜀. The first one is the linear eddy viscosity 𝑘 − 𝜀 model of Launder and Sharma as presented by Mathur and He (2013), named in this work as LS, and the second one is the 𝑘 − 𝜀 non-linear eddy viscosity model of Lien et. al (1996), named as LCL. Both models use 896 two transport equations, one for the turbulence kinetic energy 𝑘 and one for the turbulence dissipation 𝜀. Although there are numerous works in the literature regarding the modelling of low-Reynolds number airfoils, the specific turbulence models are not so widely used especially in industrial applications due to their mathematical complexity which without careful numerical adjustments may lead to major instabilities during the computations. The cubic non-linearity of the LCL model is related to the stress-strain relation between the Reynolds stresses and the strain and vorticity tensors and is given by Eq(1). For the linear models only the first linear part of equation Eq(1) is used. 𝑢𝑖𝑢𝑗̅̅ ̅̅ = 2 3 𝑘𝐼 − 2𝑓 𝜇 𝐶𝜇𝑆⏟ 𝑙𝑖𝑛𝑒𝑎𝑟 𝑝𝑎𝑟𝑡 + 𝑏1 (𝑆 2 − 1 3 𝑡𝑟𝑎𝑐𝑒𝑆 2) + 𝑏2(𝑊𝑆 − 𝑆𝑊) + 𝑏3 (𝑊 2 − 1 3 𝑡𝑟𝑎𝑐𝑒𝑊 2) −(𝛾1𝑡𝑟𝑎𝑐𝑒𝑆 2 + 𝛾2𝑡𝑟𝑎𝑐𝑒𝑊 2)𝑆 − 𝛾3 (𝑊 2𝑆 + 𝑆𝑊2 − 𝑡𝑟𝑎𝑐𝑒(𝑊2)𝑆 − 2 3 𝑡𝑟𝑎𝑐𝑒(𝑊𝑆𝑊)𝐼) −𝛾4(𝑊𝑆 2 − 𝑆2𝑊) (1) Where 𝑆 and 𝑊 are the dimensional strain and vorticity tensors while all the other constants and the damping function 𝑓𝜇 are calibrated by talking into account specific modelling calibrations regarding the turbulence equilibrium and anisotropy for each turbulence model. Another important difference of the non-linear model is the use of a strain-sensitized relation of the eddy viscosity. This relation is able to model more accurately the complex turbulent diffusion process in flow areas where laminar to turbulent boundary layer transition occurs and in the blade wake region such as the flow which is present in low Reynolds number airfoils and on SHWTs. In the general modelling approach, the eddy viscosity correlates the turbulent shear and normal stresses with the mean rate of strain and is provided by Eq(2). 𝑣𝑡 = 𝐶𝜇 ∗ 𝑓𝜇 ∗ 𝑘2 𝜀 (2) For the linear model, the coefficients are: 𝐶𝜇 = 0.09, 𝑓𝜇 = 𝑒 −3.4/(1+𝑅𝑡/50)2 (3) While for the non-linear model the coefficients are: 𝐶𝜇 = 2 3⁄ 1.25+�̅�+0.9�̅� , 𝑓 𝜇 = (1−𝑒−0.0198𝑦 ∗ )∗(𝑦∗+2𝜅/𝐶𝜇 3 4⁄ ) 𝑦∗ (4) Where �̅� and �̅� are the strain and vorticity dimensionless invariants, 𝑦∗ is a dimensionless quantity related to the wall normal distance and 𝑅𝑡 is the turbulent Reynolds number. 2.2 Numerical procedure For the numerical modelling an unstructured mesh was constructed with structured-like mesh characteristics in the near wall region having approximately 900,000 computational cells. Special care was taken for obtaining a qualitative mesh especially around the airfoil region, since this is an important adjustment for the low-Reynolds number models stability and numerical consistency. Near the airfoil boundaries the mesh inflation ensured a 𝑦+ value less than unity for approximately five to ten cells inside the viscous sublayer depending on the angle of attack (AOA) and the flow Reynolds number. This action is crucial for obtaining numerical stability for non-linear low-Reynolds turbulence models, such as the LCL, and avoid non-physical behaviour. Representative views of the overall computational mesh and the near wall region mesh are presented in Figure 2. Figure 2: Representative views of the computational mesh, overall mesh (left) and mesh inflation (right) For having a comparison with all the available experimental data, AOA ranging from -2 ° to 10 ° were used for the modelling. The inlet conditions correspond to Reynolds numbers of 200,000 and 300,000 based on the chord length and two turbulence intensities were used for these Reynolds numbers for low and high turbulence intensity cases. The inlet conditions of the modelling are summarized in Table 1 and are the same for both 897 turbulence models. The turbulence intensity decay was calibrated and compared from the domain inlet to the airfoil leading edge for having the same free stream turbulence for the LS and the LCL models. This led to the determination of the inlet turbulent length scales to the values presented in Table 1. Table 1: Inlet and operational conditions Reynolds number Inlet Velocity Inlet Turbulence intensity AOA Turbulent length scale 200,000 15 m/s 0.7 % and 7 % -2 ° to 10 ° 0.02 m and 0.06 m 300,000 10 m/s The models presented similar behaviour for the freestream decay for both Reynolds numbers and all the variations that were calculated in the aerodynamic coefficients are due to the different turbulence modelling approaches. This behaviour in the free stream region is expected since the models are based on the same modelling principles regarding turbulent kinetic energy and turbulent dissipation. The differential equations for the momentum and turbulent quantities where discretized with high order schemes and the SIMPLE algorithm was used for the pressure velocity coupling. Finally, in order to be sure that the results are grid independent, a finer grid having the double size of computational cells was used for the 10 ° AOA, which was the most computational demanding case in respect to numerical convergence, for the LS and LCL models. A small difference in the drag and lift coefficient was observed as presented in Table 2, which led to the assumption of grid independent computations for all AOA. Table 2: Aerodynamic coefficients difference for the basic and finer grid Turbulence model 𝑐𝑙 difference 𝑐𝑑 difference LS 0.2 % 3 % LCL 3. Results and discussion The results of the computations that are presented in comparison with the experimental data, refer to the aerodynamic lift and drag coefficients 𝑐𝑙 and 𝑐𝑑, which describe the main aerodynamic characteristics for the selected airfoil and are closely linked to the overall SHWT aerodynamic performance. The results for the 200,000 Reynolds number and the 0.7 % turbulence intensity are presented in Figure 4. Figure 4: Aerodynamic coefficients (1) for 200,000 Reynolds and low turbulence intensity (0.7 %) Figure 5: Aerodynamic coefficients (1) for 200,000 Reynolds and high turbulence intensity (7 %) Regarding the 𝑐𝑙 coefficient, the LCL model captures with a better accuracy the experimental values for the whole range of AOA but provides a slightly increased aerodynamic drag coefficient. As the turbulence intensity 898 increases, both models present similar behaviour as it can be seen in Figure 5. Figures 6 and 7 present the computations for the 300,000 Reynolds number and for the low and high turbulence intensities. Both models have similar behaviour especially for the 7 % turbulence intensity where the boundary layer has more turbulent characteristics and the linear model is closer to its calibration conditions. For the low turbulence intensity value of 0.7 %, which gives a boundary layer development being in the transitional region, the LCL non-linear model provides more accurate aerodynamic coefficients for all the AOA as it can be seen in Figure 6. This behaviour of the LCL model is due to the Reynolds stress non-linearity and the strain sensitized eddy viscosity and this is the reason for recommending non-linear models for demanding transitional flows, especially for low freestream turbulence intensity values and higher AOA. This observed behaviour is enhanced as the AOA increases above 6 degrees, where the aerodynamic lift coefficient distribution leaves its linear slope and flow separation due to adverse pressure gradient is present on the suction side near the airfoil trailing edge. Figure 6: Aerodynamic coefficients (1) for 300,000 Reynolds and low turbulence intensity (0.7 %) Figure 7: Aerodynamic coefficients (1) for 300,000 Reynolds and high turbulence intensity (7 %) As a general remark for the modelling results, the non-linear LCL model has the potential to compute with a better accuracy the aerodynamic coefficients, especially the lift coefficient, in transitional regions with low turbulence intensity, as the one presented in the flow around the NRL S834 airfoil. The deviation of the results is greater for the drag coefficients, since 𝑐𝑑 is smaller than 𝑐𝑙, and it is difficult to be modelled with increased accuracy. This behaviour is probably linked to the inability of the models to model all the boundary layer transitional characteristics for the selected Reynolds numbers, especially the development of turbulence near the airfoil wall where turbulence generation and dissipation are not in equilibrium. However, more experimental results are needed regarding the flow turbulence quantities for an in-depth analysis of the various mathematical modelling expressions that each turbulence model adopts. 4. Conclusions A performance assessment of two turbulence models is presented for the flow modelling of an airfoil that is used in SHWT blades. The selected models are low-Reynolds eddy viscosity models with a different mathematical expression for the Reynolds stresses. The first one is a linear eddy viscosity model with a linear expression of the Reynolds stresses (LS model) while the second one extends the linear relation to a cubic non-linear one, in respect to the strain and vorticity tensors 𝑆 and 𝑊 (LCL model). The LCL model adopts modifications in the eddy viscosity expression by adding strain-sensitized terms which model the boundary layer transitional characteristics and not only turbulent equilibrium phenomena, as the linear model does, which also computes a non-physical excess prediction in turbulence production in the transitional region. For modelling aerodynamic flows over wind turbine blades, the non-linear models are not widely used due to their numerical complexity which demands a finer and more qualitative computational domain for avoiding numerical instabilities during 899 convergence. The models assessment is focused on their ability in computing the lift and drag aerodynamic coefficients of the selected airfoil in comparison with available experimental data that are acquired by the same research group. An AOA ranging from -2 ° to 10 ° degrees was examined for two Reynolds numbers (200,000 and 300,000) in respect to the chord length with a low (0.7 %) and a high (7 %) freestream turbulence intensity. The results show that both models are able to describe with a good accuracy the aerodynamic coefficients of the airfoil for all AOA. Regarding the lift coefficient, the non-linear model seems to be more capable to capturing the experimental measurements especially for the low free stream turbulence intensity were the linear model computes more turbulence production, due to the isotropic Reynolds stresses linear expression. The LCL is the model that is recommended for modelling aerodynamic flows with transitional characteristics, such as the one that is present over wind turbine blades and aircraft wings. The use of sophisticated CFD modelling methodologies such as Hybrid RANS/LES will greatly help in resolving the complex flow characteristics near the wall regions and in the wake regions, where large eddies are dominant. As a result, more efficient wind turbine blades and turbomachinery component designs with reduced environmental and noise footprint can be thoroughly investigated. Acknowledgements This work has been implemented within the project “ADVENTUS-Advanced Small Wind Turbines” which has been financially supported by the European Regional Development Fund, Partnership Agreement for the Development Framework (2014-2020), co-funded by Greece and European Union in the framework of OPERATIONAL PROGRAMME: “Competitiveness, Entrepreneurship and Innovation 2014-2020 (EPAnEK)”, Nationwide Action: “Bilateral and Multilateral R&D Cooperations”. Authors Thomas Nikou, Zinon Vlahostergios and Dimitrios Misirlis are also affiliated to the Laboratory of Fluid Mechanics and Turbomachinery, Department of Mechanical Engineering, Aristotle University of Thessaloniki. Authors Charalampos Papadopoulos and Kyros Yakinthos are also affiliated to the UAV integrated Research Center (UAV-iRC), Center for Interdisciplinary Research and Innovation (CIRI), Aristotle University of Thessaloniki, 57001, Thessaloniki, Greece. References Balduzzi F., Bianchini A., Ferrara G., Holst D., Church B., Wegner F., Pechlivanoglou G., Nayeri C. N., Paschereit C. O., 2018, Static and dynamic analysis of a NACA 0021 airfoil section at low Reynolds numbers based on experiments and CFD, Proceedings of ASME Turbo Expo 2018 Turbomachinery Technical Conference and Exposition GT2018, June 11-15, Oslo, Norway, GT2018-75426. Di Z., Daniel R. C., Eric G. P., K. T. L., 2019, Hybrid RANS/LES Turbulence Model Applied to a Transitional Unsteady Boundary Layer on Wind Turbine Airfoil, Fluids, 4, 128. Jeong J.-H., Ha K., 2021, Numerical Investigation of Three-Dimensional and Vortical Flow Phenomena to Enhance the Power Performance of a Wind Turbine Blade, Applied Sciences, 11, 72. Li Z., Yang X., 2021, Large-eddy simulation on the similarity between wakes of wind turbines with different yaw angles, Journal of Fluid Mechanics, 921-A11, 1-24. Lin M., Sarlak C., H., 2016, A comparative study on the flow over an airfoil using transitional turbulence models, AIP Conference Proceedings, 1738, 500006, DOI: 10.1063/1.4951806. Lien F. S., Chen W.L., Leschziner M. A., 1996. Low-Reynolds-Number Eddy-Viscosity Modelling Based on Non- Linear Stress-Strain/Vorticity Relations. Engineering Turbulence Modelling and Experiments, 3, 91–100. Mathur A., He S., 2013, Performance and implementation of the Launder–Sharma low-Reynolds number turbulence model, Computers & Fluids, 79, 134-139. Muiruri P. I., Motsamai O. S., Ndeda R., 2019, A comparative study of RANS‑based turbulence models for an upscale wind turbine blade, SN Applied Sciences, 1:237. Open source field operation and manipulation-openfoam version 8.0, User Guide, 22nd July 2020, accessed 20.07.2021. Papadopoulos C., Kaparos P., Vlahostergios Z., Misirlis D., 2019, Numerical Analysis and Experimental Measurements of a Small Horizontal Wind Turbine Blade Profile for Low Reynolds Numbers, Chemical Engineering Transactions, 76, 187-192. Somers D. M., 2005, The S833, S834, and S835 Airfoils, Subcontractor Report NREL/SR-500-36340, Port Matilda, Pennsylvania. Van Treuren K.W., 2016, Small Horizontal Axis Wind Turbines: Current Status and Future Challenges, Proceedings of ASME Turbo Expo 2016, Seoul, South Korea, GT2016-57701. United Nations Industrial Development Organization – UNIDO, 2015, report: Industrial Prosumers of Renewable Energy, Contribution to Inclusive and Sustainable Industrial Development, Vienna, < www.unido.org/sites/default/files/2015-04/PROSUMERS_Energy_0.pdf> accessed 20.7.2021. 900 http://www.unido.org/sites/default/files/2015-04/PROSUMERS_Energy_0.pdf