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.

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 ( μ m)

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)

Φ u

left eigenvector

Φ v

right eigenvector

Latin Symbols
A

swept area A = H × D

c

airfoil chord (m)

CD

drag coefficient C D = F D / ( 1 / 2 ρ V 2 A )

CP

power coefficient C P = M λ / ( 1 / 2 ρ V 3 A R )

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)

m ̇ com

impinging water rate (kg/m2 s)

m ̇ evap

evaporated water rate (kg/m2 s)

m ̇ ice

iced water rate (kg/m2 s)

M

moment exerted on the rotor shaft (N m)

Q ̇ h

convective heat flux (J/s)

Re

Reynolds number

Sij

rate-of-strain tensor (1/s)

t

time (sec)

T

temperature (k)

u

streamwise component of velocity vector (m/s)

v

transversal component of velocity vector (m/s)

V

cell volume

Subscripts
a

air

d

droplet

evap

evaporation

f

water film

fus

fusion

r

relative frame

subl

sublimation

w

water

free stream

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.

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 1 × 10 8. 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.

FIG. 1.

Solution strategy of the quasi-steady solver for the icing of the wind turbine.

FIG. 1.

Solution strategy of the quasi-steady solver for the icing of the wind turbine.

Close modal

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 30 ° of rotation equal to 1 12 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.

FIG. 2.

Ice formation evolution at different azimuthal positions of the blade.

FIG. 2.

Ice formation evolution at different azimuthal positions of the blade.

Close modal

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 360 °, 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 θ = 30 ° , 90 ° , 180 ° , 270 ° , and 360 °, respectively.

ANSYS FLUENT is a well-accepted CFD simulation software capable of modeling a broad range of fluid flow problems including turbulence and heat transfer (Wang , 2021). The package is highly scalable and robust in terms of high-performance computing (HPC) for solving complex CFD simulations, such as wind turbines (Thangavelu , 2021). The unsteady Reynolds average Navier–Stokes (URANS) formulation has been used to model fluid flow around the rotating wind turbine, which is as follows:
(1)
(2)
where U i ¯ and ρ u i u j ¯ are the averaged velocity and Reynolds stress tensor, respectively. The latter can be written as
(3)
where the mean strain rate tensor is
(4)
The turbulent viscosity term, μt, can be addressed with common one- or two-equation eddy viscosity models to provide turbulence closure.
As for the LES model, filtered Navier–Stokes is utilized,
(5)
where the last term, τ ij, is the sub-grid scale (SGS) stress and is defined as follows:
(6)
There exist a couple of closure models to address SGS in terms of the local resolved flow. The selection of the appropriate turbulence closure model for RANS and LES will be discussed further on.

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 γ R e θ , k t k l ω and k ω SST, were performed at the rated tip speed ratio (TSR), and it was found that the latter, i.e., k ω 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.

Regarding the LES simulation of the turbine at the end of the icing formation as well as the clean turbine, wall-adapting local eddy-viscosity (WALE) model (Nicoud and Ducros, 1999) was implemented to capture the sub-grid scale (SGS) structures as suggested by Ma (2009) to be the most suitable modeling technique for the LES computation of the flow physics of wind turbines. The SGS is calculated as follows:
(7)
where
(8)
(9)
(10)
The WALE constant, Cw, was set to 0.325.

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 10 4 s as it corresponds to 0.24 ° 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 10 4 in each time step.

At each azimuthal position, the flow field of water drops is determined by solving the conservation of mass and momentum:
(11)
(12)
where α, u d , r, K, C D, and Fr are the droplet volume fraction, velocity in a relative frame of reference, inertial parameter, drag coefficient, and Froude number, respectively. u a , r is the velocity vector of air, which is obtained from the previous solver. The first and second terms on the right-hand side represent the drag force acting on the droplet and the Buoyancy force, respectively. However, in the case of VAWT, the gravity force is parallel to the rotating axis and can be neglected.
The collection efficiency on the wall surfaces can be determined by the following relation (Son and Kim, 2020):
(13)
Consequently, the mass flow rate of the impinging water on the surfaces is determined as follows (Son and Kim, 2020):
(14)
It was assumed that the colliding drops stick to wall surfaces, regardless of the material type. Finally, the mass rate is transferred to ICE3D to evaluate the freezing ice mass.
ICE3D can assess the ice thickness on the walls once the mass flow rate of the impinging water is determined. To do so, ICE3D solves mass and energy conservation equations. Mass conservation is written as follows (Beaugendre , 2006):
(15)
For the energy conservation (Beaugendre , 2006),
(16)
The terms on the right-hand side indicate sources of kinetic energy, evaporation, sublimation, latent heat of fusion, radiation, and convective cooling, respectively. More details of the icing solver can be found in ANSYS (2017).

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 1.03 m, and it consists of three NACA 0021 airfoils with a chord length of 0.086 m. The turbine rotates at 400 rpm with an inlet velocity of 9.1 m / s, 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.

FIG. 3.

Dimensions of the computational domain and POD box for wake analysis.

FIG. 3.

Dimensions of the computational domain and POD box for wake analysis.

Close modal

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 1 × 10 6 , 2.1 × 10 6, and 3.8 × 10 6 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 ∼ 1 × 10 7 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 6 × 10 6 , 1 × 10 7, and 1.3 × 10 7 number of elements, respectively. It was observed that the latter two yielded comparable outcomes, leading to the selection of the model with 1 × 10 7 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 5 × 10 6 corresponding to y + 1. The schematic presentation of the LES grid can be seen in Fig. 6.

FIG. 4.

LES mesh on the iced blades.

FIG. 4.

LES mesh on the iced blades.

Close modal
FIG. 5.

Contour of index quality for the LES mesh.

FIG. 5.

Contour of index quality for the LES mesh.

Close modal
FIG. 6.

Overview of the LES mesh of the iced model (a) details of the mesh around the turbine and consequent cell zones and (b) details of the mesh around the iced blade.

FIG. 6.

Overview of the LES mesh of the iced model (a) details of the mesh around the turbine and consequent cell zones and (b) details of the mesh around the iced blade.

Close modal

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.

FIG. 7.

Time history for the moment coefficient on one blade.

FIG. 7.

Time history for the moment coefficient on one blade.

Close modal
Proper orthogonal decomposition is employed to extract the coherent structures in both the spatial and temporal domain (Li , 2013). It should be noted that the POD method implemented in this study is based on the work of Taira (2017). In order for POD to be implemented, the singular value decomposition (SVD) method was employed owing to its robustness against round-off errors (Lengani , 2014). In order to perform SVD, the “matrix of snapshots” denoted as “X” needs to be formed. The rows of X include the value of the variable of interest corresponding to its position in the computational domain, and the columns correspond to the time instance at which the data have been extracted. Having formed the X matrix, SVD is applied as follows in Eq. (17):
(17)
where Φ u and Φ v T are called the left and right eigenvectors, respectively. Also, Σ is a diagonal matrix containing the eigenvalues organized in the descending order. Therefore, each element of the Σ matrix demonstrates the importance of each column of Φ u. It should be noted that the left eigenvector corresponds to the spatial coherent structures, which are in other words, the modes of the flow. On the other hand, the right eigenvector indicates the evolution of each mode in the time domain.

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.

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).

FIG. 8.

Power coefficient validation with respect to the experimental work of (Battisti , 2018).

FIG. 8.

Power coefficient validation with respect to the experimental work of (Battisti , 2018).

Close modal

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.

TABLE I.

Icing conditions.

Mach number 0.4
T static  258 K 
LWC  0.4 g/m3 
MVD  42 μ m 
AoA  0.6 ° 
Mach number 0.4
T static  258 K 
LWC  0.4 g/m3 
MVD  42 μ m 
AoA  0.6 ° 

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).

FIG. 9.

Water volume fraction around the static airfoil.

FIG. 9.

Water volume fraction around the static airfoil.

Close modal
FIG. 10.

Ice thickness on the static airfoil.

FIG. 10.

Ice thickness on the static airfoil.

Close modal
TABLE II.

Maximum ice thickness on the static airfoil.

Thickness (mm)
Experimental  3.6 
CFD simulation  3.54 
Error %  1.6 
Thickness (mm)
Experimental  3.6 
CFD simulation  3.54 
Error %  1.6 

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.

FIG. 11.

Comparison of the moment coefficient profiles of the clean and iced turbine.

FIG. 11.

Comparison of the moment coefficient profiles of the clean and iced turbine.

Close modal

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.

FIG. 12.

Comparison of wake profiles at TSR2.4 of the clean and iced turbine at the distance of 3.5D from the turbine.

FIG. 12.

Comparison of wake profiles at TSR2.4 of the clean and iced turbine at the distance of 3.5D from the turbine.

Close modal

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.

FIG. 13.

Iso-surfaces of Q-criterion colored by velocity magnitude for the entire flow domain.

FIG. 13.

Iso-surfaces of Q-criterion colored by velocity magnitude for the entire flow domain.

Close modal

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, C P = 0.14, in the case of icing corresponding TSR2.4.

FIG. 14.

Iso-surfaces of Q-criterion colored by velocity magnitude at different rotation section: (a) vortex Shedding at windward section; (b) vortex Shedding at upwind section; and (c) vortex Shedding at leeward section.

FIG. 14.

Iso-surfaces of Q-criterion colored by velocity magnitude at different rotation section: (a) vortex Shedding at windward section; (b) vortex Shedding at upwind section; and (c) vortex Shedding at leeward section.

Close modal

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.

FIG. 15.

Cumulative modal energy curve of (a) streamwise component of velocity and (b) transversal component of velocity.

FIG. 15.

Cumulative modal energy curve of (a) streamwise component of velocity and (b) transversal component of velocity.

Close modal

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.

FIG. 16.

(a) First, (b) second, (c) third, (d) fourth, and (e) fifth POD mode for the streamwise component of velocity at the design tip-speed ratio.

FIG. 16.

(a) First, (b) second, (c) third, (d) fourth, and (e) fifth POD mode for the streamwise component of velocity at the design tip-speed ratio.

Close modal
FIG. 17.

(a) First, (b) second, (c) third, (d) fourth, and (e) fifth POD mode for the transversal component of velocity at the design tip-speed ratio.

FIG. 17.

(a) First, (b) second, (c) third, (d) fourth, and (e) fifth POD mode for the transversal component of velocity at the design tip-speed ratio.

Close modal

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.

FIG. 18.

Time evolution of POD coefficients of the (a) streamwise and (b) transversal component of velocity.

FIG. 18.

Time evolution of POD coefficients of the (a) streamwise and (b) transversal component of velocity.

Close modal

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.

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.

The authors have no conflicts to disclose.

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).

The data that support the findings of this study are available within the article.

1.
Abdellah
,
B.
,
Belkheir
,
N.
, and
Khelladi
,
S.
, “
Effect of computational grid on prediction of a vertical axis wind turbine rotor using delayed detached-eddy simulations
,” in
2018 International Conference on Wind Energy and Applications in Algeria (ICWEAA)
(
IEEE
,
2018
), pp.
1
5
.
2.
Abkar
,
M.
and
Dabiri
,
J. O.
, “
Self-similarity and flow characteristics of vertical-axis wind turbine wakes: An LES study
,”
J. Turbul.
18
,
373
389
(
2017
).
3.
Addy
,
H. E.
,
Ice Accretions and Icing Effects for Modern Airfoils
(
National Aeronautics Administration, Glenn Research Center
,
2000
).
4.
ANSYS
. Ansys FENSAP-ICE user manual 18.2. (
2017
).
5.
Asim
,
T.
and
Islam
,
S. Z.
, “
Effects of damaged rotor on wake dynamics of vertical axis wind turbines
,”
Energies
14
,
7060
(
2021
).
6.
Baizhuma
,
Z.
,
Kim
,
T.
, and
Son
,
C.
, “
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
).
7.
Battisti
,
L.
,
Persico
,
G.
,
Dossena
,
V.
,
Paradiso
,
B.
,
Castelli
,
M. R.
,
Brighenti
,
A.
, and
Benini
,
E.
, “
Experimental benchmark data for h-shaped and troposkien VAWT architectures
,”
Renewable Energy
125
,
425
444
(
2018
).
8.
Beaugendre
,
H.
,
Morency
,
F.
, and
Habashi
,
W. G.
, “
FENSAP-ICE's three-dimensional in-flight ice accretion module: ICE3d
,”
J. Aircr.
40
,
239
247
(
2003
).
9.
Beaugendre
,
H.
,
Morency
,
F.
, and
Habashi
,
W. G.
, “
Development of a second generation in-flight icing simulation code
,”
J. Fluids Eng.
128
,
378
387
(
2006
).
10.
Celik
,
I. B.
,
Cehreli
,
Z. N.
, and
Yavuz
,
I.
, “
Index of resolution quality for large eddy simulations
,”
J. Fluids Eng.
127
,
949
958
(
2005
).
11.
Chen
,
L.
,
Xu
,
J.
, and
Dai
,
R.
, “
Numerical prediction of switching gurney flap effects on straight bladed VAWT power performance
,”
J. Mech. Sci. Technol.
34
,
4933
4940
(
2020
).
12.
Cozzi
,
L.
,
Gould
,
T.
,
Bouckart
,
S.
,
Crow
,
D.
,
Kim
,
T.
,
Mcglade
,
C.
,
Olejarnik
,
P.
,
Wanner
,
B.
, and
Wetzel
,
D.
, “
World energy outlook 2020
,” 2050 (
2020
), pp.
1
461
.
13.
Elkhoury
,
M.
,
Kiwata
,
T.
, and
Aoun
,
E.
, “
Experimental and numerical investigation of a three-dimensional vertical-axis wind turbine with variable-pitch
,”
J. Wind Eng. Ind. Aerodyn.
139
,
111
123
(
2015
).
14.
Fatahian
,
E.
,
Ismail
,
F.
,
Ishak
,
M. H. H.
, and
Chang
,
W. S.
, “
The role of wake splitter deflector on performance enhancement of Savonius wind turbine
,”
Phys. Fluids
34
,
095111
(
2022
).
15.
Franchina
,
N.
,
Persico
,
G.
, and
Savini
,
M.
, “
2d-3d computations of a vertical axis wind turbine flow field: Modeling issues and physical interpretations
,”
Renewable Energy
136
,
1170
1189
(
2019
).
16.
Franchina
,
N.
,
Persico
,
G.
, and
Savini
,
M.
, “
Three-dimensional unsteady aerodynamics of a h-shaped vertical axis wind turbine over the full operating range
,”
J. Wind Eng. Ind. Aerodyn.
206
,
104273
(
2020
).
17.
Fu
,
P.
and
Farzaneh
,
M.
, “
A CFD approach for modeling the rime-ice accretion process on a horizontal-axis wind turbine
,”
J. Wind Eng. Ind. Aerodyn.
98
,
181
188
(
2010
).
18.
Gao
,
L.
,
Tao
,
T.
,
Liu
,
Y.
, and
Hu
,
H.
, “
A field study of ice accretion and its effects on the power production of utility-scale wind turbines
,”
Renewable Energy
167
,
917
928
(
2021
).
19.
Habashi
,
W. G.
,
Aubé
,
M.
,
Baruzzi
,
G.
,
Morency
,
F.
,
Tran
,
P.
, and
Narramore
,
J. C.
, “
FENSAP-ICE: A fully-3d in-flight icing simulation system for aircraft, rotorcraft and UAVS
,” in
24th International Congress of the Aeronautical Sciences, International Council of the Aeronautical Sciences (ICAS)
Bonn, Germany (
2004
), pp.
1
10
.
20.
Hajisharifi
,
A.
,
Marchioli
,
C.
, and
Soldati
,
A.
, “
Particle capture by drops in turbulent flow
,”
Phys. Rev. Fluids
6
,
024303
(
2021
).
21.
Hajisharifi
,
A.
,
Marchioli
,
C.
, and
Soldati
,
A.
, “
Interface topology and evolution of particle patterns on deformable drops in turbulence
,”
J. Fluid Mech.
933
,
A41
(
2022
).
22.
Hann
,
R.
, “
UAV icing: Comparison of LEWICE and FENSAP-ICE for ice accretion and performance degradation
,” AIAA Paper No. 2018-2861,
2018
, p.
2861
.
23.
Hildebrandt
,
S.
and
Sun
,
Q.
, “
Evaluation of operational strategies on wind turbine power production during short icing events
,”
J. Wind Eng. Ind. Aerodyn.
219
,
104795
(
2021
).
24.
Hohman
,
T.
,
Martinelli
,
L.
, and
Smits
,
A.
, “
The effect of blade geometry on the structure of vertical axis wind turbine wakes
,”
J. Wind Eng. Ind. Aerodyn.
207
,
104328
(
2020
).
25.
Jha
,
P.
,
Brillembourg
,
D.
, and
Schmitz
,
S.
, “
Wind turbines under atmospheric icing conditions - ice accretion modeling, aerodynamics, and control strategies for mitigating performance degradation
,” AIAA Paper No. 2012-1287,
2012
.
26.
Kuang
,
L.
,
Lu
,
Q.
,
Huang
,
X.
,
Song
,
L.
,
Chen
,
Y.
,
Su
,
J.
,
Han
,
Z.
,
Zhou
,
D.
,
Zhao
,
Y.
,
Xu
,
Y.
, and
Liu
,
Y.
, “
Characterization of wake interference between two tandem offshore floating vertical-axis wind turbines: Effect of platform pitch motion
,”
Energy Convers. Manage.
265
,
115769
(
2022a
).
27.
Kuang
,
L.
,
Su
,
J.
,
Chen
,
Y.
,
Han
,
Z.
,
Zhou
,
D.
,
Zhang
,
K.
,
Zhao
,
Y.
, and
Bao
,
Y.
, “
Wind-capture-accelerate device for performance improvement of vertical-axis wind turbines: External diffuser system
,”
Energy
239
,
122196
(
2022b
).
28.
Kuang
,
L.
,
Zhang
,
R.
,
Su
,
J.
,
Shao
,
Y.
,
Zhang
,
K.
,
Chen
,
Y.
,
Zhang
,
Z.
,
Tu
,
Y.
,
Zhou
,
D.
,
Han
,
Z
et al, “
Systematic investigation of effect of rotor solidity on vertical-axis wind turbines: Power performance and aerodynamics analysis
,”
J. Wind Eng. Ind. Aerodyn.
233
,
105284
(
2023
).
29.
Larin
,
P.
,
Paraschivoiu
,
M.
, and
Aygun
,
C.
, “
CFD based synergistic analysis of wind turbines for roof mounted integration
,”
J. Wind Eng. Ind. Aerodyn.
156
,
1
13
(
2016
).
30.
Lengani
,
D.
,
Simoni
,
D.
,
Ubaldi
,
M.
, and
Zunino
,
P.
, “
POD analysis of the unsteady behavior of a laminar separation bubble
,”
Exp. Therm. Fluid Sci.
58
,
70
79
(
2014
).
31.
Lennie
,
M.
,
Dominin
,
S.
,
Marten
,
D.
,
Pechlivanoglou
,
G.
, and
Paschereit
,
C. O.
, “
Development of ice throw model for wind turbine simulation software QBlade
,” AIAA Paper No. 2019-1800,
2019
.
32.
Li
,
Y.
,
Tagawa
,
K.
, and
Liu
,
W.
, “
Performance effects of attachment on blade on a straight-bladed vertical axis wind turbine
,”
Curr. Appl. Phys.
10
,
S335
S338
(
2010a
).
33.
Li
,
Y.
,
Wang
,
L.
,
Tagawa
,
K.
, and
Li
,
S.
, “
Computer simulation on the icing accretions on a static straight blade used for the vertical axis wind turbine
,” in
2010 International Conference on Computer Design and Applications
(
IEEE
,
2010b
).
34.
Li
,
Y.
,
Wang
,
S.
,
Liu
,
Q.
,
Feng
,
F.
, and
Tagawa
,
K.
, “
Characteristics of ice accretions on blade of the straight-bladed vertical axis wind turbine rotating at low tip speed ratio
,”
Cold Reg. Sci. Technol.
145
,
1
13
(
2018
).
35.
Li
,
C.
,
Zhu
,
S.
,
Xu
,
Y.
, and
Xiao
,
Y.
, “
2.5d large eddy simulation of vertical axis wind turbine in consideration of high angle of attack flow
,”
Renewable Energy
51
,
317
330
(
2013
).
36.
Ma
,
J.
,
Wang
,
F.
, and
Tang
,
X.
, “
Comparison of several subgrid-scale models for large-eddy simulation of turbulent flows in water turbine
,” in
Fluid Machinery and Fluid Mechanics
(
Springer
,
Berlin Heidelberg
,
2009
), pp.
328
334
.
37.
Manatbayev
,
R.
,
Baizhuma
,
Z.
,
Bolegenova
,
S.
, and
Georgiev
,
A.
, “
Numerical simulations on static vertical axis wind turbine blade icing
,”
Renewable Energy
170
,
997
1007
(
2021
).
38.
Mat Yazik
,
M. H.
,
Chang
,
W. S.
,
Ishak
,
M. H. H.
,
Fatahian
,
E.
, and
Ismail
,
F.
, “
Effect of surface roughness and blade material on the performance of a stationary Savonius wind turbine under different operating conditions
,”
Phys. Fluids
35
,
035133
(
2023
).
39.
Mousavi
,
S. M.
and
Kamali
,
R.
, “
Mathematical modeling of the vortex shedding structure and sound pressure level of a large wind turbine tower
,”
Int. J. Appl. Mech.
12
,
2050070
(
2020
).
40.
Nakakita
,
K.
,
Nadarajah
,
S.
, and
Habashi
,
W.
, “
Toward real-time aero-icing simulation of complete aircraft via FENSAP-ICE
,”
J. Aircr.
47
,
96
109
(
2010
).
41.
Nicoud
,
F.
and
Ducros
,
F.
, “
Subgrid-scale stress modelling based on the square of the velocity gradient tensor
,”
Flow, Turbul. Combust.
62
,
183
200
(
1999
).
42.
Ou
,
M.
,
Yan
,
L.
,
Huang
,
W.
,
Li
,
S.
, and
Li
,
L.
, “
Detailed parametric investigations on drag and heat flux reduction induced by a combinational spike and opposing jet concept in hypersonic flows
,”
Int. J. Heat Mass Transfer
126
,
10
31
(
2018
).
43.
Peng
,
H.
,
Liu
,
H.
, and
Yang
,
J.
, “
A review on the wake aerodynamics of h-rotor vertical axis wind turbines
,”
Energy
232
,
121003
(
2021
).
44.
Pope
,
S. B.
,
Turbulent Flows
(
Cambridge University Press
,
2000
).
45.
Posa
,
A.
, “
Dependence of the wake recovery downstream of a vertical axis wind turbine on its dynamic solidity
,”
J. Wind Eng. Ind. Aerodyn.
202
,
104212
(
2020a
).
46.
Posa
,
A.
, “
Influence of tip speed ratio on wake features of a vertical axis wind turbine
,”
J. Wind Eng. Ind. Aerodyn.
197
,
104076
(
2020b
).
47.
Ramirez
,
J. M.
and
Saravia
,
M.
, “
Assessment of Reynolds-averaged Navier–Stokes method for modeling the startup regime of a Darrieus rotor
,”
Phys. Fluids
33
,
037125
(
2021
).
48.
Rezaeiha
,
A.
,
Montazeri
,
H.
, and
Blocken
,
B.
, “
On the accuracy of turbulence models for CFD simulations of vertical axis wind turbines
,”
Energy
180
,
838
857
(
2019
).
49.
Romanò
,
F.
,
Hajisharifi
,
A.
, and
Kuhlmann
,
H. C.
, “
Cellular flow in a partially filled rotating drum: Regular and chaotic advection
,”
J. Fluid Mech.
825
,
631
650
(
2017
).
50.
Shamsoddin
,
S.
and
Porté-Agel
,
F.
, “
Effect of aspect ratio on vertical-axis wind turbine wakes
,”
J. Fluid Mech.
889
,
R1
(
2020
).
51.
Shao
,
Y.
,
Su
,
J.
,
Tu
,
Y.
,
Kuang
,
L.
,
Han
,
Z.
,
Zhang
,
K.
, and
Zhou
,
D.
, “
Assessment of the aerodynamic benefits of collocating horizontal-and vertical-axis wind turbines in tandem using actuator line model
,”
Phys. Fluids
35
,
075115
(
2023
).
52.
Sheidani
,
A.
,
Salavatidezfouli
,
S.
, and
Schito
,
P.
, “
Study on the effect of raindrops on the dynamic stall of a NACA-0012 airfoil
,”
J. Braz. Soc. Mech. Sci. Eng.
44
,
203
(
2022
).
53.
Sheidani
,
A.
,
Salavatidezfouli
,
S.
,
Stabile
,
G.
, and
Rozza
,
G.
, “
Assessment of URANS and LES methods in predicting wake shed behind a vertical axis wind turbine
,”
J. Wind Eng. Ind. Aerodyn.
232
,
105285
(
2023
).
54.
Siddiqui
,
M. S.
,
Durrani
,
N.
, and
Akhtar
,
I.
, “
Quantification of the effects of geometric approximations on the performance of a vertical axis wind turbine
,”
Renewable Energy
74
,
661
670
(
2015
).
55.
Silva
,
J. E.
and
Danao
,
L. A. M.
, “
Varying VAWT cluster configuration and the effect on individual rotor and overall cluster performance
,”
Energies
14
,
1567
(
2021
).
56.
Son
,
C.
and
Kim
,
T.
, “
Development of an icing simulation code for rotating wind turbines
,”
J. Wind Eng. Ind. Aerodyn.
203
,
104239
(
2020
).
57.
Sotoudeh
,
F.
,
Kamali
,
R.
, and
Mousavi
,
S. M.
, “
Field tests and numerical modeling of invelox wind turbine application in low wind speed region
,”
Energy
181
,
745
759
(
2019
).
58.
Tabib
,
M.
,
Siddiqui
,
M. S.
,
Rasheed
,
A.
, and
Kvamsdal
,
T.
, “
Industrial scale turbine and associated wake development-comparison of rans based actuator line vs sliding mesh interface vs multiple reference frame method
,”
Energy Procedia
137
,
487
496
(
2017
).
59.
Taira
,
K.
,
Brunton
,
S. L.
,
Dawson
,
S. T. M.
,
Rowley
,
C. W.
,
Colonius
,
T.
,
McKeon
,
B. J.
,
Schmidt
,
O. T.
,
Gordeyev
,
S.
,
Theofilis
,
V.
, and
Ukeiley
,
L. S.
, “
Modal analysis of fluid flows: An overview
,”
AIAA J.
55
,
4013
4041
(
2017
).
60.
Thangavelu
,
S. K.
,
Chow
,
S. F.
,
Sia
,
C. C. V.
, and
Chong
,
K. H.
, “
Aeroelastic performance analysis of horizontal axis wind turbine (HAWT) swept blades
,”
Mater. Today
47
,
4965
4972
(
2021
).
61.
Villalpando
,
F.
,
Reggio
,
M.
, and
Ilinca
,
A.
, “
Prediction of ice accretion and anti-icing heating power on wind turbine blades using standard commercial software
,”
Energy
114
,
1041
1052
(
2016
).
62.
Wang
,
H.
,
Yi
,
M.
,
Zeng
,
X.
,
Zhang
,
T.
,
Luo
,
D.
, and
Zhang
,
Z.
, “
A hybrid, self-adapting drag-lift conversion wind energy harvesting system for railway turnout monitoring on the Tibetan plateau
,”
Sustainable Energy Technol. Assess.
46
,
101262
(
2021
).
63.
Xu
,
H. Y.
,
Qiao
,
C. L.
,
Yang
,
H. Q.
, and
Ye
,
Z. Y.
, “
Delayed detached eddy simulation of the wind turbine airfoil s809 for angles of attack up to 90 degrees
,”
Energy
118
,
1090
1109
(
2017
).
64.
Yagmur
,
S.
and
Kose
,
F.
, “
Numerical evolution of unsteady wake characteristics of h-type Darrieus hydrokinetic turbine for a hydro farm arrangement
,”
Appl. Ocean Res.
110
,
102582
(
2021
).
65.
Yirtici
,
O.
,
Sevine
,
T.
,
Tuncer
,
I. H.
, and
Ozgen
,
S.
, “
Wind turbine performance losses due to the ice accretion on the turbine blades
,” AIAA Paper No. 2018-4291,
2018
.
66.
Zare Chavoshi
,
M.
and
Ebrahimi
,
A.
, “
Plasma actuator effects on the flow physics of dynamic stall for a vertical axis wind turbine
,”
Phys. Fluids
34
,
075131
(
2022
).
67.
Zhang
,
T.
,
Elsakka
,
M.
,
Huang
,
W.
,
Wang
,
Z.
,
Ingham
,
D. B.
,
Ma
,
L.
, and
Pourkashanian
,
M.
, “
Winglet design for vertical axis wind turbines based on a design of experiment and CFD approach
,”
Energy Convers. Manage.
195
,
712
726
(
2019a
).
68.
Zhang
,
R.
,
Huang
,
W.
,
Yan
,
L.
,
Chen
,
Z.
, and
Moradi
,
R.
, “
Drag and heat flux reduction induced by the pulsed counterflowing jet with different waveforms on a blunt body in supersonic flows
,”
Acta Astronaut.
160
,
635
645
(
2019b
).
69.
Zhang
,
T.
,
Wang
,
Z.
,
Huang
,
W.
,
Ingham
,
D.
,
Ma
,
L.
, and
Pourkashanian
,
M.
, “
A numerical study on choosing the best configuration of the blade for vertical axis wind turbines
,”
J. Wind Eng. Ind. Aerodyn.
201
,
104162
(
2020
).
Published open access through an agreement with Universidade da Beira Interior