A new class of actuator surface models for wind turbines

4.4 MB·

by  Xiaolei Yang

WINDENERGY
WindEnerg.2017;00:1–17
DOI:10.1002/we
RESEARCHARTICLE
A new class of actuator surface models for wind turbines
XiaoleiYang1,FotisSotiropoulos1
7 1DepartmentofCivilEngineering,CollegeofEngineeringandAppliedSciences,StonyBrookUniversity,StonyBrook,NY11794,
1
USA
0
2

n
a ABSTRACT
J

2 Actuatorlinemodel hasbeenwidelyemployed inwindturbinesimulations.However, thestandardactuator linemodel
1 doesnotincludeamodelfortheturbinenacellewhichcansignificantlyimpactturbinewakecharacteristicsasshownin
theliterature(e.g.Kang, YangandSotiropoulos,Journal ofFluidMechanics 744(2014): 376-403; Violaetal.,Journal
]
ofFluidMechanics750(2014):R1;Fotietal.,PhysicalReviewFluids1(2016),044407). Anotherdisadvantageofthe
n
y standardactuatorlinemodel isthatmoregeometrical featuresofturbinebladescannot beresolvedonafinermesh.To
d alleviatethesedisadvantagesofthestandardmodel,wedevelopanewclassofactuatorsurfacemodelsforturbineblades
- andnacelletotakeintoaccount moregeometrical detailsof turbinebladesand includetheeffect of turbinenacelle. In
u
theactuator surfacemodel forblade, theaerodynamic forcescalculated usingthebladeelement methodaredistributed
l
f from the surface formed by the foil chords at different radial locations. In the actuator surface model for nacelle, the
.
s forcesaredistributedfromtheactualnacellesurfacewiththenormalforcecomponent computedinthesamewayasin
c thedirectforcingimmersedboundary methodandthetangentialforcecomponent computedusingafrictioncoefficient
i
s andareferencevelocityoftheincomingflow.Theactuatorsurfacemodelfornacelleisevaluatedbysimulatingtheflow
y over periodically placed nacelles. Boththe actuator surface simulation and thewall-resolved large-eddy simulation are
h carriedout.Thecomparisonshowsthattheactuatorsurfacemodelisabletogiveacceptableresultsespeciallyatfarwake
p
locationsonaverycoarsemesh.Thecapabilityoftheactuatorsurfacemodelinpredictingturbinewakesisassessedby
[
simulatingtheflowovertheMEXICO(ModelexperimentsinControlledConditions)turbineandthehydrokineticturbine

1 ofKang,YangandSotiropoulos(JournalofFluidMechanics744(2014):376-403).Comparisonsofthecomputedresults
v withmeasurementsshowthattheproposedactuatorsurfacemodelisabletopredictthetipvortices,turbulencestatistics
8 andmeanderingofturbinewakewithgoodaccuracy.Copyright(cid:13)c 2017JohnWiley&Sons,Ltd.
0
1 KEYWORDS
2 Turbineparameterization;Nacellemodel;Wakemeandering
0
Correspondence
.
2 FotisSotiropoulos,DepartmentofCivilEngineering,CollegeofEngineeringandAppliedSciences,StonyBrookUniversity,Stony
0 Brook,NY11794,USA
7 Emailaddress:[email protected].
1
: Received...
v
i
X
r
a
1. INTRODUCTION
Usingthecurrentstate-of-the-artsupercomputers, itisstillimpossibletosimulatewindfarmsresolvingalltherelevant
lengthscalesfromthethicknessoftheboundarylayeroverawindturbineblade(≈0.016mestimatedusingtheequation
δ/c=0.37(U c/ν)−0.2[1]withν =1.5×10−5m2/s,andairfoilchordc=1mandU =100m/satthebladetip
∞ ∞
foraMWscalewindturbinewithrotordiameterapproximately100m,tipspeedratio10andincomingwindspeed10
m/s)tothethicknessofanatmosphericboundarylayer(≈1000m[2]).Differentrotormodelsincludingtheactuatordisk
model[3,4,5],theactuatorlinemodel[6,7]andtheactuatorsurfacemodel[8]havebeendevelopedintheliteratureto
parametrizetheinteractionbetweentheturbinebladesandtheincomingflowwithoutdirectlyresolvingtheboundarylayer
flowoverturbineblades.
Intheactuatordiskmodel[9]thewindturbinerotorismodelledbyacircularpermeabledisk.Theimpactsofincoming
windonwindturbinesaretakenintoaccountbyuniformlydistributedforcesonthedisk.Intheconventionalactuatordisk
Copyright(cid:13)c 2017JohnWiley&Sons,Ltd. 1
Preparedusingweauth.cls[Version:2010/06/17v1.00]
Anewclassofactuatorsurfacemodelsforwindturbines XiaoleiYangandFotisSotiropoulos
model,theforcesareonlyappliedalongtheturbineaxialdirection,whichcanbecalculatedeitherfromone-dimensional
momentumtheorywithaprescribedinductionfactor[3]orfromtheblade-elementmomentumtheory[10].Intheworkby
WuandPorte´-Agel[11],therotationwastakenintoaccountintheactuatordiskmodel.Itwasfoundthatthepredictions
fromtheactuatordiskmodelwithrotationshowbetteragreementwiththewindtunnelmeasurementsreportedin[12].
Theactuatorlinemodelprovidesamoreaccurateapproachtoincorporatetherotationaleffect,finitebladenumbereffect
andtheeffectofnon-uniformforcedistributionintheazimuthaldirection.Inthismodel,theturbinebladeisrepresented
byarotatinglinewithdistributedforces,whicharecalculatedfromabladeelementapproach[13]combinedwithtabulated
dragandliftcoefficientsofthefoils.Theemployeddragandliftcoefficientsareusuallyfromtwo-dimensionalnumerical
simulationsorexperimentswithoutconsideringtherotationaleffects.Inordertotakeintoaccountthethree-dimensional
rotational effects, corrections on the two-dimensional force coefficients, such as those proposed in [14], are usually
employed. As discussed by Shen et al. [15], the lift needs to approach zero as the blade tip is approached. However,
theliftonthetippredictedbythebladeelementapproachisusuallynotzeroespeciallyforapitchedbladeorachambered
foil.Atip-losscorrectioninadditiontothecorrectionforthethree-dimensionaleffectisneededinordertohaveacorrect
liftdistributionintheregionnearthebladetipassuggested byShenetal.[16].Therearetwomajorlimitationsof the
standardactuatorlinemodel:(i)thelackofaneffectivenacellemodel,and(ii)thefactthatafinermeshcannotresolvemore
geometricalfeaturesoftheturbineblade.Nacelleinducedcoherentstructureswereshowntohaveasignificantimpacton
turbinewakecharacteristics,suchasvelocitydeficitinthenearwake[17,7,18]andwakemeanderingatfarwakelocations
formodelwindturbines[19,20,21]andhydrokineticturbines[17].Bothgeometry-resolving,wall-modelledlarge-eddy
simulations and actuator type large-eddy simulations were carried out in [17]. It is shown that the actuator line model
withoutanacellemodelcannotaccuratelycapturethewakemeanderingofthehydrokineticturbine,andunderpredictsthe
turbulenceintensityatfar wake locations. Thiswasattributedtothefact that thepresence of thenacelleinduces spiral
vortexbreakdown of thehubvortexcausing alarge-scaleenergeticcoherent vortexthatprecessesand expands radially
outwardtointercept andenergizethetipshearlayeratthreetofourrotordiametersdownwindoftheturbine[17].The
nacellecanberepresentedusingapermeablediskwithaspecifieddragcoefficient.However,Yangetal.[7]showedthat
such a simple nacelle model is not able to take into account the nacelle effect accurately. The computed results from
Tossasetal.[18]indicatedthatwiththeirnacelleandtowermodelthevelocitydeficitinthenearwakeispredictedmore
accurately.Theimprovementonthepredictionoftheturbulenceintensity,ontheotherhand,isverylimited.Theresultsof
Yangetal.[7]alsoshowedthatgrid-independentresultscannotbeobtainedwhenthemeshisrefinedintheactuatorline
modelsimulations.
Inorder to alleviate the above mentioned limitations of the actuator linemodel, we develop a new classof actuator
surface models for turbine blade and nacelle for taking into account the nacelle geometry and better incorporating the
geometricaleffectoftheturbineblade.Intheactuatorsurfacemodelforblade,theforcesarestillcomputedusingtheblade
elementmethodbutaredistributedonthesurfaceformedbythefoilchordsatdifferentradiallocations.Intheactuator
surfacemodelforthenacelle,thenormalforcesarecomputedbysatisfyingthenon-penetrationboundaryconditionsasin
thedirectforcingimmersedboundarymethod[22,23,24].Thetangentialforces,ontheotherhand,arecalculatedusing
africtioncoefficientandareferencevelocityoftheincomingflow.Moredetailsontheproposedactuatorsurfacemodels
arepresentedinSec.2followedbyabriefdescriptionofthegoverningequationsandflowsolverinSec.3.Somecases
forvalidatingtheproposedactuatorsurfacemodelarethenreportedinSec.4.Finally,conclusionsaredrawninSec.5.
2. ACTUATORSURFACEMODELS
Inthissection,wepresenttheproposedactuatorsurfacemodelsforthebladeandnacelle.Inbothactuatorsurfacemodels,
wehavetwosetsofindependentmeshes,i.e.,thebackgroundCartesiangridfortheflowwithitscoordinatedenotedby
x(x,y,z orx ,x ,x ),andtheLagrangiangridfollowingtheactuatorsurfaceswithitscoordinatedenotedbyX(X,
1 2 3
Y,Z orX ,X ,X ).Forbothmodels,thesmootheddiscretedeltafunctionin[24]isemployedasthekernelfunction
1 2 3
fortransferringinformationbetweenthetwomeshes.Themajordifferencebetweentheactuatorsurfacemodelforblade
andtheactuatorsurfacemodelfornacelleisthewayhowtheforcesontheactuatorsurfacesarecalculated.InSec.2.1we
presenttheforcecalculationanddistributionmethodsoftheactuatorsurfacemodelforblade.AndinSec.2.2weintroduce
theforcescalculationmethodfortheactuatorsurfacemodelfornacelle.
2.1. Actuatorsurfacemodelforblade
Intheactuatorsurfacemodelforblade,thebladegeometryisrepresentedbyasurfaceformedbythechordlinesatdifferent
radiallocationsofabladeasshowninFig.1(a).Theforcesarecalculatedinthesamewayasintheactuatorlinemodel.
Thelift(L)anddrag(D)ateachradiallocationarecalculatedby
1
L= ρC c|V |2e (1)
2 L rel L
2 WindEnerg.2017;00:1–17(cid:13)c 2017JohnWiley&Sons,Ltd.
DOI:10.1002/we
Preparedusingweauth.cls
XiaoleiYangandFotisSotiropoulos Anewclassofactuatorsurfacemodelsforwindturbines
Figure1.Aschematicoftheactuatorsurfacemodelforblade.Theliftanddragforcescalculatedusingthebladeelementmethodare
distributedovertheactuatorsurfaceformedbychordlinesofablade.
and
1
D= ρC c|V |2e , (2)
2 D rel D
whereC andC aretheliftanddragcoefficients,cisthechordlength,V istherelativeincomingvelocity,ande
L D rel L
ande aretheunitvectorsinthedirectionsofliftanddrag,respectively.TheC andC arefunctionsoftheReynolds
D L D
numberandtheangleofattack.TheangleofattackαshowninFig.1(b)iscomputedateachradiallocationby
α=φ−γ, (3)
whereφ=−tan−1(u /(u −Ωr)),andγistheangleincludingthebladetwistandthebladepitch,thelatterofwhich
x θ
isspecifiedtoagivenvalueinsteadoffromaturbinecontrolalgorithm.TherelativeincomingvelocityV ateachradial
rel
locationiscomputedby
V =u e +(u −Ωr)e , (4)
rel x x θ θ
whereΩistherotationalspeedoftherotor,e ande aretheunitvectorsintheaxialandazimuthaldirections,respectively.
x θ
The u and u are the axial and azimuthal components of the flow velocity averaged over the chord for each radial
x θ
locations,whicharecomputedusing
1
u = u(X)·e ds (5)
x c x
Zc
and
1
u = u(X)·e ds, (6)
θ c θ
Zc
where X denote the coordinates of the grid points on the actuator surfaces. Generally the grid points on the actuator
surfacesdonot coincidewithany background nodes. Weemploy asmoothed discretedeltafunction(i.e.thesmoothed
four-point cosine function) proposed by Yang et al. [24] to interpolate u(X) from the values on the background grid
nodesasfollows:
u(X)= u(x)δ (x−X)V(x), (7)
h
xX∈gx
wherexarethecoordinatesofthebackgroundgridnodes,gx isthesetofthebackgroundgridcells,V =hxhyhz (hx,
h and h are the grid spacings in the x, y and z directions, respectively) is the volume of the background grid cell,
y z
δ (x−X)= 1φ x−X φ y−Y φ z−Z is the discrete delta function, and φ is the smoothed four-point cosine
h V hx hy hz
function[24],which(cid:16)isexp(cid:17)ress(cid:16)edas (cid:17) (cid:16) (cid:17)
1 + sin(π(2|r|+1)/4) − sin(π(2|r|−1)/4), |r|≤1.5,
4 2π 2π
φ(r)= 5 − |r| − sin(π(2|r|−1)/4), 1.5≤|r|≤2.5, (8)
 8 4 2π
 0, 2.5≤|r|,

WindEnerg.2017;00:1–17(cid:13)c 2017JohnWiley&Sons,Ltd. 3
DOI:10.1002/we
Preparedusingweauth.cls
Anewclassofactuatorsurfacemodelsforwindturbines XiaoleiYangandFotisSotiropoulos
inwhichr =(x −X )/h (i=1,2,3).
i i i i
The blade rotation causes the stall delay phenomenon at the inboard sections of the blade, which increases the lift
coefficients and decreases the drag coefficients as compared with the corresponding two-dimensional airfoil data. To
accountforsuchthree-dimensionalrotationaleffect,thestalldelaymodeldevelopedbyDuandSelig[25]isemployedto
correcttheliftanddragcoefficientsfromtwo-dimensionalexperimentsornumericalsimulations.InDuandSelig’smodel,
thecorrectedliftanddragcoefficients(C andC )arecalculatedasfollows:
L,3D D,3D
C =C +f (C −C ), (9)
L,3D L,2D L L,p L,2D
and
C =C −f (C −C ), (10)
D,3D D,2D D D,2D D,0
whereC =2π(α−α ),C thetwo-dimensionaldragcoefficientatzeroangleofattack,andthecorrectionfunctions
L,p 0 D,0
f andf aredeterminedby
L D
dR
1 1.6(c/r)a−(c/r)Λ r
f = −1 , (11)
L 2π 0.1267b+(c/r)ΛdRr !
and
d R
1 1.6(c/r)a−(c/r)2Λ r
f = −1 , (12)
D 2π 0.1267b+(c/r)2dΛRr !
respectively,whereΛ=ΩR/ V2+(ΩR)2,Ristherotorradius,a,banddaretheempiricalcorrectionfactors.Inthis
w
work,a,banddareequalto1asinDuandSelig’spaper[25].
p
Non-zeroforcecanexistatthebladetipwhenthepitchangleisnonzeroorachamberedfoilisused[15].Thisisin
contradictionwiththephysicalunderstandingthattheforceshouldtendtozeroatthetipduetopressureequalizationas
discussed byShenet al.[15].Tocorrect thisnon-physical forcebehaviour, thetip-losscorrection proposed by Shenet
al.[16,15]isappliedtothedragandliftcoefficientscomputedfromEqs.(9)and(10).Withthetiplosscorrection,the
employedC andC arecalculatedas
D L
C =F C , (13)
L 1 L,3D
and
C =F C , (14)
D 1 D,3D
where
2 B(R−r)
F = cos−1 exp −g , (15)
1 π 2rsinφ
(cid:18) (cid:18) (cid:19)(cid:19)
inwhichBisthenumberofblades,giscomputedby
g=exp(−c (BΩR/U −c ))+c , (16)
1 ∞ 2 3
wherec ,c andc arethecorrectioncoefficients,whichareequalto0.125,21and0.1asin[16],respectively.
1 2 3
Aftercalculatingthelift(L)anddrag(D),theforceperunitareaontheactuatorsurfaceateachradiallocationisthen
calculatedas
f(X)=(L+D)/c. (17)
Note that the above expression essentially means that the lift and drag on the blade are uniformly distributed in the
chordwisedirection.Tocalculatetheforcesonthebackgroundmeshfortheflowfield,thecomputedforcesontheactuator
surfacesarethendistributedtothebackgroundgridnodesasfollows:
f(x)=− f(X)δ (x−X)A(X), (18)
h
XX∈gX
wheregXisthesetoftheactuatorsurfacegridcellsandAistheareaoftheactuatorsurfacegridcell.Thesamediscrete
delta function as in Eq. (7) is employed. It is noted the negative sign is because that f(x) represents the forces of the
actuatorsurfacesactingontheflow,whilef(X)denotestheforcesoftheflowactingontheactuatorsurfaces.
4 WindEnerg.2017;00:1–17(cid:13)c 2017JohnWiley&Sons,Ltd.
DOI:10.1002/we
Preparedusingweauth.cls
XiaoleiYangandFotisSotiropoulos Anewclassofactuatorsurfacemodelsforwindturbines
Figure2.Aschematicoftheactuatorsurfacemodelfornacelle.Asshownin(a),thenormalforcefn attheLagrangianpointA
iscomputedbyEq.(19)withtheestimatedvelocityu˜(X)atpointAinterpolatedfromthesurroundingCartesiangridnodesusing
Eq.(7).Themagnitudeofthetangentialforcefτ iscomputedbyEq.(21)usingtheincomingstreamwisevelocityU andafriction
coefficientgivenbyEq.(22).ThedirectionofthetangentialforceisassumedtobethesameasthetangentialvelocityatpointB,
whichisinterpolatedfromthesurroundingCartesiangridnodesusingEq.(7).Thecomputedforcesaredistributedtothesurrounding
Cartesiangridnodes(thegrayareashownin(a))usingEq.(18).
2.2. Actuatorsurfacemodelfornacelle
Resolvingtheboundarylayeroveranacelleofanutilityscalewindturbine,whichbecomesalmostimpossibleforwind
farmscalesimulations,requiresamuchfinermeshascomparedwiththemeshforresolvingthethicknessofanatmospheric
boundarylayer.Totakeintoaccountthenacelleeffectsonarelativelycoarsemesh,inthissectionwepresentanactuator
surfacemodelforparametrizingthenacelleeffectwithoutdirectlysimulatingtheboundarylayerflowsoverthenacelle.
Inthismodel,thenacellegeometryisrepresentedbytheactualsurfaceofthenacellewithdistributedforces.Asshown
inFig.2,theforceactingontheactuatorsurfaceisdecomposedintothenormalcomponentandthetangentialcomponent.
Fromaphysicalpointofview,theactuatorsurfaceofthenacellecannotactasasinkorsourceofmassnotmatterhow
theboundarylayeroverthenacelleismodelled.Basedonthisconsideration,thenormalcomponent oftheforceonthe
actuatorsurfaceiscomputedbysatisfyingthenon-penetrationconstraint.Thenon-penetrationconstraintcanbemetby
directly reconstructing the wall-normal velocity near the actuator surface as in the sharp interface immersed boundary
method[22,26,27,28].However,aphysicallyreasonableinterpolationschemeisdifficulttoimplementbecauseofthe
verycoarsegirdemployedintheactuatortypesimulations.Moreoverdirectvelocityreconstructioncannotensureasmooth
representationofthenacellegeometrybecauseoftheverycoarsegrid.Inthisworkwerepresentthenacellegeometryasa
diffusedinterfacethesameasinthedirectforcingimmersedboundarymethod[23,24,29].Inthisactuatorsurfacemodel
fornacelle,thenormalcomponentoftheforceactingonthesurfaceperunitareaiscalculatedinthesamewayasinthe
directforcingimmersedboundarymethod,whichisexpressedasfollows:
h −ud(X)+u˜(X) ·e (X)
f (X)= n e (X), (19)
n ∆t n
(cid:0) (cid:1)
whereud(X)isthedesiredvelocityonthenacellesurface,e (X)istheunitvectorinthenormaldirectionofthenacelle
n
surface,h=(h h h )1/3isthelengthscaleofthelocalbackgroundgridspacing(differentvaluesofhhavebeentested
x y z
withoutnoticinganysignificantdifferences),u˜(X)istheestimatedvelocityontheactuatorsurface,whichisinterpolated
fromthecorrespondingestimatedvelocityonthebackgroundmeshusingEq.(7)withtheu˜(x)onthebackgroundmesh
computedas
u˜(x)=un(x)+rhsn(x)∆t, (20)
where∆tisthesizeofthetimestep,theright-hand-sidetermrhsnincludestheconvection,pressuregradientanddiffusion
termscomputedfromthequantitiesofprevioustimestepn.
Thetangential force acting on the surface depends on the incoming velocity, the surface geometry and the complex
near-wallturbulence.Neitherofthelasttwocanbedirectlycapturedintheactuatortypesimulationsbyusingtheno-slip
boundaryconditionsortheshearstressboundaryconditionsforwallmodels[30,31].Inthepresentactuatorsurfacemodel,
weassumethatthetangentialforceisproportionaltothelocalincomingvelocity,andtheeffectsofsurfacegeometryand
near-wallturbulencecanbeparameterizedusingasingleparameter,thefrictioncoefficientc .Basedonthisassumption,
f
thetangentialforceactingonthesurfaceperunitareaiscomputedas
1
f (X)= c U2e (X) (21)
τ 2 f τ
WindEnerg.2017;00:1–17(cid:13)c 2017JohnWiley&Sons,Ltd. 5
DOI:10.1002/we
Preparedusingweauth.cls
Anewclassofactuatorsurfacemodelsforwindturbines XiaoleiYangandFotisSotiropoulos
wherec iscalculatedfromtheempiricalrelationproposedbyF.Schultz-Grunow[1]forturbulentboundarylayerswith
f
zeropressuregradient,whichisexpressedasfollows:
c =0.37(logRe )−2.584, (22)
f x
where Re is the Reynolds number based on the incoming velocity and the distance from the upstream edge of the
x
immersedbody.Itisnoticedthatthisexpressionisinvalidfortheregionwithalargepressuregradient.Forthepresent
nacellesimulation,however,zeropressuregradientcanberegardedasareasonableassumptionforthemostpartofthe
turbine nacelle. The framework presented in this work is applicable to more complex geometries if the corresponding
distributionofc isavailablefromexperimentsorhighfidelitysimulations.InEq.(21),U isthereferencevelocityofthe
f
incomingflow.Inthepresentwork,itisselectedasthemagnitudeoftheincomingstreamwisevelocity.Thedirectionof
thetangentialforcee (X),ontheotherhand,isdeterminedbythelocaltangentialvelocityrelativetothevelocityatthe
τ
nacellesurfaceatpointslocatedathawayfromthewall,i.e.,pointBinFig.2,asfollows:
|u(X+he (X))|
e (X)= n . (23)
τ u(X+he (X))
n
Aftercomputing theforcesontheactuator surface ofthenacelle,theforcesonthebackground mesharedistributed
fromtheactuatorsurfacemeshinthesamewayasinEq.(18).Itisnotedthattheforcesaredistributedovertheclosest
fivecellsdefinedbythewidthoftheemployedsmoothedcosinediscretedeltafunction.Thisforcedistributionwidthwill
beverysmallifwehaveasufficientlyhighgridresolutiontoresolvetheboundarylayer.Forthegridsusingintheactuator
typesimulations,however,thisforcedistributionwidthissolargethatitmaycausenon-physicalvelocityfieldnearthe
nacelleandaffectthecalculationoftherelativeincomingvelocityfortheblade.Toremedythisproblem,wesimplyletthe
relativeincomingvelocity(V inEq.(4))ofthetwoclosestradiallocationstothenacellesurfacebeequaltothatofthe
rel
thirdclosestlocationawayfromthenacelle.
3. FLOWSOLVER
The governing equations are the three-dimensional, unsteady, filtered continuity and Navier-Stokes equations in
non-orthogonal, generalized, curvilinear coordinates, which read in compact tensor notation (repeated indices imply
summation)asfollows(i,j =1,2,3):
∂Ui
J =0, (24)
∂ξi
1 ∂Ui ξi ∂ µ ∂ gjk ∂u 1 ∂ ξjp 1∂τ
= l − (Uju)+ l − ( l )− lj +f , (25)
J ∂t J ∂ξj l ρ∂ξj (cid:18) J ∂ξk(cid:19) ρ∂ξj J ρ ∂ξj l!
wherex andξiaretheCartesianandcurvilinearcoordinates,respectively,ξi =∂ξi/∂x arethetransformationmetrics,
i l l
J istheJacobianofthegeometrictransformation,u istheithcomponentofthevelocityvectorinCartesiancoordinates,
i
Ui=(ξi /J)u isthecontravariantvolumeflux,gjk =ξjξk arethecomponentsofthecontravariantmetrictensor,ρis
m m l l
thedensity,µisthedynamicviscosity,pisthepressure,f(l=1,2,3)arethebodyforcesintroducedbytheturbineblade
l
and nacelle computed using the actuator surface models presented in Sec. 2, and τ represents the anisotropic part of
ij
thesubgrid-scalestresstensor,whichismodelledbythedynamiceddyviscositysubgrid-scalestressmodel[32].Forthe
discretizationandsolversoftheaboveequations,thereaderisreferredto[33,7].
4. TESTCASES
Inthissection,weevaluatetheperformanceoftheproposedactuatorsurfacemodelsforturbinebladeandnacelle.The
nacellebodycaseinSec.4.1istoevaluatethecapabilityoftheactuatorsurfacemodelfornacelle.TheMEXICO(Model
experiments in Controlled Conditions) turbine [34, 35, 36] case in Sec. 4.2 is to show the capability of the actuator
surfacemodelsinpredictingthetipvorticesandvelocitydistributioninthenearwake.Thecapabilityofthenewactuator
surfacemodelsonthewakemeanderingpredictionisdemonstratedbythehydrokineticturbinecaseinSec.4.3.Inturbine
simulations the proposed actuator surface model is abbreviated as AS-BN when both the blade model and the nacelle
modelareemployed,andAS-Bwhenonlytheblademodelisemployed.
6 WindEnerg.2017;00:1–17(cid:13)c 2017JohnWiley&Sons,Ltd.
DOI:10.1002/we
Preparedusingweauth.cls
XiaoleiYangandFotisSotiropoulos Anewclassofactuatorsurfacemodelsforwindturbines
Figure3.Schematicofthenacellegeometrywiththetriangularsurfacemesh.Theflowisfromlefttoright.
4.1. Flowoverperiodicplacednacelles
In this section, the actuator surface model for nacelle is applied to simulate the flow over periodically placed nacelles
tomimic windfarm scenarios. The nacelle body iscomposed of a hemisphere facing upstream and a circular cylinder
in its downstream as shown in Fig. 3. Free-slip boundary condition is applied in the crosswise directions. Periodic
boundary condition is employed in the streamwise direction. The Reynolds number based on the cylinder diameter D
andthefreestreamvelocityis2000. Thesizeofthecomputational domainisL ×L ×L =30D×20D×20D,in
x y z
whichx,y andz represent thestreamwise,spanwiseandverticaldirections,respectively. .Thedownstreamendof the
nacelle body is located at the origin of the streamwise coordinate. Two grids are employed: the number of grid nodes
isN ×N ×N =502×348×348andN ×N ×N =153×80×80forthefineandcoarsegrids,respectively.
x y z x y z
Themesh isuniformintheregionaround thenacelle. Forthefinegrid, thegridspacing near the body isD/16 inthe
streamwise direction and D/66 in the other two directions; while for the coarse grid, the grid spacing is D/5 in the
streamwisedirection and D/10 in theother two directions. The employed timestep is0.01D/U and 0.1D/U for the
finegridandcoarsegridsimulations,respectively. Large-eddysimulationwiththedynamicsubgrid-scalemodel[32]is
employedfortheflowsimulation[33].Onthecoarsegrid,theactuatortypelarge-eddysimulationiscarriedoutusingthe
proposed actuator surface model for nacelleparameterization. Onthefinegrid, the so-calledgeometry-resolving, wall-
resolved large-eddy simulation is carried out with the sharp-interface immersed boundary method developed by Ge et
al.[28]forresolvingtheboundarylayerandnearwallturbulence,fromwhichthecomputedresultsisemployedtogauge
the performance of the proposed actuator surface model. It is noticed that two types of geometry-resolving large-eddy
simulationshavebeenmentioned:thegeometry-resolving,wall-modeledlarge-eddysimulationasemployedin[17]and
thegeometry-resolving,wall-resolvedlarge-eddysimulationemployedinthefinegridnacellesimulation.
WeexaminetheaccuracyoftheactuatorsurfacemodelinFigs.4and5bycomparingthemodelpredictedtime-averaged
streamwisevelocityhuiandturbulencekineticenergykprofileswiththosefromthefinegridwall-resolvedLESatdifferent
downstreamlocations.AsseeninFig.4,thevelocitydeficitcomputedfromactuatorsurfacemodelagreeswellwiththat
fromthewall-resolvedLESat1Ddownstream.At3D,5Dand7Ddownstreamlocations,however,thevelocitydeficits
computedfromtheactuatorsurfacemodelaresignificantlylarger.Forfurtherdownstreamlocations,ontheotherhand,the
velocitydeficitspredictedbytheactuatorsurfacemodelshowgoodagreement withthatcomputedbythewall-resolved
LES.Similarwiththeobservationsonthevelocityprofiles,largediscrepanciesareobservedatnearwakelocationsforthe
turbulentkineticenergyprofiles.AsshowninFig.5,notonlythemagnitudeoftheturbulencekineticenergybutalsothe
shapeofthedistributioncomputedfromtheactuatorsurfacemodelaredifferentfromthewall-resolvedsimulationresults
at near wake locations. At 1D downstream location, the peaks appear withinthe top and bottom shear layers for both
actuatorsurfacesimulationandwall-resolvedLES.At3D,5D,7Dand9Ddownstreamlocations,ontheotherhand,the
regionswithhighturbulencekineticenergyappearnearthecenterofthewakeforthewall-resolvedLES;whiletheystill
existwithintheshearlayersfortheactuatorsurfacesimulation.Itisnoticedthattheturbulencekineticenergypredicted
bythewall-resolvedLESissignificantlylargerthanthatfromtheactuatorsurfacesimulationat3Ddownstreamfromthe
body.Turbulentfluctuationshelpwakerecovery.Thisexplainswhythevelocityrecoversfasterinthenearwakeregion
for the wall-resolved LES. The actuator surface model, on the other hand, underpredicts the turbulence kinetic energy
inthenearwakeandthusunderpredictsthewakerecoveryatthoselocations.Atfarwakelocations,ontheotherhand,
theturbulencekineticenergypredictedbytheactuatorsurfacemodelisingoodagreementwiththepredictionsfromthe
wall-resolvedLES.
Tofurtherevaluatetheperformanceoftheactuatorsurfacemodel,thepowerspectraldensity(PSD)predictedbythe
actuatorsurfacesimulationiscomparedwiththatfromthewall-resolvedLESatdifferentdownstreamlocationsinFig.6.
Itisseenthelowfrequency,energy-containingpartofthePSDpredictedbytheactuatorsurfacemodelshowsanoverall
goodagreementwiththewall-resolvedLES.Boththeactuatorsurfaceandwall-resolvedsimulationspredicttheinertial
rangeofthePSDwiththeslopeof−5/3.Whiletheinertialrangecomputedbytheactuatormodelhasashorterwidth
andstartsatlowerfrequenciesinthenearwakeregion.Thismakestheenergyathighfrequenciessignificantlysmaller
fortheactuatorsurfacesimulationsatnearwakelocations.Atfarwakelocations,bothenergy-containingpartandinertial
partofthePSDpredictedbytheactuatorsurfacemodelagreewellwiththatcomputedfromthewall-resolvedLES,even
WindEnerg.2017;00:1–17(cid:13)c 2017JohnWiley&Sons,Ltd. 7
DOI:10.1002/we
Preparedusingweauth.cls
Anewclassofactuatorsurfacemodelsforwindturbines XiaoleiYangandFotisSotiropoulos
Figure4.Verticalprofilesoftime-averagedstreamwisevelocityhuiatdifferentdownstreamlocationsforthenacellecase.Symbols:
wall-resolvedLES;Solidlines:actuatorsurfacesimulation.
Figure5.Verticalprofilesofturbulencekineticenergykatdifferentdownstreamlocationsforthenacellecase.Symbols:wall-resolved
LES;Solidlines:actuatorsurfacesimulation.
thoughtheveryhighfrequency,dissipationrangeisnotresolvedbytheverycoarsegridemployedintheactuatorsurface
simulation.
8 WindEnerg.2017;00:1–17(cid:13)c 2017JohnWiley&Sons,Ltd.
DOI:10.1002/we
Preparedusingweauth.cls
XiaoleiYangandFotisSotiropoulos Anewclassofactuatorsurfacemodelsforwindturbines
Figure6.Powerspectraldensity(PSD)ofstreamwisevelocityfluctuationsatdifferentdownstreamlocationsforthenacellecase.Red
solidlines:actuatorsurfacesimulation;Greendash-dotlines:wall-resolvedLES.Theblackdashedlineshowstheslopeof−5/3.
4.2. FlowovertheMEXICOturbine
Inthissection,wesimulatetheflowovertheMEXICOturbinetoevaluatetheperformanceoftheproposedactuatorsurface
modelsonpredictingthevelocityandtipvorticesatturbinenearwakelocations.Thesizeofthecomputational domain
is L ×L ×L =4.2D×4D×4D, in which x, y and z represent the streamwise, spanwise, vertical directions,
x y z
respectively.Thenumber ofgridnodesisN ×N ×N =673×502×502.Themeshisuniformacrosstheturbine
x y z
andinthestreamwisedirectionwithgridspacingD/160,andisstretchedinthespanwiseandverticaldirectionsatregions
awayfromtheturbine.ThesizeoftimestepisT/750,whereT isrotationalperiodoftheturbinerotor.Free-slipboundary
condition is applied at the boundaries in the spanwise and vertical directions. The flow at the inlet is uniform. Three
caseswiththreedifferenttipspeedratios(TSR=ΩR/U whereΩistherotationalspeedoftherotor,Ristheradiusof
therotorandU istheincomingwindspeed), i.e.4.2,6.7and10,arecarriedout.Intheexperiment,differenttipspeed
ratioswererealizedbyadjustingtheincomingwindspeed.IntheLES,ontheotherhand,differentrotationalspeedsare
employedfordifferenttipspeedratios.TheReynoldsnumberbasedontheincomingwindspeedandtherotordiameteris
Re=6.9×106thesameforthecaseswithdifferenttipspeedratios.
Beforepresentingthecomputedresultsintheturbinenearwake,theboundcirculationsalongthebladepredictedby
theactuatorsurfacemodelarecomparedwithmeasurementsinFig.7.TheboundcirculationiscomputedbyΓ= L .
ρ|Vrel|
As seen, the location with maximum bound circulation is located near the blade root for TSR=4.2, in the middle for
TSR=6.7,whileclosetothetipforTSR=10.SuchdependenceoftheboundcirculationdistributionsonthevaluesofTSR
iswellcapturedbytheproposedactuatorsurfacemodel.Theamplitudesofthepredictedboundcirculationalsoshowan
acceptableagreementwiththemeasurementsexceptfortheTSR=10case,wherethediscrepanciesbetweenthecomputed
valuesandmeasuredonesaresomewhatlargerthantheothertwocases.
The capability of the proposed actuator surface model in predicting the tipvortices is examined in Fig. 8, in which
weshowthecomparisonofthecomputedtipvortexcirculationswiththemeasurementsforthethreecaseswithdifferent
tipspeedratios.Tocalculatethecirculationoftipvortices,thecenterofatipvortexcoreisselectedasthepositionwith
themaximumvorticitymagnitude.Thecirculationiscomputedbyintegratingthevelocityalongacirclecenteredatthe
vortex core center with radius of 0.05D. The computed tipvortex circulations show overall good agreements with the
measurements.SomediscrepanciesareobservedfortheTSR=10caseforthethirdandthefourthtipvorticesdownwind
fromtherotor.FromthecomparisonofthetipvortexcirculationsamongthethreeTSRcases,itisinterestingtonotice
thattheyareinthesameorderwhennormalizedusingU andD,althoughthedimensionalcirculationsaresignificantly
differentforthethreecases.
WindEnerg.2017;00:1–17(cid:13)c 2017JohnWiley&Sons,Ltd. 9
DOI:10.1002/we
Preparedusingweauth.cls
Anewclassofactuatorsurfacemodelsforwindturbines XiaoleiYangandFotisSotiropoulos
Figure 7.Bound circulation distributions along blade radial direction for (a) TSR=4.2, (b) TSR=6.7 and (c) TSR=10. Circles:
measurements;Redsolidline:actuatorsurfacemodel.
Figure8.Circulationoftipvorticesatdifferentdownwindlocationsfor(a)TSR=4.2,(b)TSR=6.7and(c)TSR=10.Blacklinewith
circles:measurements;Redlinewithsquares:actuatorsurfacemodel.
Thecomputedaxial(hu i)andradial(hu i)velocitiesarecomparedwithmeasurementsatdifferentradiallocationsin
a r
Fig.9.Asseen,theradialvelocityhu iincreasesastheflowapproachesturbine,anddecreasesrapidlyintheimmediate
r
downwindoftheturbineforallthefourconsideredradiallocationsforallthethreecases.Itisalsoseenthatsuchincrease
of the radial velocity issomewhat larger when TSRis higher. From these comparisons, we can see that the developed
actuatorsurfacemodelcanaccuratelypredictthedownwindvariationsofradialvelocity.Fortheaxialvelocity,asudden
decrease and increase inthe near wake of the turbine areobserved at r/D=0.52 for all the threetip-speed ratiosfor
themeasureddata.Thisis,however,becauseofthecorruptedPIVphotosattheselocationscausedbythelightreflected
ontheturbinehubasnoticedinthepaper byNilssonetal.[36].TheexperimentaldatafromthecorruptedPIVimages
willbeignoredforthecomparisons.Itisseenthatoverallgoodagreementswithmeasurementsareobtainedfortheaxial
velocitiesalthoughsomediscrepanciesareobservedfortheTSR=10case.
4.3. Flowoverahydrokineticturbine
Thecapabilityofthedevelopedactuatorsurfacemodelsinpredictingthewakemeanderingatfarwakelocationsisassessed
bysimulatingtheturbulentflowovertheaxial-flowhydrokineticturbineatSaintAnthonyFallsLaboratory,Universityof
Minnesota [37]. Kang, Yang and Sotiropoulos [17] carried out the geometry-resolving, actuator disk and actuator line
simulationsofthisturbine.Itwasfoundthatthestandardactuatorlinemodelwithoutanacellemodelunderpredictsthe
turbulenceintensityinthefarwakebecauseitcannotaccuratelycapturetheinner-outerwakeinteractionandthusthewake
meanderinginthefarwake.
The computational domain is L ×L ×L =16D×5.5D×2.3D in the streamwise, spanwise and vertical
x y z
directions, respectively, where D=0.5 m is the rotor diameter. Two different grids are employed: a coarse grid with
the grid spacing D/50 and a fine grid with grid spacing D/100 around the turbine. The number of grid nodes is
N ×N ×N =802×277×116, and 556×341×201 inthestreamwise, spanwise and vertical directions for the
x y z
coarseandfinegridsrespectively.Thetimestepis0.002D/U and0.001D/U,whereU=0.4m/sisthebulkvelocity,for
10 WindEnerg.2017;00:1–17(cid:13)c 2017JohnWiley&Sons,Ltd.
DOI:10.1002/we
Preparedusingweauth.cls