Large-Eddy Simulation of Wind Turbine Wake and Aerodynamic Performance with Actuator Line Method

Qian Yaoru，Wang TongguangJiangsu Key Laboratory of Hi-Tech Research for Wind Turbine Design，Nanjing University of Aeronautics＆Astronautics，Nanjing 210016，P.R.China(Received 24 August 2015; revised 16 December 2015; accepted 5 January 2016)



Large-Eddy Simulation of Wind Turbine Wake and Aerodynamic Performance with Actuator Line Method

Qian Yaoru，Wang Tongguang*
Jiangsu Key Laboratory of Hi-Tech Research for Wind Turbine Design，Nanjing University of Aeronautics＆Astronautics，Nanjing 210016，P.R.China
(Received 24 August 2015; revised 16 December 2015; accepted 5 January 2016)

Abstract:A hybrid method is presented to numerically investigate the wind turbine aerodynamic characteristics.The wind turbine blade is replaced by an actuator line model.Turbulence is treated using a dynamic one-equation subgrid-scale model in large eddy simulation.Detailed information on the basic characteristics of the wind turbine wake is obtained and discussed.The rotor aerodynamic performance agrees well with the measurements.The actuator line method large-eddy simulation (ALM-LES) technique demonstrates its high potential in providing accurate load prediction and high resolution of turbulent fluctuations in the wind turbine wakes and the interactions within a feasible cost.

Key words:actuator line; large-eddy simulation (LES) ; wind turbine; wake

CLC number: TK83Document code: AArticle ID: 1005-1120(2016)01-0026-11

* Corresponding author，E-mail address: tgwang@ nuaa.edu.cn.

How to cite this article: Qian Yaoru，Wang Tongguang.Large-eddy simulation of wind turbine wake and aerodynamic perform

ance with actuator line method［J］.Trans.Nanjing Univ.Aero.Astro.，2016，33(1) : 26-36.

http: / /dx.doi.org/10.16356/j.1005-1120.2016.01.026

0　Introduction

With the rapid development in the past two decades，wind energy industry has become one of the important components in renewable energy.A growing number of large scale on-land and offshore wind turbines and wind farms have been built worldwide.The size of the rotor increases to harvest more energy，leading to an increase in blade loading and stress.Due to the high turbulence levels in the various operating conditions of large scale wind turbines，accurate evaluation of the wind turbine aerodynamic performance and its wake has become a critical issue.For wind farms，the interactions from the upstream vortex systems can increase the fatigue and extreme loads on the downstream turbines.The turbulent atmosphere boundary layer (ABL) also significantly influences the prediction of the transient mechanical loading and power production.

A large amount of efforts have been made on blade loading and wake aerodynamics in the past few years，as reviewed by Sanderse et al.［1］.Most of the current evaluations of the wind turbine performance rely on the blade element momentum (BEM) method.Despite of the efficiency of the BEM theory，it lacks the ability to provide detailed information of the velocity field around the rotor.The computational fluid dynamics (CFD) method is a promising approach that can provide detailed flow field information by solving the Navier-Stokes equations on the entire computational domain，but it comes with a price of a large increase in computation cost.Reynolds-averaged Navier-Stokes method (RANS) and large-eddy simulation (LES) are the two main techniques in turbulence modeling.In comparison to the turbulence modeling based on time-averaged RANS approach，LES has attracted more attention recently with its advantages in presenting a time-dependent high quality solution，which could help us get a better understanding of the flow characteristics around the wind turbine.However，LES is currently prohibited in three-dimensional (3D) simula-tion of wind turbine aerodynamic performance and wake due to its very high computational overhead.To overcome this problem，generalized actuator approaches have been proposed for rotor modeling instead of the directly blade-resolved simulation to reduce the computational cost.LES is usually conducted with the actuator line method (ALM) in wind turbine wake simulation.In ALM，each blade is replaced by a line of force that rotates through the flow field.The reduced cost of the ALM makes LES a potential approach for the turbulence modeling in the wind turbine wake simulation.The ALM-LES method has the potential in providing the accurate load prediction and high resolution of turbulent fluctuations in the wind turbine wakes and the interactions within a feasible cost.

The current ALMs are mostly based on the work proposed by Sørensen and Shen［2］with a RANS solver.Further developments followed by Troldborg et al.［3-5］and Watters et al.［6］.Some extended work of the use of the ALM combined with LES to simulate the ABL effect was conducted by Lu and Porté-Agel［7］with turbulence modeled through a Lagrangian scale-dependent subgrid-scale (SGS) model.An ABL-LES solver has been built based on OpenFOAM［8］by the researchers at the National Renewable Energy Laboratory (NREL)，and has demonstrated its potential capabilities to model large wind farms and wake effects［9-10］.Recently，the ALM concept has been extended to actuator surface methods，such as the work of Dobrev et al.［11］，Shen et al.［12］，and Watters et al.［6］.Detailed investigations of the suitable parameters for actuator line and actuator surface methods have been conducted by Jha et al.［13］.Researches on the aerodynamic interactions under turbulent ABL in wind farms were presented by Stovall［14］，Ivanell［15］and Meyers et al.［16］.

The numerical simulations in the present work are conducted using a combination of ALM and LES techniques to explore its capabilities in aerodynamic prediction and establish the main characteristics of the wake development.This ALM-LES hybrid technique is developed here to numerically simulate the aerodynamic characteristics and wake flow fields for the NREL Phase VI rotor and the NH-1500 rotor.The development and characteristic features of the vortex system are captured using LES with a dynamic oneequation subgrid-scale model.All the numerical results of the ALM-LES have been compared against the available experimental data and the BEM method.The measured data of the Phase VI rotor were obtained from NREL Unsteady Aerodynamics Experiment［17-19］.The experimental investigations of the NH-1500 blade were performed by Xiao et al.［20-21］at China Aerodynamics Research and Development Center.The wake characters and patterns at different tip ratios have been presented.All numerical simulations presented in this paper are conducted with a steady and uniform inflow.

1　Numerical Methods

The computations of the flow field around the wind turbine are carried out using a 3 D incompressible transient solver with finite volume method to solve the incompressible Navier-Stokes equations on the non-staggered grid.Both the convective and viscous fluxes in the equations are approximated using a second-order central difference scheme.The solution is advanced in time using a second-order Crank-Nicolson implicit scheme.The pressure correction equation is solved using the pressure-implicit splitting operation (PISO) algorithm.In order to avoid the pressure decoupling，velocity fluxes at the cell faces are constructed using the Rhie-Chow interpolation approach.

1.1Governing equations

The governing equations of the wind turbine flow field are the 3 D incompressible Navier-Stokes equations with turbulence modeled by LES.Based on the theory of Kolmogrorov，the smallest scales of motion are uniform and serve mainly to drain energy from the larger scales during the cascade process.The large scales of mo-tion contain most of the energy and can be strongly affected by the boundary conditions.In LES，the contributions of the large，energy-carrying scales to momentum and energy transfer are computed directly and effects of the small scales are modeled.The variable u of the turbine flow field is decomposed into a filtered，large-scale resolved partand a sub-filtered unresolved part u' by applying a filtering operation.

where G is the filter function and Ω denotes the entire computation domain.The top-hat filter function is used.

where H (x) represents the Heaviside function.And Δ in the filter function represents the filter width，which is defined as Δ= (Δ1Δ2Δ3)1 /3with Δibeing the filter width in the i-th coordinate direction.

The governing 3D spatial filtered incompressible Navier-Stokes equations in the differential form are written as

where fεdenotes the forces on the wind turbine rotor calculated by the actuator line model and ν the kinematic viscosity coefficient.The third term on the right-hand side of the momentum Eq.(4b)，as τ=－，is the so-called sub-grid scale stress that represents the interactions between the unresolved scales of eddies.It is of great importance to set a suitable SGS stress model in LES.

1.2Sub-grid scale modeling

Turbulence is a dominator in wind turbine aerodynamics that affects both the blade performance and wake behavior.A widely used approach to compute the SGS stress is to use an eddy-viscosity model based on the hypothesis that the non-uniform component of the SGS stress tensor is locally aligned with the resolved non-uniform part of the rate of strain tensor.To capture the transient characters of the turbine flow field，the SGS stress is modeled by a dynamic one-equation model based on the eddy-viscosity concept.The one-equation SGS model uses a modeled balance equation for the SGS turbulent energy to simulate the behavior of the sub-grid scale kinetic energy KSGS.This localized dynamic KSGSequation model is proposed by Menon et al.［22］.

The transport equation for the SGS turbulent kinetic energy is written as

The sub-grid scale kinetic energy KSGSand the eddyviscosity νSGSare defined in Eqs.(6).The SGS stress tensor is modeled as

where CKis a constant with the value of 0.07 here.

The first term on the right-hand side of the transport equation ΠkSGSis the SGS dissipation of the resolved kinetic energy，andis the strain rate tensor on the grid level.ΠkSGSandare defined as

The dynamic coefficient C and the dissipation coefficient C*are computed as below.

where Lijare the dynamic Leonard stresses defined as Lij=－The explicit filtering at the test level is executed by integrating over the test cell assuming linear variation of the variables.

As been studied in the prior work［23-25］，the one-equation SGS model has the advantage of providing a more accurate time scale through the independent definition of the velocity scale when comparing with the common Smagorinsky model，especially in modeling the transitional flows orflows with large scale unsteadiness.

1.3Actuator line method

As mentioned above，the presence of the wind turbine rotor is replaced by body forces based on the actuator line model developed by Sørensen and Shen［2］.The ALM module is combined with an unsteady incompressible solver.The time-varying turbine loads are modeled by a suitable distribution of body forces along the blade determined under sectional inflow conditions and an interpolation of force coefficients from the lookup table of airfoil properties.In the present work，the body forces are distributed radially along the actuator lines.The local velocity of the blade element point is obtained using linear interpolation through the neighbor cells.Corrections for the tip loss are also applied.

The concept of ALM is shown in Fig.1 together with a cross-sectional airfoil element with the blade forces on.w denotes the angular velocity.Vxand Vθrepresent the local axial and tangential velocities at the actuator blade element point，respectively，which are interpolated from the velocities at the neighbor mesh cell center.The local relative velocity of the sectional airfoil element is determined as

The angle between Vreland the rotor plane is computed as

The local angle of attack is given by α=φ－γ，where γ denotes the local pitch angle.A tip loss correction by Glauert［26］is used.

Fig.1　Sketch of ALM for wind turbine and the forces on cross-sectional airfoil

The tip loss function is computed as

where nbis the number of the turbine blades and R the blade tip radius.For all blade elements on the actuator line，once the angle of attack and the relative velocity are determined，the lift and drag force per spanwise length can be approximated as

Tabulated lift and drag coefficients of the typical airfoils of the wind turbines are obtained from the NREL wind tunnel experiment［17-19］and the in-house NH1500 blade design report.To smooth the distribution of the rotating forces and avoid the numerical oscillations，Gaussian projection is employed to transform the line force into a volumetric force field described as

where d is the distance between the cell-centered grid point and the point at the actuator line，and ε the Gaussian width parameter to adjust the concentration of the load on the blade.The Gaussian width ε being twice the radial grid spacing Δr was suggested［4-5］to maintain the numerical stability and a good prediction of the rotor power.Another detailed research on the influence of different ε has been carried out by Jha et al.［13］.Various choices of ε in relationship with other characteristic length variables together with the performances in load prediction have been studied.In this work，Troldborg' s suggestion of ε/Δr = 2 isadopted.

2　Results and Discussion

The validation analysis of the ALM-LES solver against the wind tunnel experimental data are presented in this section.Comparison of the spanwise distributions of the normal and tangential force coefficients together with the thrust and power of the PhaseⅥrotor are provided.In the case of the NH-1500 rotor，the power coefficient Cpis compared with experimental data，while the normal and tangential force coefficients are compared with the BEM results.

The computation domains of PhaseⅥand NH-1500 rotor are both Cartesian with the size Lx×Ly× Lzbeing 8R×8R×20R，where R is the rotor radius.The inner refined region extends from－2R to 10R in the x direction.The size of the computational mesh is around 2.0×107.The computational domains are sketched in Fig.2.

Fig.2　Sketches of the computational domain and coordinates

Table 1　Parameters of phaseⅥand NH-1500 rotors

The study subject，NH-1500 turbine，is a 1.5 MW turbine developed at Nanjing University of Aeronautics and Astronautics.The rotor consists of a three-bladed rotor with the radius of 40.5 m and the rated wind speed of 10.4 m/s.The inputs of the geometry and the solver parameters of PhaseⅥand NH-1500 rotor are listed in Table 1.

2.1Aerodynamic loads

2.1.1NREL PhaseⅥrotor

The operating condition for the simulations of the NREL PhaseⅥrotor is shown in Table 2.

Table 2　Operating condition parameters for PhaseⅥrotor

At the wind speed of 7 m /s，the flow is mainly attached along the blade surface.The spanwise variations of the normal force coefficient Cnand tangential force coefficient Ctfor the NREL Phase Ⅵrotor are shown in Figs.3－4.It can be observed that both Cnand Ctfrom the simulation agree well with the experimental data except a deviation in Ctat the very inboard spanwise station.This difference may be caused by the over-predicted local angle of attack or the unsuitable 2D sectional airfoil coefficients，which leads to the overprediction of force coefficients.The over-prediction also results in a little higher torque and power as shown in Table 3.

Fig.3　Spanwise distribution of Cnfor NREL PhaseⅥrotor (V∞=7　m/s)

Fig.4　Spanwise distribution of Ctfor NREL PhaseⅥrotor (V∞=7　m/s)

Table 3　Comparison of the thrust and power for PhaseⅥ

2.1.2NH-1500 rotor

The LES-ALM solver is also applied to the NH-1500 rotor under various operating conditions as described in Table 4.The NH-1500 rotor has a pre-cone angle of 3°and a shaft tilt-up angle of 5°.

The experimental investigation on the aerodynamic performance of a 1 /16 scaled NH-1500 model was conducted in the CARDC 12 m×16 m low speed wind tunnel by Xiao et al.［20］.In Fig.5，a comparison between the computed power coefficient and the measured data is presented.It can be observed from Fig.4 that the overall power coefficient agrees well with the wind tunnel experiment.

Table 4　Simulation conditions of NH-1500 rotor

Fig.5　Comparison of the power coefficient Cpwith experiment data

Unfortunately，the experiment by Xiao et al.［20］provided neither aerodynamic force distributions along the blade nor wake data.The numerical results are compared with BEM results for now.The spanwise distribution of the normal force coefficient Cnand tangential force coefficient Ctfor the NH-1500 rotor at different wind speeds are shown in Figs.6—11.Generally speaking，the results obtained from the ALM-LES method are in good agreement with the ones from the BEM theory except that both Cnand Ctfrom ALM-LES is lower than the BEM results towards in the root and tip area.It can be reasoned that there are no interactions between the blade elements in BEM calculation.Nevertheless，the validations of blade loading with the wind tunnel experimental data need to be further investigated.

Fig.6　Spanwise distribution of Cnfor NH-1500　rotor (V∞=4　m/s)

Fig.7　Spanwise distribution of Ctfor NH-1500　rotor (V∞=4　m/s)

Fig.8　Spanwise distribution of Cnfor NH-1500　rotor (V∞=8　m/s)

Fig.9　Spanwise distribution of Ctfor NH-1500　rotor (V∞=8　m/s)

Fig.10　Spanwise distribution of Cnfor NH-1500　rotor(V∞=10　m/s)

Fig.11　Spanwise distribution of Ctfor NH-1500　rotor (V∞=10　m/s)

In order to verify whether the numerical computations can be classified as LES，the spectral characteristics of the wake are presented.Due to the homogeneous and universal character of the sub-grid scales，the energy spectrum should be proportional to k－5 /3，where k denotes the wavenumber.Figs.12(a)，(b) show the spectra obtained in the downstream wake of PhaseⅥand NH-1500.The ideal－5 /3 slope is also shown in Fig.12.It is clearly shown that the drop in the high frequency range fits the－5 /3 law which confirms that the numerical simulations could be considered a well-resolved LES.

Fig.12　Spectra of the PhaseⅥand NH-1500　wakes

2.2Wake

In contrast to RANS，LES can provide more accurate，high quality description of the wake characteristics including the main features of the wake vortex systems，the wake expansion and the wake velocity distribution.In this section，the wake development and some basic characteristics of both the PhaseⅥand NH-1500 rotors are presented.

2.2.1NREL PhaseⅥrotor

The development of the vortex system trailed from the PhaseⅥrotor in both near and far wakes at λ= 5.42，which is the smallest tip ratio among all the simulations，has been presented in Fig.13.It can be noticed that the wake downstream mainly consists of the vortices trailed from the tip and root of the rotor.It is observed from the vorticity contours presented in Fig.14 that the vortices from the tip preserve a stable circulation to about three radii downstream，and then turn into a continuous vortical sheet.The vortices break down at about 10 rotor radii downstream.The root vorticity is much less significant than the tip one and fades faster in the near wake.

2.2.2NH-1500 rotor

To better understand the vortex structures in NH-1500 rotor wake，the iso-surfaces of the second invar-iant of the velocity gradient Q in the instantaneous velocity field at different tip speed ratios are shown in Figs.15(a)—(c).Figs.16(a)—(c) show the vorticity magnitude contours of the NH-1500 wake development at different λ in the horizontal plane through the turbine center axis.

Fig.13　Vortex structure of the PhaseⅥrotor wake with iso-surface Q =0.005

Fig.14　Visualization of the vorticity contour in the Phase Ⅵrotor wake

The wake pattern of each operating condition is distinguished from one another.It can be observed that the higher the tip speed ratio is，the shorter distance the stable tip-vortex pattern preserves.At the highest tip speed ratio，λ= 10.77，the vortices shed from the rotor tip does not form into a stable circulation.The instability of the tip vortices occurs instantly downstream and the wake breaks into small-scale turbulence very fast.While at λ= 7.49，the lowest tip speed ratio，the tip vortices preserve in a smooth and stable circulation state，and then smear into a continuous vortical sheets before the wake becomes fully turbulent at around 6.5 radii downstream.For the rotor operating at λ= 9.36，the interaction occurs sooner than the wake at λ= 7.49.Before the break-down，the vortical tube appears in a halfstable state with fluctuations.

Fig.15　Vortex structure of the NH-1500　rotor wake with velocity contours and iso-surface of Q = 0.005　at different tip speed ratios

Fig.16　Vorticity contours for the NH-1500　wake development at different tip speed ratios

From the comparisons above，it can be observed that the wake can turn into turbulence even when there is no turbulence in the inflow.The wake expansion becomes more significant with the tip speed ratio increase.Besides the tip and root vortices，the simulations have captured the vortices shed from the inner sections of the rotor actuator lines，demonstrating that the ALM-LES method can simulate complicated and turbulent flow field around wind turbine.The inner section vortices become more obvious when the loadsget heavier，and thus lead to more severe interactions.It should be noticed that the jagged profile in the far wake，after 15 radii downstream，is caused by the coarse mesh in this region.

The graphs in Fig.17 show the development of the x-velocity distributions in the NH-1500 wake at various wake positions in the vertical plan through the turbine axis at different tip speed ratios.

At λ=10.77，where the rotor is the most heavily loaded in the simulations，the wake expansion effect is more apparent，which is in consistent with the discussion above.And due to the instant breakdown，the decrease in the axial velocity is so sharp that the x-velocity distribution fades into a Gaussianlike shape around two radii downstream.The more turbulent the wake gets，the faster it turns into the Gaussian-like shape distribution.

All the simulations neglect the presence of both the nacelle and the tower，which will be carried out in future work.

3　Conclusions

A numerical approach based on the LES technique in conjunction with ALM has been applied to simulate the flows around wind turbine rotors.The aerodynamic loads on the rotors have been predicted and compared with the measured data in order to validate the hybrid LES-ALM method.A detailed investigation of the wake development has been made.The spectral characteristics in the wake indicate that the small scales are well captured and the LES can be considered well resolved.

The LES-ALM method shows its ability to provide a considerably accurate prediction of the aerodynamic loads on wind turbine.The normal and tangential coefficients agree well with either the measured data or the BEM results.The wind turbine wake can be transited turbulent even without inflow turbulence.The heavier load the rotor gets，the sooner the instability occurs.The wake expansion becomes more significant with the tip speed ratio increase.tional Natural Science Foundation of China (No.11172135)，and the EU Seventh Framework Program (No.FP7-PEOPLE-2010-IRSES-269202).

Fig.17　Distributions of the mean x velocity in the wake of NH-1500　rotor at different tip speed ratios

References:

［1］SANDERSE B，VAN DER PIJL S P，KOREN B.Review of computational fluid dynamics for wind turbine wake aerodynamics［J］.Wind Energy，2011.

［2］SØRENSEN J N，SHEN W Z.Numerical modeling of wind turbine wakes［J］.ASME Journal of Fluids Engineering，2002，124: 393-399.

［3］TROLDBORG N，SØRENSEN J N，MIKKELSEN R.Actuator line simulation of wake of wind turbine operating in turbulent inflow［J］.The Science of Making Torque from Wind Journal of Physics: Conference Series 75，2007: 012063.

［4］TROLDBORG N.Actuator line modeling of wind turbine wakes［D］.Denmark: Technical University of Denmark，2008.

［5］TROLDBORG N，SØRENSEN J，MIKKELSEN R.Numerical simulations of wake characteristics of a wind turbine in uniform flow［J］.Wind Energy，2010，13: 86-99.

［6］WATTERS S C，MASSON C.Modelling of lifting-device aerodynamics using the actuator surface concept ［J］.Int J Numer Meth Fl，2010，62 (11) : 1264-1298.

［7］LU H，PORTÉ-AGEL F.Large-eddy simulation of a very large wind farm in a stable atmospheric boundary layer［J］.Physics of Fluids，2011，23: 065101.

［8］OpenFOAM Ver.2.1.x.ESI Group-OpenCFD［EB/ OL］.［2013-12-13］.http:∥openfoam.com.

［9］CHURCHFIELD M J，MORIARTY P J，VIJAYAKUMAR G，et al.Wind energy-related atmospheric boundary-layer large eddy simulation using OpenFOAM: NREL /CP-500-48905［R］.［S.l.］: NREL，2010.

［10］CHURCHFIELD M J，LEE S，MORIARTY P J，et al.A large-eddy simulation of wind-plant aerodynamics: AIAA 2012-0537［R］.Nashville，TN: AIAA，2012.

［11］DOBREV I，MASSOUH F，RAPIN M.Actuator surface hybrid model［C］∥Journal of Physics Conference Series，the Science of Making Torque from Wind.Lyngby Denmark: Technical University of Denmark，2007.

［12］SHEN W Z，ZHANG J H，SORENSEN J N.The actuator surface model: A new Navier-Stokes based model for rotor computations［J］.Sol Energy Eng Trans ASME，2009，131(1).DOI: 10.1115/1.3027502.

［13］JHA P K，CHURCHFIELD M J，MORIARTY P J.Accuracy of state-of-the-art actuator-line modeling for wind turbine wakes: AIAA 2013-0608［R］.［S.l.］: AIAA，2013.

［14］STOVALL T D，PAWLAS G，MORIARTY P J.Wind farm wake simualtions in OpenFOAM: AIAA-2010-0825 ［R］.［S.l.］: AIAA，2010.

［15］IVANELL S，MIKKELSEN R，SØRENSEN J，et al.ACD modelling of wake interaction in the horns REY wind farm［C］∥Extended Abstracts for Euromech Colloquim 508 on Wind Turbine Wakes，European Mechanics Society.Madrid，Spain:［s.n.］，2009.

［16］MEYERS J，MENEVEAU C.Large eddy simulations of large wind-turbine arrays in the atmospheric boundary layer: AIAA-2010-0827［R］.［S.l.］: AIAA，2010.

［17］HAND M M，SIMMS D A，FINGERSH L J，et al.Unsteady aerodynamics experiment phaseⅥ: Wind tunnel test configurations and available data campaigns［M］.［S.l.］: National Renewable Energy Laboratory (NREL)，2001.

［18］DUQUE E P N，BURKLUND M D，JOHNSON W.Navier-Stokes and comprehensive analysis performance predictions of the NREL phaseⅥexperiment［J］.Journal of Solar Energy Engineering，2003，125(4) : 457-467.

［19］GONZALEZ A，MUNDUATE X.Three dimensional and rotational aerodynamics on NREL phaseⅥwind turbine blade［J］.Journal of Solar Energy Engineering，2008，130(3) : 1-7.

［20］XIAO J P，CHEN L，XU B F，et al.Investigations on aerodynamic performance of a 1.5 MW wind turbine ［J］.Acta Aerodynamic Sinica，201l，29(4) : 131-136.

［21］XIAO J P，WU J，CHEN L，et al.PIV measurements of tip vortex wake structure of a wind turbine［J］.Applied Mathematics and Mechanics，2011，32(6) : 683-692.

［22］MENON S，KIM W W.High Reynolds number flow simulations using the localized dynamic subgrid-scale model: AIAA 96-0425［R］.［S.l.］: AIAA，1996.

［23］GERMANO M，PIOMELI U，MOIN P，et al.A dynamic subgrid-scale eddy viscosity model［J］.Physics of Fluids，1991，3(7) : 1760-1765.

［24］LILLY D.A proposed modification of the Germano subgrid-scale closure mthod［J］.Physics of Fluids，1992，4(3) : 633-635.

［25］DEARDORFF J W.A numerical study of three-dimensional turbulent channel flow at large Reynolds number ［J］.Fluid Mechanics，1970，41(2) : 453-480.

［26］GLAUERT H.Airplane propellers，In aerodynamic theory［M］.Durand WF (Ed.).Dover: New York，1963: 169-360.

Ms.Qian Yaoru is a Ph.D.student of Nanjing University of Aeronautics and Astronautics (NUAA).She received her B.S.degree in Aircraft Design and Engineering from NUAA in June 2011.Her recent researches focus on the numerical investigation of wind turbine aerodynamics and wake analysis.

Dr.Wang Tongguang is the Director of Jiangsu Key Laboratory for Wind Turbine Design at NUAA.He was designated as the Chief Scientist of the " 973 Program" by the Ministry of Science and Technology，China in 2007，and successfully completed the project "Fundamental Study of Large Scale Wind Turbine Aerodynamics".As the Chief Scientist again，he is now in charge of another project of the "973 Program"—"Key Mechanical Issues and Design of Large Scale Wind Turbines".Prof.Wang graduated from NUAA with B.S.degree and M.S.degree in 1983 and 1988，respectively，and he received his Ph.D.degree from the University of Glasgow in 1999 and has been a postdoctoral fellow at the University of Glasgow from October 1999 to August 2001.During this period，he was invited to join in the Wind Turbine Unsteady Aerodynamics Experiment -Blind Comparison，organized by the National Renewable Energy Laboratory (NREL) of USA.His paper "An examination of key aerodynamic modeling issues raised by the NREL Blind Comparison" was awarded The Best Paper Prize by ASME/AIAA in 2002.

(Executive Editor: Zhang Tong)