-
Notifications
You must be signed in to change notification settings - Fork 67
/
actuator_turbine_model.f90
1843 lines (1449 loc) · 72.7 KB
/
actuator_turbine_model.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
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Written by:
!!
!! Luis 'Tony' Martinez <[email protected]> (Johns Hopkins University)
!!
!! Copyright (C) 2012-2013, Johns Hopkins University
!!
!! This file is part of The Actuator Turbine Model Library.
!!
!! The Actuator Turbine Model is free software: you can redistribute it
!! and/or modify it under the terms of the GNU General Public License as
!! published by the Free Software Foundation, either version 3 of the
!! License, or (at your option) any later version.
!!
!! The Actuator Turbine Model is distributed in the hope that it will be
!! useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with Foobar. If not, see <http://www.gnu.org/licenses/>.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*******************************************************************************
module actuator_turbine_model
!*******************************************************************************
! This module has the subroutines to provide all calculations for use in the
! actuator turbine model (ATM)
! Imported modules
use atm_base ! Include basic types and precission of real numbers
use atm_input_util ! Utilities to read input files
implicit none
! Declare everything private except for subroutine which will be used
private
public :: atm_initialize, numberOfTurbines, &
atm_computeBladeForce, atm_update, &
vector_add, vector_divide, vector_mag, distance, &
atm_output, atm_process_output, &
atm_initialize_output, atm_computeNacelleForce, atm_write_restart, &
atm_compute_cl_correction
! The very crucial parameter pi
real(rprec), parameter :: pi=acos(-1._rprec)
! These are used to do unit conversions
real(rprec) :: degRad = pi/180._rprec ! Degrees to radians conversion
real(rprec) :: rpmRadSec = pi/30._rprec ! Set the revolutions/min to radians/s
logical :: pastFirstTimeStep ! Establishes if we are at the first time step
! Subroutines for the actuator turbine model
! All suboroutines names start with (atm_)
contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_initialize()
! This subroutine initializes the ATM. It calls the subroutines in
! atm_input_util to read the input data and creates the initial geometry
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
integer :: i
logical :: file_exists
pastFirstTimeStep=.false. ! The first time step not reached yet
write(*,*) 'Reading Actuator Turbine Model Input...'
call read_input_conf() ! Read input data
write(*,*) 'Done Reading Actuator Turbine Model Input'
do i = 1,numberOfTurbines
inquire(file = "./turbineOutput/"//trim(turbineArray(i) % turbineName)// &
"/actuatorPoints", exist=file_exists)
! Creates the ATM points defining the geometry
call atm_create_points(i)
! This will create the first yaw alignment
turbineArray(i) % deltaNacYaw = turbineArray(i) % nacYaw
call atm_yawNacelle(i)
if (file_exists .eqv. .true.) then
write(*,*) 'Reading bladePoints from Previous Simulation'
call atm_read_actuator_points(i)
endif
call atm_calculate_variables(i) ! Calculates variables depending on input
inquire(file = "./turbineOutput/"//trim(turbineArray(i) % turbineName)// &
"/restart", exist=file_exists)
if (file_exists .eqv. .true.) then
write(*,*) 'Reading Turbine Properties from Previous Simulation'
call atm_read_restart(i)
endif
end do
pastFirstTimeStep=.true. ! Past the first time step
end subroutine atm_initialize
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_read_actuator_points(i)
! This subroutine reads the location of the actuator points
! It is used if the simulation wants to start from a previous simulation
! without having to start the turbine from the original location
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
integer, intent(in) :: i ! Indicates the turbine number
integer :: j, m, n, q
j=turbineArray(i) % turbineTypeID ! The turbine type ID
open(unit=1, file="./turbineOutput/"//trim(turbineArray(i) % turbineName)// &
"/actuatorPoints", action='read')
do m=1, turbineModel(j) % numBl
do n=1, turbineArray(i) % numAnnulusSections
do q=1, turbineArray(i) % numBladePoints
read(1,*) turbineArray(i) % bladePoints(m,n,q,:)
enddo
enddo
enddo
close(1)
end subroutine atm_read_actuator_points
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_read_restart(i)
! This subroutine reads the rotor speed
! It is used if the simulation wants to start from a previous simulation
! without having to start the turbine from the original omega
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
integer, intent(in) :: i ! Indicates the turbine number
! Open the file at the last line (append)
open( unit=1, file="./turbineOutput/"//trim(turbineArray(i) % turbineName)// &
"/restart", action='read') !, position='append')
! Bring the pointer to the last line
!~ backspace 1
! Read past the first line
read(1,*)
! Read the restart variables
read(1,*) turbineArray(i) % rotSpeed
read(1,*) turbineArray(i) % torqueGen
read(1,*) turbineArray(i) % torqueRotor
read(1,*) turbineArray(i) % u_infinity
read(1,*) turbineArray(i) % induction_a
read(1,*) turbineArray(i) % PitchControlAngle
read(1,*) turbineArray(i) % IntSpeedError
read(1,*) turbineArray(i) % nacYaw
read(1,*) turbineArray(i) % rotorApex
read(1,*) turbineArray(i) % uvShaft
close(1)
write(*,*) ' RotSpeed Value from previous simulation is ', &
turbineArray(i) % rotSpeed
write(*,*) ' torqueGen Value from previous simulation is ', &
turbineArray(i) % torqueGen
write(*,*) ' torqueRotor Value from previous simulation is ', &
turbineArray(i) % torqueRotor
write(*,*) ' PitchControlAngle Value from previous simulation is ', &
turbineArray(i) % PitchControlAngle
write(*,*) ' IntSpeedError Value from previous simulation is ', &
turbineArray(i) % IntSpeedError
write(*,*) ' Yaw Value from previous simulation is ', &
turbineArray(i) % nacYaw
write(*,*) ' Rotor Apex Value from previous simulation is ', &
turbineArray(i) % rotorApex
write(*,*) ' uvShaft Value from previous simulation is ', &
turbineArray(i) % uvShaft
end subroutine atm_read_restart
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_write_restart(i)
! This subroutine reads the rotor speed
! It is used if the simulation wants to start from a previous simulation
! without having to start the turbine from the original omega
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
integer, intent(in) :: i ! Indicates the turbine number
integer :: pointsFile=787 ! File to write the actuator points
integer :: restartFile=21 ! File to write restart data
integer j, m,n,q ! counters
! Open the file
open( unit=restartFile, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/restart", status="replace")
write(restartFile,*) 'RotSpeed ', 'torqueGen ', 'torqueRotor ', 'u_infinity ', &
'induction_a ', 'PitchControlAngle ', 'IntSpeedError ', &
'nacYaw ', 'rotorApex ', 'uvShaft'
! Store the rotSpeed value
write(restartFile,*) turbineArray(i) % rotSpeed
write(restartFile,*) turbineArray(i) % torqueGen
write(restartFile,*) turbineArray(i) % torqueRotor
write(restartFile,*) turbineArray(i) % u_infinity
write(restartFile,*) turbineArray(i) % induction_a
write(restartFile,*) turbineArray(i) % PitchControlAngle
write(restartFile,*) turbineArray(i) % IntSpeedError
write(restartFile,*) turbineArray(i) % nacYaw
write(restartFile,*) turbineArray(i) % rotorApex
write(restartFile,*) turbineArray(i) % uvShaft
close(restartFile)
! Write the actuator points at every time-step regardless
j=turbineArray(i) % turbineTypeID ! The turbine type ID
open(unit=pointsFile, status="replace", file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/actuatorPoints")
do m=1, turbineModel(j) % numBl
do n=1, turbineArray(i) % numAnnulusSections
do q=1, turbineArray(i) % numBladePoints
! A new file will be created each time-step with the proper
! location of the blades
write(pointsFile,*) turbineArray(i) % bladePoints(m,n,q,:)
enddo
enddo
enddo
close(pointsFile)
end subroutine atm_write_restart
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_initialize_output()
! This subroutine initializes the output files for the ATM
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
logical :: file_exists
integer :: i
! Write to the screen output start
call atm_print_initialize()
do i = 1,numberOfTurbines
inquire(file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName),EXIST=file_exists)
if (file_exists .eqv. .false.) then
! Create turbineOutput directory
call system("mkdir -vp turbineOutput/"// &
trim(turbineArray(i) % turbineName))
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/power")
write(1,*) 'time PowerRotor powerGen '
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/thrust")
write(1,*) 'time thrust '
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/RotSpeed")
write(1,*) 'time RotSpeed'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/Yaw")
write(1,*) 'time deltaNacYaw NacYaw'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/lift")
write(1,*) 'turbineNumber bladeNumber '
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/drag")
write(1,*) 'turbineNumber bladeNumber '
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/Cl")
write(1,*) 'turbineNumber bladeNumber Cl'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/Cd")
write(1,*) 'turbineNumber bladeNumber Cd'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/alpha")
write(1,*) 'turbineNumber bladeNumber alpha'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/Vrel")
write(1,*) 'turbineNumber bladeNumber Vrel'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/Vaxial")
write(1,*) 'turbineNumber bladeNumber Vaxial'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/Vtangential")
write(1,*) 'turbineNumber bladeNumber Vtangential'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/tangentialForce")
write(1,*) 'turbineNumber bladeNumber tangentialForce'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/axialForce")
write(1,*) 'turbineNumber bladeNumber axialForce'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/nacelle")
write(1,*) 'time Velocity-no-correction Velocity-w-correction'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/uy_LES")
write(1,*) 'turbineNumber bladeNumber uy_LES'
close(1)
open(unit=1, file="./turbineOutput/"// &
trim(turbineArray(i) % turbineName)//"/uy_opt")
write(1,*) 'turbineNumber bladeNumber uy_opt'
close(1)
endif
enddo
end subroutine atm_initialize_output
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_create_points(i)
! This subroutine generate the set of blade points for each turbine
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
integer, intent(in) :: i ! Indicates the turbine number
integer :: j ! Indicates the turbine type
integer :: m ! Indicates the blade point number
integer :: n ! Indicates number of actuator section
integer :: k
real(rprec), dimension (3) :: root ! Location of rotor apex
real(rprec) :: beta ! Difference between coning angle and shaft tilt
real(rprec) :: dist ! Distance from each actuator point
integer, pointer :: numBladePoints
integer, pointer :: numBl
integer, pointer :: numAnnulusSections
real(rprec), pointer :: NacYaw
real(rprec), pointer :: db(:)
real(rprec), pointer :: bladePoints(:,:,:,:)
real(rprec), pointer :: bladeRadius(:,:,:)
real(rprec), pointer :: azimuth
real(rprec), pointer :: rotSpeed
real(rprec), pointer :: ShftTilt
real(rprec), pointer :: PreCone
real(rprec), pointer :: towerShaftIntersect(:)
real(rprec), pointer :: baseLocation(:)
real(rprec), pointer :: TowerHt
real(rprec), pointer :: Twr2Shft
real(rprec), pointer :: rotorApex(:)
real(rprec), pointer :: OverHang
real(rprec), pointer :: UndSling
real(rprec), pointer :: uvShaftDir
real(rprec), pointer :: uvShaft(:)
real(rprec), pointer :: uvTower(:)
real(rprec), pointer :: TipRad
real(rprec), pointer :: HubRad
real(rprec), pointer :: annulusSectionAngle
real(rprec), pointer :: solidity(:,:,:)
! Identifies the turbineModel being used
j=turbineArray(i) % turbineTypeID ! The type of turbine (eg. NREL5MW)
! Variables to be used locally. They are stored in local variables within the
! subroutine for easier code following. The values are then passed to the
! proper type
numBladePoints => turbineArray(i) % numBladePoints
numBl=>turbineModel(j) % numBl
numAnnulusSections=>turbineArray(i) % numAnnulusSections
! Allocate variables depending on specific turbine properties and general
! turbine model properties
allocate(turbineArray(i) % db(numBladePoints))
allocate(turbineArray(i) % bladePoints(numBl, numAnnulusSections, &
numBladePoints,3))
allocate(turbineArray(i) % bladeRadius(numBl,numAnnulusSections,numBladePoints))
allocate(turbineArray(i) % solidity(numBl,numAnnulusSections,numBladePoints))
! Assign Pointers turbineArray denpendent (i)
db=>turbineArray(i) % db
bladePoints=>turbineArray(i) % bladePoints
bladeRadius=>turbineArray(i) % bladeRadius
solidity=>turbineArray(i) % solidity
azimuth=>turbineArray(i) % azimuth
rotSpeed=>turbineArray(i) % rotSpeed
towerShaftIntersect=>turbineArray(i) % towerShaftIntersect
baseLocation=>turbineArray(i) % baseLocation
uvShaft=>turbineArray(i) % uvShaft
uvTower=>turbineArray(i) % uvTower
rotorApex=>turbineArray(i) % rotorApex
uvShaftDir=>turbineArray(i) % uvShaftDir
nacYaw=>turbineArray(i) % nacYaw
numAnnulusSections=>turbineArray(i) % numAnnulusSections
annulusSectionAngle=>turbineArray(i) % annulusSectionAngle
! Assign Pointers turbineModel (j)
ShftTilt=>turbineModel(j) % ShftTilt
preCone=>turbineModel(j) % preCone
TowerHt=>turbineModel(j) % TowerHt
Twr2Shft=> turbineModel(j) % Twr2Shft
OverHang=>turbineModel(j) % OverHang
UndSling=>turbineModel(j) % UndSling
TipRad=>turbineModel(j) % TipRad
HubRad=>turbineModel(j) % HubRad
PreCone=>turbineModel(j) %PreCone
!!-- Do all proper conversions for the required variables
! Convert nacelle yaw from compass directions to the standard convention
!~ call atm_compassToStandard(nacYaw)
! The nacelle Yaw is set to 0 deg in the streamwise direction
write(*,*) 'NacYaw is ', nacYaw
! Turbine specific
azimuth = degRad * azimuth
rotSpeed = rpmRadSec * rotSpeed
nacYaw = degRad * nacYaw
! Turbine model specific
shftTilt = degRad * shftTilt
preCone = degRad * preCone
! Calculate tower shaft intersection and rotor apex locations. (The i-index is
! at the turbine array level for each turbine and the j-index is for each type
! of turbine--if all turbines are the same, j- is always 0.) The rotor apex is
! not yet rotated for initial yaw that is done below.
towerShaftIntersect = turbineArray(i) % baseLocation
towerShaftIntersect(3) = towerShaftIntersect(3) + TowerHt + Twr2Shft
rotorApex = towerShaftIntersect
rotorApex(1) = rotorApex(1) + (OverHang + UndSling) * cos(ShftTilt)
rotorApex(3) = rotorApex(3) + (OverHang + UndSling) * sin(ShftTilt)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Create Nacelle Point
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
turbineArray(i) % nacelleLocation = rotorApex
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Create the first set of actuator points !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Define the vector along the shaft pointing in the direction of the wind
uvShaftDir = OverHang / abs( OverHang )
! Define the vector along the shaft pointing in the direction of the wind
uvShaft = vector_add(rotorApex , - towerShaftIntersect)
uvShaft = vector_divide(uvShaft, vector_mag(uvShaft))
uvShaft = vector_multiply( uvShaft,uvShaftDir)
! Define vector aligned with the tower pointing from the ground to the nacelle
uvTower = vector_add(towerShaftIntersect, - baseLocation)
uvTower = vector_divide( uvTower, vector_mag(uvTower))
! Define thickness of each blade section
do k=1, numBladePoints
db(k) = (TipRad - HubRad)/(numBladePoints)
enddo
! This creates the first set of points
do k=1, numBl
root = rotorApex
beta = PreCone - ShftTilt
root(1)= root(1) + HubRad*sin(beta)
root(3)= root(3) + HubRad*cos(beta)
!~ dist = HubRad
dist = 0.
! Number of blade points for the first annular section
do m=1, numBladePoints
dist = dist + 0.5*(db(m))
bladePoints(k,1,m,1) = root(1) + dist*sin(beta)
bladePoints(k,1,m,2) = root(2)
bladePoints(k,1,m,3) = root(3) + dist*cos(beta)
do n=1,numAnnulusSections
!~ bladeRadius(k,n,m) = dist
bladeRadius(k,n,m) = dist + HubRad
solidity(k,n,m)=1./numAnnulusSections
enddo
dist = dist + 0.5*db(m)
enddo
! If there are more than one blade create the points of other blades by
! rotating the points of the first blade
if (k > 1) then
do m=1, numBladePoints
bladePoints(k,1,m,:)=rotatePoint(bladePoints(k,1,m,:), rotorApex, &
uvShaft,(360.0/NumBl)*(k-1)*degRad)
enddo
endif
! Rotate points for all the annular sections
if (numAnnulusSections .lt. 2) cycle ! Cycle if only one section (ALM)
do n=2, numAnnulusSections
do m=1, numBladePoints
bladePoints(k,n,m,:) = &
rotatePoint(bladePoints(k,1,m,:), rotorApex, &
uvShaft,(annulusSectionAngle/(numAnnulusSections))*(n-1.)*degRad)
enddo
enddo
enddo
! Apply the first rotation
turbineArray(i) % deltaAzimuth = azimuth
call atm_rotateBlades(i)
end subroutine atm_create_points
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_update(i, dt, time)
! This subroutine updates the model each time-step
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
integer, intent(in) :: i ! Turbine number
real(rprec), intent(in) :: dt ! Time step
real(rprec), intent(in) :: time ! Simulation time
! Rotate the blades
call atm_computeRotorSpeed(i,dt)
call atm_rotateBlades(i)
call atm_control_yaw(i, time)
!~ if(pastFirstTimeStep) then
! Compute the lift correction for this case
!~ call atm_compute_cl_correction(i)
!~ endif
end subroutine atm_update
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_control_yaw(i, time)
! This subroutine updates the model each time-step
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
integer, intent(in) :: i ! Turbine number
real(rprec), intent(in) :: time ! Simulation time
integer :: j ! Turbine Type ID
! Identifies the turbineModel being used
j=turbineArray(i) % turbineTypeID ! The type of turbine (eg. NREL5MW)
! Will calculate the yaw angle and yaw the nacelle (from degrees to radians)
if ( turbineModel(j) % YawControllerType == "timeYawTable" ) then
turbineArray(i) % deltaNacYaw = interpolate(time, &
turbineModel(j) % yaw_time(:), turbineModel(j) % yaw_angle(:)) * degRad - &
turbineArray(i) % NacYaw
! Yaw only if angle is greater than given tolerance
if (abs(turbineArray(i) % deltaNacYaw) > 0.00000001) then
call atm_yawNacelle(i)
endif
!~ write(*,*) 'Delta Yaw is', turbineArray(i) % deltaNacYaw/degRad
!~ write(*,*) 'Nacelle Yaw is', turbineArray(i) % NacYaw/degRad
endif
end subroutine atm_control_yaw
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_computeRotorSpeed(i,dt)
! This subroutine rotates the turbine blades
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
integer, intent(in) :: i ! Turbine number
real(rprec), intent(in) :: dt ! time step
integer :: j ! Turbine type
!real(rprec) :: deltaAzimuth ! Angle of rotation
! Pointers to turbineArray(i)
real(rprec), pointer :: rotSpeed, torqueGen, torqueRotor, fluidDensity
! Pointers to turbineModel(j)
real(rprec), pointer :: GBRatio, CutInGenSpeed, RatedGenSpeed
real(rprec), pointer :: Region2StartGenSpeed, Region2EndGenSpeed
real(rprec), pointer :: CutInGenTorque,RateLimitGenTorque,RatedGenTorque
real(rprec), pointer :: KGen,TorqueControllerRelax, DriveTrainIner
! Other variables to be used
real(rprec) :: torqueGenOld, genSpeed, dGenSpeed, Region2StartGenTorque
real(rprec) :: torqueSlope, Region2EndGenTorque
! Pitch Controller values
real(rprec) :: GK, KP, KI, SpeedError
real(rprec), pointer :: IntSpeedError, PitchControlAngle
j=turbineArray(i) % turbineTypeID
rotSpeed=>turbineArray(i) % rotSpeed
torqueGen=>turbineArray(i) % torqueGen
torqueRotor => turbineArray(i) % torqueRotor
fluidDensity => turbineArray(i) % fluidDensity
GBRatio => turbineModel(j) % GBRatio
CutInGenSpeed => turbineModel(j) % CutInGenSpeed
CutInGenTorque => turbineModel(j) % CutInGenTorque
Region2StartGenSpeed => turbineModel(j) % Region2StartGenSpeed
KGen => turbineModel(j) % KGen
RatedGenSpeed => turbineModel(j) % RatedGenSpeed
Region2EndGenSpeed => turbineModel(j) % Region2EndGenSpeed
RatedGenTorque => turbineModel(j) % RatedGenTorque
RateLimitGenTorque => turbineModel(j) % RateLimitGenTorque
TorqueControllerRelax => turbineModel(j) % TorqueControllerRelax
DriveTrainIner => turbineModel(j) % DriveTrainIner
IntSpeedError => turbineArray(i) % IntSpeedError
PitchControlAngle => turbineArray(i) % PitchControlAngle
! No torque controller option
if (turbineModel(j) % TorqueControllerType == "none") then
elseif (turbineModel(j) % TorqueControllerType == "fiveRegion") then
! Get the generator speed.
genSpeed = (rotSpeed/rpmRadSec)*GBRatio
! Save the generator torque from the last time step.
torqueGenOld = torqueGen
! Region 1.
if (genSpeed < CutInGenSpeed) then
torqueGen = CutInGenTorque
! Region 1-1/2.
elseif ((genSpeed >= CutInGenSpeed) .and. &
(genSpeed < Region2StartGenSpeed)) then
dGenSpeed = genSpeed - CutInGenSpeed
Region2StartGenTorque = KGen * Region2StartGenSpeed * &
Region2StartGenSpeed
torqueSlope = (Region2StartGenTorque - CutInGenTorque) / &
( Region2StartGenSpeed - CutInGenSpeed )
torqueGen = CutInGenTorque + torqueSlope*dGenSpeed
! Region 2.
elseif ((genSpeed >= Region2StartGenSpeed) .and. &
(genSpeed < Region2EndGenSpeed)) then
torqueGen = KGen * genSpeed * genSpeed
! Region 2-1/2.
elseif ((genSpeed >= Region2EndGenSpeed) .and. &
(genSpeed < RatedGenSpeed)) then
dGenSpeed = genSpeed - Region2EndGenSpeed
Region2EndGenTorque = KGen * Region2EndGenSpeed * &
Region2EndGenSpeed
torqueSlope = (RatedGenTorque - Region2EndGenTorque) / &
( RatedGenSpeed - Region2EndGenSpeed )
torqueGen = Region2EndGenTorque + torqueSlope*dGenSpeed
! Region 3.
elseif (genSpeed >= RatedGenSpeed) then
torqueGen = RatedGenTorque
endif
! Limit the change in generator torque if after first time step
! (otherwise it slowly ramps up from its zero initialized value--we
! want it to instantly be at its desired value on the first time
! step, but smoothly vary from there).
if ((abs((torqueGen - torqueGenOld)/dt) > RateLimitGenTorque) &
.and. (pastFirstTimeStep)) then
if (torqueGen > torqueGenOld) then
torqueGen = torqueGenOld + (RateLimitGenTorque * dt);
elseif (torqueGen <= torqueGenOld) then
torqueGen = torqueGenOld - (RateLimitGenTorque * dt);
endif
endif
! Update the rotor speed.
rotSpeed = rotSpeed + TorqueControllerRelax * (dt/DriveTrainIner) * &
(torqueRotor*fluidDensity - GBRatio*torqueGen)
if (turbineModel(j) % PitchControllerType == "none") then
! Limit the rotor speed to be positive and such that the generator
!does not turn faster than rated.
rotSpeed = max(0.0,rotSpeed)
rotSpeed = min(rotSpeed,(RatedGenSpeed*rpmRadSec)/GBRatio)
endif
! Torque control for fixed tip speed ratio
! Note that this current method does NOT support Coning in the rotor
elseif (turbineModel(j) % TorqueControllerType == "fixedTSR") then
if (pastFirstTimeStep) then
! Integrate the velocity along all actuator points
call atm_integrate_u(i)
! Match the rotor speed to a given TSR
rotSpeed = turbineArray(i) % u_infinity_mean * &
turbineArray(i) % TSR / turbineModel(j) % tipRad
! Important to get rid of negative values
rotSpeed = max(0.0,rotSpeed)
endif
endif
! Pitch controllers (If there's no pitch controller, then don't do anything)
if (turbineModel(j) % PitchControllerType == "gainScheduledPI") then
! Get the generator speed.
genSpeed = (rotSpeed/rpmRadSec)*GBRatio
! Calculate the gain
GK = 1.0/(1.0 + PitchControlAngle/turbineModel(j) % PitchControlAngleK)
! Calculate the Proportional and Integral terms
KP = GK*turbineModel(j) % PitchControlKP0
KI = GK*turbineModel(j) % PitchControlKI0
! Get speed error (generator in rpm) and update integral
! Integral is saturated to not push the angle beyond its limits
SpeedError = genSpeed - RatedGenSpeed
!write(*,*) 'Speed Error is: ', speedError
IntSpeedError = IntSpeedError + SpeedError*dt
IntSpeedError = min( max(IntSpeedError, &
turbineModel(j) % PitchControlAngleMin/KI), &
turbineModel(j) % PitchControlAngleMax/KI)
! Apply PI controller and saturate
PitchControlAngle = KP*SpeedError + KI*IntSpeedError
PitchControlAngle = min( max( PitchControlAngle, &
turbinemodel(j) % PitchControlAngleMin), &
turbineModel(j) % PitchControlAngleMax)
endif
! Compute the change in blade position at new rotor speed.
turbineArray(i) % deltaAzimuth = rotSpeed * dt
end subroutine atm_computeRotorSpeed
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_rotateBlades(i)
! This subroutine rotates the turbine blades
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
implicit none
integer, intent(in) :: i ! Turbine number
!real(rprec), intent(in) :: dt ! time step
integer :: j ! Turbine type
integer :: m, n, q ! Counters tu be used in do loops
real(rprec) :: deltaAzimuth, deltaAzimuthI ! Angle of rotation
real(rprec), pointer :: rotorApex(:)
real(rprec), pointer :: rotSpeed
real(rprec), pointer :: uvShaft(:)
real(rprec), pointer :: azimuth
j=turbineArray(i) % turbineTypeID
! Variables which are used by pointers
rotorApex=> turbineArray(i) % rotorApex
rotSpeed=>turbineArray(i) % rotSpeed
uvShaft=>turbineArray(i) % uvShaft
azimuth=>turbineArray(i) % azimuth
! Angle of rotation
deltaAzimuth = turbineArray(i) % deltaAzimuth
! Check the rotation direction first and set the local delta azimuth
! variable accordingly.
if (turbineArray(i) % rotationDir == "cw") then
deltaAzimuthI = deltaAzimuth
else if (turbineArray(i) % rotationDir == "ccw") then
deltaAzimuthI = -deltaAzimuth
endif
! Loop through all the points and rotate them accordingly
do q=1, turbineArray(i) % numBladePoints
do n=1, turbineArray(i) % numAnnulusSections
do m=1, turbineModel(j) % numBl
turbineArray(i) % bladePoints(m,n,q,:)=rotatePoint( &
turbineArray(i) % bladePoints(m,n,q,:), rotorApex, uvShaft, &
deltaAzimuthI)
enddo
enddo
enddo
if (pastFirstTimeStep) then
azimuth = azimuth + deltaAzimuth;
if (azimuth .ge. 2.0 * pi) then
azimuth =azimuth - 2.0 *pi;
endif
endif
end subroutine atm_rotateBlades
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine atm_compute_cl_correction(i)
! This subroutine compute the correction for the Cl
! It needs to be run only once in the initialization
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
integer, intent(in) :: i ! Turbine number
integer :: j ! Turbine type
integer :: m, n, q, k ! Counters tu be used in do loops
real(rprec) :: z_k, z_q ! The radial distance from the tip
real(rprec) :: eps_opt ! The optimal epsilon
real(rprec) :: eps_s ! The epsilon value in the simulation
real(rprec) :: dG ! The finite difference G
real(rprec) :: uy_LES, uy_opt ! Perturbation velocity
real(rprec) :: f ! The relaxation factor
real(rprec) :: cd_eff ! The effective cd coefficient
real(rprec) :: dir_vec_tan(3) ! The direction vector tangential
real(rprec) :: dir_vec_s(3) ! The direction vector streamwise
real(rprec) :: vec_kq(3) ! The vector between points k and q
real(rprec), pointer :: uy_opt_vec(:,:,:,:)
real(rprec), pointer :: uy_LES_vec(:,:,:,:)
real(rprec), pointer :: ux_LES_vec(:,:,:,:)
real(rprec), pointer :: Uinf_vec(:,:,:,:)
real(rprec), pointer :: windVectors(:,:,:,:)
real(rprec), pointer :: cl(:,:,:)
real(rprec), pointer :: cd(:,:,:)
real(rprec), pointer :: Vmag(:,:,:)
real(rprec), pointer :: chord(:,:,:)
! The wind vector
uy_opt_vec => turbineArray(i) % uy_opt_vec
uy_LES_vec => turbineArray(i) % uy_LES_vec
Uinf_vec => turbineArray(i) % Uinf_vec
windVectors => turbineArray(i) % windVectors
cl => turbineArray(i) % cl
cd => turbineArray(i) % cd
Vmag => turbineArray(i) % Vmag
chord => turbineArray(i) % chord
ux_LES_vec => turbineArray(i) % ux_LES_vec
! Simulation epsilon
eps_s = turbineArray(i) % epsilon
! The turbine type number
j=turbineArray(i) % turbineTypeID
! The relaxation factor (between zero and one)
f = 0.1
! Determine the inflow velocity U infinity relative to the blade
do q=1, turbineArray(i) % numBladePoints
do n=1, turbineArray(i) % numAnnulusSections
do m=1, turbineModel(j) % numBl
! Compute Uinf from the previous time by subtracting the
! induced velocity from the optimal solution and the
! drag from the used epsilon
Uinf_vec(m,n,q,:) = windVectors(m,n,q,:) - uy_opt_vec(m,n,q,:)
! Determine the drag vector. This only acts locally on the
! actuator point
dir_vec_s = vector_divide(Uinf_vec(m,n,q,:), &
vector_mag(Uinf_vec(m,n,q,:)))
! The effective drag coefficient
cd_eff = cd(m,n,q) + cl(m,n,q) * vector_mag(uy_opt_vec(m,n,q,:)) &
/ Vmag(m,n,q)
! Compute the induced velocity
ux_LES_vec(m,n,q,:) = cd_eff * chord(m,n,q) / eps_s * 1. &
/ (4. * pi**0.5) * Vmag(m,n,q) * dir_vec_s(:)
! Add the induced velocity from the epsilon LES
!~ Uinf_vec(m,n,q,:) = Uinf_vec(m,n,q,:) + ux_LES_vec(m,n,q,:)
! Save the wind vectors as U inf to compute the forces in the
! perpendicular and parallel directions
!~ windVectors(m,n,q,:) = Uinf_vec(m,n,q,:)
!~ Vmag(m,n,q) = vector_mag(Uinf_vec(m,n,q,:))
enddo
enddo
enddo
! First compute the function G from the previous time
do q=1, turbineArray(i) % numBladePoints
do n=1, turbineArray(i) % numAnnulusSections
do m=1, turbineModel(j) % numBl
! Compute G
turbineArray(i) % G(m,n,q) = 1./2. * &
turbineArray(i) % Cl(m,n,q) * &
turbineArray(i) % chord(m,n,q) * &
turbineArray(i) % Vmag(m,n,q)**2
! The optimal epsilon
turbineArray(i) % epsilon_opt(m,n,q) = &
turbineArray(i) % chord(m,n,q) * &
turbineArray(i) % optimalEpsilonChord
enddo
enddo
enddo
! Now compute the difference dG
do m=1, turbineModel(j) % numBl
do n=1, turbineArray(i) % numAnnulusSections
do q=1, turbineArray(i) % numBladePoints
! Notice this is the G form the previous time, not Gb
! Compute the difference in G using finite difference
! The tip and root are a special case
if (q.eq.1) then
turbineArray(i) % dG(m,n,q) = turbineArray(i) % G(m,n,1)
elseif (q.eq.turbineArray(i) % numBladePoints) then
turbineArray(i) % dG(m,n,q) = &
- turbineArray(i) % G(m,n,turbineArray(i) % numBladePoints)
else
! The central finite difference in G
turbineArray(i) % dG(m,n,q) = (turbineArray(i) % G(m,n,q+1) - &
turbineArray(i) % G(m,n,q-1)) / 2.
!~ turbineArray(i) % dG(m,n,q) = (turbineArray(i) % G(m,n,q+1) - &
!~ turbineArray(i) % G(m,n,q))
endif
enddo
enddo
enddo
! Now compute the difference dG
do m=1, turbineModel(j) % numBl
do n=1, turbineArray(i) % numAnnulusSections
! Now compute the induced velocity at every point
do q=1, turbineArray(i) % numBladePoints
turbineArray(i) % uy_LES_vec(m,n,q,1:3) = 0.
turbineArray(i) % uy_opt_vec(m,n,q,1:3) = 0.
! The z location at the point to compute the correction
z_q = turbineArray(i) % bladeRadius(m,n,q)
! The perturbation velocity
uy_LES = 0.
uy_opt = 0.
! Store the optimal epsilon in this variable
eps_opt = turbineArray(i) % epsilon_opt(m,n,q)
! Loop through all point and compute the perturbation velocity
do k=1, turbineArray(i) % numBladePoints
! Only evaluate if not at the actuator point
! At the actuator point the function is zero
if (k /= q) then
! The change in circulation
dG = turbineArray(i) % dG(m,n,k)
! The z location at the point to compute the correction
z_k = turbineArray(i) % bladeRadius(m,n,k)
! The vector in the z direction
vec_kq(:) = (/ 0._rprec, 0._rprec, 1._rprec /)
! Compute the vector which is perpendicular to Uinf
! This vector is used to project the induced velocity
dir_vec_tan = cross_product(Uinf_vec(m,n,k,:), vec_kq(:))
dir_vec_tan = vector_divide(dir_vec_tan , &
vector_mag(dir_vec_tan))
! Perturbation velocity from the previous time