Merge pull request #110 from WISDEM/run_olaf

Run olaf

examples/01_aeroelasticse/OpenFAST_models/IEA-15-240-RWT/IEA-15-240-RWT/IEA-15-240-RWT_BeamDyn.dat

@@ -75,7 +75,7 @@ True          RotStates       - Orient states in the rotating frame during linea
 ---------------------- MESH PARAMETER ------------------------------------------
           5   order_elem     - Order of interpolation (basis) function (-)
 ---------------------- MATERIAL PARAMETER --------------------------------------
-"IEA 15MW Offshore Reference Turbine, with taped chord tip design_BeamDyn_Blade.dat"    BldFile - Name of file containing properties for blade (quoted string)
+"../IEA-15-240-RWT/IEA-15-240-RWT_BeamDyn_blade.dat"    BldFile - Name of file containing properties for blade (quoted string)
 ---------------------- PITCH ACTUATOR PARAMETERS -------------------------------
 False         UsePitchAct - Whether a pitch actuator should be used (flag)
         200   PitchJ      - Pitch actuator inertia (kg-m^2) [used only when UsePitchAct is true]

examples/01_aeroelasticse/run_OLAF.py

@@ -58,7 +58,7 @@
 case_inputs[("ElastoDyn","BlPitch1")]   = {'vals': pitch_init, 'group': 1}
 case_inputs[("ElastoDyn","BlPitch2")]   = case_inputs[("ElastoDyn","BlPitch1")]
 case_inputs[("ElastoDyn","BlPitch3")]   = case_inputs[("ElastoDyn","BlPitch1")]
-dt_wanted, nNWPanel, nFWPanel, nFWPanelFree = OLAFParams(omega_init)
+dt_wanted, tMax, nNWPanel, nFWPanel, nFWPanelFree = OLAFParams(omega_init)
 dt_olaf = np.zeros_like(dt_wanted)
 dt = case_inputs[("Fst","DT")]["vals"]
 n_dt = dt_wanted / dt

weis/aeroelasticse/CaseGen_IEC.py

@@ -53,6 +53,7 @@ def __init__(self):
         self.TStart                      = 30.
         self.uniqueSeeds                 = False
         self.uniqueWaveSeeds             = False
+        self.grid_center_over_hh         = 1.      # Ratio between turbsim grid center and hub height
 
         self.debug_level                 = 2
         self.parallel_windfile_gen       = False
@@ -134,6 +135,7 @@ def execute(self, case_inputs={}):
             iecwind.dir_change       = self.transient_dir_change
             iecwind.shear_orient     = self.transient_shear_orientation
             iecwind.z_hub            = self.z_hub
+            iecwind.z_grid_center    = self.grid_center_over_hh * iecwind.z_hub
             iecwind.D                = self.D
             iecwind.PLExp            = alpha
 

weis/aeroelasticse/FAST_reader.py

@@ -428,8 +428,9 @@ def execute(self):
         if self.fst_vt['Fst']['CompMooring'] == 3: # MoorDyn implimented
             self.read_MoorDyn()
 
-        if self.fst_vt['Fst']['CompElast'] == 2: # BeamDyn read assumes all 3 blades are the same
-            self.read_BeamDyn()
+        bd_file = os.path.join(self.FAST_directory, self.fst_vt['Fst']['BDBldFile(1)'])
+        if os.path.isfile(bd_file):
+            self.read_BeamDyn(bd_file)
 
     def read_MainInput(self):
         # Main FAST v8.16-v8.17 Input File
@@ -705,12 +706,9 @@ def read_ElastoDyn(self):
 
         f.close()
 
-    def read_BeamDyn(self):
+    def read_BeamDyn(self, bd_file):
         # BeamDyn Input File
-
-        bd_file = os.path.join(self.FAST_directory, self.fst_vt['Fst']['BDBldFile(1)'])
         f = open(bd_file)
-
         f.readline()
         f.readline()
         f.readline()
@@ -1126,10 +1124,9 @@ def read_AeroDyn15(self):
         self.read_AeroDyn15Blade()
         self.read_AeroDyn15Polar()
         self.read_AeroDyn15Coord()
-        if self.fst_vt['AeroDyn15']['WakeMod'] == 3:
-            if self.fst_vt['AeroDyn15']['AFAeroMod'] == 2:
-                raise Exception('OLAF is called with unsteady airfoil aerodynamics, but OLAF currently only supports AFAeroMod == 1')
-            self.read_AeroDyn15OLAF()
+        olaf_filename = os.path.join(self.FAST_directory, self.fst_vt['AeroDyn15']['OLAFInputFileName'])
+        if os.path.isfile(olaf_filename):
+            self.read_AeroDyn15OLAF(olaf_filename)
 
     def read_AeroDyn15Blade(self):
         # AeroDyn v5.00 Blade Definition File
@@ -1274,10 +1271,9 @@ def read_AeroDyn15Coord(self):
 
                 f.close()
 
-    def read_AeroDyn15OLAF(self):
+    def read_AeroDyn15OLAF(self, olaf_filename):
 
         self.fst_vt['AeroDyn15']['OLAF'] = {}
-        olaf_filename = os.path.join(self.FAST_directory, self.fst_vt['AeroDyn15']['OLAFInputFileName'])
         f = open(olaf_filename)
         f.readline()
         f.readline()
@@ -1323,7 +1319,6 @@ def read_AeroDyn15OLAF(self):
         f.readline()
         f.close()
 
-
     def read_ServoDyn(self):
         # ServoDyn v1.05 Input File
         # Currently no differences between FASTv8.16 and OpenFAST.

weis/aeroelasticse/FAST_writer.py

@@ -878,6 +878,8 @@ def write_AeroDyn15(self):
             self.write_AeroDyn15Coord()
 
         if self.fst_vt['AeroDyn15']['WakeMod'] == 3:
+            if self.fst_vt['AeroDyn15']['AFAeroMod'] == 2:
+                raise Exception('OLAF is called with unsteady airfoil aerodynamics, but OLAF currently only supports AFAeroMod == 1')
             self.write_OLAF()
 
         # Generate AeroDyn v15.03 input file

weis/aeroelasticse/openmdao_openfast.py

@@ -628,7 +628,7 @@ def update_FAST_model(self, fst_vt, inputs, discrete_inputs):
 
         # Update OLAF
         if fst_vt['AeroDyn15']['WakeMod'] == 3:
-            _, nNWPanel, nFWPanel, nFWPanelFree = OLAFParams(fst_vt['ElastoDyn']['RotSpeed'])
+            _, _, nNWPanel, nFWPanel, nFWPanelFree = OLAFParams(fst_vt['ElastoDyn']['RotSpeed'])
             fst_vt['AeroDyn15']['OLAF']['nNWPanel'] = nNWPanel
             fst_vt['AeroDyn15']['OLAF']['nFWPanel'] = nFWPanel
             fst_vt['AeroDyn15']['OLAF']['nFWPanelFree'] = nFWPanelFree
@@ -1152,11 +1152,14 @@ def DLC_creation_IEC(self, inputs, discrete_inputs, fst_vt, powercurve=False):
         iec.init_cond[("ElastoDyn","BlPitch2")]        = iec.init_cond[("ElastoDyn","BlPitch1")]
         iec.init_cond[("ElastoDyn","BlPitch3")]        = iec.init_cond[("ElastoDyn","BlPitch1")]
 
-        # Set DT according to OLAF guidelines
+        # If running OLAF...
         if fst_vt['AeroDyn15']['WakeMod'] == 3:
-            dt_wanted, _, _, _ = OLAFParams(inputs['Omega_init'])
+            # Set DT according to OLAF guidelines
+            dt_wanted, _, _, _, _ = OLAFParams(inputs['Omega_init'])
             iec.init_cond[("Fst","DT")]        = {'U':inputs['U_init']}
             iec.init_cond[("Fst","DT")]['val'] = dt_wanted
+            # Raise the center of the grid 50% above hub height because the wake will expand
+            iec.grid_center_over_hh = 1.5
 
         # Todo: need a way to handle Metocean conditions for Offshore
         # if offshore:
@@ -1273,7 +1276,7 @@ def DLC_creation_powercurve(self, inputs, discrete_inputs, fst_vt):
 
                 # Set DT according to OLAF guidelines
                 if fst_vt['AeroDyn15']['WakeMod'] == 3:
-                    dt_wanted, _, _, _ = OLAFParams(inputs['Omega_init'])
+                    dt_wanted, _, _, _, _ = OLAFParams(inputs['Omega_init'])
                     case_inputs[("Fst","DT")]               = {'vals':dt_wanted, 'group':1}
 
                 case_list, case_name = CaseGen_General(case_inputs, self.FAST_runDirectory, self.FAST_namingOut + '_powercurve')
@@ -1994,36 +1997,49 @@ def writeCpsurfaces(self, inputs):
 def RayleighCDF(x, xbar=10.):
     return 1.0 - np.exp(-np.pi/4.0*(x/xbar)**2)
 
-def OLAFParams(omega_rpm, deltaPsiDeg=6, nNWrot=2, nFWrot=10, nFWrotFree=3, nPerAzimuth=None):
+def OLAFParams(omega_rpm, deltaPsiDeg=6, nNWrot=2, nFWrot=10, nFWrotFree=3, nPerRot=None, totalRot=None, show=False):
     """
     Computes recommended time step and wake length based on the rotational speed in RPM
-
+ 
     INPUTS:
      - omega_rpm: rotational speed in RPM
      - deltaPsiDeg : azimuthal discretization in deg
      - nNWrot : number of near wake rotations
      - nFWrot : total number of far wake rotations
      - nFWrotFree : number of far wake rotations that are free
-
-        deltaPsiDeg  -  nPerAzimuth
-             5            72
-             6            60
-             7            51.5
-             8            45
+ 
+        deltaPsiDeg  -  nPerRot
+             5            72   
+             6            60   
+             7            51.5 
+             8            45   
     """
+    omega_rpm = np.asarray(omega_rpm)
     omega = omega_rpm*2*np.pi/60
     T = 2*np.pi/omega
-    if nPerAzimuth is not None:
-        dt_wanted    = np.around(T/nPerAzimuth,2)
+    if nPerRot is not None:
+        dt_wanted    = np.around(T/nPerRot,4)
     else:
-        dt_wanted    = np.around(deltaPsiDeg/(6*omega_rpm), 4)
-        nPerAzimuth = int(2*np.pi /(deltaPsiDeg*np.pi/180))
-
-    nNWPanel     = nNWrot*nPerAzimuth
-    nFWPanel     = nFWrot*nPerAzimuth
-    nFWPanelFree = nFWrotFree*nPerAzimuth
-
-    return dt_wanted, nNWPanel, nFWPanel, nFWPanelFree
+        dt_wanted    = np.around(deltaPsiDeg/(6*omega_rpm),4)
+        nPerRot = int(2*np.pi /(deltaPsiDeg*np.pi/180))
+
+    nNWPanel     = nNWrot*nPerRot
+    nFWPanel     = nFWrot*nPerRot
+    nFWPanelFree = nFWrotFree*nPerRot
+
+    if totalRot is None:
+        totalRot = (nNWrot + nFWrot)*3 # going three-times through the entire wake
+
+    tMax = dt_wanted*nPerRot*totalRot
+
+    if show:
+        print(dt_wanted              , '  dt')
+        print(int      (nNWPanel    ), '  nNWPanel          ({} rotations)'.format(nNWrot))
+        print(int      (nFWPanel    ), '  FarWakeLength     ({} rotations)'.format(nFWrot))
+        print(int      (nFWPanelFree), '  FreeFarWakeLength ({} rotations)'.format(nFWrotFree))
+        print(tMax              , '  Tmax ({} rotations)'.format(totalRot))
+
+    return dt_wanted, tMax, nNWPanel, nFWPanel, nFWPanelFree
 
 class ModesElastoDyn(ExplicitComponent):
     """

weis/aeroelasticse/pyIECWind.py

@@ -432,6 +432,7 @@ def __init__(self):
         self.seed             = np.random.uniform(1, 1e8)
         self.Turbulence_Class = 'B'  # IEC Turbulance Class
         self.z_hub            = 90.  # wind turbine hub height (m)
+        self.z_grid_center    = 90   # height of grid center from the ground (m)
         self.D                = 126. # rotor diameter (m)
         self.PLExp            = 0.2
         self.AnalysisTime     = 720.
@@ -443,9 +444,9 @@ def setup(self):
         turbsim_vt.runtime_options.RandSeed1  = self.seed
         turbsim_vt.runtime_options.WrADTWR    = False
         turbsim_vt.tmspecs.AnalysisTime       = self.AnalysisTime
-        turbsim_vt.tmspecs.HubHt              = self.z_hub # wind grid centered at hub height
-        turbsim_vt.tmspecs.GridHeight         = (self.z_hub - 1.) * 2.0 # wind grid stops 1 meter above the ground
-        turbsim_vt.tmspecs.GridWidth          = (self.z_hub - 1.) * 2.0 # squared wind grid
+        turbsim_vt.tmspecs.HubHt              = self.z_grid_center # wind grid may be centered at hub height, but it may also go higher. RefHt stays at hub height
+        turbsim_vt.tmspecs.GridHeight         = (self.z_grid_center - 1.) * 2.0 # wind grid stops 1 meter above the ground
+        turbsim_vt.tmspecs.GridWidth          = (self.z_grid_center - 1.) * 2.0 # squared wind grid
         turbsim_vt.tmspecs.NumGrid_Z          = 25
         turbsim_vt.tmspecs.NumGrid_Y          = 25
         turbsim_vt.tmspecs.HFlowAng           = 0.0