-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathFAST_Lin.f90
2949 lines (2238 loc) · 177 KB
/
FAST_Lin.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
MODULE FAST_Lin_Subs
USE NWTC_Library
CONTAINS
!=======================================================================
SUBROUTINE CalcSteady
! CalcSteady is used to march in time until a periodic steady state
! solution (if rotor is spinning) or static equilibrium (if rotor
! is parked) is found via a simple tolerance check. Once found
! the periodic solution is interpolated to find the periodic
! operating point displacements, velocities, and accelerations of
! each DOF at each of the equally-spaced azimuth steps.
! FAST Modules:
USE Blades
USE DOFs
USE DriveTrain
USE Features
USE InitCond
USE Linear
USE Modes
USE NacelleYaw
USE SimCont
USE Tower
USE TurbConf
USE TurbCont
! AeroDyn Modules:
USE AeroDyn
USE FASTsubs !Solver
IMPLICIT NONE
! Local variables.
REAL(ReKi) :: AvgSpeed ! Average rotor speed over one Period.
REAL(ReKi), ALLOCATABLE :: CBEStart (:,:,:) ! The value of CBE upon entering this routine.
REAL(ReKi), ALLOCATABLE :: CBFStart (:,:,:) ! The value of CBF upon entering this routine.
REAL(ReKi) :: CTFAStart (2,2) ! The value of CTFA upon entering this routine.
REAL(ReKi) :: CTSSStart (2,2) ! The value of CTSS upon entering this routine.
REAL(ReKi) :: DampFact ! Damping factor that scales with the velocity 2-norm for artificially increasing the damping to aid in convergence.
REAL(ReKi) :: DTTorDmpStart ! The value of DTTorDmp upon entering this routine.
REAL(ReKi) :: Fract ! Ratio (fraction) used during an interpolation calculation.
REAL(ReKi) :: HHWndVec (3) ! Current hub-height wind vector in the AeroDyn coordinate system.
REAL(ReKi) :: HHWndVecS (3) ! Hub-height wind vector in the AeroDyn coordinate system at the start of the period.
REAL(ReKi) :: KBlPitch ! Gain (in rad) for adjusting the rotor collective blade pitch from the relative rotor speed error (used with TrimCase = 3).
REAL(ReKi) :: KDamp ! Gain for artificially increasing the damping for the blades, tower, and drivetrain to aid in convergence (used during all steady state solution calculations).
REAL(ReKi) :: KGenTrq ! Gain (in N-m) for adjusting the electrical generator torque from the relative rotor speed error (used with TrimCase = 2).
REAL(ReKi) :: KYawPos ! Gain (in rad) for adjusting the nacelle yaw command from the relative rotor speed error (used with TrimCase = 1).
REAL(ReKi) :: MinDif ! The minimum value of Qop4MinusQ4 which is greater than or equal to zero.
REAL(ReKi), ALLOCATABLE :: QDStart (:) ! Velocities saved from the start of the period.
REAL(ReKi), ALLOCATABLE :: QDWeight (:) ! Weights used to scale each individual velocities when computing the velocity 2-norm, AbsQDNorm.
REAL(ReKi), ALLOCATABLE :: QStart (:) ! Displacements saved from the start of the period.
REAL(ReKi), ALLOCATABLE :: QWeight (:) ! Weights used to scale each individual displacements when computing the displacement 2-norm, AbsQNorm.
REAL(ReKi) :: RotAzimDif ! The difference between the RotAzimop's and the current rotor azimuth angles.
REAL(ReKi), ALLOCATABLE :: RotAzimOp (:) ! The rotor azimuth angles that we will linearize about.
REAL(ReKi) :: SpeedErr ! The relative rotor speed error for the current iteration.
INTEGER(4) :: I ! Loops through all DOFs
INTEGER(4) :: K ! Loops through blades.
INTEGER(4) :: L ! Generic index.
INTEGER(4) :: LStart ! The index of the first azimuth step of Qop greater than or equal to the current azimuth step in Q(DOF_DrTr).
INTEGER :: Sttus ! Status returned by an attempted allocation.
CHARACTER(10) :: AbsQDNStr ! String containing the current value of AbsQDNorm.
CHARACTER(10) :: AbsQNStr ! String containing the current value of AbsQNorm.
CHARACTER( 9) :: AvgSpdStr ! String containing the current value of AvgSpeed.
CHARACTER(10) :: GenTrqStr ! String containing the current value of GenTrq.
CHARACTER( 4) :: IterStr ! String containing the current value of Iteration.
CHARACTER(10) :: PitchStr ! String containing the current value of BlPitch(1).
CHARACTER(10) :: YawPosStr ! String containing the current value of the command (demand) yaw angle.
CHARACTER( 7) :: ZTimeStr ! String containing the current value of ZTime.
! Specify the gains used during the trim analysis:
! (THE GAINS SPECIFIED BELOW WERE FOUND BY TRIAL AND ERROR THROUGH THE USE
! OF SEVERAL DIFFERENT WIND TURBINE MODELS. EXPERIENCED USERS OF FAST MAY
! WISH TO PLAY AROUND WITH THESE GAINS IN ORDER TO IMPROVE CONVERGENCE FOR
! THEIR ANALYSIS.)
KDamp = 0.0 ! Gain for artificially increasing the damping for the blades, tower, and drivetrain to aid in convergence (used during all steady state solution calculations). !JASON: This didn't improve the convergence times much and, in fact, made the solution more unstable. Thus, I disabled this by setting the gain to zero.
KYawPos = 0.1 ! Gain (in rad) for adjusting the nacelle yaw command from the relative rotor speed error (used with TrimCase = 3).
KGenTrq = 1.5*AvgNrmTpRd**3/GBRatio ! Gain (in N-m) for adjusting the electrical generator torque from the relative rotor speed error (used with TrimCase = 2). The rated rotor torque of a typical wind turbine in N-m is roughly 15*Radius^3 where the radius is in meters. The 1.5 is 1/10th of the factor of 15. The 1/GBRatio is used to cast the torque onto the HSS (i.e., generator) side of the gearbox. (Nm)
KBlPitch = 0.1 ! Gain (in rad) for adjusting the rotor collective blade pitch from the relative rotor speed error (used with TrimCase = 3).
! If trimming electrical generator torque, let's override the generator
! torque models calculated in SUBROUTINE DrvTrTrq(). The trimmed torque
! value is determined with VSContrl = 9999.
! Also, let's specify a starting "guess" value for the generator torque:
IF ( ( GenDOF ) .AND. ( TrimCase == 2 ) ) THEN ! We will be trimming generator torque
VSContrl = 9999 ! Override the selected VSContrl and GenModel input options
GenTrq = 0.0 ! Use a defult generator torque of zero--I found that if I guessed high on the initial generator torque, the results were very unstable; but if I started at zero, the solution always converged.
! NOTE: This method of specifying the default generator torque using the intial
! aerodynamic torque (converted to the generator end of the drivetrain),
! did not work:
! ! This default generator torque is determined by calling RtHS once (with
! ! initial conditions) and estimating the generator torque from the
! ! aerodynamic torque given in vector, MomLPRott.
! ! NOTE: Using this method requires the USE of MODULE RtHndSid() and
! ! declaration of global function, DotProd:
! GenTrq = 0.0 ! Define a dummy generator torque
! QT = Q (:,IC(1)) ! Use initial conditions
! QDT = QD(:,IC(1)) ! for the DOFs
! CALL RtHS ! Call the dynamics routine once
! GenTrq = DotProd( MomLPRott, e1 )*GBoxEffFac/GBRatio ! Convert the aerodynamic torque to the HSS-side
ENDIF
! Allocate some arrays.
ALLOCATE ( CBEStart(NumBl,1,1) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the CBEStart array.' )
ENDIF
ALLOCATE ( CBFStart(NumBl,2,2) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the CBFStart array.' )
ENDIF
ALLOCATE ( QWeight(NDOF) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the QWeight array.' )
ENDIF
ALLOCATE ( QDWeight(NDOF) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the QDWeight array.' )
ENDIF
ALLOCATE ( QStart(NDOF) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the QStart array.' )
ENDIF
ALLOCATE ( QDStart(NDOF) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the QDStart array.' )
ENDIF
ALLOCATE ( RotAzimop( NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the RotAzimop array.' )
ENDIF
! Initialize the displacement weightings to unity:
QWeight = 1.0
! Use weights that nondimensionalize the displacement DOFs for the flexible
! members (i.e., blade and tower flexibility). These weights are chosen
! so that the 2-norms are computed using the slopes (dimensionless angles)
! of the tips of the flexible members, instead of the actual dimensional
! translational displacements of the flexible members. This ensures that
! each DOF in the 2-norm calculation has the same units (radians).
QWeight(DOF_TFA1) = SHP( 1.0, TwrFlexL, TwFAM1Sh(:), 1 ) ! The slope of the 1st tower fore-aft mode shape at the tower-top.
QWeight(DOF_TSS1) = SHP( 1.0, TwrFlexL, TwSSM1Sh(:), 1 ) ! The slope of the 1st tower side-to-side mode shape at the tower-top.
QWeight(DOF_TFA2) = SHP( 1.0, TwrFlexL, TwFAM2Sh(:), 1 ) ! The slope of the 2nd tower fore-aft mode shape at the tower-top.
QWeight(DOF_TSS2) = SHP( 1.0, TwrFlexL, TwSSM2Sh(:), 1 ) ! The slope of the 2nd tower side-to-side mode shape at the tower-top.
DO K = 1,NumBl ! Loop through all blades
QWeight(DOF_BF(K,1)) = SHP( 1.0, BldFlexL, BldFl1Sh(:,K), 1 ) ! The slope of the 1st blade flapwise mode shape at the blade tip.
QWeight(DOF_BF(K,2)) = SHP( 1.0, BldFlexL, BldFl2Sh(:,K), 1 ) ! The slope of the 2nd blade flapwise mode shape at the blade tip.
QWeight(DOF_BE(K,1)) = SHP( 1.0, BldFlexL, BldEdgSh(:,K), 1 ) ! The slope of the 1st blade edgewise mode shape at the blade tip.
ENDDO ! K - Blades
! Like the flexible members, use weights that nondimensionalize the
! translational displacement DOFs of the platform. These weights are
! chosen so that the 2-norms are computed using the slopes (dimensionless
! angles) at a reference depth between the displaced and undisplaced
! platform position, instead of the actual dimensional translational
! displacements of the platform. This ensures that each DOF in the 2-norm
! calculation has the same units (radians).
QWeight(DOF_Sg ) = 0.01 ! Assume 100m
QWeight(DOF_Sw ) = 0.01 ! reference
QWeight(DOF_Hv ) = 0.01 ! depth
! Use the same weightings for velocity. The resulting units of the velocity
! 2-norm will be rad/sec
QDWeight = QWeight
!JASON: WE SHOULD ADD ARTIFICIAL DAMPING IN OTHER DOFs, LIKE THE FURLING, TEETER, AND PLATFORM MOTIONS.
! Let's store the original damping values of the blades, tower, and
! drivetrain:
DTTorDmpStart = DTTorDmp
DO I = 1,2 ! Loop through all tower DOFs in one direction
DO L = 1,2 ! Loop through all tower DOFs in one direction
CTFAStart( I,L) = CTFA( I,L)
CTSSStart( I,L) = CTSS( I,L)
ENDDO ! L - All tower DOFs in one direction
ENDDO ! I - All tower DOFs in one direction
DO K = 1,NumBl ! Loop through all blades
DO I = 1,2 ! Loop through all flap DOFs
DO L = 1,2 ! Loop through all flap DOFs
CBFStart(K,I,L) = CBF(K,I,L)
ENDDO ! L - All flap DOFs
ENDDO ! I - All flap DOFs
CBEStart (K,1,1) = CBE(K,1,1)
ENDDO ! K - All blades
! Get the current hub-height wind speed:
HHWndVec(:) = AD_GetUndisturbedWind( ZTime, (/REAL(0.0, ReKi), REAL(0.0, ReKi), FASTHH /), Sttus )
! Begin the search for a steady state solution:
! Inform the users of what we are about to do:
CALL WrScr1( ' Beginning iteration to find a steady state solution of type:' )
IF ( .NOT. GenDOF ) THEN ! Constant speed case
IF ( RotSpeed == 0.0 ) THEN
CALL WrScr ( ' Static equilibrium (RotSpeed = 0.0)' )
ELSE
CALL WrScr ( ' Constant speed (GenDOF = False)' )
ENDIF
ELSEIF ( TrimCase == 1 ) THEN ! Trimming nacelle yaw
CALL WrScr ( ' Trimmed nacelle yaw (TrimCase = 1)' )
ELSEIF ( TrimCase == 2 ) THEN ! Trimming generator torque
CALL WrScr ( ' Trimmed generator torque (TrimCase = 2)' )
ELSEIF ( TrimCase == 3 ) THEN ! Trimming collective blade pitch
CALL WrScr ( ' Trimmed collective blade pitch (TrimCase = 3)' )
ENDIF
CALL WrScr1( ' Avg Rotor Nacelle Generator Blade Dsplcmnt Velocity' )
CALL WrScr ( ' Iter Time Speed Yaw Demand Torque Pitch 2-norm 2-norm' )
CALL WrScr ( ' Nmbr (sec) (rpm) (deg) (kN-m) (deg) (rad) (rad/s)' )
CALL WrScr ( ' -----------------------------------------------------------------------------' )
DO
! Increment the Iteration number:
Iteration = Iteration + 1
! Reinitialize the average rotor speed for this Period:
AvgSpeed = 0.0
! Save the current values of the DOFs and their 1st time derivatives:
QStart = Q (:,IC(1))
QDStart = QD(:,IC(1))
! If we will have reached TMax during this iteration or when calculating
! Qop, QDop, or QD2op later, quit trying to find a steady state solution:
IF ( ( Step + 2*NStep )*DT > TMax ) THEN
CALL WrScr1( ' The solution does not appear to converge after '//TRIM( Flt2LStr( TMax ) )//' seconds!' )
CALL WrScr1( ' Try increasing the total run time, TMax, increasing system damping values,'// &
' or increasing the convergence tolerances, DispTol and/or VelTol.' )
CALL ProgAbort ( ' The linearized system matrices were not formed.' )
ENDIF
! Loop through NStep time steps (one complete period for periodic systems):
DO L = 1,NStep
! Save the current hub-height wind velocity vector:
HHWndVecS = HHWndVec
! Call predictor-corrector routine:
CALL Solver
! Make sure the rotor azimuth is not greater or equal to 360 degrees:
IF ( ( Q(DOF_GeAz,IC(1)) + Q(DOF_DrTr,IC(1)) ) >= TwoPi ) THEN
Q(DOF_GeAz,IC(1)) = Q(DOF_GeAz,IC(1)) - TwoPi
ENDIF
! Advance time:
Step = Step + 1
ZTime = Step*DT
! NOTE: There is no sense CALLing SUBROUTINES CalcOuts(), PredictNoise(),
! WrOutput(), and SimStatus() here as in SUBROUTINE TimeMarch().
! Make sure the wind hasn't changed. If so, Abort since we can't find a
! periodic steady state solution with time varying winds.
HHWndVec(:) = AD_GetUndisturbedWind( ZTime, (/REAL(0.0, ReKi), REAL(0.0, ReKi), FASTHH /), Sttus )
IF ( ( HHWndVecS(1) /= HHWndVec(1) ) .OR. &
( HHWndVecS(2) /= HHWndVec(2) ) .OR. &
( HHWndVecS(3) /= HHWndVec(3) ) ) THEN
CALL WrScr1( ' Time varying winds discovered!' )
CALL WrScr ( ' A steady state solution can''t be found using time varying winds.' )
CALL ProgAbort ( ' Use a steady hub-height wind input instead.' )
ENDIF
! Compute the average rotor speed for this Period:
AvgSpeed = AvgSpeed + ( QD(DOF_GeAz,IC(1)) + QD(DOF_DrTr,IC(1)) )/NStep
ENDDO
! Find the absolute displacement and velocity 2-norms:
AbsQNorm = 0.0
AbsQDNorm = 0.0
DO I = 1,NDOF ! Loop through all DOFs
AbsQNorm = AbsQNorm + ( ( Q (I,IC(1)) - QStart (I) )*QWeight (I) )**2
AbsQDNorm = AbsQDNorm + ( ( QD(I,IC(1)) - QDStart(I) )*QDWeight(I) )**2
ENDDO ! I - All DOFs
AbsQNorm = SQRT( AbsQNorm )
AbsQDNorm = SQRT( AbsQDNorm )
! Inform the users of the current status:
WRITE(IterStr ,'(I4)' ) Iteration
WRITE(ZTimeStr ,'(F7.2)' ) ZTime
WRITE(AvgSpdStr ,'(F9.5)' ) AvgSpeed *RPS2RPM
IF ( YawDOF ) THEN ! Yaw DOF is currently enabled (using FAST's built-in actuator) - echo out the current command yaw angle, YawNeut.
WRITE(YawPosStr,'(F10.5)') YawNeut *R2D
ELSE ! Yaw DOF is currently disabled (no built-in actuator) - echo out the current yaw angle, Q(DOF_Yaw,IC(1)).
WRITE(YawPosStr,'(F10.5)') Q(DOF_Yaw,IC(1))*R2D
ENDIF
WRITE(GenTrqStr ,'(F10.5)') 0.001*GenTrq
WRITE(PitchStr ,'(F10.5)') BlPitch(1) *R2D ! only output the pitch of blade 1
WRITE(AbsQNStr ,'(F10.7)') AbsQNorm
WRITE(AbsQDNStr ,'(F10.7)') AbsQDNorm
CALL WrScr ( ' '//IterStr//' '//ZTimeStr//' '//AvgSpdStr//' '//YawPosStr// &
' '//GenTrqStr//' '//PitchStr//' '//AbsQNStr//' '//AbsQDNStr )
! If displacement and velocity convergence tolerances have been met, quit
! iterating (EXIT this DO...LOOP), else continue on...:
IF ( ( AbsQNorm <= DispTol ) .AND. ( AbsQDNorm <= VelTol ) ) EXIT
! Let's add artificial damping to the blade modes, tower modes, and
! torsional compliance of the drivetrain to aid in convergence. Let's
! have this artificial damping scale with the velocity 2-norm (because
! we need more damping the more velocity error we have) and the
! corresponding stiffnesses. Don't add damping to the yaw, furl, or
! teeter hinges, or to the variable speed generator (like SymDyn):
DampFact = KDamp*AbsQDNorm ! Damping factor that scales with the velocity 2-norm, which is used in the following equations
DTTorDmp = DTTorDmpStart + DampFact*DTTorSpr
DO I = 1,2 ! Loop through all tower DOFs in one direction
DO L = 1,2 ! Loop through all tower DOFs in one direction
CTFA( I,L) = CTFAStart( I,L) + DampFact*KTFA( I,L)
CTSS( I,L) = CTSSStart( I,L) + DampFact*KTSS( I,L)
ENDDO ! L - All tower DOFs in one direction
ENDDO ! I - All tower DOFs in one direction
DO K = 1,NumBl ! Loop through all blades
DO I = 1,2 ! Loop through all flap DOFs
DO L = 1,2 ! Loop through all flap DOFs
CBF(K,I,L) = CBFStart(K,I,L) + DampFact*KBF(K,I,L)
ENDDO ! L - All flap DOFs
ENDDO ! I - All flap DOFs
CBE (K,1,1) = CBEStart(K,1,1) + DampFact*KBE(K,1,1)
ENDDO ! K - All blades
! Trim the selected control input only if necessary (i.e., only when GenDOF
! is True):
IF ( GenDOF ) THEN
! Compute the relative rotor speed error:
SpeedErr = ( AvgSpeed - RotSpeed ) / RotSpeed ! AvgSpeed is the average speed over the current iteration; RotSpeed is the desired speed over one period (which can't be zero since GenDOF is True)
! Trim the selected control input using feedback from the relative rotor
! error:
SELECT CASE ( TrimCase ) ! Which control input are we trimming?
CASE ( 1 ) ! Nacelle yaw
! If the yaw DOF is enabled, the trimmed command yaw angle becomes the
! neutral yaw angle in FAST's built-in second-order actuator model defined
! by inputs YawSpr and YawDamp. If the yaw DOF is disabled (no yaw DOF),
! then the trimmed command yaw angle becomes the actual yaw angle (no
! built-in actuator). The commanded rates are always left zero during
! trim since the trim solution must be steady. Also, the SIGN of the yaw
! control gain is determined by the sign of the commanded yaw value since
! the rotor torque versus yaw angle curve is virtually symmetric about yaw
! equals zero:
IF ( YawDOF ) THEN ! Yaw DOF is currently enabled (use FAST's built-in actuator).
YawNeut = YawNeut + SIGN( KYawPos, YawNeut )*SpeedErr
! Leave the neutral yaw rate, YawRateNeut, zero here.
ELSE ! Yaw DOF is currently disabled (no built-in actuator).
Q(DOF_Yaw,IC(1)) = Q(DOF_Yaw,IC(1)) + SIGN( KYawPos, Q(DOF_Yaw,IC(1)) )*SpeedErr ! Update the current saved value used in routine Solver()
! Leave the current yaw rate, QD(DOF_Yaw,IC(1)), zero here.
ENDIF
CASE ( 2 ) ! Electrical generator torque
GenTrq = GenTrq + KGenTrq *SpeedErr
CASE ( 3 ) ! Rotor collective blade pitch
BlPitch = BlPitch + KBlPitch*SpeedErr
ENDSELECT
ENDIF
ENDDO
! We're done!
! Inform the users of this great news!
CALL WrScr1( ' Steady state solution found in '//TRIM(Int2LStr(Iteration))//' iterations!' )
! Let's remove the artificial damping in the blade, tower, and drivetrain by
! restoring the original values:
DTTorDmp = DTTorDmpStart
DO I = 1,2 ! Loop through all tower DOFs in one direction
DO L = 1,2 ! Loop through all tower DOFs in one direction
CTFA( I,L) = CTFAStart( I,L)
CTSS( I,L) = CTSSStart( I,L)
ENDDO ! L - All tower DOFs in one direction
ENDDO ! I - All tower DOFs in one direction
DO K = 1,NumBl ! Loop through all blades
DO I = 1,2 ! Loop through all flap DOFs
DO L = 1,2 ! Loop through all flap DOFs
CBF(K,I,L) = CBFStart(K,I,L)
ENDDO ! L - All flap DOFs
ENDDO ! I - All flap DOFs
CBE (K,1,1) = CBEStart(K,1,1)
ENDDO ! K - All blades
! Now we have to calculate the operating values of the displacement (Qop),
! velocities (QDop), and accelerations (QD2op) at each azimuth step
! (determined by NAzimStep):
! Save the periodic steady state operating point values:
IF ( RotSpeed == 0.0 ) THEN ! Rotor is parked, therefore save the current states.
Qop (:,1) = Q (:,IC(1))
QDop (:,1) = QD (:,IC(1))
QD2op(:,1) = QD2(:,IC(1))
ELSE ! Rotor is spinning, therefore save the states every NAzimStep equally-spaced azimuth steps.
! NOTE: all of the code in this section assumes that the rotor azimuth
! orientation will monotonically increase throughout one period. If
! this condition is not met, the rotor must be oscillating so severely
! that it moves forward then backward slightly, then forward then backward
! slightly, etc... (i.e., the whole 2 step forward, 1 step backward dance
! move). If this is so, it is impossible to find unique periodic
! operating states of each DOF since there will are multiple operating
! states at each azimuth step. To ensure the monotonically increasing
! condition, the zero to TwoPi constraint is not enforced during this
! algorithm and the rotor azimth is checked at each time step to make sure
! it has monotonically increased. To ensure the monotonically increasing
! condition, both the rotor speed and shaft stiffness should be set to
! reasonable values.
! Inform the users of what we are about to do:
CALL WrScr1( ' Interpolating to find the operating point values of each DOF.' )
! Store the rotor azimuth locations, which define the periodic operating
! points. The first azimuth location is always the initial azimuth
! orientation (QAzimInit). Make sure the azimuth angles are between 2*Pi
! (inclusive) and 4*Pi (exclusive):
DO L = 1,NAzimStep ! Loop through all equally-spaced azimuth steps
RotAzimop(L) = QAzimInit + ( TwoPi/NAzimStep )*( L - 1 ) + TwoPi ! 2*Pi <= RotAzimop(L) < 4*Pi (for L = 1,2,...,NAzimStep)
IF ( RotAzimop(L) >= 2.0*TwoPi ) RotAzimop(L) = RotAzimop(L) - TwoPi !
ENDDO ! L - Equally-spaced azimuth steps
! Find the index, LStart, of the 1st azimuth step of RotAzimop greater
! than or equal to the current rotor azimuth location, which is
! ( Q(DOF_GeAz,IC(1)) + Q(DOF_DrTr,IC(1)) ). At the same time, make
! RotAzimop increase monotonically from its lowest value at LStart to its
! highest value at ( LStart - 1 ):
LStart = 1
MinDif = TwoPi
DO L = 1,NAzimStep ! Loop through all equally-spaced azimuth steps
RotAzimDif = RotAzimop(L) - ( Q(DOF_GeAz,IC(1)) + Q(DOF_DrTr,IC(1)) ) ! 0 < RotAzimDif < 4*Pi
IF ( RotAzimDif >= TwoPi ) THEN
RotAzimop(L) = RotAzimop(L) - TwoPi ! current rotor azimuth angle <= RotAzimop(L) < 4*Pi (L = 1,2,...,NAzimStep)
RotAzimDif = RotAzimDif - TwoPi ! 0 <= RotAzimDif < 2*Pi
ENDIF
IF ( RotAzimDif < MinDif ) THEN
LStart = L
MinDif = RotAzimDif
ENDIF
ENDDO ! L - Equally-spaced azimuth steps
! Initialize L to LStart:
L = LStart
! Loop through one complete period until we've calculated the periodic
! operating point displacements and velocities of each DOF at each of the
! equally-spaced azimuth steps:
DO
! Integrate in time:
DO
! If the condition "previous rotor azimuth <= RotAzimop(L) < current rotor
! azimuth" is met, then interpolate to find the periodic operating point
! displacements and velocities of all DOFs at the current rotor azimuth
! step, RotAzimop(L), then increment L to the next azimuth step (outside
! of this DO...LOOP) [where previous rotor azimuth = ( Q(DOF_GeAz,IC(2))
! + Q(DOF_DrTr,IC(2)) ) and current rotor azimuth = ( Q(DOF_GeAz,IC(1))
! + Q(DOF_DrTr,IC(1)) ) in the above notation]. Else, integrate in time
! until the condition is met:
IF ( ( RotAzimop(L) >= ( Q(DOF_GeAz,IC(2)) + Q(DOF_DrTr,IC(2)) ) ) .AND. &
( RotAzimop(L) < ( Q(DOF_GeAz,IC(1)) + Q(DOF_DrTr,IC(1)) ) ) ) EXIT
! Save the current hub-height wind velocity vector:
HHWndVecS = HHWndVec
! Call predictor-corrector routine:
CALL Solver
! NOTE: we do not want to enforce the condition that the rotor azimuth be
! between zero and TwoPi here as in SUBROUTINE TimeMarch(), since we need
! the rotor azimuth to monotonically increase for this algorithm to work.
! Check the monotonically increasing condition on the rotor azimuth DOF:
IF ( ( Q(DOF_GeAz,IC(1)) + Q(DOF_DrTr,IC(1)) ) <= ( Q(DOF_GeAz,IC(2)) + Q(DOF_DrTr,IC(2)) ) ) THEN
CALL WrScr1( ' The rotor azimuth DOF does not monotically increase in SUBROUTINE CalcSteady()'// &
' when finding the periodic operating points!' )
CALL WrScr1( ' Make sure the rotor speed is fast enough and/or the drivetrain is stiff enough'// &
' to ensure that the rotor azimuth does not double back on itself.' )
CALL ProgAbort ( ' The linearized system matrices were not formed.' )
ENDIF
! Advance time:
Step = Step + 1
ZTime = Step*DT
! NOTE: There is no sense CALLing SUBROUTINES CalcOuts(), PredictNoise(),
! WrOutput(), and SimStatus() here as in SUBROUTINE TimeMarch().
! Make sure the wind hasn't changed. If so, ProgAbort since we can't find a
! periodic steady state solution with time varying winds.
HHWndVec(:) = AD_GetUndisturbedWind( ZTime, (/REAL(0.0, ReKi), REAL(0.0, ReKi), FASTHH /), Sttus )
IF ( ( HHWndVecS(1) /= HHWndVec(1) ) .OR. &
( HHWndVecS(2) /= HHWndVec(2) ) .OR. &
( HHWndVecS(3) /= HHWndVec(3) ) ) THEN
CALL WrScr1( ' Time varying winds discovered!' )
CALL WrScr ( ' A steady state solution can''t be found using time varying winds.' )
CALL ProgAbort ( ' Use a steady hub-height wind input instead.' )
ENDIF
ENDDO
! Stop integrating in time, at least for the moment.
! Here is the interpolation:
Fract = ( RotAzimop(L) - ( Q(DOF_GeAz,IC(2)) + Q(DOF_DrTr,IC(2)) ) ) &
/ ( ( Q(DOF_GeAz,IC(1)) + Q(DOF_DrTr,IC(1)) ) - ( Q(DOF_GeAz,IC(2)) + Q(DOF_DrTr,IC(2)) ) )
DO I = 1,NDOF ! Loop through all DOFs
Qop (I,L) = Q (I,IC(2)) + Fract*( Q (I,IC(1)) - Q (I,IC(2)) )
QDop (I,L) = QD (I,IC(2)) + Fract*( QD (I,IC(1)) - QD (I,IC(2)) )
QD2op(I,L) = QD2(I,IC(2)) + Fract*( QD2(I,IC(1)) - QD2(I,IC(2)) )
ENDDO ! I - All DOFs
! Here is the increment of L:
L = L + 1
IF ( L > NAzimStep ) L = 1
! Here is the check to see if we have looped around all of the
! equally-spaced azimuth steps. If we have, exit this loop since we're
! done:
IF ( L == LStart ) EXIT
ENDDO
! We're done!
! Make sure the rotor azimuths are not greater or equal to 360 degrees:
DO L = 1,NAzimStep ! Loop through all equally-spaced azimuth steps
IF ( ( Qop(DOF_GeAz,L) + Qop(DOF_DrTr,L) ) >= TwoPi ) THEN
Qop(DOF_GeAz,L) = Qop(DOF_GeAz,L) - TwoPi
ENDIF
ENDDO ! L - Equally-spaced azimuth steps
ENDIF
! deallocate arrays
IF ( ALLOCATED( CBEStart ) ) DEALLOCATE( CBEStart )
IF ( ALLOCATED( CBFStart ) ) DEALLOCATE( CBFStart )
IF ( ALLOCATED( QDStart ) ) DEALLOCATE( QDStart )
IF ( ALLOCATED( QDWeight ) ) DEALLOCATE( QDWeight )
IF ( ALLOCATED( QStart ) ) DEALLOCATE( QStart )
IF ( ALLOCATED( QWeight ) ) DEALLOCATE( QWeight )
IF ( ALLOCATED( RotAzimOp ) ) DEALLOCATE( RotAzimOp )
RETURN
END SUBROUTINE CalcSteady
!=======================================================================
SUBROUTINE Linearize
! Linearize is used to perturb each displacement and velocity DOF
! about its periodic operating point at each of equally-spaced
! azimuth steps to compute the linearized state matrices.
! FAST Modules:
USE Blades
USE DOFs
USE DriveTrain
USE Features
USE General
USE InitCond
USE Linear
USE Modes
USE NacelleYaw
USE Output
USE RtHndSid
USE Tower
USE TurbConf
USE TurbCont
USE FASTsubs !RtHS(), CalcOuts()
USE AeroElem, ONLY: ADIntrfaceOptions
! AeroDyn Modules:
USE AeroDyn
IMPLICIT NONE
! Local variables:
REAL(ReKi), ALLOCATABLE :: AMat (:,:,:) ! Periodic state matrix (includes only active DOFs)--i.e. [A].
REAL(ReKi), ALLOCATABLE :: BdMat (:,:,:) ! Periodic input disturbance matrix (includes only active DOFs)--i.e. [Bd].
REAL(ReKi), ALLOCATABLE :: BlPitchop(:) ! Steady operating point blade pitch angles.
REAL(ReKi), ALLOCATABLE :: BMat (:,:,:) ! Periodic input matrix (includes only active DOFs)--i.e. [B].
REAL(ReKi), ALLOCATABLE :: CMat (:,:,:) ! Periodic output matrix (includes only active DOFs)--i.e. [C].
REAL(ReKi), ALLOCATABLE :: DdMat (:,:,:) ! Periodic drct. trans. disturb. matrix (includes only active DOFs)--i.e. [Dd].
REAL(ReKi), ALLOCATABLE :: DMat (:,:,:) ! Periodic direct transmission matrix (includes only active DOFs)--i.e. [D].
REAL(ReKi), ALLOCATABLE :: DampMat (:,:,:) ! Periodic damping matrix (includes only active DOFs).
REAL(ReKi) :: DelBlPtch ! Magnitude of pertubation in blade pitch.
REAL(ReKi) :: DelDELTA ! Magnitude of pertubation in horizontal wind direction.
REAL(ReKi) :: DelGenTq ! Magnitude of pertubation in generator torque.
REAL(ReKi) :: DelHSHR ! Magnitude of pertubation in horizontal shear coefficient.
REAL(ReKi), ALLOCATABLE :: DelQ (:) ! Pertubations of displacements about the periodic steady state operating points.
REAL(ReKi), ALLOCATABLE :: DelQD (:) ! Pertubations of velocities about the periodic steady state operating points.
REAL(ReKi) :: DELTAop ! Steady operating point horizontal wind direction.
REAL(ReKi) :: DelV ! Magnitude of pertubation in horizontal hub height wind speed (same as SymDyn).
REAL(ReKi) :: DelVGUST ! Magnitude of pertubation in horizontal hub height wind gust (same as SymDyn).
REAL(ReKi) :: DelVLINSHR ! Magnitude of pertubation in vertical (linear) shear coefficient (same as SymDyn).
REAL(ReKi) :: DelVSHR ! Magnitude of pertubation in vertical shear coefficient (same as SymDyn).
REAL(ReKi) :: DelVZ ! Magnitude of pertubation in vertical wind speed (same as SymDyn).
REAL(ReKi) :: DelYawPos ! Mangitude of pertubation in nacelle yaw angle.
REAL(ReKi) :: DelYawRate ! Mangitude of pertubation in nacelle yaw rate.
REAL(ReKi), ALLOCATABLE :: FdMat (:,:,:) ! Periodic input disturbance matrix (includes only active DOFs)--i.e. [Fd].
REAL(ReKi), ALLOCATABLE :: FMat (:,:,:) ! Periodic input matrix (includes only active DOFs)--i.e. [F].
REAL(ReKi) :: GenTrqop ! Steady operating point generator torque.
REAL(ReKi) :: HSHRop ! Steady operating point horizontal shear coefficient
REAL(ReKi) :: LinPerturbations(7) ! A vector to send the HH wind speed perturbations to the wind module (replace this method!)
REAL(ReKi), ALLOCATABLE :: MassMat (:,:,:) ! Periodic mass matrix (includes only active DOFs).
REAL(ReKi), ALLOCATABLE :: StffMat (:,:,:) ! Periodic stiffness matrix (includes only active DOFs).
REAL(ReKi) :: VGUSTop ! Steady operating point horizontal hub height wind gust
REAL(ReKi) :: VLINSHRop ! Steady operating point vertical (linear) shear coefficient
REAL(ReKi) :: Vop ! Steady operating point horizontal hub height wind speed
REAL(ReKi) :: VSHRop ! Steady operating point vertical shear coefficient
REAL(ReKi) :: VZop ! Steady operating point vertical wind speed
REAL(ReKi) :: YawPosop ! Steady operating point command yaw angle (position).
REAL(ReKi) :: YawRateop ! Steady operating point command yaw rate.
REAL(ReKi), ALLOCATABLE :: Yop (:,:) ! Periodic steady state operating output measurements.
INTEGER :: ErrStat ! Determines if an error was encountered in the Inflow Wind Module
INTEGER(4) :: I ! Generic index.
INTEGER(4) :: I1 ! Loops through all active (enabled) DOFs (rows).
INTEGER(4) :: I2 ! Loops through all active (enabled) DOFs (cols).
INTEGER(4) :: I3 ! Loops through all active (enabled) DOFs (rows & cols).
!bjj: SAVE these variables (or should they be initialized below?)
INTEGER(4) :: IndxCPtch = 0 ! Column index of the rotor collective blade pitch control input in the [B ] and [D ] matrices.
INTEGER(4) :: IndxDELTA = 0 ! Column index of the horizontal wind direction disturbance in the [Bd] and [Dd] matrices.
INTEGER(4) :: IndxGenTq = 0 ! Column index of the generator torque control input in the [B ] and [D ] matrices.
INTEGER(4) :: IndxHSHR = 0 ! Column index of the horizontal wind shear disturbance in the [Bd] and [Dd] matrices.
INTEGER(4) :: IndxIPtch(3) = 0 ! Column indices of the individual blade pitch control inputs in the [B ] and [D ] matrices.
INTEGER(4) :: IndxV = 0 ! Column index of the horizontal hub-height wind speed disturbance in the [Bd] and [Dd] matrices.
INTEGER(4) :: IndxVLSHR = 0 ! Column index of the vertical wind shear (linear) disturbance in the [Bd] and [Dd] matrices.
INTEGER(4) :: IndxVGUST = 0 ! Column index of the horizontal hub-height wind gust disturbance in the [Bd] and [Dd] matrices.
INTEGER(4) :: IndxVSHR = 0 ! Column index of the vertical wind shear disturbance in the [Bd] and [Dd] matrices.
INTEGER(4) :: IndxVZ = 0 ! Column index of the vertical wind speed disturbance in the [Bd] and [Dd] matrices.
INTEGER(4) :: IndxYawPos = 0 ! Column index of the nacelle yaw (position) control input in the [B ] and [D ] matrices.
INTEGER(4) :: IndxYawRate = 0 ! Column index of the nacelle rate control input in the [B ] and [D ] matrices.
INTEGER(4), PARAMETER :: IterRtHS = 2 ! Number of times to iterate on RtHS() when pertubing DOFs and linearizing FAST. As Karl Stol pointed out when developing SymDyn, the iteration for the induction factor in AeroDyn starts from the previous states of the induction factors for each element. Thus, he has found that he gets better convergence when he calls AeroDyn multiple times for each time step. He currently implements 4 loops when he linearizes his model. However, I found that only 2 was necessary.
INTEGER(4) :: K ! Loops through blades.
INTEGER(4) :: L ! Generic index.
INTEGER(4) :: Sttus ! Status returned by an attempted allocation.
CHARACTER( 8) :: FmtHead = '(//,A,/)' ! Format for outputting headings.
CHARACTER( 3) :: FmtText = '(A)' ! Format for outputting pure text.
CHARACTER(200) :: Frmt ! A string to hold a format specifier.
CHARACTER(200) :: Frmt1 ! A string to hold a format specifier.
CHARACTER(200) :: Frmt2 ! A string to hold a format specifier.
! Allocate some arrays:
ALLOCATE ( DelQ (NDOF) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the DelQ array.' )
ENDIF
ALLOCATE ( DelQD(NDOF) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the DelQD array.' )
ENDIF
ALLOCATE ( AMat(2*NActvDOF,2*NActvDOF,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the AMat array.' )
ENDIF
ALLOCATE ( BMat(2*NActvDOF,NInputs,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the BMat array.' )
ENDIF
ALLOCATE ( BdMat(2*NActvDOF,NDisturbs,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the BdMat array.' )
ENDIF
ALLOCATE ( CMat(NumOuts,2*NActvDOF,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the CMat array.' )
ENDIF
ALLOCATE ( DMat(NumOuts,NInputs,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the DMat array.' )
ENDIF
ALLOCATE ( DdMat(NumOuts,NDisturbs,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the DdMat array.' )
ENDIF
ALLOCATE ( BlPitchop(NumBl) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the BlPitchop array.' )
ENDIF
ALLOCATE ( Yop(NumOuts,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the Yop array.' )
ENDIF
IF ( MdlOrder == 2 ) THEN ! 2nd order model
ALLOCATE ( MassMat(NActvDOF,NActvDOF,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the MassMat array.' )
ENDIF
ALLOCATE ( DampMat(NActvDOF,NActvDOF,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the DampMat array.' )
ENDIF
ALLOCATE ( StffMat(NActvDOF,NActvDOF,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the StffMat array.' )
ENDIF
ALLOCATE ( FMat(NActvDOF,NInputs,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the FMat array.' )
ENDIF
ALLOCATE ( FdMat(NActvDOF,NDisturbs,NAzimStep) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort ( ' Error allocating memory for the FdMat array.' )
ENDIF
ENDIF
! Input the magnitudes of the displacement and velocity pertubations:
! Initialize the displacement pertubations to 2 degrees (same as SymDyn):
DelQ = 2.0*D2R
! For displacement DOFs of the flexible members (i.e., blade and tower
! flexibility), dimensionalize these pertubations using the reciprical of
! the weights used to nondimensionalize the dislacement norm. Also,
! change the 2 degrees to 0.2 degrees, to make the magnitudes of the
! pertubations more reasonable.
DelQ(DOF_TFA1) = 0.2*D2R/SHP( 1.0, TwrFlexL, TwFAM1Sh(:), 1 )
DelQ(DOF_TSS1) = 0.2*D2R/SHP( 1.0, TwrFlexL, TwSSM1Sh(:), 1 )
DelQ(DOF_TFA2) = 0.2*D2R/SHP( 1.0, TwrFlexL, TwFAM2Sh(:), 1 )
DelQ(DOF_TSS2) = 0.2*D2R/SHP( 1.0, TwrFlexL, TwSSM2Sh(:), 1 )
DO K = 1,NumBl ! Loop through all blades
DelQ(DOF_BF(K,1)) = 0.2*D2R/SHP( 1.0, BldFlexL, BldFl1Sh(:,K), 1 )
DelQ(DOF_BF(K,2)) = 0.2*D2R/SHP( 1.0, BldFlexL, BldFl2Sh(:,K), 1 )
DelQ(DOF_BE(K,1)) = 0.2*D2R/SHP( 1.0, BldFlexL, BldEdgSh(:,K), 1 )
ENDDO ! K - Blades
! Like the flexible members, for the translational platform DOFs,
! dimensionalize these pertubations using the reciprical of the weights
! used to nondimensionalize the dislacement norm. Also, change the 2
! degrees to 0.2 degrees, to make the magnitudes of the pertubations
! more reasonable.
DelQ(DOF_Sg ) = 0.2*D2R/0.01 ! Assume 100m
DelQ(DOF_Sw ) = 0.2*D2R/0.01 ! reference
DelQ(DOF_Hv ) = 0.2*D2R/0.01 ! depth
! Use the same pertubations for velocity:
DelQD = DelQ
! Input the magnitudes of control input pertubations:
DelYawPos = 2.0 *D2R ! Pertubation in nacelle yaw angle (same as for DOFs) (rad)
DelYawRate = 2.0 *D2R ! Pertubation in nacelle yaw rate (same as for DOF velocities) (rad/s)
DelGenTq = ( 0.15*AvgNrmTpRd**3 )/GBRatio ! Pertubation in generator torque. The rated rotor torque of a typical wind turbine in N-m is roughly 15*Radius^3 where the radius is in meters. The 0.15 is 1/100th of the factor of 15. The 1/GBRatio is used to cast the torque onto the HSS (i.e., generator) side of the gearbox. (Nm)
DelBlPtch = 0.2 *D2R ! Pertubation in blade pitch (same as SymDyn) (rad)
DelV = 0.1 ! Pertubation in horizontal hub height wind speed (same as SymDyn) (m/s)
DelDELTA = 2.0 *D2R ! Pertubation in horizontal wind direction (rad)
DelVZ = 0.1 ! Pertubation in vertical wind speed (same as SymDyn) (m/s)
DelHSHR = 0.01 ! Pertubation in horizontal shear coefficient (same as SymDyn)
DelVSHR = 0.01 ! Pertubation in vertical shear coefficient (same as SymDyn)
DelVLINSHR = 0.01 ! Pertubation in vertical (linear) shear coefficient (same as SymDyn)
DelVGUST = 0.1 ! Pertubation in horizontal hub height wind gust (same as SymDyn) (m/s)
! Let's save the current (operating point) values of the inputs and
! disturbances:
IF ( YawDOF ) THEN ! Yaw DOF is currently enabled (using FAST's built-in actuator) - use the current command yaw angle, YawNeut, and command yaw rate, YawRateNeut.
YawPosop = YawNeut ! Yaw angle (position)
YawRateop = YawRateNeut ! Yaw rate
ELSE ! Yaw DOF is currently disabled (no built-in actuator) - use the current (operating point) yaw angle, Qop(DOF_Yaw,1), and yaw rate, QD(DOF_Yaw,1)--this is constant for all azimuth steps; therefore, use the first one which is always available and not dependent on input NAzimStep.
YawPosop = Qop (DOF_Yaw,1) ! Yaw angle (position)
YawRateop = QDop(DOF_Yaw,1) ! Yaw rate
ENDIF
! No need to store the current generator torque since we perturb this input slightly differently
BlPitchop = BlPitch ! Blade pitch
GenTrqop = GenTrq ! Generator torque
! Open the FAST linear file (.lin) and give it a heading:
CALL OpenFOutFile ( UnLn, TRIM( RootName )//'.lin' )
WRITE (UnLn,'(/,A)' ) 'This linearized model file was generated by '//TRIM(ProgName)//' '//TRIM( ProgVer )// &
' on '//CurDate()//' at '//CurTime()//'.'
WRITE (UnLn,FmtText ) 'The aerodynamic calculations were made by '//TRIM(AD_Prog%Name)//' '//TRIM(AD_Prog%Ver)//'.'
WRITE (UnLn,'(/,1X,A)') TRIM( FTitle )
! Output some useful information:
WRITE (UnLn,FmtHead) 'Some Useful Information:'
IF ( .NOT. CalcStdy ) THEN ! No steady state solution found (linearized about initial conditions)
WRITE(UnLn,FmtText ) ' Type of steady state solution found None (linearized about initial conditions)'
ELSEIF ( .NOT. GenDOF ) THEN ! Constant speed case
IF ( RotSpeed == 0.0 ) THEN
WRITE(UnLn,FmtText ) ' Type of steady state solution found Static equilibrium (RotSpeed = 0.0)'