/
olaf.py
211 lines (179 loc) · 9.39 KB
/
olaf.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
"""
Tools to work with OLAF the vortex code implemented in openfast
"""
import numpy as np
def OLAFParams(omega_rpm, U0, R, a=0.3, aScale=1.2,
deltaPsiDeg=6, nPerRot=None,
targetFreeWakeLengthD=1,
targetWakeLengthD=4.,
nNWrot=8,
nNWrotFree=1,
nFWrot=0, nFWrotFree=0,
verbose=True, dt_glue_code=None):
"""
Computes recommended time step and wake length for OLAF based on:
INPUTS:
- omega_rpm: rotational speed [RPM]
- U0: mean wind speed [m/s]
- R: rotor radius [m]
OPTIONS FOR TIME STEP:
- either:
- deltaPsiDeg : target azimuthal discretization [deg]
or
- nPerRot : number of time step per rotations.
deltaPsiDeg - nPerRot
5 72
6 60
7 51.5
8 45
- dt_glue_code: glue code time step. If provided, the time step of OLAF will be approximated
such that it is a multiple of the glue-code time step.
OPTIONS FOR WAKE LENGTH:
- a: average axial induction factor at the rotor [-]
- aScale: scaling factor to estimate induction, such that the wake convection velocity is:
Uc=U0(1-aScale*a)
- targetWakeLengthD: target wake length in diameter [D]
- nNWrot : minimum number of near wake rotations
- nFWrot : minimum number of far wake rotations
- nFWrotFree : minimum number of far wake rotations (free panels)
"""
def myprint(*args, **kwargs):
if verbose:
print(*args, **kwargs)
# Rotational period
omega = omega_rpm*2*np.pi/60
T = 2*np.pi/omega
# Convection velocity
Uc = U0 * (1-aScale*a)
# Desired time step
if nPerRot is not None:
dt_wanted = np.around(T/nPerRot,5)
deltaPsiDeg = np.around(omega*dt_wanted*180/np.pi ,2)
else:
dt_wanted = np.around(deltaPsiDeg/(6*omega_rpm),5)
nPerRot = int(2*np.pi /(deltaPsiDeg*np.pi/180))
# Adapting desired time step based on glue code time step
if dt_glue_code is not None:
dt_rounded = round(dt_wanted/dt_glue_code)*dt_glue_code
deltaPsiDeg2 = np.around(omega*dt_rounded *180/np.pi ,2)
myprint('>>> To satisfy glue-code dt:')
myprint(' Rounding dt from {} to {}'.format(dt_wanted, dt_rounded ))
myprint(' Changing dpsi from {} to {}'.format(deltaPsiDeg, deltaPsiDeg2))
dt_fvw = dt_rounded
deltaPsiDeg = deltaPsiDeg2
nPerRot = int(2*np.pi /(deltaPsiDeg*np.pi/180))
else:
dt_fvw = dt_wanted
# Useful functions
n2L = lambda n: (n * dt_fvw * Uc)/(2*R) # convert number of panels to distance
n2R = lambda n: n * dt_fvw / T # convert number of panels to number of rotations
n2t = lambda n: n * dt_fvw # convert number of panels to time
# All Wake (AW) panels - Wake length from mean wind speed
targetWakeLength = targetWakeLengthD * 2 * R
nAWPanels_FromU0 = int(targetWakeLength / (Uc*dt_fvw))
# Free near wake panels (based on distance)
targetFreeWakeLength = targetFreeWakeLengthD * 2 * R
nNWPanelsFree_FromU0 = int(targetFreeWakeLength / (Uc*dt_fvw))
nNWPanelsFree_FromRot = int(nNWrotFree*nPerRot)
# Far wake (FW) panels, always from number of rotations
nFWPanels = int(nFWrot*nPerRot)
nFWPanelsFree = int(nFWrotFree*nPerRot)
# Wake length from rotational speed and number of rotations
nAWPanels_FromRot = int(nNWrot*nPerRot) # Total number of panels NW+FW
# Below we chose between criteria on number of rotation or donwstream distance
# This can be adapted/improved
myprint('Number of panels (NW free) from wind speed and distance:{:15d}'.format(nNWPanelsFree_FromU0))
myprint('Number of panels (NW free) from number of rotations :{:15d}'.format(nNWPanelsFree_FromRot))
myprint('Number of panels (NW+FW) from wind speed and distance:{:15d}'.format(nAWPanels_FromU0))
myprint('Number of panels (NW+FW) from number of rotations :{:15d}'.format(nAWPanels_FromRot))
myprint('Number of panels (NW+FW) from average between two :{:15d}'.format(int((nAWPanels_FromRot+nAWPanels_FromU0)/2)))
if nAWPanels_FromRot>nAWPanels_FromU0:
# Criteria based on rotation wins:
myprint('[INFO] Using number of rotations to setup number of panels')
nAWPanels = nAWPanels_FromRot # Total number of panels NW+FW
else:
myprint('[INFO] Using wind speed and distance to setup number of panels')
# Wake distance wins, we keep the nFW from rot but increase nNW
nAWPanels = nAWPanels_FromU0 # Total number of panels NW+FW
if nNWPanelsFree_FromRot>nNWPanelsFree_FromU0:
# Criteria based on rotation wins:
myprint('[INFO] Using number of rotations to setup number of free panels')
nNWPanelsFree = nNWPanelsFree_FromRot
else:
myprint('[INFO] Using wind speed and distance to setup number of free panels')
# Wake distance wins, we keep the nFW from rot but increase nNW
nNWPanelsFree = nNWPanelsFree_FromU0
nNWPanels = nAWPanels - nFWPanels # nNW = All-Far Wake
# See "free" near wake
if nNWPanelsFree>nNWPanels:
nNWPanelsFree=nNWPanels
myprint('[INFO] Capping number of free NW panels to max.')
if nNWPanelsFree<nNWPanels and nFWPanelsFree>0:
nFWPanelsFree=0
myprint('[INFO] Setting number of Free FW panels to zero because a frozen near wake is used')
# Transient time (twice the time to develop the full wake extent)
# This is the minimum recommended time before convergence of the wake is expected
# (might be quite long)
tMin = 2 * dt_fvw*nAWPanels
if verbose:
myprint('')
myprint('{:15.2f} Transient time ({:5.1f} rot)'.format(tMin, tMin/T))
myprint('{:15d} nAWPanels ({:5.1f} rot, {:5.1f}D)'.format(nAWPanels, n2R(nAWPanels), n2L(nAWPanels)))
myprint('')
myprint('OLAF INPUT FILE:')
myprint('----------------------- GENERAL OPTIONS ---------------------')
myprint('{:15.6f} DTFVW (delta psi = {:5.1f}deg)'.format(dt_fvw, deltaPsiDeg))
myprint('--------------- WAKE EXTENT AND DISCRETIZATION --------------')
myprint('{:15d} nNWPanels ({:5.1f} rot, {:5.1f}D)'.format(nNWPanels , n2R(nNWPanels ), n2L(nNWPanels )))
myprint('{:15d} nNWPanelsFree ({:5.1f} rot, {:5.1f}D {:5.1f})'.format(nNWPanelsFree, n2R(nNWPanelsFree), n2L(nNWPanelsFree), n2t(nNWPanelsFree)))
myprint('{:15d} nFWPanels ({:5.1f} rot, {:5.1f}D)'.format(nFWPanels , n2R(nFWPanels ), n2L(nFWPanels )))
myprint('{:15d} nFWPanelsFree ({:5.1f} rot, {:5.1f}D)'.format(nFWPanelsFree, n2R(nFWPanelsFree), n2L(nFWPanelsFree)))
return dt_fvw, tMin, nNWPanels, nNWPanelsFree, nFWPanels, nFWPanelsFree
def OLAFParamsRPM(omega_rpm, deltaPsiDeg=6, nNWrot=2, nFWrot=10, nFWrotFree=3, nPerRot=None, verbose=True, dt_glue_code=None, totalRot=None):
"""
Computes recommended time step and wake length based on the rotational speed in RPM
INPUTS:
- omega_rpm: rotational speed in RPM
- deltaPsiDeg : azimuthal discretization in deg
- nNWrot : number of near wake rotations
- nFWrot : total number of far wake rotations
- nFWrotFree : number of far wake rotations that are free
deltaPsiDeg - nPerRot
5 72
6 60
7 51.5
8 45
"""
omega_rpm = np.asarray(omega_rpm)
omega = omega_rpm*2*np.pi/60
T = 2*np.pi/omega
if nPerRot is not None:
dt_wanted = np.around(T/nPerRot,5)
deltaPsiDeg = np.around(omega*dt_wanted*180/np.pi ,2)
else:
dt_wanted = np.around(deltaPsiDeg/(6*omega_rpm),5)
nPerRot = int(2*np.pi /(deltaPsiDeg*np.pi/180))
if dt_glue_code is not None:
dt_rounded = round(dt_wanted/dt_glue_code)*dt_glue_code
deltaPsiDeg2 = np.around(omega*dt_rounded *180/np.pi ,2)
print('>>> To satisfy glue-code dt:')
print(' Rounding dt from {} to {}'.format(dt_wanted, dt_rounded ))
print(' Changing dpsi from {} to {}'.format(deltaPsiDeg, deltaPsiDeg2))
dt_wanted = dt_rounded
deltaPsiDeg = deltaPsiDeg2
nPerRot = int(2*np.pi /(deltaPsiDeg*np.pi/180))
nNWPanel = int(nNWrot*nPerRot)
nFWPanel = int(nFWrot*nPerRot)
nFWPanelFree = int(nFWrotFree*nPerRot)
if totalRot is None:
totalRot = (nNWrot + nFWrot)*3 # going three-times through the entire wake
tMax = dt_wanted*nPerRot*totalRot
if verbose:
print(dt_wanted , ' dt')
print(int (nNWPanel ), ' nNWPanel ({} rotations)'.format(nNWrot))
print(int (nFWPanel ), ' FarWakeLength ({} rotations)'.format(nFWrot))
print(int (nFWPanelFree), ' FreeFarWakeLength ({} rotations)'.format(nFWrotFree))
print(tMax , ' Tmax ({} rotations)'.format(totalRot))
return dt_wanted, tMax, nNWPanel, nFWPanel, nFWPanelFree
if __name__ == '__main__':
OLAFParams(omega_rpm = 4.87558, U0=8, R=63, deltaPsiDeg=6, verbose=True)