/
openmdao_openfast.py
2678 lines (2295 loc) · 165 KB
/
openmdao_openfast.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
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
import numpy as np
import pandas as pd
import os, shutil, sys, platform
import copy
import glob
from pathlib import Path
from scipy.interpolate import PchipInterpolator
from openmdao.api import ExplicitComponent
from wisdem.commonse.mpi_tools import MPI
from wisdem.commonse import NFREQ
from wisdem.commonse.cylinder_member import get_nfull
import wisdem.commonse.utilities as util
from wisdem.rotorse.rotor_power import eval_unsteady
from weis.aeroelasticse.FAST_writer import InputWriter_OpenFAST
from weis.aeroelasticse.FAST_reader import InputReader_OpenFAST
import weis.aeroelasticse.runFAST_pywrapper as fastwrap
from weis.aeroelasticse.FAST_post import FAST_IO_timeseries
from wisdem.floatingse.floating_frame import NULL, NNODES_MAX, NELEM_MAX
from weis.dlc_driver.dlc_generator import DLCGenerator
from weis.aeroelasticse.CaseGen_General import CaseGen_General
from functools import partial
from pCrunch import PowerProduction
from weis.aeroelasticse.LinearFAST import LinearFAST
from weis.control.LinearModel import LinearTurbineModel, LinearControlModel
from weis.aeroelasticse import FileTools
from weis.aeroelasticse.turbsim_file import TurbSimFile
from weis.aeroelasticse.turbsim_util import generate_wind_files
from weis.aeroelasticse.utils import OLAFParams
from ROSCO_toolbox import control_interface as ROSCO_ci
from pCrunch.io import OpenFASTOutput
from pCrunch import LoadsAnalysis, PowerProduction, FatigueParams
from weis.control.dtqp_wrapper import dtqp_wrapper
from weis.aeroelasticse.StC_defaults import default_StC_vt
from weis.aeroelasticse.CaseGen_General import case_naming
weis_dir = os.path.dirname(os.path.dirname(os.path.dirname(__file__)))
import pickle
if MPI:
from mpi4py import MPI
if platform.system() == 'Windows':
lib_ext = '.dll'
elif platform.system() == 'Darwin':
lib_ext = '.dylib'
else:
lib_ext = '.so'
weis_dir = os.path.dirname(os.path.dirname(os.path.dirname(__file__)))
def make_coarse_grid(s_grid, diam):
s_coarse = [s_grid[0]]
slope = np.diff(diam) / np.diff(s_grid)
for k in range(slope.size-1):
if np.abs(slope[k]-slope[k+1]) > 1e-2:
s_coarse.append(s_grid[k+1])
s_coarse.append(s_grid[-1])
return np.array(s_coarse)
class FASTLoadCases(ExplicitComponent):
def initialize(self):
self.options.declare('modeling_options')
self.options.declare('opt_options')
def setup(self):
modopt = self.options['modeling_options']
rotorse_options = modopt['WISDEM']['RotorSE']
mat_init_options = modopt['materials']
self.n_blades = modopt['assembly']['number_of_blades']
self.n_span = n_span = rotorse_options['n_span']
self.n_pc = n_pc = rotorse_options['n_pc']
# Environmental Conditions needed regardless of where model comes from
self.add_input('V_cutin', val=0.0, units='m/s', desc='Minimum wind speed where turbine operates (cut-in)')
self.add_input('V_cutout', val=0.0, units='m/s', desc='Maximum wind speed where turbine operates (cut-out)')
self.add_input('Vrated', val=0.0, units='m/s', desc='rated wind speed')
self.add_input('hub_height', val=0.0, units='m', desc='hub height')
self.add_discrete_input('turbulence_class', val='A', desc='IEC turbulence class')
self.add_discrete_input('turbine_class', val='I', desc='IEC turbine class')
self.add_input('Rtip', val=0.0, units='m', desc='dimensional radius of tip')
self.add_input('shearExp', val=0.0, desc='shear exponent')
if not self.options['modeling_options']['Level3']['from_openfast']:
self.n_pitch = n_pitch = rotorse_options['n_pitch_perf_surfaces']
self.n_tsr = n_tsr = rotorse_options['n_tsr_perf_surfaces']
self.n_U = n_U = rotorse_options['n_U_perf_surfaces']
self.n_mat = n_mat = mat_init_options['n_mat']
self.n_layers = n_layers = rotorse_options['n_layers']
self.n_xy = n_xy = rotorse_options['n_xy'] # Number of coordinate points to describe the airfoil geometry
self.n_aoa = n_aoa = rotorse_options['n_aoa']# Number of angle of attacks
self.n_Re = n_Re = rotorse_options['n_Re'] # Number of Reynolds, so far hard set at 1
self.n_tab = n_tab = rotorse_options['n_tab']# Number of tabulated data. For distributed aerodynamic control this could be > 1
self.te_ss_var = rotorse_options['te_ss']
self.te_ps_var = rotorse_options['te_ps']
self.spar_cap_ss_var = rotorse_options['spar_cap_ss']
self.spar_cap_ps_var = rotorse_options['spar_cap_ps']
n_freq_blade = int(rotorse_options['n_freq']/2)
n_pc = int(rotorse_options['n_pc'])
self.n_xy = n_xy = rotorse_options['n_xy'] # Number of coordinate points to describe the airfoil geometry
self.n_aoa = n_aoa = rotorse_options['n_aoa']# Number of angle of attacks
self.n_Re = n_Re = rotorse_options['n_Re'] # Number of Reynolds, so far hard set at 1
self.n_tab = n_tab = rotorse_options['n_tab']# Number of tabulated data. For distributed aerodynamic control this could be > 1
self.te_ss_var = rotorse_options['te_ss']
self.te_ps_var = rotorse_options['te_ps']
self.spar_cap_ss_var = rotorse_options['spar_cap_ss']
self.spar_cap_ps_var = rotorse_options['spar_cap_ps']
# ElastoDyn Inputs
# Assuming the blade modal damping to be unchanged. Cannot directly solve from the Rayleigh Damping without making assumptions. J.Jonkman recommends 2-3% https://wind.nrel.gov/forum/wind/viewtopic.php?t=522
self.add_input('r', val=np.zeros(n_span), units='m', desc='radial positions. r[0] should be the hub location \
while r[-1] should be the blade tip. Any number \
of locations can be specified between these in ascending order.')
self.add_input('le_location', val=np.zeros(n_span), desc='Leading-edge positions from a reference blade axis (usually blade pitch axis). Locations are normalized by the local chord length. Positive in -x direction for airfoil-aligned coordinate system')
self.add_input('beam:Tw_iner', val=np.zeros(n_span), units='m', desc='y-distance to elastic center from point about which above structural properties are computed')
self.add_input('beam:rhoA', val=np.zeros(n_span), units='kg/m', desc='mass per unit length')
self.add_input('beam:EIyy', val=np.zeros(n_span), units='N*m**2', desc='flatwise stiffness (bending about y-direction of airfoil aligned coordinate system)')
self.add_input('beam:EIxx', val=np.zeros(n_span), units='N*m**2', desc='edgewise stiffness (bending about :ref:`x-direction of airfoil aligned coordinate system <blade_airfoil_coord>`)')
self.add_input('x_tc', val=np.zeros(n_span), units='m', desc='x-distance to the neutral axis (torsion center)')
self.add_input('y_tc', val=np.zeros(n_span), units='m', desc='y-distance to the neutral axis (torsion center)')
self.add_input('flap_mode_shapes', val=np.zeros((n_freq_blade,5)), desc='6-degree polynomial coefficients of mode shapes in the flap direction (x^2..x^6, no linear or constant term)')
self.add_input('edge_mode_shapes', val=np.zeros((n_freq_blade,5)), desc='6-degree polynomial coefficients of mode shapes in the edge direction (x^2..x^6, no linear or constant term)')
self.add_input('gearbox_efficiency', val=1.0, desc='Gearbox efficiency')
self.add_input('gearbox_ratio', val=1.0, desc='Gearbox ratio')
self.add_input('platform_displacement', val=1.0, desc='Volumetric platform displacement', units='m**3')
# ServoDyn Inputs
self.add_input('generator_efficiency', val=1.0, desc='Generator efficiency')
self.add_input('max_pitch_rate', val=0.0, units='deg/s', desc='Maximum allowed blade pitch rate')
# StC or TMD inputs; structural control and tuned mass dampers
# tower properties
n_height_tow = modopt['WISDEM']['TowerSE']['n_height']
n_full_tow = get_nfull(n_height_tow, nref=modopt['WISDEM']['TowerSE']['n_refine'])
n_freq_tower = int(NFREQ/2)
self.add_input('fore_aft_modes', val=np.zeros((n_freq_tower,5)), desc='6-degree polynomial coefficients of mode shapes in the flap direction (x^2..x^6, no linear or constant term)')
self.add_input('side_side_modes', val=np.zeros((n_freq_tower,5)), desc='6-degree polynomial coefficients of mode shapes in the edge direction (x^2..x^6, no linear or constant term)')
self.add_input('mass_den', val=np.zeros(n_height_tow-1), units='kg/m', desc='sectional mass per unit length')
self.add_input('foreaft_stff', val=np.zeros(n_height_tow-1), units='N*m**2', desc='sectional fore-aft bending stiffness per unit length about the Y_E elastic axis')
self.add_input('sideside_stff', val=np.zeros(n_height_tow-1), units='N*m**2', desc='sectional side-side bending stiffness per unit length about the Y_E elastic axis')
self.add_input('tor_stff', val=np.zeros(n_height_tow-1), units='N*m**2', desc='torsional stiffness per unit length about the Y_E elastic axis')
self.add_input('tor_freq', val=0.0, units='Hz', desc='First tower torsional frequency')
self.add_input('tower_outer_diameter', val=np.zeros(n_height_tow), units='m', desc='cylinder diameter at corresponding locations')
self.add_input('tower_z', val=np.zeros(n_height_tow), units='m', desc='z-coordinates of tower and monopile used in TowerSE')
self.add_input('tower_z_full', val=np.zeros(n_full_tow), units='m', desc='z-coordinates of tower and monopile used in TowerSE')
self.add_input('tower_height', val=0.0, units='m', desc='tower height from the tower base')
self.add_input('tower_base_height', val=0.0, units='m', desc='tower base height from the ground or mean sea level')
self.add_input('tower_cd', val=np.zeros(n_height_tow), desc='drag coefficients along tower height at corresponding locations')
self.add_input("tower_I_base", np.zeros(6), units="kg*m**2", desc="tower moments of inertia at the tower base")
# These next ones are needed for SubDyn
n_height_mon = n_full_mon = 0
if modopt['flags']['offshore']:
self.add_input('transition_piece_mass', val=0.0, units='kg')
self.add_input('transition_piece_I', val=np.zeros(3), units='kg*m**2')
if modopt['flags']['monopile']:
n_height_mon = modopt['WISDEM']['FixedBottomSE']['n_height']
n_full_mon = get_nfull(n_height_mon, nref=modopt['WISDEM']['FixedBottomSE']['n_refine'])
self.add_input('monopile_z', val=np.zeros(n_height_mon), units='m', desc='z-coordinates of tower and monopile used in TowerSE')
self.add_input('monopile_z_full', val=np.zeros(n_full_mon), units='m', desc='z-coordinates of tower and monopile used in TowerSE')
self.add_input('monopile_outer_diameter', val=np.zeros(n_height_mon), units='m', desc='cylinder diameter at corresponding locations')
self.add_input('monopile_wall_thickness', val=np.zeros(n_height_mon-1), units='m')
self.add_input('monopile_E', val=np.zeros(n_height_mon-1), units='Pa')
self.add_input('monopile_G', val=np.zeros(n_height_mon-1), units='Pa')
self.add_input('monopile_rho', val=np.zeros(n_height_mon-1), units='kg/m**3')
self.add_input('gravity_foundation_mass', val=0.0, units='kg')
self.add_input('gravity_foundation_I', val=np.zeros(3), units='kg*m**2')
monlen = max(0, n_height_mon-1)
monlen_full = max(0, n_full_mon-1)
# DriveSE quantities
self.add_input('hub_system_cm', val=np.zeros(3), units='m', desc='center of mass of the hub relative to tower to in yaw-aligned c.s.')
self.add_input('hub_system_I', val=np.zeros(6), units='kg*m**2', desc='mass moments of Inertia of hub [Ixx, Iyy, Izz, Ixy, Ixz, Iyz] around its center of mass in yaw-aligned c.s.')
self.add_input('hub_system_mass', val=0.0, units='kg', desc='mass of hub system')
self.add_input('above_yaw_mass', val=0.0, units='kg', desc='Mass of the nacelle above the yaw system')
self.add_input('yaw_mass', val=0.0, units='kg', desc='Mass of yaw system')
self.add_input('rna_I_TT', val=np.zeros(6), units='kg*m**2', desc=' moments of Inertia for the rna [Ixx, Iyy, Izz, Ixy, Ixz, Iyz] about the tower top')
self.add_input('nacelle_cm', val=np.zeros(3), units='m', desc='Center of mass of the component in [x,y,z] for an arbitrary coordinate system')
self.add_input('nacelle_I_TT', val=np.zeros(6), units='kg*m**2', desc=' moments of Inertia for the nacelle [Ixx, Iyy, Izz, Ixy, Ixz, Iyz] about the tower top')
self.add_input('distance_tt_hub', val=0.0, units='m', desc='Vertical distance from tower top plane to hub flange')
self.add_input('twr2shaft', val=0.0, units='m', desc='Vertical distance from tower top plane to shaft start')
self.add_input('GenIner', val=0.0, units='kg*m**2', desc='Moments of inertia for the generator about high speed shaft')
self.add_input('drivetrain_spring_constant', val=0.0, units='N*m/rad', desc='Moments of inertia for the generator about high speed shaft')
self.add_input("drivetrain_damping_coefficient", 0.0, units="N*m*s/rad", desc='Equivalent damping coefficient for the drivetrain system')
# AeroDyn Inputs
self.add_input('ref_axis_blade', val=np.zeros((n_span,3)),units='m', desc='2D array of the coordinates (x,y,z) of the blade reference axis, defined along blade span. The coordinate system is the one of BeamDyn: it is placed at blade root with x pointing the suction side of the blade, y pointing the trailing edge and z along the blade span. A standard configuration will have negative x values (prebend), if swept positive y values, and positive z values.')
self.add_input('chord', val=np.zeros(n_span), units='m', desc='chord at airfoil locations')
self.add_input('theta', val=np.zeros(n_span), units='deg', desc='twist at airfoil locations')
self.add_input('rthick', val=np.zeros(n_span), desc='relative thickness of airfoil distribution')
self.add_input('ac', val=np.zeros(n_span), desc='aerodynamic center of airfoil distribution')
self.add_input('pitch_axis', val=np.zeros(n_span), desc='1D array of the chordwise position of the pitch axis (0-LE, 1-TE), defined along blade span.')
self.add_input('Rhub', val=0.0, units='m', desc='dimensional radius of hub')
self.add_input('airfoils_cl', val=np.zeros((n_span, n_aoa, n_Re, n_tab)), desc='lift coefficients, spanwise')
self.add_input('airfoils_cd', val=np.zeros((n_span, n_aoa, n_Re, n_tab)), desc='drag coefficients, spanwise')
self.add_input('airfoils_cm', val=np.zeros((n_span, n_aoa, n_Re, n_tab)), desc='moment coefficients, spanwise')
self.add_input('airfoils_aoa', val=np.zeros((n_aoa)), units='deg', desc='angle of attack grid for polars')
self.add_input('airfoils_Re', val=np.zeros((n_Re)), desc='Reynolds numbers of polars')
self.add_input('airfoils_Ctrl', val=np.zeros((n_span, n_Re, n_tab)), units='deg',desc='Airfoil control paremeter (i.e. flap angle)')
# Airfoil coordinates
self.add_input('coord_xy_interp', val=np.zeros((n_span, n_xy, 2)), desc='3D array of the non-dimensional x and y airfoil coordinates of the airfoils interpolated along span for n_span stations. The leading edge is place at x=0 and y=0.')
# Floating platform inputs
self.add_input("transition_node", np.zeros(3), units="m")
self.add_input("platform_nodes", NULL * np.ones((NNODES_MAX, 3)), units="m")
self.add_input("platform_elem_n1", NULL * np.ones(NELEM_MAX, dtype=np.int_))
self.add_input("platform_elem_n2", NULL * np.ones(NELEM_MAX, dtype=np.int_))
self.add_input("platform_elem_D", NULL * np.ones(NELEM_MAX), units="m")
self.add_input("platform_elem_t", NULL * np.ones(NELEM_MAX), units="m")
self.add_input("platform_elem_rho", NULL * np.ones(NELEM_MAX), units="kg/m**3")
self.add_input("platform_elem_E", NULL * np.ones(NELEM_MAX), units="Pa")
self.add_input("platform_elem_G", NULL * np.ones(NELEM_MAX), units="Pa")
self.add_discrete_input("platform_elem_memid", [0]*NELEM_MAX)
self.add_input("platform_total_center_of_mass", np.zeros(3), units="m")
self.add_input("platform_mass", 0.0, units="kg")
self.add_input("platform_I_total", np.zeros(6), units="kg*m**2")
if modopt['flags']["floating"]:
n_member = modopt["floating"]["members"]["n_members"]
for k in range(n_member):
n_height_mem = modopt["floating"]["members"]["n_height"][k]
self.add_input(f"member{k}:joint1", np.zeros(3), units="m")
self.add_input(f"member{k}:joint2", np.zeros(3), units="m")
self.add_input(f"member{k}:s", np.zeros(n_height_mem))
self.add_input(f"member{k}:s_ghost1", 0.0)
self.add_input(f"member{k}:s_ghost2", 0.0)
self.add_input(f"member{k}:outer_diameter", np.zeros(n_height_mem), units="m")
self.add_input(f"member{k}:wall_thickness", np.zeros(n_height_mem-1), units="m")
# Turbine level inputs
self.add_discrete_input('rotor_orientation',val='upwind', desc='Rotor orientation, either upwind or downwind.')
self.add_input('control_ratedPower', val=0., units='W', desc='machine power rating')
self.add_input('control_maxOmega', val=0.0, units='rpm', desc='maximum allowed rotor rotation speed')
self.add_input('control_maxTS', val=0.0, units='m/s', desc='maximum allowed blade tip speed')
self.add_input('cone', val=0.0, units='deg', desc='Cone angle of the rotor. It defines the angle between the rotor plane and the blade pitch axis. A standard machine has positive values.')
self.add_input('tilt', val=0.0, units='deg', desc='Nacelle uptilt angle. A standard machine has positive values.')
self.add_input('overhang', val=0.0, units='m', desc='Horizontal distance from tower top to hub center.')
# Initial conditions
self.add_input('U', val=np.zeros(n_pc), units='m/s', desc='wind speeds')
self.add_input('Omega', val=np.zeros(n_pc), units='rpm', desc='rotation speeds to run')
self.add_input('pitch', val=np.zeros(n_pc), units='deg', desc='pitch angles to run')
# Cp-Ct-Cq surfaces
self.add_input('Cp_aero_table', val=np.zeros((n_tsr, n_pitch, n_U)), desc='Table of aero power coefficient')
self.add_input('Ct_aero_table', val=np.zeros((n_tsr, n_pitch, n_U)), desc='Table of aero thrust coefficient')
self.add_input('Cq_aero_table', val=np.zeros((n_tsr, n_pitch, n_U)), desc='Table of aero torque coefficient')
self.add_input('pitch_vector', val=np.zeros(n_pitch), units='deg', desc='Pitch vector used')
self.add_input('tsr_vector', val=np.zeros(n_tsr), desc='TSR vector used')
self.add_input('U_vector', val=np.zeros(n_U), units='m/s', desc='Wind speed vector used')
# Environmental conditions
self.add_input('V_R25', val=0.0, units='m/s', desc='region 2.5 transition wind speed')
self.add_input('Vgust', val=0.0, units='m/s', desc='gust wind speed')
self.add_input('V_extreme1', val=0.0, units='m/s', desc='IEC extreme wind speed at hub height for a 1-year retunr period')
self.add_input('V_extreme50', val=0.0, units='m/s', desc='IEC extreme wind speed at hub height for a 50-year retunr period')
self.add_input('V_mean_iec', val=0.0, units='m/s', desc='IEC mean wind for turbulence class')
self.add_input('rho', val=0.0, units='kg/m**3', desc='density of air')
self.add_input('mu', val=0.0, units='kg/(m*s)', desc='dynamic viscosity of air')
self.add_input('speed_sound_air', val=340., units='m/s', desc='Speed of sound in air.')
self.add_input(
"water_depth", val=0.0, units="m", desc="Water depth for analysis. Values > 0 mean offshore"
)
self.add_input('rho_water', val=0.0, units='kg/m**3', desc='density of water')
self.add_input('mu_water', val=0.0, units='kg/(m*s)', desc='dynamic viscosity of water')
self.add_input('beta_wave', val=0.0, units='deg', desc='Incident wave propagation heading direction')
self.add_input('Hsig_wave', val=0.0, units='m', desc='Significant wave height of incident waves')
self.add_input('Tsig_wave', val=0.0, units='s', desc='Peak-spectral period of incident waves')
# Blade composite material properties (used for fatigue analysis)
self.add_input('gamma_f', val=1.35, desc='safety factor on loads')
self.add_input('gamma_m', val=1.1, desc='safety factor on materials')
self.add_input('E', val=np.zeros([n_mat, 3]), units='Pa', desc='2D array of the Youngs moduli of the materials. Each row represents a material, the three columns represent E11, E22 and E33.')
self.add_input('Xt', val=np.zeros([n_mat, 3]), units='Pa', desc='2D array of the Ultimate Tensile Strength (UTS) of the materials. Each row represents a material, the three columns represent Xt12, Xt13 and Xt23.')
self.add_input('Xc', val=np.zeros([n_mat, 3]), units='Pa', desc='2D array of the Ultimate Compressive Strength (UCS) of the materials. Each row represents a material, the three columns represent Xc12, Xc13 and Xc23.')
self.add_input('m', val=np.zeros([n_mat]), desc='2D array of the S-N fatigue slope exponent for the materials')
# Blade composit layup info (used for fatigue analysis)
self.add_input('sc_ss_mats', val=np.zeros((n_span, n_mat)), desc="spar cap, suction side, boolean of materials in each composite layer spanwise, passed as floats for differentiablity, used for Fatigue Analysis")
self.add_input('sc_ps_mats', val=np.zeros((n_span, n_mat)), desc="spar cap, pressure side, boolean of materials in each composite layer spanwise, passed as floats for differentiablity, used for Fatigue Analysis")
self.add_input('te_ss_mats', val=np.zeros((n_span, n_mat)), desc="trailing edge reinforcement, suction side, boolean of materials in each composite layer spanwise, passed as floats for differentiablity, used for Fatigue Analysis")
self.add_input('te_ps_mats', val=np.zeros((n_span, n_mat)), desc="trailing edge reinforcement, pressure side, boolean of materials in each composite layer spanwise, passed as floats for differentiablity, used for Fatigue Analysis")
self.add_discrete_input('definition_layer', val=np.zeros(n_layers), desc='1D array of flags identifying how layers are specified in the yaml. 1) all around (skin, paint, ) 2) offset+rotation twist+width (spar caps) 3) offset+user defined rotation+width 4) midpoint TE+width (TE reinf) 5) midpoint LE+width (LE reinf) 6) layer position fixed to other layer (core fillers) 7) start and width 8) end and width 9) start and end nd 10) web layer')
# self.add_discrete_input('layer_name', val=n_layers * [''], desc='1D array of the names of the layers modeled in the blade structure.')
# self.add_discrete_input('layer_web', val=n_layers * [''], desc='1D array of the names of the webs the layer is associated to. If the layer is on the outer profile this entry can simply stay empty.')
# self.add_discrete_input('layer_mat', val=n_layers * [''], desc='1D array of the names of the materials of each layer modeled in the blade structure.')
self.layer_name = rotorse_options['layer_name']
# MoorDyn inputs
mooropt = modopt["mooring"]
if self.options["modeling_options"]["flags"]["mooring"]:
n_nodes = mooropt["n_nodes"]
n_lines = mooropt["n_lines"]
self.add_input("line_diameter", val=np.zeros(n_lines), units="m")
self.add_input("line_mass_density", val=np.zeros(n_lines), units="kg/m")
self.add_input("line_stiffness", val=np.zeros(n_lines), units="N")
self.add_input("line_transverse_added_mass", val=np.zeros(n_lines), units="kg/m")
self.add_input("line_tangential_added_mass", val=np.zeros(n_lines), units="kg/m")
self.add_input("line_transverse_drag", val=np.zeros(n_lines))
self.add_input("line_tangential_drag", val=np.zeros(n_lines))
self.add_input("nodes_location_full", val=np.zeros((n_nodes, 3)), units="m")
self.add_input("nodes_mass", val=np.zeros(n_nodes), units="kg")
self.add_input("nodes_volume", val=np.zeros(n_nodes), units="m**3")
self.add_input("nodes_added_mass", val=np.zeros(n_nodes))
self.add_input("nodes_drag_area", val=np.zeros(n_nodes), units="m**2")
self.add_input("unstretched_length", val=np.zeros(n_lines), units="m")
self.add_discrete_input("node_names", val=[""] * n_nodes)
# Inputs required for fatigue processing
self.add_input('lifetime', val=25.0, units='yr', desc='Turbine design lifetime')
self.add_input('blade_sparU_wohlerexp', val=1.0, desc='Blade root Wohler exponent, m, in S/N curve S=A*N^-(1/m)')
self.add_input('blade_sparU_wohlerA', val=1.0, units="Pa", desc='Blade root parameter, A, in S/N curve S=A*N^-(1/m)')
self.add_input('blade_sparU_ultstress', val=1.0, units="Pa", desc='Blade root ultimate stress for material')
self.add_input('blade_sparL_wohlerexp', val=1.0, desc='Blade root Wohler exponent, m, in S/N curve S=A*N^-(1/m)')
self.add_input('blade_sparL_wohlerA', val=1.0, units="Pa", desc='Blade root parameter, A, in S/N curve S=A*N^-(1/m)')
self.add_input('blade_sparL_ultstress', val=1.0, units="Pa", desc='Blade root ultimate stress for material')
self.add_input('blade_teU_wohlerexp', val=1.0, desc='Blade root Wohler exponent, m, in S/N curve S=A*N^-(1/m)')
self.add_input('blade_teU_wohlerA', val=1.0, units="Pa", desc='Blade root parameter, A, in S/N curve S=A*N^-(1/m)')
self.add_input('blade_teU_ultstress', val=1.0, units="Pa", desc='Blade root ultimate stress for material')
self.add_input('blade_teL_wohlerexp', val=1.0, desc='Blade root Wohler exponent, m, in S/N curve S=A*N^-(1/m)')
self.add_input('blade_teL_wohlerA', val=1.0, units="Pa", desc='Blade root parameter, A, in S/N curve S=A*N^-(1/m)')
self.add_input('blade_teL_ultstress', val=1.0, units="Pa", desc='Blade root ultimate stress for material')
self.add_input('blade_root_sparU_load2stress', val=np.ones(6), units="m**2", desc='Blade root upper spar cap coefficient between axial load and stress S=C^T [Fx-z;Mx-z]')
self.add_input('blade_root_sparL_load2stress', val=np.ones(6), units="m**2", desc='Blade root lower spar cap coefficient between axial load and stress S=C^T [Fx-z;Mx-z]')
self.add_input('blade_maxc_teU_load2stress', val=np.ones(6), units="m**2", desc='Blade max chord upper trailing edge coefficient between axial load and stress S=C^T [Fx-z;Mx-z]')
self.add_input('blade_maxc_teL_load2stress', val=np.ones(6), units="m**2", desc='Blade max chord lower trailing edge coefficient between axial load and stress S=C^T [Fx-z;Mx-z]')
self.add_input('lss_wohlerexp', val=1.0, desc='Low speed shaft Wohler exponent, m, in S/N curve S=A*N^-(1/m)')
self.add_input('lss_wohlerA', val=1.0, desc='Low speed shaft parameter, A, in S/N curve S=A*N^-(1/m)')
self.add_input('lss_ultstress', val=1.0, units="Pa", desc='Low speed shaft Ultimate stress for material')
self.add_input('lss_axial_load2stress', val=np.ones(6), units="m**2", desc='Low speed shaft coefficient between axial load and stress S=C^T [Fx-z;Mx-z]')
self.add_input('lss_shear_load2stress', val=np.ones(6), units="m**2", desc='Low speed shaft coefficient between shear load and stress S=C^T [Fx-z;Mx-z]')
self.add_input('tower_wohlerexp', val=np.ones(n_height_tow-1), desc='Tower Wohler exponent, m, in S/N curve S=A*N^-(1/m)')
self.add_input('tower_wohlerA', val=np.ones(n_height_tow-1), desc='Tower parameter, A, in S/N curve S=A*N^-(1/m)')
self.add_input('tower_ultstress', val=np.ones(n_height_tow-1), units="Pa", desc='Tower ultimate stress for material')
self.add_input('tower_axial_load2stress', val=np.ones([n_height_tow-1,6]), units="m**2", desc='Tower coefficient between axial load and stress S=C^T [Fx-z;Mx-z]')
self.add_input('tower_shear_load2stress', val=np.ones([n_height_tow-1,6]), units="m**2", desc='Tower coefficient between shear load and stress S=C^T [Fx-z;Mx-z]')
self.add_input('monopile_wohlerexp', val=np.ones(monlen), desc='Tower Wohler exponent, m, in S/N curve S=A*N^-(1/m)')
self.add_input('monopile_wohlerA', val=np.ones(monlen), desc='Tower parameter, A, in S/N curve S=A*N^-(1/m)')
self.add_input('monopile_ultstress', val=np.ones(monlen), units="Pa", desc='Tower ultimate stress for material')
self.add_input('monopile_axial_load2stress', val=np.ones([monlen,6]), units="m**2", desc='Tower coefficient between axial load and stress S=C^T [Fx-z;Mx-z]')
self.add_input('monopile_shear_load2stress', val=np.ones([monlen,6]), units="m**2", desc='Tower coefficient between shear load and stress S=C^T [Fx-z;Mx-z]')
# TMD params
if self.options['modeling_options']['flags']['TMDs']:
n_TMDs = self.options['modeling_options']['TMDs']['n_TMDs']
self.add_input('TMD_mass', val=np.zeros(n_TMDs), units='kg', desc='TMD Mass')
self.add_input('TMD_stiffness', val=np.zeros(n_TMDs), units='N/m', desc='TMD Stiffnes')
self.add_input('TMD_damping', val=np.zeros(n_TMDs), units='N/(m/s)', desc='TMD Damping')
# DLC options
n_ws_dlc11 = modopt['DLC_driver']['n_ws_dlc11']
# OpenFAST options
OFmgmt = modopt['General']['openfast_configuration']
self.model_only = OFmgmt['model_only']
FAST_directory_base = OFmgmt['OF_run_dir']
# If the path is relative, make it an absolute path to current working directory
if not os.path.isabs(FAST_directory_base):
FAST_directory_base = os.path.join(os.getcwd(), FAST_directory_base)
# Flag to clear OpenFAST run folder. Use it only if disk space is an issue
self.clean_FAST_directory = False
self.FAST_InputFile = OFmgmt['OF_run_fst']
# File naming changes whether in MPI or not
if MPI:
rank = MPI.COMM_WORLD.Get_rank()
self.FAST_runDirectory = os.path.join(FAST_directory_base,'rank_%000d'%int(rank))
self.FAST_namingOut = self.FAST_InputFile+'_%000d'%int(rank)
else:
self.FAST_runDirectory = FAST_directory_base
self.FAST_namingOut = self.FAST_InputFile
self.wind_directory = os.path.join(self.FAST_runDirectory, 'wind')
if not os.path.exists(self.FAST_runDirectory):
os.makedirs(self.FAST_runDirectory, exist_ok=True)
if not os.path.exists(self.wind_directory):
os.mkdir(self.wind_directory)
# Number of cores used outside of MPI. If larger than 1, the multiprocessing module is called
self.cores = OFmgmt['cores']
self.case = {}
self.channels = {}
self.mpi_run = False
if 'mpi_run' in OFmgmt.keys():
self.mpi_run = OFmgmt['mpi_run']
if self.mpi_run:
self.mpi_comm_map_down = OFmgmt['mpi_comm_map_down']
# User-defined FAST library/executable
if OFmgmt['FAST_exe'] != 'none':
if os.path.isabs(OFmgmt['FAST_exe']):
self.FAST_exe = OFmgmt['FAST_exe']
else:
self.FAST_exe = os.path.join(os.path.dirname(self.options['modeling_options']['fname_input_modeling']),
OFmgmt['FAST_exe'])
else:
self.FAST_exe = 'none'
if OFmgmt['FAST_lib'] != 'none':
if os.path.isabs(OFmgmt['FAST_lib']):
self.FAST_lib = OFmgmt['FAST_lib']
else:
self.FAST_lib = os.path.join(os.path.dirname(self.options['modeling_options']['fname_input_modeling']),
OFmgmt['FAST_lib'])
else:
self.FAST_lib = 'none'
# Rotor power outputs
self.add_output('V_out', val=np.zeros(n_ws_dlc11), units='m/s', desc='wind speed vector from the OF simulations')
self.add_output('P_out', val=np.zeros(n_ws_dlc11), units='W', desc='rotor electrical power')
self.add_output('Cp_out', val=np.zeros(n_ws_dlc11), desc='rotor aero power coefficient')
self.add_output('Omega_out', val=np.zeros(n_ws_dlc11), units='rpm', desc='rotation speeds to run')
self.add_output('pitch_out', val=np.zeros(n_ws_dlc11), units='deg', desc='pitch angles to run')
self.add_output('AEP', val=0.0, units='kW*h', desc='annual energy production reconstructed from the openfast simulations')
self.add_output('My_std', val=0.0, units='N*m', desc='standard deviation of blade root flap bending moment in out-of-plane direction')
self.add_output('flp1_std', val=0.0, units='deg', desc='standard deviation of trailing-edge flap angle')
self.add_output('rated_V', val=0.0, units='m/s', desc='rated wind speed')
self.add_output('rated_Omega', val=0.0, units='rpm', desc='rotor rotation speed at rated')
self.add_output('rated_pitch', val=0.0, units='deg', desc='pitch setting at rated')
self.add_output('rated_T', val=0.0, units='N', desc='rotor aerodynamic thrust at rated')
self.add_output('rated_Q', val=0.0, units='N*m', desc='rotor aerodynamic torque at rated')
self.add_output('loads_r', val=np.zeros(n_span), units='m', desc='radial positions along blade going toward tip')
self.add_output('loads_Px', val=np.zeros(n_span), units='N/m', desc='distributed loads in blade-aligned x-direction')
self.add_output('loads_Py', val=np.zeros(n_span), units='N/m', desc='distributed loads in blade-aligned y-direction')
self.add_output('loads_Pz', val=np.zeros(n_span), units='N/m', desc='distributed loads in blade-aligned z-direction')
self.add_output('loads_Omega', val=0.0, units='rpm', desc='rotor rotation speed')
self.add_output('loads_pitch', val=0.0, units='deg', desc='pitch angle')
self.add_output('loads_azimuth', val=0.0, units='deg', desc='azimuthal angle')
# Control outputs
self.add_output('rotor_overspeed', val=0.0, desc='Maximum percent overspeed of the rotor during all OpenFAST simulations') # is this over a set of sims?
self.add_output('max_nac_accel', val=0.0, units='m/s**2', desc='Maximum nacelle acceleration magnitude all OpenFAST simulations') # is this over a set of sims?
self.add_output('avg_pitch_travel', val=0.0, units='deg/s', desc='Average pitch travel') # is this over a set of sims?
self.add_output('pitch_duty_cycle', val=0.0, units='deg/s', desc='Average pitch travel') # is this over a set of sims?
# Blade outputs
self.add_output('max_TipDxc', val=0.0, units='m', desc='Maximum of channel TipDxc, i.e. out of plane tip deflection. For upwind rotors, the max value is tower the tower')
self.add_output('max_RootMyb', val=0.0, units='kN*m', desc='Maximum of the signals RootMyb1, RootMyb2, ... across all n blades representing the maximum blade root flapwise moment')
self.add_output('max_RootMyc', val=0.0, units='kN*m', desc='Maximum of the signals RootMyb1, RootMyb2, ... across all n blades representing the maximum blade root out of plane moment')
self.add_output('max_RootMzb', val=0.0, units='kN*m', desc='Maximum of the signals RootMzb1, RootMzb2, ... across all n blades representing the maximum blade root torsional moment')
self.add_output('DEL_RootMyb', val=0.0, units='kN*m', desc='damage equivalent load of blade root flap bending moment in out-of-plane direction')
self.add_output('max_aoa', val=np.zeros(n_span), units='deg', desc='maxima of the angles of attack distributed along blade span')
self.add_output('std_aoa', val=np.zeros(n_span), units='deg', desc='standard deviation of the angles of attack distributed along blade span')
self.add_output('mean_aoa', val=np.zeros(n_span), units='deg', desc='mean of the angles of attack distributed along blade span')
# Blade loads corresponding to maximum blade tip deflection
self.add_output('blade_maxTD_Mx', val=np.zeros(n_span), units='kN*m', desc='distributed moment around blade-aligned x-axis corresponding to maximum blade tip deflection')
self.add_output('blade_maxTD_My', val=np.zeros(n_span), units='kN*m', desc='distributed moment around blade-aligned y-axis corresponding to maximum blade tip deflection')
self.add_output('blade_maxTD_Fz', val=np.zeros(n_span), units='kN', desc='distributed force in blade-aligned z-direction corresponding to maximum blade tip deflection')
# Hub outputs
self.add_output('hub_Fxyz', val=np.zeros(3), units='kN', desc = 'Maximum hub forces in the non rotating frame')
self.add_output('hub_Mxyz', val=np.zeros(3), units='kN*m', desc = 'Maximum hub moments in the non rotating frame')
self.add_output('max_TwrBsMyt',val=0.0, units='kN*m', desc='maximum of tower base bending moment in fore-aft direction')
self.add_output('max_TwrBsMyt_ratio',val=0.0, desc='ratio of maximum of tower base bending moment in fore-aft direction to maximum allowable bending moment')
self.add_output('DEL_TwrBsMyt',val=0.0, units='kN*m', desc='damage equivalent load of tower base bending moment in fore-aft direction')
self.add_output('DEL_TwrBsMyt_ratio',val=0.0, desc='ratio of damage equivalent load of tower base bending moment in fore-aft direction to maximum allowable bending moment')
# Tower outputs
if not self.options['modeling_options']['Level3']['from_openfast']:
self.add_output('tower_maxMy_Fx', val=np.zeros(n_full_tow-1), units='kN', desc='distributed force in tower-aligned x-direction corresponding to maximum fore-aft moment at tower base')
self.add_output('tower_maxMy_Fy', val=np.zeros(n_full_tow-1), units='kN', desc='distributed force in tower-aligned y-direction corresponding to maximum fore-aft moment at tower base')
self.add_output('tower_maxMy_Fz', val=np.zeros(n_full_tow-1), units='kN', desc='distributed force in tower-aligned z-direction corresponding to maximum fore-aft moment at tower base')
self.add_output('tower_maxMy_Mx', val=np.zeros(n_full_tow-1), units='kN*m', desc='distributed moment around tower-aligned x-axis corresponding to maximum fore-aft moment at tower base')
self.add_output('tower_maxMy_My', val=np.zeros(n_full_tow-1), units='kN*m', desc='distributed moment around tower-aligned x-axis corresponding to maximum fore-aft moment at tower base')
self.add_output('tower_maxMy_Mz', val=np.zeros(n_full_tow-1), units='kN*m', desc='distributed moment around tower-aligned x-axis corresponding to maximum fore-aft moment at tower base')
# Monopile outputs
self.add_output('max_M1N1MKye',val=0.0, units='kN*m', desc='maximum of My moment of member 1 at node 1 (base of the monopile)')
self.add_output('monopile_maxMy_Fx', val=np.zeros(monlen_full), units='kN', desc='distributed force in monopile-aligned x-direction corresponding to max_M1N1MKye')
self.add_output('monopile_maxMy_Fy', val=np.zeros(monlen_full), units='kN', desc='distributed force in monopile-aligned y-direction corresponding to max_M1N1MKye')
self.add_output('monopile_maxMy_Fz', val=np.zeros(monlen_full), units='kN', desc='distributed force in monopile-aligned z-direction corresponding to max_M1N1MKye')
self.add_output('monopile_maxMy_Mx', val=np.zeros(monlen_full), units='kN*m', desc='distributed moment around tower-aligned x-axis corresponding to max_M1N1MKye')
self.add_output('monopile_maxMy_My', val=np.zeros(monlen_full), units='kN*m', desc='distributed moment around tower-aligned x-axis corresponding to max_M1N1MKye')
self.add_output('monopile_maxMy_Mz', val=np.zeros(monlen_full), units='kN*m', desc='distributed moment around tower-aligned x-axis corresponding to max_M1N1MKye')
# Floating outputs
self.add_output('Max_PtfmPitch', val=0.0, desc='Maximum platform pitch angle over a set of OpenFAST simulations')
self.add_output('Std_PtfmPitch', val=0.0, units='deg', desc='standard deviation of platform pitch angle')
self.add_output('Max_Offset', val=0.0, units='m', desc='Maximum distance in surge/sway direction')
# Fatigue output
self.add_output('damage_blade_root_sparU', val=0.0, desc="Miner's rule cumulative damage to upper spar cap at blade root")
self.add_output('damage_blade_root_sparL', val=0.0, desc="Miner's rule cumulative damage to lower spar cap at blade root")
self.add_output('damage_blade_maxc_teU', val=0.0, desc="Miner's rule cumulative damage to upper trailing edge at blade max chord")
self.add_output('damage_blade_maxc_teL', val=0.0, desc="Miner's rule cumulative damage to lower trailing edge at blade max chord")
self.add_output('damage_lss', val=0.0, desc="Miner's rule cumulative damage to low speed shaft at hub attachment")
self.add_output('damage_tower_base', val=0.0, desc="Miner's rule cumulative damage at tower base")
self.add_output('damage_monopile_base', val=0.0, desc="Miner's rule cumulative damage at monopile base")
# Simulation output
self.add_output('openfast_failed', val=0.0, desc="Numerical value for whether any openfast runs failed. 0 if false, 2 if true")
# Open loop to closed loop error
if self.options['modeling_options']['OL2CL']['flag']:
self.add_output('OL2CL_pitch', val=0.0, desc="Open loop to closed loop avarege error")
self.add_discrete_output('fst_vt_out', val={})
self.add_discrete_output('ts_out_dir', val={})
# Iteration counter for openfast calls. Initialize at -1 so 0 after first call
self.of_inumber = -1
self.sim_idx = -1
if modopt['Level2']['flag']:
if MPI:
rank = MPI.COMM_WORLD.Get_rank()
lin_pkl_dir = os.path.join(self.options['opt_options']['general']['folder_output'], 'lin', 'rank_{}'.format(rank))
if not os.path.exists(lin_pkl_dir):
os.makedirs(lin_pkl_dir, exist_ok=True)
self.lin_pkl_file_name = os.path.join(lin_pkl_dir, 'ABCD_matrices.pkl')
else:
lin_pkl_dir = os.path.join(self.options['opt_options']['general']['folder_output'], 'lin')
self.lin_pkl_file_name = os.path.join(lin_pkl_dir, 'ABCD_matrices.pkl')
self.ABCD_list = []
path = '.'.join(self.lin_pkl_file_name.split('.')[:-1])
os.makedirs(path, exist_ok=True)
with open(self.lin_pkl_file_name, 'wb') as handle:
pickle.dump(self.ABCD_list, handle)
self.lin_idx = 0
def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
modopt = self.options['modeling_options']
#print(impl.world_comm().rank, 'Rotor_fast','start')
sys.stdout.flush()
if modopt['Level2']['flag']:
self.sim_idx += 1
ABCD = {
'sim_idx' : self.sim_idx,
'A' : None,
'B' : None,
'C' : None,
'D' : None,
'x_ops':None,
'u_ops':None,
'y_ops':None,
'u_h':None,
'omega_rpm' : None,
'DescCntrlInpt' : None,
'DescStates' : None,
'DescOutput' : None,
'StateDerivOrder' : None,
'ind_fast_inps' : None,
'ind_fast_outs' : None,
'AEP':None
}
with open(self.lin_pkl_file_name, 'rb') as handle:
ABCD_list = pickle.load(handle)
ABCD_list.append(ABCD)
self.ABCD_list = ABCD_list
with open(self.lin_pkl_file_name, 'wb') as handle:
pickle.dump(ABCD_list, handle)
fst_vt = self.init_FAST_model()
if not modopt['Level3']['from_openfast']:
fst_vt = self.update_FAST_model(fst_vt, inputs, discrete_inputs)
else:
fast_reader = InputReader_OpenFAST()
fast_reader.FAST_InputFile = modopt['Level3']['openfast_file'] # FAST input file (ext=.fst)
if os.path.isabs(modopt['Level3']['openfast_dir']):
fast_reader.FAST_directory = modopt['Level3']['openfast_dir'] # Path to fst directory files
else:
fast_reader.FAST_directory = os.path.join(weis_dir, modopt['Level3']['openfast_dir'])
fast_reader.path2dll = modopt['General']['openfast_configuration']['path2dll'] # Path to dll file
fast_reader.execute()
fst_vt = fast_reader.fst_vt
fst_vt = self.load_FAST_model_opts(fst_vt)
# Fix TwrTI: WEIS modeling options have it as a single value...
if not isinstance(fst_vt['AeroDyn15']['TwrTI'],list):
fst_vt['AeroDyn15']['TwrTI'] = [fst_vt['AeroDyn15']['TwrTI']] * len(fst_vt['AeroDyn15']['TwrElev'])
# Fix AddF0: Should be a n x 1 array (list of lists):
fst_vt['HydroDyn']['AddF0'] = [[F0] for F0 in fst_vt['HydroDyn']['AddF0']]
if modopt['ROSCO']['flag']:
fst_vt['DISCON_in'] = modopt['General']['openfast_configuration']['fst_vt']['DISCON_in']
if self.model_only == True:
# Write input OF files, but do not run OF
self.write_FAST(fst_vt, discrete_outputs)
else:
# Write OF model and run
summary_stats, extreme_table, DELs, Damage, case_list, case_name, chan_time, dlc_generator = self.run_FAST(inputs, discrete_inputs, fst_vt)
# Set up linear turbine model
if modopt['Level2']['flag']:
try:
LinearTurbine = LinearTurbineModel(
self.FAST_runDirectory,
self.lin_case_name,
nlin=modopt['Level2']['linearization']['NLinTimes'],
reduceControls=True
)
except FileNotFoundError as e:
print('FileNotFoundError: {} {}'.format(e.strerror, e.filename))
return
# Save linearizations
print('Saving ABCD matrices!')
ABCD = {
'sim_idx' : self.sim_idx,
'A' : LinearTurbine.A_ops,
'B' : LinearTurbine.B_ops,
'C' : LinearTurbine.C_ops,
'D' : LinearTurbine.D_ops,
'x_ops':LinearTurbine.x_ops,
'u_ops':LinearTurbine.u_ops,
'y_ops':LinearTurbine.y_ops,
'u_h':LinearTurbine.u_h,
'omega_rpm' : LinearTurbine.omega_rpm,
'DescCntrlInpt' : LinearTurbine.DescCntrlInpt,
'DescStates' : LinearTurbine.DescStates,
'DescOutput' : LinearTurbine.DescOutput,
'StateDerivOrder' : LinearTurbine.StateDerivOrder,
'ind_fast_inps' : LinearTurbine.ind_fast_inps,
'ind_fast_outs' : LinearTurbine.ind_fast_outs,
}
with open(self.lin_pkl_file_name, 'rb') as handle:
ABCD_list = pickle.load(handle)
ABCD_list[self.sim_idx] = ABCD
with open(self.lin_pkl_file_name, 'wb') as handle:
pickle.dump(ABCD_list, handle)
lin_files = glob.glob(os.path.join(self.FAST_runDirectory, '*.lin'))
dest = os.path.join(self.FAST_runDirectory, f'copied_lin_files_{self.lin_idx}')
Path(dest).mkdir(parents=True, exist_ok=True)
for file in lin_files:
shutil.copy2(file, dest)
self.lin_idx += 1
# Shorten output names from linearization output to one like level3 openfast output
# This depends on how openfast sets up the linearization output names and may break if that is changed
OutList = [out_name.split()[1][:-1] for out_name in LinearTurbine.DescOutput]
OutOps = {}
for i_out, out in enumerate(OutList):
OutOps[out] = LinearTurbine.y_ops[i_out,:]
# save to yaml, might want in analysis outputs
FileTools.save_yaml(
self.FAST_runDirectory,
'OutOps.yaml',OutOps)
# Set up Level 2 disturbance (simulation or DTQP)
if modopt['Level2']['simulation']['flag'] or modopt['Level2']['DTQP']['flag']:
# Extract disturbance(s)
level2_disturbance = []
for case in case_list:
ts_file = TurbSimFile(case[('InflowWind','FileName_BTS')])
ts_file.compute_rot_avg(fst_vt['ElastoDyn']['TipRad'])
u_h = ts_file['rot_avg'][0,:]
tt = ts_file['t']
level2_disturbance.append({'Time':tt, 'Wind': u_h})
# Run linear simulation:
# Get case list, wind inputs should have already been generated
if modopt['Level2']['simulation']['flag']:
if modopt['Level2']['DTQP']['flag']:
raise Exception('Only DTQP or simulation flag can be set to true in Level2 modeling options')
# This is going to use the last discon_in file of the linearization set as the simulation file
# Currently fine because openfast is executed (or not executed if overwrite=False) after the file writing
if 'DLL_InFile' in self.fst_vt['ServoDyn']: # if using file inputs
discon_in_file = os.path.join(self.FAST_runDirectory, self.fst_vt['ServoDyn']['DLL_InFile'])
else: # if using fst_vt inputs from openfast_openmdao
discon_in_file = os.path.join(self.FAST_runDirectory, self.lin_case_name[0] + '_DISCON.IN')
lib_name = os.path.join(os.path.dirname(os.path.realpath(__file__)),'../../local/lib/libdiscon'+lib_ext)
ss = {}
et = {}
dl = {}
dam = {}
ct = []
for i_dist, dist in enumerate(level2_disturbance):
sim_name = 'l2_sim_{}'.format(i_dist)
controller_int = ROSCO_ci.ControllerInterface(
lib_name,
param_filename=discon_in_file,
DT=1/80, # modelling input?
sim_name = os.path.join(self.FAST_runDirectory,sim_name)
)
l2_out, _, P_op = LinearTurbine.solve(dist,Plot=False,controller=controller_int)
output = OpenFASTOutput.from_dict(l2_out, sim_name, magnitude_channels=self.magnitude_channels)
_name, _ss, _et, _dl, _dam = self.la._process_output(output)
ss[_name] = _ss
et[_name] = _et
dl[_name] = _dl
dam[_name] = _dam
ct.append(l2_out)
output.df.to_pickle(os.path.join(self.FAST_runDirectory,sim_name+'.p'))
summary_stats, extreme_table, DELs, Damage = self.la.post_process(ss, et, dl, dam)
# Overwrite timeseries with simulated data instead of saved linearization timeseries
chan_time = ct
elif modopt['Level2']['DTQP']['flag']:
summary_stats, extreme_table, DELs, Damage = dtqp_wrapper(
LinearTurbine,
level2_disturbance,
self.options['opt_options'],
self.options['modeling_options'],
self.fst_vt,
self.la,
self.magnitude_channels,
self.FAST_runDirectory
)
# TODO: pull chan_time out of here
# Post process regardless of level
self.post_process(summary_stats, extreme_table, DELs, Damage, case_list, dlc_generator, chan_time, inputs, discrete_inputs, outputs, discrete_outputs)
# Save AEP value to linear pickle file
if modopt['Level2']['flag']:
with open(self.lin_pkl_file_name, 'rb') as handle:
ABCD_list = pickle.load(handle)
ABCD_list[self.sim_idx]['AEP'] = outputs['AEP']
with open(self.lin_pkl_file_name, 'wb') as handle:
pickle.dump(ABCD_list, handle)
# delete run directory. not recommended for most cases, use for large parallelization problems where disk storage will otherwise fill up
if self.clean_FAST_directory:
try:
shutil.rmtree(self.FAST_runDirectory)
except:
print('Failed to delete directory: %s'%self.FAST_runDirectory)
def init_FAST_model(self):
modopt = self.options['modeling_options']
fst_vt = modopt['General']['openfast_configuration']['fst_vt']
# Main .fst file`
fst_vt['Fst'] = {}
fst_vt['ElastoDyn'] = {}
fst_vt['ElastoDynBlade'] = {}
fst_vt['ElastoDynTower'] = {}
fst_vt['AeroDyn15'] = {}
fst_vt['AeroDynBlade'] = {}
fst_vt['ServoDyn'] = {}
fst_vt['InflowWind'] = {}
fst_vt['SubDyn'] = {}
fst_vt['HydroDyn'] = {}
fst_vt['MoorDyn'] = {}
fst_vt['MAP'] = {}
# List of structural controllers
fst_vt['TStC'] = {}; fst_vt['TStC'] = []
fst_vt['SStC'] = {}; fst_vt['SStC'] = []
fst_vt = self.load_FAST_model_opts(fst_vt)
return fst_vt
def load_FAST_model_opts(self,fst_vt):
modeling_options = self.options['modeling_options']
for key in modeling_options['Level3']['simulation']:
fst_vt['Fst'][key] = modeling_options['Level3']['simulation'][key]
for key in modeling_options['Level3']['ElastoDyn']:
fst_vt['ElastoDyn'][key] = modeling_options['Level3']['ElastoDyn'][key]
for key in modeling_options['Level3']['ElastoDynBlade']:
fst_vt['ElastoDynBlade'][key] = modeling_options['Level3']['ElastoDynBlade'][key]
for key in modeling_options['Level3']['ElastoDynTower']:
fst_vt['ElastoDynTower'][key] = modeling_options['Level3']['ElastoDynTower'][key]
for key in modeling_options['Level3']['AeroDyn']:
fst_vt['AeroDyn15'][key] = copy.copy(modeling_options['Level3']['AeroDyn'][key])
for key in modeling_options['Level3']['InflowWind']:
fst_vt['InflowWind'][key] = modeling_options['Level3']['InflowWind'][key]
for key in modeling_options['Level3']['ServoDyn']:
fst_vt['ServoDyn'][key] = modeling_options['Level3']['ServoDyn'][key]
for key in modeling_options['Level3']['SubDyn']:
fst_vt['SubDyn'][key] = modeling_options['Level3']['SubDyn'][key]
for key in modeling_options['Level3']['HydroDyn']:
fst_vt['HydroDyn'][key] = modeling_options['Level3']['HydroDyn'][key]
for key in modeling_options['Level3']['MoorDyn']:
fst_vt['MoorDyn'][key] = modeling_options['Level3']['MoorDyn'][key]
for key1 in modeling_options['Level3']['outlist']:
for key2 in modeling_options['Level3']['outlist'][key1]:
fst_vt['outlist'][key1][key2] = modeling_options['Level3']['outlist'][key1][key2]
fst_vt['ServoDyn']['DLL_FileName'] = modeling_options['General']['openfast_configuration']['path2dll']
if fst_vt['AeroDyn15']['IndToler'] == 0.:
fst_vt['AeroDyn15']['IndToler'] = 'Default'
if fst_vt['AeroDyn15']['DTAero'] == 0.:
fst_vt['AeroDyn15']['DTAero'] = 'Default'
if fst_vt['AeroDyn15']['OLAF']['DTfvw'] == 0.:
fst_vt['AeroDyn15']['OLAF']['DTfvw'] = 'Default'
if fst_vt['ElastoDyn']['DT'] == 0.:
fst_vt['ElastoDyn']['DT'] = 'Default'
return fst_vt
def update_FAST_model(self, fst_vt, inputs, discrete_inputs):
modopt = self.options['modeling_options']
# Update fst_vt nested dictionary with data coming from WISDEM
# Update ElastoDyn
fst_vt['ElastoDyn']['NumBl'] = self.n_blades
fst_vt['ElastoDyn']['TipRad'] = inputs['Rtip'][0]
fst_vt['ElastoDyn']['HubRad'] = inputs['Rhub'][0]
if discrete_inputs['rotor_orientation'] == 'upwind':
k = -1.
else:
k = 1
fst_vt['ElastoDyn']['PreCone(1)'] = k*inputs['cone'][0]
fst_vt['ElastoDyn']['PreCone(2)'] = k*inputs['cone'][0]
fst_vt['ElastoDyn']['PreCone(3)'] = k*inputs['cone'][0]
fst_vt['ElastoDyn']['ShftTilt'] = k*inputs['tilt'][0]
fst_vt['ElastoDyn']['OverHang'] = k*inputs['overhang'][0] / np.cos(np.deg2rad(inputs['tilt'][0])) # OpenFAST defines the overhang tilted (!)
fst_vt['ElastoDyn']['GBoxEff'] = inputs['gearbox_efficiency'][0] * 100.
fst_vt['ElastoDyn']['GBRatio'] = inputs['gearbox_ratio'][0]
# Update ServoDyn
fst_vt['ServoDyn']['GenEff'] = float(inputs['generator_efficiency']/inputs['gearbox_efficiency']) * 100.
fst_vt['ServoDyn']['PitManRat(1)'] = float(inputs['max_pitch_rate'])
fst_vt['ServoDyn']['PitManRat(2)'] = float(inputs['max_pitch_rate'])
fst_vt['ServoDyn']['PitManRat(3)'] = float(inputs['max_pitch_rate'])
# Update ServoDyn
fst_vt['ServoDyn']['GenEff'] = float(inputs['generator_efficiency']/inputs['gearbox_efficiency']) * 100.
fst_vt['ServoDyn']['PitManRat(1)'] = float(inputs['max_pitch_rate'])
fst_vt['ServoDyn']['PitManRat(2)'] = float(inputs['max_pitch_rate'])
fst_vt['ServoDyn']['PitManRat(3)'] = float(inputs['max_pitch_rate'])
# Masses and inertias from DriveSE
fst_vt['ElastoDyn']['HubMass'] = inputs['hub_system_mass'][0]
fst_vt['ElastoDyn']['HubIner'] = inputs['hub_system_I'][0]
fst_vt['ElastoDyn']['HubCM'] = inputs['hub_system_cm'][0] # k*inputs['overhang'][0] - inputs['hub_system_cm'][0], but we need to solve the circular dependency in DriveSE first
fst_vt['ElastoDyn']['NacMass'] = inputs['above_yaw_mass'][0]
fst_vt['ElastoDyn']['YawBrMass'] = inputs['yaw_mass'][0]
# Advice from R. Bergua, add 1/3 the tower yaw inertia here because it is activated as a lumped modal mass
fst_vt['ElastoDyn']['NacYIner'] = inputs['nacelle_I_TT'][2] + inputs['tower_I_base'][2]/3.0
fst_vt['ElastoDyn']['NacCMxn'] = -k*inputs['nacelle_cm'][0]
fst_vt['ElastoDyn']['NacCMyn'] = inputs['nacelle_cm'][1]
fst_vt['ElastoDyn']['NacCMzn'] = inputs['nacelle_cm'][2]
tower_top_height = float(inputs['hub_height']) - float(inputs['distance_tt_hub']) # Height of tower above ground level [onshore] or MSL [offshore] (meters)
# The Twr2Shft is just the difference between hub height, tower top height, and sin(tilt)*overhang
fst_vt['ElastoDyn']['Twr2Shft'] = float(inputs['hub_height']) - tower_top_height - abs(fst_vt['ElastoDyn']['OverHang'])*np.sin(np.deg2rad(inputs['tilt'][0]))
fst_vt['ElastoDyn']['GenIner'] = float(inputs['GenIner'])
# Mass and inertia inputs
fst_vt['ElastoDyn']['TipMass(1)'] = 0.
fst_vt['ElastoDyn']['TipMass(2)'] = 0.
fst_vt['ElastoDyn']['TipMass(3)'] = 0.
tower_base_height = max(float(inputs['tower_base_height']), float(inputs["platform_total_center_of_mass"][2]))
fst_vt['ElastoDyn']['TowerBsHt'] = tower_base_height # Height of tower base above ground level [onshore] or MSL [offshore] (meters)
fst_vt['ElastoDyn']['TowerHt'] = tower_top_height
# TODO: There is some confusion on PtfmRefzt
# DZ: based on the openfast r-tests:
# if this is floating, the z ref. point is 0. Is this the reference that platform_total_center_of_mass is relative to?
# if fixed bottom, it's the tower base height.
if modopt['flags']['floating']:
fst_vt['ElastoDyn']['PtfmMass'] = float(inputs["platform_mass"])
fst_vt['ElastoDyn']['PtfmRIner'] = float(inputs["platform_I_total"][0])
fst_vt['ElastoDyn']['PtfmPIner'] = float(inputs["platform_I_total"][1])
fst_vt['ElastoDyn']['PtfmYIner'] = float(inputs["platform_I_total"][2])
fst_vt['ElastoDyn']['PtfmCMxt'] = float(inputs["platform_total_center_of_mass"][0])
fst_vt['ElastoDyn']['PtfmCMyt'] = float(inputs["platform_total_center_of_mass"][1])
fst_vt['ElastoDyn']['PtfmCMzt'] = float(inputs["platform_total_center_of_mass"][2])
fst_vt['ElastoDyn']['PtfmRefzt'] = 0. # Vertical distance from the ground level [onshore] or MSL [offshore] to the platform reference point (meters)
else:
# Ptfm* can capture the transition piece for fixed-bottom, but we are doing that in subdyn, so only worry about getting height right
fst_vt['ElastoDyn']['PtfmMass'] = 0.
fst_vt['ElastoDyn']['PtfmRIner'] = 0.
fst_vt['ElastoDyn']['PtfmPIner'] = 0.
# Advice from R. Bergua- Use a dummy quantity (at least 1e4) here when have fixed-bottom support
fst_vt['ElastoDyn']['PtfmYIner'] = 1e8 if modopt['flags']['offshore'] else 0.0
fst_vt['ElastoDyn']['PtfmCMxt'] = 0.
fst_vt['ElastoDyn']['PtfmCMyt'] = 0.
fst_vt['ElastoDyn']['PtfmCMzt'] = float(inputs['tower_base_height'])
fst_vt['ElastoDyn']['PtfmRefzt'] = tower_base_height # Vertical distance from the ground level [onshore] or MSL [offshore] to the platform reference point (meters)
# Drivetrain inputs
fst_vt['ElastoDyn']['DTTorSpr'] = float(inputs['drivetrain_spring_constant'])
fst_vt['ElastoDyn']['DTTorDmp'] = float(inputs['drivetrain_damping_coefficient'])
# Update Inflowwind
fst_vt['InflowWind']['RefHt'] = float(inputs['hub_height'])
fst_vt['InflowWind']['RefHt_Uni'] = float(inputs['hub_height'])
fst_vt['InflowWind']['PLexp'] = float(inputs['shearExp'])
if fst_vt['InflowWind']['NWindVel'] == 1:
fst_vt['InflowWind']['WindVxiList'] = 0.
fst_vt['InflowWind']['WindVyiList'] = 0.
fst_vt['InflowWind']['WindVziList'] = float(inputs['hub_height'])
else:
raise Exception('The code only supports InflowWind NWindVel == 1')
# Update AeroDyn Tower Input File starting one station above ground to avoid error because the wind grid hits the ground
twr_elev = inputs['tower_z']
twr_d = inputs['tower_outer_diameter']
twr_index = np.argmin(abs(twr_elev - np.maximum(1.0, tower_base_height)))
cd_index = 0
if twr_elev[twr_index] <= 1.:
twr_index += 1
cd_index += 1
fst_vt['AeroDyn15']['NumTwrNds'] = len(twr_elev[twr_index:])
fst_vt['AeroDyn15']['TwrElev'] = twr_elev[twr_index:]
fst_vt['AeroDyn15']['TwrDiam'] = twr_d[twr_index:]
fst_vt['AeroDyn15']['TwrCd'] = inputs['tower_cd'][cd_index:]
fst_vt['AeroDyn15']['TwrTI'] = np.ones(len(twr_elev[twr_index:])) * fst_vt['AeroDyn15']['TwrTI']
fst_vt['AeroDyn15']['tau1_const'] = 0.24 * float(inputs['Rtip']) # estimated using a=0.3 and U0=7.5
z_tow = twr_elev
z_sec, _ = util.nodal2sectional(z_tow)
sec_loc = (z_sec - z_sec[0]) / (z_sec[-1] - z_sec[0])
fst_vt['ElastoDynTower']['NTwInpSt'] = len(sec_loc)
fst_vt['ElastoDynTower']['HtFract'] = sec_loc
fst_vt['ElastoDynTower']['TMassDen'] = inputs['mass_den']
fst_vt['ElastoDynTower']['TwFAStif'] = inputs['foreaft_stff']
fst_vt['ElastoDynTower']['TwSSStif'] = inputs['sideside_stff']
fst_vt['ElastoDynTower']['TwFAM1Sh'] = inputs['fore_aft_modes'][0, :] / sum(inputs['fore_aft_modes'][0, :])
fst_vt['ElastoDynTower']['TwFAM2Sh'] = inputs['fore_aft_modes'][1, :] / sum(inputs['fore_aft_modes'][1, :])
fst_vt['ElastoDynTower']['TwSSM1Sh'] = inputs['side_side_modes'][0, :] / sum(inputs['side_side_modes'][0, :])
fst_vt['ElastoDynTower']['TwSSM2Sh'] = inputs['side_side_modes'][1, :] / sum(inputs['side_side_modes'][1, :])
# Calculate yaw stiffness of tower (springs in series) and use in servodyn as yaw spring constant
k_tow_tor = inputs['tor_stff'] / np.diff(inputs['tower_z'])
k_tow_tor = 1.0/np.sum(1.0/k_tow_tor)
# R. Bergua's suggestion to set the stiffness to the tower torsional stiffness and the
# damping to the frequency of the first tower torsional mode- easier than getting the yaw inertia right
damp_ratio = 0.01
f_torsion = float(inputs['tor_freq'])
fst_vt['ServoDyn']['YawSpr'] = k_tow_tor
if f_torsion > 0.0:
fst_vt['ServoDyn']['YawDamp'] = damp_ratio * k_tow_tor / np.pi / f_torsion
else:
fst_vt['ServoDyn']['YawDamp'] = 2 * damp_ratio * np.sqrt(k_tow_tor * inputs['rna_I_TT'][2])
# Update ElastoDyn Blade Input File
fst_vt['ElastoDynBlade']['NBlInpSt'] = len(inputs['r'])
fst_vt['ElastoDynBlade']['BlFract'] = (inputs['r']-inputs['Rhub'])/(inputs['Rtip']-inputs['Rhub'])
fst_vt['ElastoDynBlade']['BlFract'][0] = 0.
fst_vt['ElastoDynBlade']['BlFract'][-1]= 1.
fst_vt['ElastoDynBlade']['PitchAxis'] = inputs['le_location']
fst_vt['ElastoDynBlade']['StrcTwst'] = inputs['theta'] # to do: structural twist is not nessessarily (nor likely to be) the same as aero twist
fst_vt['ElastoDynBlade']['BMassDen'] = inputs['beam:rhoA']
fst_vt['ElastoDynBlade']['FlpStff'] = inputs['beam:EIyy']
fst_vt['ElastoDynBlade']['EdgStff'] = inputs['beam:EIxx']
fst_vt['ElastoDynBlade']['BldFl1Sh'] = np.zeros(5)
fst_vt['ElastoDynBlade']['BldFl2Sh'] = np.zeros(5)
fst_vt['ElastoDynBlade']['BldEdgSh'] = np.zeros(5)
for i in range(5):
fst_vt['ElastoDynBlade']['BldFl1Sh'][i] = inputs['flap_mode_shapes'][0,i] / sum(inputs['flap_mode_shapes'][0,:])
fst_vt['ElastoDynBlade']['BldFl2Sh'][i] = inputs['flap_mode_shapes'][1,i] / sum(inputs['flap_mode_shapes'][1,:])
fst_vt['ElastoDynBlade']['BldEdgSh'][i] = inputs['edge_mode_shapes'][0,i] / sum(inputs['edge_mode_shapes'][0,:])