-
Notifications
You must be signed in to change notification settings - Fork 12
/
dynstall_legacy.f90
393 lines (315 loc) · 11.7 KB
/
dynstall_legacy.f90
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
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
module dynstall_legacy
use airfoils
implicit none
! Leishman-Beddoes model parameters
type LB_Type
logical :: StallFlag = .false.
real(mytype) :: dp
real(mytype) :: dF
real(mytype) :: dCNv
real(mytype) :: sLEv
integer :: LESepState
real(mytype) :: CLRef
real(mytype) :: CLRefLE
real(mytype) :: CLA
real(mytype) :: CLCritP
real(mytype) :: CLCritN
integer :: CLRateFlag
real(mytype) :: Fstat
real(mytype) :: F
real(mytype) :: cv
real(mytype) :: dcv
real(mytype) :: CLRef_Last
real(mytype) :: CLRefLE_Last
real(mytype) :: Fstat_Last
real(mytype) :: cv_Last
! Additional LB diagnostic output
integer, dimension(9) :: LB_LogicOutputs
integer :: Logic_W(9)=[1,1,1,1,1,3,2,1,1]
end type LB_Type
Type(LB_Type) :: lb
contains
subroutine dystl_init_LB(lb,dynstallfile)
implicit none
type(LB_Type) :: lb
character :: dynstallfile*80
real(mytype) :: CLcritp, CLcritn, CLalpha
NAMELIST/LBParam/CLcritp,CLcritn,CLalpha
lb%StallFlag = .true.
lb%dp=0.0
lb%dF=0.0
lb%dCNv=0.0
lb%LESepState=0.0
lb%sLEv=0.0
lb%CLRef_Last=0.0
lb%CLRefLE_Last=0.0
lb%Fstat_Last=1.0
lb%cv_Last=0.0
lb%LB_LogicOutputs(:)=0
open(30,file=dynstallfile)
read(30,nml=LBParam)
close(30)
lb%CLCritP=CLcritp
lb%CLCritN=CLcritn
lb%CLA=CLalpha
return
end subroutine dystl_init_LB
subroutine LB_EvalIdealCL(AOA,AOA0,CLa,RefFlag,CLID)
! AOA inputs in radians
! AOA0 is zero lift AOA
! RefFlag defines whether to output reference CL or ideal CL
! CLa is reference lift slope (per radian) to be used for reference CL (ideal CLa is 2*pi)
implicit none
real(mytype) :: AOA, AOA0, CLa
integer :: RefFlag
real(mytype) :: CLID
real(mytype) :: IDS, aID, d1, CLaI, aIDc
aID=AOA-AOA0
call Force180(aID)
! reflect function across axis for abs(aID)>90
IDS=1
if (aID>pi/2.0) then
aID=pi-aID
IDS=-1.0
else if (aID<-pi/2.0) then
aID=-pi-aID
IDS=-1.0
end if
! If RefFlag is 1, output reference CL, otherwise round off ideal CL at high AOA
if (RefFlag==1) then
CLID=IDS*CLa*aID
else
! round off the ideal CL after cutoff AOA
aIDc=30.0*conrad
d1=1.8
CLaI=2.0*pi
if (abs(aID)<aIDc) then
CLID=IDS*CLaI*aID
else
CLID=IDS*(CLaI*(aIDc-1.0/d1*sin(d1*aIDc))+CLaI/d1*sin(d1*aID))
end if
end if
end subroutine LB_EvalIdealCL
subroutine Force180(a)
implicit none
real(mytype) :: a
! alpha in radians
if (a>pi) then
a=a-2.0*pi
elseif (a<-pi) then
a=a+2.0*pi
endif
end subroutine Force180
SUBROUTINE LB_UpdateStates(lb,airfoil,Re,ds)
implicit none
type(LB_Type),intent(inout) :: lb
type(AirfoilType),intent(in) :: airfoil
real(mytype), intent(in) :: Re, ds
integer :: i, nei, j, IsBE
real(mytype) :: Tf,TfRef,Tp,TvRef, Tv, TvL
! Set model parameters. All of these are potentially a function of Mach
! and are set to low mach values...
Tp=4.0 ! time constant on LE pressure response to change in CL
TfRef=3.0 ! time constant on TE separation point travel
TvRef=6.3 ! time constant on LE vortex lift indicial function
TvL=11.0 ! Characteristic LE vortex travel time
! Eval LE separation state
! Note: CLCrit is critical (ideal) CL value for LE separation. This is
! approximately equal to the CL that would exist at the angle of attack at
! max CL if the CL curve had remained linear
if (lb%LESepState==0 .AND. (lb%CLRefLE>lb%CLCritP .OR. lb%CLRefLE<lb%CLCritN)) then
! In LE separation state
lb%LESepState=1
lb%sLEv=0 ! reset leading edge vortex time counter
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(5)=1
else if (lb%LESepState==1 .AND. (lb%CLRefLE<lb%CLCritP .AND. lb%CLRefLE>lb%CLCritN)) then
! Out of LE separation state
lb%LESepState=0
lb%sLEv=0 ! reset leading edge vortex time counter
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(5)=2
end if
! Set time constants based on LE separation state and TE separation point
! location. Different time constants for abs(CL) increasing vs. decreasing
if (lb%LESepState==1) then
if (lb%sLEv<TvL) then
if (lb%CLRateFlag>0) then
!Tf=3.0*TfRef ! original
Tf=4.0*TfRef
Tv=TvRef
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(8)=1
else
Tf=1.0/2.0*TfRef
Tv=1.0/2.0*TvRef
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(8)=2
end if
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(7)=1
else if (lb%sLEv<2.0*TvL) then
if (lb%CLRateFlag>0) then
! orig
!Tf=1.0/3.0*TfRef
!Tv=1.0/4.0*TvRef
Tf=2.0*TfRef
Tv=TvRef
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(8)=3
else
Tf=1.0/2.0*TfRef
Tv=1.0/2.0*TvRef
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(8)=4
end if
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(7)=2
else
! orig
!Tf=4.0*TfRef
!Tv=0.9*TvRef
Tf=TfRef
Tv=TvRef
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(7)=3
end if
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(6)=1
else
if (lb%F>0.7) then
Tf=TfRef
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(7)=4
else
Tf=2*TfRef
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(7)=5
end if
Tv=TvRef
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(6)=2
end if
! update LE vortex time counter if in LE separation state
if (lb%LESepState==1) then
lb%sLEv=lb%sLEv+ds
! Set logic state flags (for model diagnosis output)
lb%LB_LogicOutputs(9)=1
end if
! Update states, first order lag equations, exponential recursion form (midpoint rule version)
lb%dp=lb%dp*exp(-ds/Tp)+(lb%CLRef-lb%CLRef_Last)*exp(-ds/(2*Tp))
lb%dF=lb%dF*exp(-ds/Tf)+(lb%Fstat-lb%Fstat_Last)*exp(-ds/(2*Tf))
lb%dCNv=lb%dCNv*exp(-ds/Tv)+lb%dcv*exp(-ds/(2*Tv))
! update lagged values
lb%CLRef_Last=lb%CLRef
lb%CLRefLE_Last=lb%CLRefLE
lb%Fstat_Last=lb%Fstat
lb%cv_Last=lb%cv
End SUBROUTINE LB_UpdateStates
SUBROUTINE LB_LogicChecksum(LBCheck)
implicit none
integer :: LBCheck
integer :: Loop
! Calculates a checksum for the logic states in the LB model
LBCheck=0
do Loop=1,9
LBCheck=LBCheck+lb%LB_LogicOutputs(Loop)*lb%Logic_W(Loop)
end do
End SUBROUTINE LB_LogicChecksum
subroutine LB_DynStall(airfoil,lb,CLstat,CDstat,alphaL,alpha5,Re,CL,CD)
! GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
! Routine that computes the Leishmann-Beddoes dynamic stall model
! with incompressible reductionand returns corrected values for
! CL and CD having taken into account the dynamic stall effects
! GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
implicit none
type(AirfoilType) :: airfoil ! Airfoil structure
type(LB_Type) :: lb ! Leishmann-Beddoes model structure
real(mytype) :: CLstat, CDstat, alphaL, alpha5, Re, CL, CD
real(mytype) :: AOA0, CLID, Trans, dCLRefLE, dAOARefLE, AOARefLE, CLstatF, C, C1, CLIDF
real(mytype) :: CLRatio, CLsep, CLF, dCDF, KD, CLa, NOF, dCLv, dCDv, acut, CLCritP, CLCritN
! Airfoil data
AOA0=airfoil%alzer
! Model constants
KD=0.1 ! Trailing Edge separation drag factor
! Evaluate the ideal CL curve at current AOA
call LB_EvalIdealCL(alphaL,AOA0,lb%CLa,1,lb%CLRef)
call LB_EvalIdealCL(alphaL,AOA0,lb%CLa,0,CLID)
! calc lagged ideal CL for comparison with critical LE separation CL
Trans=(cos(alphaL-AOA0))**2 ! fair effect to zero at 90 deg. AOA...
dCLRefLE=Trans*lb%dp ! dp is lagged CLRef change
dAOARefLE=dCLRefLE/CLa
! define reference LE CL and AOA
lb%CLRefLE=lb%CLRef-dCLRefLE
if (lb%CLRefLE*(lb%CLRefLE-lb%CLRefLE_Last) > 0) then
lb%CLRateFlag=1
else
lb%CLRateFlag=0
end if
AOARefLE=alphaL-dAOARefLE
Call Force180(AOARefLE)
! calc effective static TE separation point using effective LE AOA
Call EvalStaticCoeff(Re,AOARefLE*condeg,CLstatF,C,C1,airfoil)
Call LB_EvalIdealCL(AOARefLE,AOA0,CLa,0,CLIDF)
if (abs(CLIDF)<0.001) then
CLRatio=999
else
CLRatio=CLstatF/CLIDF;
end if
if (CLRatio > 0.25) then
lb%Fstat=min((sqrt(4.0*CLRatio)-1.0)**2,1.0)
! Test logic
lb%LB_LogicOutputs(1)=1
else
lb%Fstat=0
! Test logic
lb%LB_LogicOutputs(1)=2
end if
! calc lagged Fstat to represent dynamic TE separation point
lb%F=lb%Fstat-lb%dF
! force limits on lagged F (needed due to discretization error...)
lb%F=min(max(lb%F,0.0),1.0)
! Calc dynamic CL due to TE separation as fairing between fully attached and fully separated predictions from the Kirchoff approximation at current AOA
if (abs(CLID)<0.001) then
CLRatio=999
else
CLRatio=CLstat/CLID
end if
if (CLRatio > 1.0) then
CLID=CLstat
! Test logic
lb%LB_LogicOutputs(2)=1
end if
if (CLRatio > 0.25) then
CLsep=CLID/4.0
! Test logic
lb%LB_LogicOutputs(3)=1
else
CLsep=CLstat
! Test logic
lb%LB_LogicOutputs(3)=2
end if
CLF=CLsep+CLID*0.25*(lb%F+2.0*sqrt(lb%F))
dCDF=KD*(CLstat-CLF)*sign(1.0d0,CLstat)
! LE vortex lift component, dCNv is a lagged change in the added normal force due
! to LE vortex shedding. Assumed to affect lift coeff as an added circulation...
dCLv=lb%dCNv*cos(alpha5)
dCDv=lb%dCNv*sin(alpha5)
! vortex feed is given by the rate at which lift (circulation) is being shed due to dynamic separation. Lift component due to separation is defined by the
! difference between the ideal lift and the lift including dynamic separation effects.
lb%cv=CLID-CLF
lb%dcv=lb%cv-lb%cv_Last
! If the sign of dcv is opposite the reference LE CL, set to zero to disallow negative vorticity from shedding from the leading edge. Also, limit the model
! at AOA>acut or if the magnitude of the reference CL is decreasing...
acut=50.0*conrad
if (sign(1.0d0,lb%dcv*lb%CLRefLE)<0 .OR. abs(alphaL-AOA0)>acut .OR. lb%CLRateFlag<0) then
lb%dcv=0.0
! Test logic
lb%LB_LogicOutputs(4)=1
end if
! Total lift and drag
CL=CLF+dCLv
CD=CDstat+dCDF+dCDv
return
end subroutine LB_DynStall
end module dynstall_legacy