To shed light on the effect of the icing phenomenon on the vertical-axis wind turbine (VAWT) wake characteristics, we present a high-fidelity computational fluid dynamics simulation of the flow field of H-Darrieus turbine under the icing conditions. To address continuous geometry alteration due to the icing and predefined motion of the VAWT, a pseudo-steady approach proposed by Baizhuma et al. [“Numerical method to predict ice accretion shapes and performance penalties for rotating vertical axis wind turbines under icing conditions,” J. Wind Eng. Ind. Aerodyn. 216, 104708 (2021)] was implemented, which enables the utilization of appropriate approaches for handling turbine rotation and turbulence prediction for each solver. Proper orthogonal decomposition (POD) was utilized to perform a deep analysis of the wake and aerodynamics of the wind turbine for the clean and iced turbines with large eddy simulation turbulence method. Icing causes the leading edge vortex and trailing edge vortex to separate faster than the clean case resulting in a steeper drop in the power coefficient. As for POD modes, those of the streamwise component of velocity illustrated more difference in the amount of modal energy especially at the first modes proving that the icing phenomenon mainly affects the vortex shedding of the flow structures with larger energy and size. The modes of the transversal component of velocity of the clean and iced cases demonstrated more similarity in essence, which could also be understood from the accumulated energy curve.
NOMENCLATURE
Acronyms
- AoA
-
angle of attack (degree)
- CFD
-
computational fluid dynamics
- DNS
-
direct numerical simulation
- HPC
-
high performance computing
- IEA
-
international energy agency
- LES
-
large eddy simulation
- LEV
-
leading edge vortex
- LWC
-
liquid water content (g/m3)
- MRF
-
moving reference frame
- MVD
-
median volume droplet diameter ( )
- PIV
-
particle image velocimetry
- POD
-
proper orthogonal decomposition
- RANS
-
Reynolds-averaged Navier–Stokes
- SGS
-
sub-grid scales
- SMM
-
sliding mesh method
- SVD
-
singular value decomposition
- TEV
-
trailing edge vortex
- TSR (=λ)
-
tip speed ratio
- URANS
-
unsteady Reynolds average Navier Stokes
- VAWT
-
vertical axis wind turbine
- WALE
-
wall adapted local eddy viscosity
Greek Letters
- α
-
water volume fraction
- β
-
water collection efficiency
-
azimuthal step (=30°)
- θ
-
azimuthal angle
- κ
-
von Kármán constant (-)
- μ
-
dynamic viscosity (Pa s)
- μt
-
sub-grid scale eddy viscosity (Pa s)
- ν
-
kinematic viscosity (m2/s)
- Ω
-
rotational velocity (rad/s)
- ρ
-
density (kg/m3)
- Σ
-
eigenvectors
- τ
-
wall shear stress (Pa)
-
left eigenvector
-
right eigenvector
Latin Symbols
- A
-
swept area
- c
-
airfoil chord (m)
- CD
-
drag coefficient
- CP
-
power coefficient
- D =(2R)
-
rotor diameter (m)
- Fr
-
Froude number
- g
-
gravitational acceleration (m/s2)
- H
-
turbine span (m)
- hf
-
water film height (m)
- L
-
latent heat (J/kg)
- Ls
-
mixing length of sub-grid scales (m)
-
impinging water rate (kg/m2 s)
-
evaporated water rate (kg/m2 s)
-
iced water rate (kg/m2 s)
- M
-
moment exerted on the rotor shaft (N m)
-
convective heat flux (J/s)
- Re
-
Reynolds number
- Sij
-
rate-of-strain tensor (1/s)
- t
-
time (sec)
- T
-
temperature (k)
-
streamwise component of velocity vector (m/s)
-
transversal component of velocity vector (m/s)
- V
-
cell volume
Subscripts
I. INTRODUCTION
The legislative requirements regarding environmental protection, especially those concerning greenhouse gas reduction, have led to wide interest and investments in renewable energies. In this regard, wind energy has proven to be one of the most popular forms of renewable resources (Zhang , 2019b; 2020; Posa, 2020a; Franchina , 2020; Hohman , 2020; Larin , 2016; and Posa, 2020b). To illustrate, the International Energy Agency (IEA) estimated that wind energy will gain a significant share of nearly 40% of the total energy production market by 2050 (Cozzi , 2020). As a whole, wind energy can be harvested by two main scenarios known as horizontal axis wind turbine (HAWT) and vertical axis wind turbine (VAWT). Owing to the considerably higher efficiency and power output, HAWTs have absorbed more research interest and investments from industries. However, thanks to their exceptional characteristics, VAWTs have drawn the researchers' attention in recent years (Ramirez and Saravia, 2021; Fatahian , 2022; Zare Chavoshi and Ebrahimi, 2022; Kuang , 2022b; Sheidani , 2023; Shao , 2023; Mat Yazik , 2023; and Kuang , 2023). Less sensitivity to the wind direction and turbulence, easier installation and lower noise level can be considered as the main advantages of VAWTs over HAWTs making them a suitable solution to wind energy extraction in urban areas.
Despite the pros and cons mentioned above, both types of wind turbines are prone to a sharp reduction in efficiency when exposed to severe low-temperature weather conditions leading to a phenomenon known as “ice accretion” (Gao , 2021). The importance of the blade icing could be better understood by considering the fact that wind turbines deliver the best performance in cold climates owing to the lower air density. In this regard, there have been several studies investigating the effect of ice formation on the aerodynamic performance of wind turbines. It should be noted that HAWTs have drawn almost all parts of the research activities in this field for being the primary source of wind energy production (Yirtici , 2018; Lennie , 2019; Fu and Farzaneh, 2010; Jha , 2012; and Villalpando , 2016). As for the studies conducted on the icing of VAWTs blades, there are challenges, which are needed to be dealt with which have gained research attention in the last decade. In one of the most primitive studies, Li (2010a) attached mud to the leading edge of a blade to experimentally investigate the effect of ice formation on the VAWT performance. Obviously, due to the random nature of mud attachment, such methods for the investigation of the icing effect cannot be totally trusted. Later, Li (2018) investigated the effect of ice accretion on the performance of a scaled VAWT and extracted the iced blade profile for different minutes of being exposed to icing conditions at different Tip Speed Ratios (TSR). As one of the main findings, it was proved that the blade airfoil profile loses its symmetrical profile as the exposure time increases. Nevertheless, experimental methods on the prototypes cannot still be considered highly accurate in this matter due to the fact that the non-dimensional numbers for icing are yet unknown to be matched (Baizhuma , 2021). As a result, in order to address this issue, the numerical simulation of the icing phenomenon can be implemented to achieve higher accuracy.
The very first numerical studies (Beaugendre , 2003; Li , 2010b) regarding the ice accretion phenomenon were performed on static airfoils at different angles of attack aiming to represent the VAWT blade azimuthal positions. However, the numerical simulation of static airfoils fails to completely capture the flow characteristics past a VAWT blade as the effect of rotation, which leads to dynamic stall at high relative angles of attack (Franchina , 2019; Sheidani , 2022), is not taken into consideration. It should be noted that the numerical simulation of ice accretion on rotating zones is associated with some difficulties concerning the continual alteration in the blade geometry, which needs to be updated after the formation of each layer of ice. In this regard, Manatbayev (2021) proposed a novel method, which combines the moving reference frame (MRF) in FENSAP-ICE and the sliding mesh model (SSM) in ANSYS FLUENT to perform the simulation of ice accretion on VAWT blades. The results demonstrated good agreement with the experimental data, which is going to be explicated later in this paper.
Regardless of the cardinal importance of ice accretion, the study of the wake shed behind a VAWT has risen into importance in recent years due to the fact they are being used in more compact and denser areas. On this subject, Silva and Danao (2021) conducted a numerical study to find the most optimum position of VAWTs in a tandem configuration. It was shown that the performance of the wind turbines increases with the rise in the angle with the streamwise component of velocity since the effect of wake interaction on the downstream VAWT becomes less noticeable indicating the importance of this matter. Kuang (2022a) employed improved delayed detached eddy simulation (IDDES) to investigate the effect of an upstream floating VAWT on a downstream one. It was found that pitch motion decreases the velocity deficit in the core of the wake region. It should be noted that the wake behind a single VAWT has been studied in several studies owing to the cardinal importance of this phenomenon. For instance, Abkar and Dabiri (2017) employed the actuator line (AL) and large eddy simulation (LES) methods to study the wake of a VAWT. Scale-resolving simulation models, including DNS, LES, and DES, are proven to be able to capture important flow features (Hajisharifi , 2021; 2022). The presence of self-similarity in the central part of the wake region was reported in this work. Shamsoddin and Porté-Agel (2020) studied the effect of the aspect ratio of a VAWT on the wake behind the rotor. It was reported that while the near wake region mostly keeps the main shape of the rotor in the near wake, it transmutes to an elliptical geometry in the far wake. Because of the importance of the wake shed behind a VAWT, there have been several analytical models proposed to describe the wake characteristics. In this regard, Peng (2021) drew a comparison between the analytical and high-fidelity methods, such as PIV and computational fluid dynamics (CFD), to assess the performance of the former.
Apart from the most prevalent techniques to investigate the wake characteristics, namely, pure numerical and experimental studies, proper orthogonal decomposition (POD) method can be used to reveal the coherent structures of the flow. For instance, Chen (2020) employed RANS methods to study the effect of gurney flaps on the performance improvement of an H-shaped VAWT. The POD analysis on the wind turbine wake revealed that the utilization of gurney flaps leads to the limitation of vortex strength after the fifth mode. In another study, Asim and Islam (2021) numerically investigated the effect of damage caused to the rotor of a Savonius VAWT on the wake behind the wind turbine. It was reported that the pressure modes are more noticeably affected due to the damaged rotor compared to other variables.
In addition to the scarce number of numerical studies conducted so far on ice accretion of a VAWT, to the best of our knowledge, there is no study on the wake characteristics behind a VAWT in icing conditions. Accordingly, this study is dedicated to the assessment of icing effects on the wake region of VAWT by utilizing a high-fidelity turbulence model. Therefore, the primary objective of this work is to deeply investigate how the icing condition affects the performance and, more importantly, the wake region behind the rotor of a VAWT as the mentioned goals, especially the latter, have remained unnoticed in the literature. Accordingly, at first, a validation study on a clean VAWT with respect to the experimental work of Battisti (2018) using LES will be performed. In the next step, to validate the icing solver, i.e., FENSAP-ICE, the numerical simulation for the ice accretion compared to the experimental results of Addy (2000) has been carried out. Following the validation steps, the POD analysis of the wake region based on the LES simulation of clean and iced turbines is presented followed by a conclusion on the results.
II. NUMERICAL METHODOLOGY
A. Solution strategy
Unlike typical fluid mechanic problems in which the icing takes place on a static surface, wind turbines undergo ice accretion during rotary motion. A certain number of studies tried to model icing by considering fixed blades at different angles of attack (AoA) (Manatbayev , 2021). However, as the blade rotates, the ice covers the entire surface, and thus, this assumption may be hardly applicable to the situation.
There exist two distinct methods to handle the rotational motion of a wind turbine numerically, i.e., moving reference frame (MRF) and sliding mesh method (SMM). As for the former, the rotating motion of the body is addressed by adding Coriolis and centrifugal forces as source terms into the momentum equation without any need to change the position of the body in time. More precisely, the MRF method is generally accompanied by a steady solver. As for the icing of wind turbines, the time-averaged solution obtained by the MRF method can hardly capture the phenomena caused by the unsteady motion of blades, such as the blade-wake interaction and dynamic stall (Baizhuma , 2021). Conversely, the SMM method considers the full revolution of the body, and consequently, there will be a need for a transient solver. Previous research reveals that the SMM can better capture the complex 3-D nature of the wind turbine wake region as compared to the MRF method (Tabib , 2017).
FENSAP-ICE (finite element Navier–Stokes analysis package) is a state-of-the-art icing code capable of 3-D simulations for a large variety of aerodynamic applications (Habashi , 2004). It is well established in the literature that the package offers a higher degree of fidelity compared to the other icing software (Hann, 2018). The code includes several modules, which aim to capture each step of the icing process (Hildebrandt and Sun, 2021):
-
FENSAP: the basic CFD solver based on the FEM for modeling a wide range of aerodynamic single-phase flows.
-
DROP3D: an Eulerian–Eulerian multiphase flow solver for modeling droplets. This module accounts for various interphase phenomena, including breakup, splashing, bounce, and impingement.
-
ICE3D: a 3-D module used to compute ice accretion using the frictional forces and heat fluxes from the viscous flow solution provided by FENSAP and the water volume fraction provided by DROP3D (Nakakita , 2010).
FENSAP-ICE only provides the MRF method and, therefore, has limited accuracy in modeling ice accretion for aerodynamically transient problems involving rotation. To overcome this shortcoming for the current study, the quasi-steady methodology proposed by Baizhuma (2021) has been implemented in which coupling of different solvers and numerical methods are utilized to reach the optimum accuracy. According to Fig. 1, this methodology consists of four main steps in different software. In the first step, the single-phase turbulent flow field around the rotating wind turbine is solved with the SMM solver in ANSYS FLUENT. The simulation continues until the flow field reaches a statistically stationary state, typically four revolution cycles. Afterward, the single-phase results are imported into the DROP3D to simulate the droplet flow field around the blades. Eulerian frame of reference is utilized for modeling droplets in the domain while the rotation of the blades is addressed by the MRF method. The steady simulation continues until the residuals of the equations reach . ICE3D is utilized as the final solver to account for the ice accretion on the wall surfaces, including the shaft and blades. The solver makes use of the air and droplet flow fields acquired by FLUENT and DROP3D, respectively.
The mentioned process efficiently makes use of MRF and SMM methods to accurately model the unsteady phenomena of the VAWTs where a wide range of AoAs happens along with the blade-to-blade and wake-to-blade interactions. However, if the ice accretion shape is being obtained by only one-shot, i.e., single azimuthal position, the flow characteristics of other positions will be ignored, and consequently, effects induced by the blade rotation will not be precisely captured. To address this issue, the authors proposed the utilization of several azimuthal angles or simply the multi-shot method. To this end, the ice accretion process is applied every of rotation equal to of a complete revolution. Thus, the ice accretion process is going to be carried out 12 times, continuing from the shape and position of the last simulation.
As shown in Fig. 1, the solution procedure starts with solving the flow field around the clean turbine. This simulation, which involves the SMM model, continues until the stationary flow is achieved, corresponding to at least four revolution cycles, and moreover, the blades are at the position of the next azimuthal angle. The dry air solution is then considered as the initial condition of the DROP3D solver to calculate the droplet collection efficiency on the wall boundaries. Then, the obtained air and water flow fields are imported to the ICE3D to predict the ice accretion at that azimuthal position. Here, it was considered that the turbine was under icing conditions for 30 min. It should be noted that the DROP3D and ICE3D utilize the MRF method to account for the turbine rotation. This corresponds to the ice shape for the first azimuthal angle shown in Fig. 2.
Afterward, once the iced shape is known, the new geometry is imported to ICEM CFD to generate a hexahedral grid around the turbine. The same procedure is done iteratively to obtain ice shape for the next azimuthal angle. The final ice accretion shape of the wind turbine, i.e., azimuthal angle of , is achieved by applying the intermediate ice shape process twelve times. Figure 2 shows the ice accretion through multi-shot steps for the first, third, sixth, ninth, and twelfth azimuthal angles equivalent to and , respectively.
B. Dry flow field
C. Turbulence modeling techniques
The selection of turbulence model is of great importance in fluid flow problems (Romanò , 2017; Hajisharifi , 2021; 2022). Available icing solvers are only able to capture turbulence with one- or two-equation RANS models. However, as the blades during one period of rotation undergo a wide range of Reynolds numbers from very low to high, a suitable RANS model needs to be found. To this end, several fully turbulent and transient models, including and SST, were performed at the rated tip speed ratio (TSR), and it was found that the latter, i.e., SST yields the most accurate results, similar to the findings of Rezaeiha (2019), and hence, was utilized during the ice formation simulations. In addition, all the spatial and temporal terms in Eq. (2) were discretized by the second-order method and SIMPLE algorithm was employed to address the pressure–velocity coupling.
Additionally, the convection term in Eq. (5) was discretized by the central differencing (CD) method as suggested by Li (2013) where similar to the current study, a 2.5D LES of a VAWT had been implemented. Also, the time step size was set to be as it corresponds to of rotation leading to a CFL number less than unity. It should be noted that the convergence criteria for all the variables have been set to in each time step.
D. Droplet flow field
E. Icing solver
F. Problem description
This study is aimed to perform the LES turbulence model on the clean and iced VAWT at the rated condition (TSR = 2.4). The turbine has dimensions similar to that of the experimental work of Battisti (2018) and our previous work (Sheidani , 2023). The rotor diameter of the VAWT is equal to , and it consists of three NACA 0021 airfoils with a chord length of . The turbine rotates at with an inlet velocity of , corresponding to a tip speed ratio of 2.4. Also, this study is intended to perform POD analysis and investigate the wakes shed behind the rotor as shown in the box presented in Fig. 3 where the dimensions used in the computational domains have been designated. The inlet and outlet boundaries are, respectively, on the left and right sides of the domain, and the upper and lower boundaries are considered to be of pressure outlet type.
G. Grid generation and convergence
The necessity of the full-scale simulation of the turbine, including blades and support structure (support arms and central hub), is dependent on the type of the wind turbine (Sotoudeh , 2019; Mousavi and Kamali, 2020). Siddiqui (2015) concluded that for this specific type of VAWT, i.e., H-Darrieus wind turbine, effects of the upper and lower tips are negligible. Hence, the employment of a 3D model composed of a portion of the blades is sufficiently accurate. To this end, 3D meshes were created by extrusion of a 2D mesh into the normal direction by means of one chord length, as recommended by Xu (2017) and Li (2013), while considering the periodic boundary condition for two sides of the extrusion.
The pseudo-steady procedure of the icing formation on the turbine blades requires multiple airflow simulations in which SMM was implemented. Hence, utilization of the LES turbulence model for these intermediate runs requires severe computational resources. Additionally, as mentioned earlier, common icing solvers only incorporate medium-fidelity turbulence models, i.e., RANS. To perform a grid independency study for the RANS section, three sets of meshes with , and hexahedral elements were created, respectively. There observed less than a 1% difference between the results of two later meshes, and hence, the second mesh was selected to perform the pseudo-steady icing procedure.
Mesh generation is a crucial step in the CFD simulation (Ou , 2018; Zhang , 2019a). As for the LES mesh of the clean and iced simulations, the models consisted of ∼ hexahedral elements with a total of 28 layers in the spanwise direction. The number of mesh is resulted after conducting a mesh sensitivity test on models with about , and number of elements, respectively. It was observed that the latter two yielded comparable outcomes, leading to the selection of the model with elements as the final mesh. For the clean turbine model, a fully structured grid has been employed. However, in the case of the icing, the utilization of a structured grid would result in less mesh quality due to the complex shape of the iced region of the blade, as shown in Fig. 4. Thus, an unstructured hexahedral grid was used, as it was shown to be of good accuracy as well (Elkhoury , 2015; Abdellah , 2018). The minimum and average orthogonal qualities of the icing grid were 0.5 and 0.9, respectively, which is regarded as a good quality for LES simulations (Yagmur and Kose, 2021). The minimum resolution percentage based on the index quality criterion proposed by Celik (2005) was found to be 81% (Fig. 5). According to (Pope, 2000; Celik , 2005), the value above 80% corresponds to an adequate resolved portion of the turbulent kinetic energy for the LES simulation indicating an appropriate mesh density. Finally, the first element height for wall boundaries was set to corresponding to . The schematic presentation of the LES grid can be seen in Fig. 6.
Wind turbines exhibit a transient phenomenon during their operation. CFD simulation should be continued until a statistically steady solution is obtained, i.e., characterized by repetitive flow patterns over time. Hence, monitors of critical variables have been defined during the simulation. Figure 7 demonstrates the evolution of the moment of a blade throughout the simulation time. As it is shown, a repeating behavior is emerged in the moment after six cycles, signifying the achievement of a steady solution.
H. Proper orthogonal decomposition method
To derive POD, snapshots of the velocity in the x- and y-direction were taken at each time step. After the completion of one cycle, i.e., 1500 time steps, those velocity snapshots were imported into MATLAB. A code has been written in order to derive the temporal and spatial modes of the velocity.
I. Validation study
The utilization of the quasi-steady ice accretion method has been successfully validated by Baizhuma (2021) for modeling VAWTs under icing conditions. Therefore, it is exclusively necessary to get assured of the accuracy of discrete solvers. As for the clean solver, i.e., ANSYS FLUENT, the RANS and LES simulations were performed at the peak power TSR 2.4 and the two extremes of off-design points, i.e., TSRs 3.3 and 1.5. Following this, the results were extracted and the moment coefficient vs TSR was plotted as shown in Fig. 8. The results are in good agreement with that of the experimental study of Battisti (2018) and the RANS numerical results of Franchina (2019, 2020).
For validating the multiphase and icing solvers, a static case of ice accretion on a commercial transport type airfoil representing the horizontal tailplane section of an airplane has been simulated. Geometrical details including the airfoil data points can be found in Addy (2000). Inlet boundary conditions are presented in Table I. It should be noted that the aforementioned data of water droplets, i.e., the liquid water content (LWC) and median volume diameter (MVD), was utilized for the ice accretion on the VAWT as well.
Mach number . | 0.4 . |
---|---|
T static | 258 K |
LWC | 0.4 g/m3 |
MVD | |
AoA |
Mach number . | 0.4 . |
---|---|
T static | 258 K |
LWC | 0.4 g/m3 |
MVD | |
AoA |
The simulation has been carried out using FENSAP-ICE modules. The water volume fraction and ice thickness on the airfoil are shown in Figs. 9 and 10. The maximum total thickness of ice on the airfoil after 2 min is compared for the experimental and numerical models, which indicate an error of less than 2% (Table II).
III. RESULTS AND DISCUSSION
In Fig. 11, the total moment coefficients of the wind turbine for the clean and iced cases with respect to the azimuthal angle have been presented. Overall, the icing phenomenon results in a drop in the amount of moment coefficient and consequently the performance of the wind turbine. In addition, icing causes a phase shift in the trend of the iced case as the peaks are hit at an angle in advance of the clean case. This is due to the fact that the icing phenomenon leads to an augmentation in the blade dimensions, especially near the leading edge. As a result, the maximum velocity increase on the suction side, for the iced case, is met at a point prior to that of the clean case. Moreover, the augmented leading edge causes the flow to undergo a higher amount of rotation compared to the clean case leading to a forwarded separation, which will be discussed further on completely (according to Fig. 14). Furthermore, it is also observed in Fig. 11 that the most important advantage of Darrieus over Savonius wind turbines, i.e., the nonexistence of adverse moment for power production, has clearly lost due to the ice accretion on the blades. In fact, in the Savonius wind turbines or drag tubes, as a result of the geometry of the blades, power production is impeded when one of the blades sweeps the leeward section. However, owing to the utilization of the airfoils in Darrieus wind turbines, power production is not hindered even at the most critical sections. Nevertheless, as seen in Fig. 11 the icing condition causes the moment coefficient drops to considerable negative amounts.
On the other hand, it is evident that the plummet of the moment coefficient of the iced case occurs at an azimuthal angle after that of the clean case. This observation along with considering the fact that the moment coefficient of the clean case at its lowest demonstrates a sharper variation in value proves that the icing condition leads to a more extended period of adverse power production. In addition, despite the similar trend between both cases, the moment coefficient for the iced simulation is associated with more oscillations due to the deformed shape of the blade causing a more intense vortex shedding, indicating that the fatigue load will be more critical in the case of icing.
In Fig. 12, the velocity profile in the wake region for the clean and iced cases has been presented. In general, the velocity deficit for the clean case is more significant compared to that of the icing condition, which shows the higher efficiency of the clean case. However, a similar trend between the rotor ends in both cases exists except that for the iced case, this region is shifted mostly toward the windward section. As result, the performance of the wind turbine in the icing condition drops, especially in the leeward section. This is due to the fact that as in the leeward section, the blade sees a higher angle of attack compared to other rotor sections, the distorted shape of the airfoil causes more separations and, therefore, the performance at this section will be lower even compared to that of the clean case. Moreover, for the same reason, it is evident in Fig. 12 that the area being laterally away from the rotor has been more affected by the iced blades in the vicinity of the leeward section.
Iso-surfaces of the Q-criterion are commonly used to visualize the unsteady 3D flow structures of the fluid domain. To gain an initial insight into the wake and the propagation of the structures emanating from the rotor, in Fig. 13, the Q-criterion iso-surface of the clean and iced cases at the same blade positions have been presented. Overall, as seen in for the iced case, a more extended area in the wake region is affected by the rotor vortex shedding compared to the clean case. Moreover, the clean rotor causes the flow to undergo more expansion behind the rotor which clearly demonstrates the better performance and higher power coefficient of this case. This observation is also more pronounced near the windward section of both cases owing to the low angle of attack experienced by the blade at this phase of rotation. This proves that the distorted geometry of the blade even causes the performance to drop at the sections of the lowest angles of attack. In this regard, considering the structures stemming from the leeward section, the lower expansion downstream of this phase rotation is also observed, however, at a smaller amount compared to the windward section due to the generally lower performance of the blades at this section. It should be noted that in order to have a better understanding of the effect of the iced blade on the wake region, a POD analysis of the region will be presented.
Furthermore, it is evident that the rotor is highly affected by the icing condition. On the whole, the lower energy extraction is clear due to the fewer structures being present in the iced contour of Fig. 13. As expected, the flow structures, including vortices shape and sizes downstream of each blade, do not demonstrate the same behavior. For instance, it appears that the blade in the downwind section of the iced case sees the flow separation at a point prior to that of the clean case. Therefore, the icing displaces the point of separation on the blade, which deeply affects the performance of the turbine.
In order to gain a better understanding of the icing effects on the flow separation and blade performance, the flow around the blades at different positions has been presented in Fig. 14. As it is evident in Figs. 14(a-1) and 14(a-2) for the clean and iced cases, respectively, the iced case is prone to a deeper separation even at the windward section. To illustrate, as shown in Fig. 14(a-2), the separation occurs at a point fairly closer to the leading edge compared to the clean case [Fig. 14(a-1)]. In addition, it is evident that the separation becomes significantly deeper after almost half of the chord length compared to the clean case. As a result, it is concluded that the distorted geometry of the blade causes a drop in the performance of the wind turbine even at the windward section where the blades undergo the lowest relative angle of attack. In Fig. 14(b), the Q-criterion of the blade at the upwind section demonstrates the separation of the leading edge vortex (LEV), which positively contributed to the stall procrastination and lift augmentation. However, no trace of LEV for the iced blade at the same position in the upwind section is observed. This proves that the shape of the iced blade impedes the formation and longer presence of the LEVs. The vortical structures formed in the leeward have been presented in Figs. 14(c-1) and 14(c-2) for the clean and iced cases, respectively. The vortical structures suggest that like the other sections, the flow around the iced blade separates at a point closer to the leading edge compared to the clean case. Hence, it is proven that one of the main effects of icing on the aerodynamics of the wind turbine is the advanced separation with respect to the clean aerodynamics behavior. It should be noted that this observation applies to the trailing edge vortex (TEV) whose presence is detrimental to the lift generation. To illustrate, as demonstrated in Fig. 14(c-1) the TEV is located nearer to the blade compared to the iced case, indicating that even the detrimental vortices separate faster in the iced case. As a consequence of the mentioned phenomena, there observed around 50% of the power coefficient reduction, , in the case of icing corresponding TSR2.4.
The cumulative modal energy of the clean and iced cases has been presented in Fig. 15 for the streamwise and transversal components of velocity. Overall, the results reveal that the main effect of the icing is on the streamwise component of velocity as the drop in the modal energy is more distinguishable especially, at the first modes. This offers more distinct vortical structures of flow features that could be expected for the modes of the streamwise component of velocity. Also, the iced case suggests a faster growth rate in the cumulative modal energy compared to the clean case, illustrating that the main flow structures could be found in the initial mode numbers namely before the tenth mode for both the components of velocity. Therefore, it is found that the distortion caused in the geometry caused by the icing phenomenon leads to structures larger in size and energy amount. However, as seen in Fig. 15(b) cumulative modal energy of the transversal component of velocity is less affected by the distorted geometry of the blade suggesting that the mixing phenomenon mechanisms being the main characteristics of modes of the transversal component of velocity will act more or less the same for both the clean and iced cases.
In Figs. 16 and 17, the spatial modes for streamwise and transversal components of velocity have been presented. As seen in Figs. 16(a-1) and 16(a-2), the wake region collapse point has receded by almost one diameter, indicating that the flow momentum in the icing condition is dissipated due to the flow separation stemming from the geometry alteration of the blades. In addition, the wake region in the iced case is considerably less uniform compared to the clean condition proving the lower efficiency of the icing condition. As expected from the accumulation curves, moving on to the next modes, the discrepancies between the spatial modes begin to be more apparent. In this regard, in spite of the fact that both the clean and iced simulations represent the entrainment at the point of the wake region collapse along with those emanating from the rotor, the structures for the iced case are mostly concentrated near the downwind section, while those of the clean case is more marked near the leeward. In addition, the structures corresponding to the wake region collapse appear to be denser for the clean case. Moving onto the third mode, the decay of the structures from the previous mode for both cases is apparent as expected. More importantly, it is observed in the third mode of the clean case presented in Fig. 16(c-1) that the structures have formed more symmetrically around the downwind section compared to the iced case presented in Fig. 16(c-2) where the structures could be generally found around the leeward section. This is attributed to the fact that the distorted aerodynamics of the blades in the iced case lead to a more intense vortex shedding. With the increase in the mode number, the decay of the structures in the preceding mode will occur. To illustrate, moving on to the third and the fourth modes, the structures in the second mode are observed, however, at the same location of the wake region with a lower energy content, suggesting that both the clean and iced cases follow the same trend in this regard. Furthermore, as seen in the right column of Fig. 16, representing the POD modes of the iced case, the symmetrical rolled-up structures, which come into existence due to the difference in the velocity between the wake region and the far field areas, are significantly less pronounced under the icing condition compared to the clean case. This is due to the fact that in the iced case, the aerodynamics of the VAWT is quite distorted and the efficiency drops sharply. As a result, the rotor extracts energy from the flow to a considerably lower extent leading to less difference in the streamwise velocity between the wake region and the far field area. This trend could be better seen in the fifth mode of the streamwise component of velocity presented in Fig. 16(e) for the clean and iced cases, respectively.
The POD modes of the vertical component of velocity are presented in Fig. 17. Beginning with the first modes, as shown in Figs. 17(a-1) and 17(a-2) for the clean and iced cases respectively, first, it is evident that both cases demonstrate the structures due to the wake region collapse. In addition, as seen in Fig. 17(a-2), the wake region contains more structures as a result of the entrainment before the wake region collapse point. This stems from what was previously explicated for the first streamwise component of velocity POD mode of the iced case, where the nonuniform wake region indicated the lower efficiency of the wind turbine. As a result, a more extensive part of the wake region is filled with structures representing the mixing process unlike that of the clean case. The second mode of the transversal component of velocity as shown in Figs. 17(b-1) and 17(b-2) once again proves the absence of the symmetrical rolled-up under the icing condition. As shown in Fig. 17(b-1), the flow structures mostly appear around the dashed lines being an indicator of the rotor diameter, while for the iced case, the structures are concentrated in the inner part of the wake region, which is due to the wake shed from the rotor as a result of low aerodynamics efficiency. However, it should be noted that the decay of the structures at the point of the wake region collapse could be regarded as the only similarity between the second modes of the transversal component of velocity. In addition, it is evident in Fig. 17(b-1) that the mixing structures on the verge of the wake region are more pronounced at the leeward since in clean conditions this part of the rotation is associated with the highest angle of attack and, therefore, flow momentum drop leading to more difference in velocity between the inner and outer parts of the wake region. However, with the increase in the mode number, the structures on the lower part of the wake region begin to decay and those on the upper part, i.e., leeward section, become more notable, which could be clearly seen in Figs. 17(c-1) and 17(c-2) for the clean case. However, as also observed for the modes of the streamwise component of velocity of the iced case, the flow undergoes the largest amount of momentum drop on the leeward section since the most distorted part of the blade, which is the leading edge, interacts most at the highest relative angle of attack with the wind compared to the windward where this interaction is quite negligible leading to more distinct structures at the windward section contrary to what was observed for the clean case. Moreover, it is clear in Figs. 17(c-2) and 17(d-2) that unlike the third and fourth modes of the clean case presented in Figs. 17(c-1) and 17(d-1), respectively, the structures on the edge of the wake region are not the most prominent as those emanating from the rotor are equally marked in the third and fourth modes of the iced case as a result of the lower aerodynamic efficiency. Moving on to the fifth mode, almost the same behavior is observed, however, as the structures emerging on the windward in the previous modes decay and, therefore, the both upper and lower parts of the wake region are pronounced to the same extent. It should be noted this trend of transition could be clearly seen from the second to the fifth mode of the clean case. In addition, as for the iced case, the same structures are mainly visible in the third and the fourth modes and their presence can hardly be traced in the consecutive modes.
The temporal evolution of the streamwise and transversal POD modes of the clean and the iced cases are presented in Fig. 18. As for the modes of the streamwise component of velocity presented in Fig. 18(a) for the clean and iced cases, respectively, it is observed that in both conditions, the first mode evolution is constant in time and, therefore, as expected they both prove to be non-convective. In addition, with the increase in the mode number, the frequency of the modal evolution increases for both cases with quite a similar trend. It should be noted that the modal evolution of transversal components of velocity follows the same trend. As seen in Fig. 18(b), the frequency of the modal evolution increases at the same rate, though quite symmetrically. Moreover, contrary to the streamwise component of velocity, the evolution of the first mode of the transversal component of velocity is not constant in time proving to be convective. This observation is completely in agreement with that of the spatial mode discussed previously. Therefore, considering the results of the temporal evolution of the POD modes, it is concluded that the icing condition preserves the frequency growth with the increase in the mode number through the temporal domain. It is also found that the icing condition alters the amplitude of the temporal evolution to a much lower extent with respect to the turbine's working condition, which could be regarded as different TSRs. In other words, the effects of icing could be mainly traced to an enormous phase shift in the temporal evolution of the POD modes for both of the velocity components, while a change in the working condition primarily affects the amplitude more significantly.
IV. CONCLUSION
In this study, the effect of the icing phenomenon on the POD modes of the wake shed behind the rotor of a VAWT was investigated. The main challenge that had to be overcome was to perform the numerical coupling of the continuous alteration in the geometry due to the icing and the predefined mesh motion. In this regard, the validated method of Baizhuma (2021) was implemented.
As for the practice-related effects of icing, it was found that the first exposure of the turbine to half an hour of icing condition results in a sharp drop in the power coefficient, which was calculated to be almost 50%. In addition, it was demonstrated that the icing phenomenon impairs the main characteristics of the Darrieus wind turbine as around the extremes of the windward and leeward sections, the moment coefficient drops to a negative amount.
The POD results showed that the icing leads to a significant drop in the power coefficient. Also, the results revealed that the iced blade causes the LEV and TEV to separate faster than the clean case resulting in a steeper growth and drop in the moment coefficient. As for the POD modes, those of the streamwise component of velocity proved that the first mode of the iced case is quite nonuniform which is an indicator of the efficiency of the rotor in energy extraction. Moreover, it was found that while in the clean case, the symmetrical rolled structures are vividly apparent on the borders of the wake region, the same structures are mostly absent in the iced case as a result of the lower efficiency of the wind turbine and therefore less difference in velocity between the wake and the outer region. However, the modes of the transversal component of velocity of the clean and iced case demonstrated more similarity in essence, which could also be understood from the accumulated energy curve where the modal energy of the transversal component bore a marked resemblance.
ACKNOWLEDGMENTS
This work was supported by Portuguese national funds by FCT—Foundation for Science and Technology, I.P., within the C-MAST—No. UIDB/00151/2020.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Armin Sheidani: Conceptualization (equal); Data curation (equal). Sajad Salavatidezfouli: Conceptualization (equal); Data curation (equal). Giovanni Stabile: Methodology (equal); Writing – original draft (equal). Mostafa Barzegar Gerdroodbary: Methodology (equal). Gianluigi Rozza: Supervision (equal).
DATA AVAILABILITY
The data that support the findings of this study are available within the article.