CFD modeling of varying complexity for aerodynamic analysis of H-vertical axis wind turbines

Renewable Energy 145 (2020) 2658e2670

Contents lists available at ScienceDirect

Renewable Energy journal homepage: www.elsevier.com/locate/renene

CFD modeling of varying complexity for aerodynamic analysis of H-vertical axis wind turbines Jiao He, Xin Jin*, Shuangyi Xie, Le Cao, Yaming Wang, Yifan Lin, Ning Wang College of Mechanical Engineering, Chongqing University, Chongqing, 400044, PR China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 18 May 2018 Received in revised form 16 May 2019 Accepted 28 July 2019 Available online 7 August 2019

Computational fluid mechanics (CFD) is considered as an efficient approach for studying aerodynamic characteristics of vertical axis wind turbines (VAWTs). Currently, 2D Unsteady Reynolds-Averaged Naviere Stokes (URANS) is widely applied, although previous researches revealed its limit accuracy in the aerodynamic analysis. This paper investigates the accuracy and feasibility of various CFD modeling techniques, namely 2D URANS, 2.5D URANS, 2.5D large eddy simulations (LES), 3D URANS and 3D LES, in the aerodynamic study of VAWTs through a comparison with the wind tunnel results. Compared with the URANS method, the LES approach can provide more accurate prediction on the aerodynamic performance for VAWTs operating at the dynamic stall. The significant improved simulation results by 2.5D LES imply that the neglect of tip vortices may not be the major mechanism causing the over prediction in 2D and 2.5D URANS. 2.5D LES can be regarded as a promising and efficient approach to investigate the aerodynamic behaviors of VAWTs, considering the compromise between the accuracy and computational cost among 2.5D LES, 3D LES and 3D URANS. Furthermore, considering the huge amount of time consumed by CFD simulations, a hybrid meta-model is therefore proposed to predict the power coefficient of VAWTs. The prediction results show that the accuracy of the hybrid meta-model satisfies the requirements, and the calculation time is also reduced. © 2019 Elsevier Ltd. All rights reserved.

Keywords: H-type VAWT CFD Aerodynamic performance Meta-model

1. Introduction As a clean and renewable energy, wind power has attracted more and more attention all over the world. A large amount of clean electricity provided by wind, which has great ecological and social benefits, can replace the conventional thermal power generation. The most typical device to convert wind energy into electricity is the wind turbine, which generally is divided into horizontal axis wind turbines (HAWTs) and vertical axis wind turbines (VAWTs). Compared with HAWTs, VAWTs have many advantages [1]: (1) no yawing devices; (2) less noise; (3) relatively easy installment and maintenance, and so on. Therefore, many researches have been focused on the VAWTs in recent years. Currently, computational fluid mechanics (CFD) method, which can simulate the complex flow field, is regarded as an important tool for studying the aerodynamic performance of VAWTs. Reviewing literatures suggest that CFD methods can be divided into

* Corresponding author. E-mail address: [email protected] (X. Jin). https://doi.org/10.1016/j.renene.2019.07.132 0960-1481/© 2019 Elsevier Ltd. All rights reserved.

2D CFD, 2.5D CFD and 3D CFD. Due to the constraint of computing resources and time, the 2D Unsteady Reynolds-Averaged NavierStokes (URANS) is widely used to simulate the aerodynamic performance of VAWTs. Nevertheless, compared with the experimental results, the 2D URANS approach generally over predicts the power coefficients (Cp) of VAWTs, although this approach can basically reflect the change trend of the power coefficient [2e5]. Many studies consider the difference is due to the neglect of the blade tip vortices in the VAWTs when the 2D URANS method is used. 2.5D CFD simulations, which are different from the conventional 2D and 3D CFD simulations, are investigated to solve the limitation of the 2D CFD model. 2.5D CFD model is different from 2D CFD model in that a certain height of blade is modeled with a period condition at the top and bottom boundaries. Li et al. [6] simulated the flow field around airfoil by large eddy based on a 2.5D model. Travin et al. [7], Bertagnolio et al. [8] and IM and Zha [9] analyzed the airfoil stall through a 2.5D CFD model. For the 3D CFD simulation, it is still desirable to estimate the aerodynamics of VAWTs due to its high accuracy, although it needs greatly high computational cost. In Ref. [5], it was found that the maximum torque obtained by the 2D SST k-u model is larger than that by the

J. He et al. / Renewable Energy 145 (2020) 2658e2670

Nomenclature VAWTs URANS LES HAWT CFD

l r Ud

U q

4 W Fn Ft

vertical axis wind turbines unsteady Reynolds-averaged Naviere Stokes large eddy simulations horizontal axis wind turbines computational fluid mechanics tip speed ratio rotor radius inflow wind speed rotational speed of the rotor azimuthal angle attack angle relative wind speed normal force tangential force

3D model. The reason for this deviation is due to the influence of the finite aspect ratio and traverse drag in the 3D CFD model. McLaren [2] found that the correction factors of the inlet wind speed and the blade height can be used in the 2D URANS model to account for the 3D flow influences. These researches indicates that the reasonable estimation for the aerodynamic performance of VAWTs could be realized when the finite aspect ratio of blades is appropriately simulated in the 3D CFD model. Nevertheless, several key factors, such as the aerodynamic characteristics when the blade is in stall at a low blade tip speed ratio (TSR) and the abilities of different turbulence models, have not attracted enough attention in most cases, although it is well recognized that these factors may have significant influences on the prediction of the power coefficient of VAWTs. VAWTs frequently operates at a low blade tip speed ratio in their practical application, resulting in the angle of attack (AOA) is larger than the stall angle. The study by Rom [10] showed that the flow over the blade surface at a high TSR is unsteady, highly separated and three dimensional. In Ref. [11], it is found that some factors, such as the selection of turbulence models and domain dimensions, are important to accurately estimate the flow field characteristics at low TSRs. However, these factors are not paid considerable cautions in previous simulations of the aerodynamics of VAWTs. For example, the URANS models frequently used in the analysis of the aerodynamic performance of VAWTs cannot achieve accurate prediction for the low-TSR operating case, in which the large-scale and complex turbulence structures exist in the separated flow. On the other hand, the large eddy simulation (LES), although considered to be a more powerful and advanced technique, is rarely used in the previous studies of VAWTs because it needs a large amount of computational cost. Meanwhile, the shortcoming of the 2D CFD approach is apt to be ignored. For example, the flow diversities in the vertical direction cannot be represented in the 2D CFD models. Therefore, an efficient and accurate simulation scheme for the VAWTs is still being sought. Based on the above analysis, a systematic investigation into CFD modeling is performed in the present paper to study the accuracy and feasibility of the 2D URANS, 2.5D URANS, 3D URANS, 2.5D LES and 3D LES in estimating the aerodynamic characteristics of VAWTs. Furthermore, considering the high computational cost of CFD approaches, particularly of the sophisticated 3D CFD model, it has great significance to find a relatively accurate and efficient way to predict the aerodynamic performance of VAWTs. Therefore, a

Cn Ct Cp H

r T d c b yi yi y MSE R2 RE

2659

normal force coefficient tangential force coefficient power coefficient blade height air density mechanical torque rotor diameter blade chord length simulated data predicted data the average of simulation data mean square error multiple correlation coefficient relative error

prediction meta-model is also investigated in this paper. In view of this, this paper is arranged as follows: In Section 2, the operation principle of the H-type VAWT is introduced; in Section 3, three CFD modeling schemes, namely the 2D CFD, 2.5D CFD and 3D CFD approaches, are investigated in detail, including the determination of the computational domain and boundary, mesh grid independence study, simulation setup, and so on; in Section 4, a test platform of the H-type VAWT is established to evaluate different CFD modeling techniques, the comparison between simulation and the experimental data are discussed, and 2.5D LES is chosen as a promising method; in Section 5, in order to achieve a relatively accurate and efficient approach to predict the power coefficient of the VAWTs, a hybrid meta-model is proposed based on two independent meta-models. Moreover, the effectiveness of the hybrid meta-model is evaluated by a variety of evaluation indices; in Section 6, the conclusions are provided. 2. Operation principle of an H-type VAWT Fig. 1 shows the operation schematic of a single blade of an Htype VAWT. The tip speed ratio (TSR), which is a non-unit value, can

Wt

Wn

D

FC

W

Ud

V

r F L

Fig. 1. The operation schematic of a single blade.

FN

2660

J. He et al. / Renewable Energy 145 (2020) 2658e2670

be defined as [12].

TSR; l ¼

rU Ud

(1)

where r is the rotor radius, Ud is the inflow wind speed, and U is the rotational speed of the rotor. Continuous change of the AOA causes variation of the lift and drag forces acting on the blade. The AOA is defined as follows [12]:

 sin q cos q þ l

 4 ¼ arctan

(2)

where q represents a azimuthal angle of the blade. According to the geometrical relationship shown in Fig. 1, the expression of the relative wind speed, W, can be expressed by:

W ¼ Ud

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðcos q þ lÞ2 þ sin2 q

(3)

The resultant aerodynamic force acting on the blade can be decomposed into two components:the normal force Fn and the tangential force Ft. The two forces are perpendicular and parallel to the blade chord, respectively. The normal force coefficient Cn and the tangential force coefficient Ct are represented as [2]:

Cn ¼

Ct ¼

Fn 0:5rcHU 2d Ft 0:5rcHU 2d

(4)

(5)

For a straight-bladed H-type VAWT, the power coefficient Cp can be used to illustrate its performance [13e17]:

TU Cp ¼ . 1 2rð2RHÞrU 3d

(6)

where H represents the blade height; r represents the air density; T represents the mechanical torque. According to Eq. (2), the change of the AOA is related to the azimuthal angle and the TSR. When the H-type VAWTs operate at a low TSR corresponding to a high AOA, the blades generally experience stall.The results in Ref. [18] demonstrated that the dynamic stall will happen when the VAWTs operated at the TSR less or equal than four, which was confirmed through the test in Ref. [19]. Especially, when the H-type VAWT is in the start-up stage, the TSR is very low, causing the AOA far larger than the stall angle. Accordingly, when evaluating the performance of the VAWTs, the flow analysis for the low TSR is necessary.

and rotating domains are set to interface to ensure the data exchange. The division of the computational domain in this plane is also applicable to the 2.5D and 3D CFD models. 3.1.2. 2.5D CFD model The 2.5D CFD model is obtained by extending the 2D model, to a certain length in the spanwise direction but not the actual blade length. This method is between 2D and 3D approaches. In the 2.5D modeling scheme, the blade height is 3 times the blade chord length c, and the boundary conditions at the top and bottom of the computational domain are set to translational periodic [22], as shown in Fig. 2 (b). 3.1.3. 3D CFD model In the 3D CFD model, the blade height is modeled as the actual value. The height of the calculation domain is set as 8d, and the boundary conditions at the top and bottom of the computational domain are set to symmetry, as shown in Fig. 2 (c). 3.2. Mesh generation In this paper, the CFD pre-processing tool ICEM is used to mesh the computational domain. The sliding mesh technique is adopted to ensure the data exchange between two adjacent flow fields [23]. 3.2.1. 2D CFD model According to the findings of Versteeg and Malalasekera [24] and White [25], the boundary layers can be divided into four parts, namely the viscous sub-layer (yþ  5), the buffer layer (5 < yþ  30), the log-law region (<30 yþ  500), and the outer layer. When generating the mesh, the viscous sub layer near the blade surface should be paid much attention. In order to directly resolve the laminar sub-layer, the first mesh close to the blade surface is refined to ensure yþ less than 1, and the growth rate of mesh size is set to 1.2, as shown in Fig. 3. This mesh setting is also applied to the 2.5D and 3D CFD models. Moreover, the size of grids on the interface boundary of the adjacent fluid domain should be approximately consistent to ensure the continuity, accuracy and fast convergence of CFD calculations [26]. The generated grids of the 2D rotating and blade domains are shown in Fig. 4(a). 3.2.2. 2.5D CFD model The computational grids of the 2.5D model are obtained by extruding the 2D meshes to the 3c height, as shown in Fig. 4(b).

3. CFD modeling techniques

3.2.3. 3D CFD model The grids of 3D CFD model are obtained by stretching the 2D meshes to the actual blade length. The total meshes of the 3D CFD model are shown in Fig. 4 (c).

3.1. Computational domain and boundary condition

3.3. Mesh independency study

3.1.1. 2D CFD model Fig. 2 shows the computational domains for different CFD schemes. The rotating area of a rotor is divided into three independent blade subdomains (Fig. 2 (a)). The data exchange between adjacent fields is realized through an interface boundary condition [20]. According to the proposal of Ferreira et al. [21], the inlet and outlet boundary conditions are, respectively, placed at 10d (d is the rotor diameter) upwind and 15d downwind away from the rotor center to allow full development of the wake. The width of the flow field is set to 10d, and the boundary conditions on both sides are set to symmetry. Meanwhile, boundary conditions between the static

In order to choose the appropriate mesh number to satisfy the calculational accuracy and cost, a mesh independency study is performed. CFD simulations are carried out by using three sets of grid resolutions, i.e. G1, G2 and G3, as shown in Table 1. Fig. 5 shows the mesh independence test for different CFD models. It can be seen that the power coefficients obtained by G2 and G3 are very close. Therefore, G2 is selected to perform the following simulations, so that not only the computation accuracy is ensured, but also the additional calculation cost is not caused. Additionally, the meshes for the LES model is enough sufficient for the URANS model, thus leading to smaller power coefficient difference for the URANS.

J. He et al. / Renewable Energy 145 (2020) 2658e2670

Wind tunnel flow field domain

10d

2661

Symmetry Free-slip wall

Pressure outlet

Velocity inlet Rotating domain Interface

Symmetry

Blade subdomain

10d

15d (a) 2D model

%ODGH WRS translational periodic

+ c ,QOHW c

ERWWRP translational periodic

2XWOHW

(b) 2.5D model

%ODGH + G

WRS

Symmetry

,QOHW 2XWOHW

G ERWWRP

Symmetry

(c) 3D model Fig. 2. Computational domains for different CFD schemes.

Fig. 3. Meshes around the blade.

3.4. Simulation setup 3.4.1. Turbulence model The key in CFD simulation is to solve the turbulence generated by fluid flow. Currently, the turbulence numerical models used in engineering applications include the URANS and large eddy models. The SST k-u model based on the URANS outperforms in the aerodynamic performance prediction because it can effectively mix the robust k-u model in the near wall range and convert to the standard k-3 model in the far flow filed to resolve the complex flow with the adverse gradient [27e29]. SST k-u model is a two-equation eddy-viscosity model. The utilization of a k-u formulation in the inner parts of the boundary layer makes the model directly usable all the way down to the wall by the viscous sub-layer; therefore the SST k-u model can be adopted as a Low-Re turbulence model without any additional damping functions. The findings in Ref. [30] showed that the results simulated by the SST k-u model were in good agreement with the experimental data while the numerical stability, reliability, and flexibility were kept in performing the near wall resolution. Similar conclusions were also obtained in Ref. [29], suggesting that the results obtained by the SST k-u approach are closest to the test results.

On the contrary, the conservation equations in the LES approach are averaged in space not in time. The LES model can directly solve the large and boundary-dependent eddy by the governing equation, and can consider the effects of small and more homogenous eddy through a sub-grid model [31], which costs considerable computing resources. Nevertheless, because the LES model can maintain the unsteady large-scale coherent structures, it can be applied to a wide scope of turbulent flows. Some studies have indicated that the LES model is promising in the simulation of the unsteady separated flows [32e34]. The LES model can directly solve the large and boundary-dependent eddy by the governing equation, and can consider the effects of small and more homogenous eddies through a sub-grid model, which costs considerable computing resources. Nevertheless, because the LES model can maintain the unsteady large-scale coherent structures, it can be applied to a wide scope of turbulent flows. The Smagorinsky-Lilly model overcomes the limitation of the original Smagorinsky subgrid model [35] and is selected in the present work. Therefore, the SST k-u and LES models are selected in this work to evaluate their influences on the performance of the H-type VAWT. Because the 2D LES approach exhibits an obvious inconsistent performance under various situations, the LES method is only examined in 2.5D and 3D simulations, and the SST k-u approach is tested in 2D, 2.5D and 3D CFD simulations. 3.4.2. Solver setting The SIMPLE algorithm is used to solve the pressure velocity coupling for accelerating the convergence of CFD models. Due to the poor initial guesses, the first order upwind discrete scheme is chosen in the first revolution to achieve the properly developed fluid flow. The second-order and the third-order upwind discrete schemes are then adopted for the URANS and LES, respectively. In

2662

J. He et al. / Renewable Energy 145 (2020) 2658e2670

(a) 2D CFD model

(b) 2.5D CFD model

(b) 2.5D CFD model

(c) 3D CFD model Fig. 4. Meshes for different CFD models.

Table 1 Cell number versus different CFD models. CFD model

2D model 2.5D model 3D model

Cell number G1

G2

G3

132373 3196100 5843560

159805 4794150 7990250

181188 6392200 10015960

addition, the double precision segregated solver with an implicit method is selected to resolve the discretized algebraic equations. Moreover, the time step is an important factor that affects the accuracy and efficiency of CFD numerical simulations [36]. The independence of the time step is studied applying 0.25 , 0.5 , 1 and 2 rotor rotation per time step for the LES and URANS models, as shown in Fig. 6. The results show that there is no great

improvement between 0.25 and 0.5 , so 0.5 is chosen to enhance the calculation efficiency for the two models. In Refs. [37e39], the time step was set as 1 or 2 , indicating that the time step of 0.5 in our research is considered appropriate.

4. Experiment setup and result discussions 4.1. Experiment setup To validate the accuracy of CFD models, a torque test platform was built at College of Mechanical Engineering of Chongqing University. Fig. 7 shows the experimental wind tunnel, testing platform and control device. The experimental wind tunnel with a length of 6 m and a cross section of 4.5  4.5m consists of nine 80 kW fans. The honeycomb with a diameter of 0.9 cm and a thickness of 2.7 cm

J. He et al. / Renewable Energy 145 (2020) 2658e2670

0.45

2663

0.35 G1

G2

G3

G1

G2

G3

Power coefficient

Power coefficient

0.4

0.35

0.3

0.3

0.25

0.25

0.2

0.2 2D

2.5D

3D

2.5D

(a) URANS

3D (b) LES

0.38

0.4

0.36

0.38

Poweer coefficient

Poweer coefficient

Fig. 5. Mesh independence test for different CFD models.

0.34 0.32 0.3

0.36 0.34 0.32 0.3

0.28

0.28

0.26 0

0.5

1 1.5 Time step (e)

2

2.5

0

(a) LES

0.5

1 1.5 Time step (e)

2

2.5

(b) URANS

Fig. 6. Mean power coefficient for one revolution against different time steps.

(a) Wind tunnel

(b) H-type VAWT

(c)Operation platform

Fig. 7. Experimental platform.

is adopted to reduce the flow fluctuations. The wind speed scope is 2.0e12.5 m/s. In this work, a small vertical axis wind turbine with three blades of NACA0021 airfoil is established. The rotor diameter is 1030 mm, and the blade height h and chord length c are 1456.4 mm and 85.8 mm, respectively. The connecting location of the pitching shaft is almost at 1/4 of chord from the leading edge. The H-type VAWT is mounted in the downwind area with relatively uniform wind field, and 6.2 m away from the vane axial fans. The wind speed is measured by a digital anemometer (precision is±5%). The speed of the rotor, which is connected with a torque sensor used to measure the rotational speed and torque, is controlled by a magnetic powder load, as shown in Fig. 8. The wind speed is recorded using a non-contact digital tachometer (the minimum count is 1/rpm).

4.2. Results and discussions 4.2.1. Performance comparison Fig. 9 shows the mean power coefficients at various TSRs and turbulence models for one revolution. The mean values of the experimental Cp for different TSRs are also shown in the figure. It can be found that all CFD modeling techniques can approximately replicate the change trend of the power coefficient, namely the Cp increases with the increase of the TSR until reaches the maximum, and then decreases with the increase of the TSR. 2D and 2.5D URANS approaches overestimate Cp, which is in agreement with the findings in Refs. [2e5]. 2.5D LES, 3D LES and 3D URANS can provide a good estimation of the power coefficient. However, for the low TSR cases, 2.5D LES and 3D LES show better prediction than the 3D URANS method. For the high TSRs, 2.5D LES, which although

2664

J. He et al. / Renewable Energy 145 (2020) 2658e2670

Wind speed sensor

z

only shows fair agreement with the experimental data, still exhibits much better performance than the 2D and 2.5D URANS approaches. It needs to be noted that considering the less modeling complexity of 2.5D LES than that of 3D LES, this finding indicates that 2.5D LES can be used as a promising and efficient approach to predict the power coefficient of the VAWT, especially applicable for the low TSR case corresponding to the start-up stage of VAWTs. Some studies have revealed the limitation of 2D URANS models in simulating the aerodynamic performance of VAWTs. However, the limitation of 2D URANS approaches is usually thought to that the blade tip vortices are not simulated [2,5,40,41]. Because none of three approaches, namely the 2D URANS, 2.5D URANS and 2.5D LES, can properly simulate the blade tip vortices, the overestimation of the power coefficient of the H-type VAWT by the 2D URANS and 2.5D URANS methods is probably due to the inherent limitation of the URANS approach in simulating the lift force beyond stall [6]. Scheurich et al. [42] also found the limited influence of tip vortices on the wind turbine performance by the vorticity transport model. Therefore, although the results obtained by 2.5D URANS are slightly improved than those by 2D URANS, the 2.5D URANS is also limited due to the insufficient ability to resolve the flow containing large separations for the low TSRs. Fig. 10 shows the changes of the tangential force coefficient (Ct) and normal force coefficient (Cn) of a single blade with the variation of azimuthal angles for various CFD schemes when TSR equals 1.5 and 2. Cn and Ct are calculated by Eq. (4) and Eq. (5), respectively. Chandramouli et al. [43] found that the forces experienced by blades were similar except for the existing phase differences among the forces. Therefore, only the tangential and normal force coefficients of a single blade are plotted in Fig. 10. The change of Ct can be confirmed by observing the flow behaviors around the blade at various azimuthal angles (shown in Fig. 11), which will be discussed next. Similar to findings in other studies [2e4,6], 2D and 2.5D

y

Rotor

U

x

o

Artificial wind source

Wind tunnel

Converter

Torque sensor

Controller

Magnetic powder load Fig. 8. Schematic diagram of the wind tunnel test.

Experiment 3D URANS 3D LES 2.5D URANS 2.5D LES 2D URANS

0.3 0.2 0.1 0.0 1.4

1.6

1.8

2.0

2.2

2.4

2.6

Tip speed ratio

1.5

2D SST 2.5D SST 3D SST 2.5D LES 3D LES

1.2 0.9

Tangential force coefficient

Tangential force coefficient

Fig. 9. Comparison of experimental and simulated Cp at different TSRs.

0.6 0.3 0.0 -0.3 0

45

90

135

180

225

270

315

2D SST 2.5D SST 3D SST 2.5D LES 3D LES

5 4 3 2 1 0 -1 0

360

45

90

135

180

225

270

315

360

Azimuthal angle ( )

Azimuthal angle ( )

(a) TSR=2 2D SST 2.5D SST 3D SST 2.5D LES 3D LES

1.2 0.9 0.6

Tangential force coefficient

Tangential force coefficient

Power coefficient

0.4

0.3 0.0 -0.3 -0.6 0

45

90

135

180

225

270

315

360

2D SST 2.5D SST 3D SST 2.5D LES 3D LES

6 4 2 0 -2 0

Azimuthal angle ( )

45

90

135

180

225

270

Azimuthal angle ( )

(b) TSR=1.5 Fig. 10. Tangential and normal force coefficients of a single blade against the azimuthal angle for various CFD schemes.

315

360

J. He et al. / Renewable Energy 145 (2020) 2658e2670

URANS in this work exhibit larger peak amplitude of the tangential force coefficient, resulting in significantly overestimation of the power coefficient, as shown in Fig. 9. This is essentially because of the overestimation of the lift forces for the low TSR cases by 2D and 2.5D URANS approaches. The maximum Ct obtained by the LES approach appears earlier than that by the URANS method, which is because the flow separation happens at a smaller azimuthal angle for the LES approach, as shown in Fig. 11(a). In addition, it also can be found that the azimuthal angle corresponding to the peak value of Ct decreases with the decrease of the TSR, which is because that stall is more likely to occur at the low TSR corresponding to a large AOA. For the normal force coefficient Cn, the differences in all CFD methods are relatively small. 4.2.2. Flow mechanism In order further insight into the flow development of the VAWT, the analysis of vorticity and pressure is needed. Fig. 11 shows the vorticity and pressure contours around the mid-plane of a single blade at different azimuthal angles, i.e., q ¼ 0 , 30 , 60 , 90 , 120 , 150 , 180 , for TSR ¼ 1.5. In Fig. 11 (a), the shedding vortices in counter-clockwise directions are represented by red color. Because the flow field in the mid-plane is almost undisturbed by the tip vortices, it can be seen from Fig. 11 that when the same URANS model is applied to the 2D, 2.5D and 3D CFD schemes respectively, differences of the vorticity and pressure field characteristics are relatively small. However, the flow characteristics of the LES model are quite different from those of the URANS models. For URANS approaches, the flow attaches to the blade surface when the blade operates between azimuthal angles of 0 and approximately 60 .

2665

Then, flow separation starts to happen for URANS models, whereas for the LES models, the flow separation occurs at a smaller azimuthal angle (approximately 30 ), resulting in different azimuthal angles corresponding to the maximum tangential force coefficient (Ct), as shown in Fig. 10. Afterwards, the vortices shed from the blade leading edge area and gradually flow through the blade surface, resulting in the change of pressure around the blade, as shown in Fig. 11 (b). The lift force therefore decreases significantly, resulting in the great reduction of Ct. It can also be seen from Figs. 11(b) and Fig. 10 that the larger the pressure difference between the two surfaces of the blade, the larger the tangential force coefficient. Moreover, significant discrepancy in the vorticity and pressure fields simulated by the URANS and LES approaches is observed, which can be used to explain the various peak values of the tangential force coefficient estimated by the two approaches. LES approach generates the obvious shedding vortices from the blade leading edge, and then the vortices develop through the blade surface. This was also found in the particle image velocimetry results in Ref. [44]. Just as the name of URANS illustrates, this approach uses the Reynolds time-averaged treatment to simulate the turbulence, where the Reynolds stress term, treated as the energy transferred in the random turbulence variation, is unresolved but simulated by the time-averaged term [6]. Therefore, it is unsurprising that the URANS approach cannot capture the dynamic stall vortices at a low TSR, thus leading to the over prediction of the output power. To further investigate the development of three-dimensional flows, the Q criterion approach is used to identify and visualize

Fig. 11. Comparison of vorticity and pressure contours for 2D Vs 2.5D Vs 3D CFD schemes at the mid-plane.

2666

J. He et al. / Renewable Energy 145 (2020) 2658e2670

the complex vortex structure [45], where Q denotes the second invariant of velocity gradient tensor. Fig. 12 shows a transient vortical field from 2.5D URANS, 2.5D LES, 3D URANS and 3D LES for TSR ¼ 1.5. It can be seen that the vortex structure simulated by the URANS approach is less sophisticated than that by the LES method due to the inherent treatment of the Navier-Stokes equation. When the dynamic stall occurs, the shedding vortices have an important effect on the power output of the VAWT. As illustrated in Fig. 12, the separated columnar vortices in 2.5D URANS indicate a significant pressure variation but a relatively uniform pressure profile along the blade spanwise direction, which probably can be used to explain why the aerodynamic forces predicted by 2D and 2.5D

Blade 1 =0°

The separated columnar vortex

URANS are close when stall happens. Compared with the URANS simulation, the LES approach could capture the essential structures of the three-dimensional flow in the stall region. In addition, the large-scale vortices have a great effect on the transfer of the momentum and energy in the turbulence flow, which also affect the pressure distribution, resulting in different aerodynamic forces [46], as shown in Fig. 10. Thus, the LES model has stronger ability to represent the three-dimensional flow structures than URANS. Comparing Compared to 3D URANS, 3D LES can better capture the details of separated vortices, resulting in more accurate simulation results than 3D URANS for the low TSRs. However, it needs to be noted that considering the less modeling complexity of the 2.5D

Blade 1 =0°

Blade 1 =30°

The separated columnar vortex

Blade 1 =30°

Blade 3 =270°

Blade 3 =270° The separated columnar vortex

The separated columnar vortex

Blade 3 =240°

The separated columnar vortex

The separated columnar vortex

Blade 3 =240°

Blade 2 =120°

Blade 2 =120° Blade 2 =150°

Blade 2 =150°

Blade 3 =330°

The separated columnar vortex

Blade 3 =330° Blade 3 =300°

The separated columnar vortex

Blade 3 =300°

The separated columnar vortex

The separated columnar vortex

The separated columnar vortex

Blade 1 =60°

Blade 1 =60°

Blade1 =90°

Blade1 =90°

Blade 2 =210°

Blade 2 =180°

Blade 2 =210°

Blade 2 =180°

(a) 2.5D SST k-

(b) 2.5D LES

Blade 1 =30°

Blade 1 =0°

Blade 1 =0°

The separated columnar vortex

Blade 3 =270°

Blade 1 =30°

The separated columnar vortex Blade 3 =270°

The separated columnar vortex

The separated columnar vortex

Blade 3 =240°

Blade 3 =240°

Blade 2 =120°

Blade 2 =120°

Blade 2 =150°

Blade 2 =150°

Blade 3 =330°

Blade 3 =300°

Blade 3 =330°

Blade 3 =300°

The separated columnar vortex

The separated columnar vortex

The separated columnar vortex

The separated columnar vortex

Blade 1 =60° Blade1 =90°

Blade 1 =60° Blade1 =90° Blade 2 =180°

Blade 2 =180°

Blade 2 =210°

Blade 2 =210°

(c) 3D SST k-

(d) 3D LES

Fig. 12. Comparison of iso-surfaces of Q-criterion for 2.5D Vs 3D CFD models (colored by the pressure).

J. He et al. / Renewable Energy 145 (2020) 2658e2670

Testing and validation data set

Simulation model

Training data set

Fitting

2667

Test and validation

Ploynomial Regression metamodel

Hybrid metamodel M

Residual data set Fitting Neural network metamodel Fig. 13. Scheme of composition meta-model.

Table 2 Training data set. TSR

V ¼ 5 m/s

V ¼ 5.5 m/s

V ¼ 6 m/s

V ¼ 6.5 m/s

V ¼ 7 m/s

V ¼ 7.5 m/s

V ¼ 8 m/s

V ¼ 8.5 m/s

V ¼ 9 m/s

0.01639 0.10153 0.14092 0.20336 0.22271 0.14914 0.11369 0.05946

0.01745 0.10994 0.15172 0.21071 0.22534 0.15339 0.12018 0.06216

0.01821 0.12293 0.16798 0.22423 0.24512 0.15708 0.12771 0.07397

0.02037 0.13205 0.17173 0.22932 0.24871 0.16917 0.13672 0.07836

0.02092 0.13891 0.19516 0.26125 0.25651 0.17331 0.14092 0.08528

0.02196 0.14524 0.20318 0.27152 0.26119 0.17379 0.14365 0.09186

0.02282 0.15644 0.20949 0.28072 0.26866 0.17775 0.15251 0.09509

0.02372 0.15791 0.22167 0.28839 0.28697 0.19786 0.16093 0.10319

Cp 1.4 2.0 2.3 2.5 2.8 3.1 3.6 4.0

0.01479 0.09396 0.13377 0.19752 0.21189 0.1421 0.10876 0.04956

Table 3 Validating and test data set. Data classification

TSR

V ¼ 5 m/s

V ¼ 5.5 m/s

V ¼ 6 m/s

V ¼ 6.5 m/s

V ¼ 7 m/s

V ¼ 7.5 m/s

V ¼ 8 m/s

V ¼ 8.5 m/s

V ¼ 9 m/s

0.04355 0.23449 0.14603 0.10997

0.04503 0.23562 0.14938 0.11634

0.04777 0.25811 0.15377 0.12382

0.05262 0.26106 0.16639 0.13237

0.05418 0.27009 0.16965 0.13466

0.05842 0.27749 0.16856 0.14013

0.06172 0.28588 0.17401 0.14755

0.06545 0.30331 0.19369 0.15656

Cp Validation and Test data

1.7 2.7 3.2 3.7

0.04046 0.22393 0.13754 0.10308

model than that of the 3D model, 2.5D LES is selected as a promising and efficient approach to analyze the aerodynamic forces of VAWTs, especially applicable for the start-up case of VAWTs.

adequately solved by a standard curve fitting algorithm, a novelty hybrid meta-model is therefore constructed.

5. Performance prediction based on a hybrid meta-model for the VAWT

5.1. Hybrid meta-model establishment

Based on the analysis in Section 4, it can be known that CFD simulations cost a huge amount of computing resources. Especially, when CFD simulations with a large number of cells and a significantly complex turbulence model is performed, the simulation work will become extremely difficult. In addition, the extra CFD simulations of VAWTs are needed if the wind speed and tip speed ratio change, resulting in the increase of workload and cost. Therefore, a prediction model will be proposed in this section to efficiently and cost-effectively estimate the performance of the VAWT. However, considering the strong nonlinearity among the wind speed, tip speed ratio and power coefficient, which cannot be

The basic ideology of the hybrid meta-model is: firstly, several appropriate meta-models are chosen according to the application purpose; then, the meta-models are combined by one certain method; finally, the hybrid meta-model is tested and verified. In this work, the specific steps are as follows: (1) Run the simulation model and get two sets of input-output data. One set includes the training data, and the other set includes the verification and test data. (2) The training data set is used to train and fit the polynomial regression meta-model, and then the residual data set can be achieved through comparing the predicted data obtained by

J. He et al. / Renewable Energy 145 (2020) 2658e2670

6%

6%

4%

4%

5%

2%

2%

4%

0% -2%

0% -2%

-4%

-4%

-6%

-6% 1.7

2.7

3.2

Error rate

6%

Error rate

Error rate

2668

0% 2.7

(a) V=5m/s

3.7

1.7

(a) V=5.5m/s

4%

4%

4%

2%

2%

2%

-2%

0% -2%

-4%

-4%

-6%

-6% 2.7

3.2

Error rate

6%

0%

3.7

2.7

3.2

3.7

1.7

(a) V=6.5m/s

(a) V=7m/s

4%

4%

4%

2%

Error rate

Error rate

6%

0% -2% -4% -6%

2.7

3.2

TSR (a) V=8m/s

3.7

2% 0% -2% -4%

1.7

3.7

3.2

(a) V=7.5m/s

6%

-4%

3.7

TSR

6%

1.7

2.7

TSR

-2%

3.2

0% -2% -6%

1.7

0%

3.7

-4%

TSR

2%

3.2

(a) V=6m/s

6%

1.7

2.7

TSR

6%

Error rate

Error rate

3.2

TSR

TSR

Error rate

2% 1%

1.7

3.7

3%

2.7

3.2

3.7

1.7

TSR (a) V=8.5m/s

2.7

TSR (a) V=9m/s

7E-05 6E-05

4E-05

R2

MSE

5E-05

3E-05 2E-05 1E-05 Hybrid meta-model

Neural network meta-model

Polynomial regression meta-model

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 Hybrid meta-model

Neural network meta-model

Polynomial regression meta-model

Maximum absolute value of RE (%)

Fig. 14. Result comparison of simulation and prediction data for test and validation data set.

7 6 5 4 3 2 1 0 Hybrid meta-model

Neural network meta-model

Polynomial regression meta-model

Fig. 15. Performance comparison between the hybrid meta-model and single meta-model.

the polynomial regression meta-model with the training data. (3) Neural network meta-model is trained and fitted using the residual data set. (4) The hybrid meta-model is evaluated by the test and verification data set for satisfying the given requirement.

n P 2

R ¼1 

(8) ðyi  yÞ2

i¼1

RE ¼ 5.2. Prediction analysis In order to evaluate the prediction accuracy of the hybrid metamodel, three statistics, i.e. mean square error (MSE), multiple correlation coefficient (R2) and relative error (RE), are used. The specific expressions are defined as follows: n 1X 2 ðyi  b yiÞ n

ðyi  b y i Þ2

yis the average of simulation data.

The specific process is shown in Fig. 13.

MSE ¼

i¼1 n P

(7)

i¼1

where, b y i and yi are the simulated and predicted data, respectively.

b y i  yi  100% yi

(9)

The simulation data obtained by 2.5D LES are used as the training, validating and test data of the hybrid meta-model, as shown in Table 2 and Table 3. The input parameters are the tip speed ratio and the wind speed, and the outputs are the power coefficients. Fig. 14 shows the comparison of the simulated and predicted values for the verification and test data set. It can be seen that the maximum relative error is about 6%. To illustrate the significant prediction performance of the designed hybrid meta-model, three statistics obtained by the hybrid meta-model, neural network meta-model and polynomial regression meta-model are showed in

J. He et al. / Renewable Energy 145 (2020) 2658e2670

2669

Fig. 15. It can be found that the performance of the hybrid metamodel is the best for the three evaluation indices, indicating the rationality of the proposed hybrid meta-model.

(approval No.: 51005255). The authors appreciate the anonymous reviewers for their valuable comments, which are helpful to improve this paper.

6. Conclusions

References

In this paper, the accuracy and feasibility of 2D, 2.5D and 3D CFD modeling approaches in the aerodynamic performance of the Htype vertical axis wind turbine are investigated. The results are then compared with the experiment data. The comparisons show that the 2D URANS approach largely overestimates the power coefficient of the wind turbine because the 2D Navier-Stokes equation cannot accurately represent the characteristics of vorticity diffusion. For 2.5D URANS, this approach can generate slightly improved simulation results than 2D URANS, but is still limited by the ability of URANS to simulate the large flow separations at a low tip speed ratio. By contrast, 2.5D LES, 3D LES and 3D URANS can generate much improved predictions. Previous studies have found the shortcoming of the 2D URANS approach in estimating the aerodynamic characteristics of the vertical axis wind turbine. Moreover, these studies usually attribute this shortcoming to the neglect of the blade tip vortices in 2D CFD simulations. Nevertheless, the significantly improved estimation in 2.5D LES approach indicates that this influence factor might not be the dominant mechanism that leads to the over prediction of the power coefficient obtained by using 2D URANS approaches. A detailed investigation into the flow field of the wind turbine demonstrates that URANS methods overestimate the tangential force coefficient and delay the dynamic stall compared with LES models. By contrast, LES approach can well reproduce the dynamic flow behaviors. The inherent deficiency of the URANS approach can be used to explain the over prediction of the power coefficient. Therefore, the predictive capability of 2.5D LES is better than that of both 2D and 2.5D URANS models, especially for the low TSR cases corresponding to the startup stage of the wind turbine. Additionally, considering the complex modeling process and a large number of grids of 3D LES, although 2.5D LES can only provide fair simulation results at the high TSR, this approach can still be regarded as an efficient and promising approach to study the aerodynamic performance of the vertical axis wind turbine, especially for the investigation into the self-starting performance of the turbine. Furthermore, it should be pointed out that CFD simulation time will be effectively reduced with the development of the hardware and parallel computing ability of computers. However, for the small research institutes or underfunded research institutes, they still encounter great difficulties in dealing with CFD simulations with a huge amount of grids or significantly complex turbulence models. Therefore, this paper proposes a prediction model based on the hybrid meta-model to estimate the power coefficient of the vertical axis wind turbines. Considering the strong nonlinearity among the blade tip speed ratio, wind speed and power coefficient, two independent meta-models are innovatively combined to form a hybrid meta-model, which overcomes the shortage of a single meta-model in the strong nonlinear fitting. The results show that the hybrid meta-model proposed in this paper not only has powerful ability in the highly nonlinear fitting, but also greatly shortens the simulation time. In addition, the hybrid-metal model constructed in this work not only can be used for the prediction of the power coefficient of vertical axis wind turbines, more importantly, but also can provide an engineering experience for estimating other complex and time-consuming simulations. Acknowledgment This paper is sponsored by Natural Science Foundation of China

[1] M. Islam, D.S.K. Ting, A. Fartaj, Aerodynamic models for Darrieus-type straight-bladed vertical axis wind turbines, Renew Sust Energ 12 (4) (2008) 1087e1109. Rev. [2] K.W. McLaren, A Numerical and Experimental Study of Unsteady Loading of High Solidity Vertical axis Wind Turbines, McMaster University, 2011. [3] M.R. Castelli, A. Englaro, E. Benini, The Darrieus wind turbine: proposal for a new performance prediction model based on CFD, Energy 36 (8) (2011) 4919e4934. [4] R. Howell, N. Qin, J. Edwards, N. Durrani, Wind tunnel and numerical study of a small vertical axis wind turbine, Renew. Energy 35 (2) (2010) 412e422. [5] M. Mukinovi, G. Brenner, A. Rahimi, Analysis of vertical axis wind turbines, in numerical & experimental fluid mechanics VII, NNFM 112 (2010) 587e594. [6] C. Li, S. Zhu, Y.L. Xu, et al., 2.5D large eddy simulation of vertical axis wind turbine in consideration of high angle of attack flow, Renew. Energy 51 (2013) 317e330. [7] A. Travin, M. Shur, M. Strelets, P.R. Spalart, Physical and numerical upgrades in the detached-eddy simulation of complex turbulent flows, Adv. LES Complex. Flows (2004) 239e254. [8] F. Bertagnolio, N.N. Sørensen, J. Johansen, Profile Catalogue for Airfoil Sections Based on 3D Computations. Risø-R-1581(EN). Roskilde, Risø National Laboratory, Denmark, 2006. [9] ZhaG. IMH, Delayed detached eddy simulation of a stall flow over NACA0012 airfoil using high order scheme, in: 49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2011. Orlando, Florida. [10] J. Rom, High Angle of Attack Aerodynamics, Springer, New York, 1992. [11] R.M. Cummings, J.R. Forsythe, S.A. Morton, K.D. Squires, Computational challenges in high angle of attack flow prediction, Prog. Aerosp. Sci. 39 (5) (2003) 369e384. [12] Q.a. Li, T. Maeda, Y. Kamada, J. Murata, K. Furukawa, M. Yamamoto, The influence of flow field and aerodynamic forces on a straight-bladed vertical axis wind turbine, Energy 111 (2016) 260e271. [13] P. Marsh, D. Ranmuthugala, I. Penesis, G. Thomas, The influence of turbulence model and two and three-dimensional domain selection on the simulated performance characteristics of vertical axis tidal turbines, Renew. Energy 105 (2017) 106e116. [14] A. Orlandi, M. Collu, S. Zanforlin, A. Shires, 3D URANS analysis of a vertical axis wind turbine in skewed flows, J. Wind Eng. Ind. Aerodyn. 147 (2015) 77e84. [15] Y.-T. Lee, H.-C. Lim, Numerical study of the aerodynamic performance of a 500W Darrieus-type vertical-axis wind turbine, Renew. Energy 83 (2015) 407e415. [16] E. Sobhani, M. Ghaffari, M.J. Maghrebi, Numerical investigation of dimple effects on darrieus vertical axis wind turbine, Energy 133 (2017) 231e241. [17] Q.a. Li, T. Maeda, Y. Kamada, J. Murata, T. Kawabata, K. Shimizu, T. Ogasawara, A. Nakai, T. Kasuya, Wind tunnel and numerical study of a straight-bladed vertical axis wind turbine in three-dimensional analysis (Part I: for predicting aerodynamic loads and performance), Energy 106 (2016) 443e452. [18] A. Laneville, P. Vittecoq, Dynamic stall: the case of the vertical axis wind turbine, J. Sol. Energy Eng. 108 (1986) 141. [19] Robert E. Akins, Dale E. Berg, W Tait Cyrus, Measurements and Calculations of Aerodynamic Torques for a Vertical-axis Wind Turbine, Sandia National Laboratories, 1987. [20] Fluent Boundary Conditions. Fluent Inc..[EB/OL], 2016, 04-05]. Giahi MH, Dehkordi AJ, https://lyle.smu.edu/me/7337/fluent_bcs.pdf. [21] Ferreira C J S, Bijl H, Bussel G V. Simulating dynamic stall in a 2D VAWT: modeling strategy, verification and validation with particle image velocimetry data. J. Phys. Conf. Ser., 2007(75). [22] F. Bertagnolio, N.N. Sørensen, J. Johansen, Profile Catalogue for Airfoil Sections Based on 3D Computations. Risø-R-1581(EN). Roskilde, Risø National Laboratory, Denmark, 2006. [23] F. Trivellato, M. Raciti Castelli, On the CouranteFriedrichseLewy criterion of rotating grids in 2D vertical-axis wind turbine analysis, Renew. Energy 62 (2014) 53e62. [24] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics, Pearson Education Limited, England, UK, 2007, 978-0-13-127498-3. [25] F.M. White, Fluid Mechanics, The McGraw-Hill Companies, Inc., New York, 2009. [26] Fluent Inc., Fluent User's Manual, 2015. [27] H.F. Lam, H.Y. Peng, Study of wake characteristics of a vertical axis wind turbine by two- and three-dimensional computational fluid dynamics simulations, Renew. Energy 90 (2016) 386e398. [28] G. Bedon, S. De Betta, E. Benini, A computational assessment of the aerodynamic performance of a tilted Darrieus wind turbine, J. Wind Eng. Ind. Aerodyn. 145 (Supplement C) (2015) 263e269. [29] A.M. Chowdhury, H. Akimoto, Y. Hara, Comparative CFD analysis of Vertical Axis Wind Turbine in upright and tilted configuration, Renew. Energy 85 (2016) 327e337.

2670

J. He et al. / Renewable Energy 145 (2020) 2658e2670

[30] F. Balduzzi, A. Bianchini, R. Maleci, G. Ferrara, L. Ferrari, Critical issues in the CFD simulation of Darrieus wind turbines, Renew. Energy 85 (2016) 419e435. [31] J. Smagorinsky, General circulation experiments with the primitive equations. I. The basic experiment, Mon. Weather Rev. 91 (3) (1963) 99e164. [32] R.M. Cummings, J.R. Forsythe, S.A. Morton, K.D. Squires, Computational challenges in high angle of attack flow prediction, Prog. Aerosp. Sci. 39 (5) (2003) 369e384. [33] H. Gao, H. Hu, Z.J. Wang, Computational study of unsteady flows around dragonfly and smooth airfoils at low Reynolds numbers, in: 46th AIAA Aerospace Sciences Meeting and Exhibit, Nevada, Reno, 2008. [34] C. Li, Q.S. Li, S.H. Huang, J.Y. Fu, Y.Q. Xiao, Large eddy simulation of wind loads on a long-span spatial lattice roof, Wind Struct. 13 (1) (2010) 57e82. [35] M. Germano, U. Piomelli, P. Moin, W.H. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A Fluid Dyn. 3 (1991) 1760. [36] D.K. Lilly, A proposed modification of the germano-subgrid-scale closure method, Phys. Fluids A Fluid Dyn. 4 (3) (1992) 633e5. [37] K.H. Wong, W.T. Chong, S.C. Poh, 3D CFD simulation and parametric study of a flat plate deflector for vertical axis wind turbine, Renew. Energy 129 (2018) 32e55. [38] B. Shahizare, N. Nik-Ghazali, W.T. Chong, et al., Novel investigation of the different Omni-direction-guide-vane angles effects on the urban vertical axis wind turbine output power via three-dimensional numerical simulation, Energy Convers. Manag. 117 (2016) 206e217.

[39] A. Subramanian, S.A. Yogesh, H. Sivanandan, et al., Effect of airfoil and solidity on performance of small scale vertical axis wind turbine using three dimensional CFD model, Energy 133 (2017) 179e190. [40] M.H. Mohamed, Performance investigation of H-rotor Darrieus turbine with new airfoil shapes, Energy 47 (2012) 522d530. [41] R. Lanzafame, S. Mauro, M. Messina, 2D CFD modeling of H-Darrieus wind turbines using a transition turbulence model, Energy Procedia 45 (2014) 131e140. [42] F. Scheurich, T.M. Fletcher, R.E. Brown, The influence of blade curvature and helical blade twist on the performance of a vertical-axis wind turbine, in: 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2010 (Orlando, Florida, USA). [43] S. Chandramouli, T. Premsai, P. Prithviraj, V. Mugundhan, R.K. Velamati, Numerical analysis of effect of pitch angle on a small scale vertical Axis wind turbine, Int. J. Renew. Energy Resour. 4 (4) (2014) 935. [44] C. Sim~ ao Ferreira, G. van Kuik, G. van Bussel, F. Scarano, Visualization by PIV of dynamic stall on a vertical axis wind turbine, Exp. Fluid 46 (1) (2009) 97e108. [45] J. Joeng, F. Hussain, On the identification of a vortex, J. Fluid Mech. 285 (1995) 69e94. [46] H. Zhong, P. Du, F. Tang, L. Wang, Lagrangian dynamic large-eddy simulation of wind turbine near wakes combined with an actuator line method, Appl. Energy 144 (2015) 224e233.