Merge pull request #113 from IEAWindSystems/perfdata

Improved rotor performance data

WISDEM/generateTables.py

@@ -76,6 +76,57 @@ def vabs_load(fname):
     return A, r
 
 
+def read_rotor_performance_yaml(fname):
+    data = load_yaml(fname)
+
+    cases = data['cases']
+    n_ws = len(cases)
+
+    col_rename = {'wind_speed':'Wind Speed [m/s]',
+                  'rotor_speed':'Rotor Speed [rpm]',
+                  'pitch':'Pitch [deg]',
+                  'tip_speed_ratio':'TSR',
+                  'mechanical_power':'Mech Power [kW]',
+                  'electrical_power':'Elec Power [kW]',
+                  'generator_torque':'Generator Torque [kNm]',
+                  'rotor_thrust':'Rotor Thrust [kN]',
+                  'rotor_torque':'Rotor Torque [kNm]'}
+
+    # Read each run configuration
+    configuration_keys = list(data["cases"][0]["configuration"].keys())
+    col_labels = configuration_keys[:]
+    config_data = np.zeros((n_ws, len(configuration_keys)))
+    for ik, k in enumerate(configuration_keys):
+        config_data[:,ik] = [c["configuration"][k] for c in cases]
+
+    # Read steady state integrated values
+    integrated_keys = list(data["cases"][0]["outputs"]["integrated"].keys())
+    col_labels += integrated_keys
+    integrated_data = np.zeros((n_ws, len(integrated_keys)))
+    for ik, k in enumerate(integrated_keys):
+        coeff = 1.0
+        if k.lower().find('power') >= 0 or k.lower().find('torque') >= 0:
+            coeff = 1e-3
+        integrated_data[:,ik] = [coeff * c["outputs"]["integrated"][k] for c in cases]
+
+    # Create DataFrame for data so far
+    rotDF = pd.DataFrame(np.hstack((config_data, integrated_data)), columns=col_labels)
+    rotDF.rename(columns = col_rename, inplace=True)
+
+    # Read steady state distributed data along span
+    rotDF_dist = []
+    distributed_keys = list(data["cases"][0]["outputs"]["distributed"].keys())
+    for ic, c in enumerate(cases):
+        igrid = c["outputs"]["distributed"][distributed_keys[0]]['grid']
+        distributed_data = np.zeros((len(igrid), 1+len(distributed_keys)))
+        distributed_data[:,0] = igrid
+        for ik, k in enumerate(distributed_keys):
+            distributed_data[:,ik+1] = c["outputs"]["distributed"][k]['values']
+        rotDF_dist.append( pd.DataFrame(distributed_data, columns=['Spanwise grid']+distributed_keys) )
+
+    return rotDF, rotDF_dist
+
+
 class RWT_Tabular(object):
     def __init__(self, finput, towDF=None, rotDF=None, layerDF=None, nacDF=None, overview=None):
 
@@ -88,7 +139,7 @@ def __init__(self, finput, towDF=None, rotDF=None, layerDF=None, nacDF=None, ove
         froot, _      = os.path.splitext( fname )
         self.fout     = folder_output + os.sep + froot + '_tabular.xlsx'
 
-        # If provided, store blade, tower, and rotor data
+        # If provided, store blade, tower, and rotor data from WISDEM
         self.towDF   = towDF
         self.rotDF   = rotDF
         self.layerDF = layerDF
@@ -812,49 +863,63 @@ def write_materials(self):
 
 
     def write_rotor_performance(self):
-        ws = self.wb.create_sheet(title = 'Rotor Performance')
-
-        # Use OpenFAST output if it is available
-        foutput = '..'+os.sep+'OpenFAST'+os.sep+'outputs'+os.sep+'IEA-22-280-RWT_steady.yaml'
-        if os.path.exists(foutput):
-            fastout = load_yaml(foutput)
-
-            rotmat = np.c_[fastout['Wind1VelX']['mean'],
-                           fastout['BldPitch1']['mean'],
-                           1e-3*np.array(fastout['GenPwr']['mean']),
-                           fastout['RotSpeed']['mean'],
-                           120.0*np.array(fastout['RotSpeed']['mean'])*2*np.pi/60.0,
-                           1e-3*np.array(fastout['RotThrust']['mean']),
-                           1e-3*np.array(fastout['RotTorq']['mean']),
-                           1e-3*np.array(fastout['GenTq']['mean']),
-                           fastout['RtAeroCp']['mean'],
-                           fastout['RtAeroCt']['mean']]
-            cols = ['Wind [m/s]',
-                    'Pitch [deg]',
-                    'Power [MW]',
-                    'Rotor Speed [rpm]',
-                    'Tip Speed [m/s]',
-                    'Rotor Thrust [MN]',
-                    'Rotor Torque [MNm]',
-                    'Generator Torque [MNm]',
-                    'Rotor Cp [-]',
-                    'Rotor Ct [-]']
-            myDF = pd.DataFrame(rotmat, columns=cols)
+        if self.rotDF is not None:
+            # WISDEM Rotor Performance
+            ws = self.wb.create_sheet(title = 'Rotor Performance - WISDEM')
+
+            # Write to the file
+            for r in dataframe_to_rows(self.rotDF, index=False, header=True):
+                ws.append(r)
+
+            # Header row style formatting
+            for cell in ws["1:1"]:
+                cell.style = 'Headline 2'
+
+        def write_yaml_data(rotDF, rotDF_dist, fname):
+            # Rotor Performance
+            ws = self.wb.create_sheet(title = f'Rotor Performance - {fname}')
+
+            # Write to the file
+            rotDF, rotDF_dist = read_rotor_performance_yaml(fname_hawc2)
+            for r in dataframe_to_rows(rotDF, index=False, header=True):
+                ws.append(r)
 
-        elif not self.rotDF is None:
-            # Use WISDEM output
-            myDF = self.rotDF
+            # Header row style formatting
+            for cell in ws["1:1"]:
+                cell.style = 'Headline 2'
 
-        else:
-            return
+            # Write distributed data
+            row_start = 3+len(rotDF)
+            for iws in range(len(rotDF_dist)):
+                # Reynolds number label
+                labstr = 'Wind Speed [m/s]'
+                ws.cell(row=row_start, column=1, value=labstr)
+                ws.cell(row=row_start, column=2, value=rotDF[labstr].iloc[iws])
 
-        # Write to the file
-        for r in dataframe_to_rows(myDF, index=False, header=True):
-            ws.append(r)
+                # Append to worksheet
+                for r in dataframe_to_rows(rotDF_dist[iws], index=False, header=True):
+                    ws.append(r)
+
+                # Header row style formatting
+                for cell in ws[str(row_start)+':'+str(row_start)]:
+                    cell.style = 'Headline 1'
+                for cell in ws[str(row_start+1)+':'+str(row_start+1)]:
+                    cell.style = 'Headline 2'
+
+                # Prep for next polar
+                row_start += len(rotDF_dist[iws]) + 4
+
+        # Grab OF/HAWC2 data
+        fname_hawc2 = os.path.join('..','outputs','01_steady_states','HAWC2','iea-22-280-rwt-steady-states-hawc2.yaml')
+        fname_of    = os.path.join('..','outputs','01_steady_states','OpenFAST','iea-22-280-rwt-steady-states-of.yaml')
 
-        # Header row style formatting
-        for cell in ws["1:1"]:
-            cell.style = 'Headline 2'
+        if os.path.exists(fname_hawc2):
+            rotDF_hawc2, rotDF_dist_hawc2 = read_rotor_performance_yaml(fname_hawc2)
+            write_yaml_data(rotDF_hawc2, rotDF_dist_hawc2, 'HAWC2')
+
+        if os.path.exists(fname_of):
+            rotDF_of, rotDF_dist_of = read_rotor_performance_yaml(fname_of)
+            write_yaml_data(rotDF_of, rotDF_dist_of, 'OpenFAST')
 
 
     def write_nacelle(self):

WISDEM/run_model.py

@@ -3,7 +3,6 @@
 from wisdem import run_wisdem
 import wisdem.postprocessing.compare_designs as compare_designs
 import wisdem.postprocessing.wisdem_get as getter
-import wisdem.commonse.utilities as util
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
@@ -16,72 +15,20 @@
 fname_analysis_options = os.path.join(thisdir, "analysis_options.yaml")
 folder_output = os.path.join(thisdir, "outputs")
 
-def run_15mw(fname_wt_input):
+def run_22mw(fname_wt_input):
     float_flag = fname_wt_input.find('Float') >= 0
 
     # Run WISDEM
     prob, modeling_options, analysis_options = run_wisdem(fname_wt_input, fname_modeling_options, fname_analysis_options)
 
     # Produce standard plots
-    compare_designs.run([prob], ['IEA Wind 15-MW'], modeling_options, analysis_options)
+    compare_designs.run([prob], ['IEA Wind 22-MW'], modeling_options, analysis_options)
 
     # Tabular output: Blade Shape
-    blade_shape  = np.c_[prob.get_val('blade.outer_shape_bem.s'),
-                         prob.get_val('blade.outer_shape_bem.ref_axis','m')[:,2],
-                         prob.get_val('blade.outer_shape_bem.chord','m'),
-                         prob.get_val('blade.outer_shape_bem.twist', 'deg'),
-                         prob.get_val('blade.interp_airfoils.r_thick_interp')*100,
-                         prob.get_val('blade.outer_shape_bem.pitch_axis')*100,
-                         prob.get_val('blade.outer_shape_bem.ref_axis','m')[:,0],
-                         prob.get_val('blade.outer_shape_bem.ref_axis','m')[:,1],
-                         ]
-    blade_shape_col = ['Blade Span','Rotor Coordinate [m]',
-                       'Chord [m]', 'Twist [deg]',
-                       'Relative Thickness [%]', 'Pitch Axis Chord Location [%]',
-                       'Prebend [m]', 'Sweep [m]']
-    bladeDF = pd.DataFrame(data=blade_shape, columns=blade_shape_col)
+    bladeDF = getter.get_blade_shape(prob)
 
     # Tabular output: Blade Stiffness
-    blade_stiff  = np.c_[prob.get_val('rotorse.r','m'),
-                         prob.get_val('rotorse.A','m**2'),
-                         prob.get_val('rotorse.EA','N'),
-                         prob.get_val('rotorse.EIxx','N*m**2'),
-                         prob.get_val('rotorse.EIyy','N*m**2'),
-                         prob.get_val('rotorse.EIxy','N*m**2'),
-                         prob.get_val('rotorse.GJ','N*m**2'),
-                         prob.get_val('rotorse.rhoA','kg/m'),
-                         prob.get_val('rotorse.rhoJ','kg*m'),
-                         prob.get_val('rotorse.x_ec','mm'),
-                         prob.get_val('rotorse.y_ec','mm'),
-                         prob.get_val('rotorse.re.x_tc','mm'),
-                         prob.get_val('rotorse.re.y_tc','mm'),
-                         prob.get_val('rotorse.re.x_sc','mm'),
-                         prob.get_val('rotorse.re.y_sc','mm'),
-                         prob.get_val('rotorse.re.x_cg','mm'),
-                         prob.get_val('rotorse.re.y_cg','mm'),
-                         prob.get_val('rotorse.re.precomp.flap_iner','kg/m'),
-                         prob.get_val('rotorse.re.precomp.edge_iner','kg/m')]
-    blade_stiff_col = ['Blade Span [m]',
-                       'Cross-sectional area [m^2]',
-                       'Axial stiffness [N]',
-                       'Edgewise stiffness [Nm^2]',
-                       'Flapwise stiffness [Nm^2]',
-                       'Flap-edge coupled stiffness [Nm^2]',
-                       'Torsional stiffness [Nm^2]',
-                       'Mass density [kg/m]',
-                       'Polar moment of inertia density [kg*m]',
-                       'X-distance to elastic center [mm]',
-                       'Y-distance to elastic center [mm]',
-                       'X-distance to tension center [mm]',
-                       'Y-distance to tension center [mm]',
-                       'X-distance to shear center [mm]',
-                       'Y-distance to shear center [mm]',
-                       'X-distance to mass center [mm]',
-                       'Y-distance to mass center [mm]',
-                       'Section flap inertia [kg/m]',
-                       'Section edge inertia [kg/m]',
-                       ]
-    bladeStiffDF = pd.DataFrame(data=blade_stiff, columns=blade_stiff_col)
+    bladeStiffDF = getter.get_blade_elasticity(prob)
 
     # Blade internal laminate layer details
     layerDF = []
@@ -97,101 +44,13 @@ def run_15mw(fname_wt_input):
         layerDF.append( pd.DataFrame(data=ilay, columns=layer_cols) )
 
     # Tabular output: Rotor Performance
-    rotor_perf = np.c_[prob.get_val("rotorse.rp.powercurve.V",'m/s'),
-                       prob.get_val("rotorse.rp.powercurve.pitch",'deg'),
-                       prob.get_val("rotorse.rp.powercurve.P",'MW'),
-                       prob.get_val("rotorse.rp.powercurve.Cp"),
-                       prob.get_val("rotorse.rp.powercurve.Cp_aero"),
-                       prob.get_val("rotorse.rp.powercurve.Omega",'rpm'),
-                       prob.get_val("rotorse.rp.powercurve.Omega",'rad/s')*0.5*prob["configuration.rotor_diameter_user"],
-                       prob.get_val("rotorse.rp.powercurve.T",'MN'),
-                       prob.get_val("rotorse.rp.powercurve.Ct_aero"),
-                       prob.get_val("rotorse.rp.powercurve.Q",'MN*m'),
-                       prob.get_val("rotorse.rp.powercurve.Cq_aero"),
-                       prob.get_val("rotorse.rp.powercurve.M",'MN*m'),
-                       prob.get_val("rotorse.rp.powercurve.Cm_aero"),
-                       ]
-    rotor_perf_col = ['Wind [m/s]','Pitch [deg]',
-                      'Power [MW]','Power Coefficient [-]','Aero Power Coefficient [-]',
-                      'Rotor Speed [rpm]','Tip Speed [m/s]',
-                      'Thrust [MN]','Thrust Coefficient [-]',
-                      'Torque [MNm]','Torque Coefficient [-]',
-                      'Blade Moment [MNm]','Blade Moment Coefficient [-]',
-                      ]
-    perfDF = pd.DataFrame(data=rotor_perf, columns=rotor_perf_col)
+    perfDF = getter.get_rotor_performance(prob)
 
     # Nacelle mass properties tabular
-    # Columns are ['Mass', 'CoM_x', 'CoM_y', 'CoM_z',
-    #              'MoI_cm_xx', 'MoI_cm_yy', 'MoI_cm_zz', 'MoI_cm_xy', 'MoI_cm_xz', 'MoI_cm_yz',
-    #              'MoI_TT_xx', 'MoI_TT_yy', 'MoI_TT_zz', 'MoI_TT_xy', 'MoI_TT_xz', 'MoI_TT_yz']
-    nacDF = prob.model.wt.wt_rna.drivese.nac._mass_table
-    hub_cm = prob["drivese.hub_system_cm"][0]
-    L_drive = prob["drivese.L_drive"][0]
-    tilt = prob.get_val('nacelle.uptilt', 'rad')[0]
-    shaft0 = prob["drivese.shaft_start"]
-    Cup = -1.0
-    hub_cm = R = shaft0 + (L_drive + hub_cm) * np.array([Cup * np.cos(tilt), 0.0, np.sin(tilt)])
-    hub_mass = prob['drivese.hub_system_mass']
-    hub_I = prob["drivese.hub_system_I"]
-    hub_I_TT = util.rotateI(hub_I, -Cup * tilt, axis="y")
-    hub_I_TT = util.unassembleI( util.assembleI(hub_I_TT) +
-                                 hub_mass * (np.dot(R, R) * np.eye(3) - np.outer(R, R)) )
-    blades_mass = prob['drivese.blades_mass']
-    blades_I = prob["drivese.blades_I"]
-    blades_I_TT = util.rotateI(blades_I, -Cup * tilt, axis="y")
-    blades_I_TT = util.unassembleI( util.assembleI(blades_I_TT) +
-                                    blades_mass * (np.dot(R, R) * np.eye(3) - np.outer(R, R)) )
-    rna_mass = prob['drivese.rna_mass']
-    rna_cm = R = prob['drivese.rna_cm']
-    rna_I_TT = prob['drivese.rna_I_TT']
-    rna_I = util.unassembleI( util.assembleI(rna_I_TT) +
-                                    rna_mass * (np.dot(R, R) * np.eye(3) + np.outer(R, R)) )
-    nacDF.loc['Blades'] = np.r_[blades_mass, hub_cm, blades_I, blades_I_TT].tolist()
-    nacDF.loc['Hub_System'] = np.r_[hub_mass, hub_cm, hub_I, hub_I_TT].tolist()
-    nacDF.loc['RNA'] = np.r_[rna_mass, rna_cm, rna_I, rna_I_TT].tolist()
+    nacDF = getter.get_nacelle_mass(prob)
 
     # Tabular output: Tower
-    water_depth = prob['env.water_depth']
-    h_trans = getter.get_transition_height(prob)
-    htow = getter.get_zpts(prob) #np.cumsum(np.r_[0.0, prob['towerse.tower_section_height']]) + prob['towerse.z_start']
-    t = getter.get_tower_thickness(prob)
-    towdata = np.c_[htow,
-                    getter.get_tower_diameter(prob),
-                    np.r_[t[0], t]]
-    rowadd = []
-    for k in range(towdata.shape[0]):
-        if k==0: continue
-        if k+1 < towdata.shape[0]:
-            rowadd.append([towdata[k,0]+1e-3, towdata[k,1], towdata[k+1,2]])
-    towdata = np.vstack((towdata, rowadd))
-    towdata[:,-1] *= 1e3
-    towdata = np.round( towdata[towdata[:,0].argsort(),], 3)
-    colstr = ['Height [m]','OD [m]', 'Thickness [mm]']
-    towDF = pd.DataFrame(data=towdata, columns=colstr)
-    mycomments = ['']*towdata.shape[0]
-    if not float_flag:
-        #breakpoint()
-        mycomments[0] = 'Monopile start'
-        mycomments[np.where(towdata[:,0] == -water_depth)[0][0]] = 'Mud line'
-        mycomments[np.where(towdata[:,0] == 0.0)[0][0]] = 'Water line'
-    idx_tow = np.where(towdata[:,0] == h_trans)[0][0]
-    mycomments[idx_tow] = 'Tower start'
-    mycomments[-1] = 'Tower top'
-    towDF['Location'] = mycomments
-    towDF = towDF[['Location']+colstr]
-    A = 0.25*np.pi*(towDF['OD [m]']**2 - (towDF['OD [m]']-2*1e-3*towDF['Thickness [mm]'])**2)
-    I = (1/64.)*np.pi*(towDF['OD [m]']**4 - (towDF['OD [m]']-2*1e-3*towDF['Thickness [mm]'])**4)
-    outfitting = np.zeros( len(A) )
-    if not float_flag:
-        outfitting[:idx_tow] = prob['fixedse.outfitting_factor_in']
-    outfitting[idx_tow:] = prob['towerse.outfitting_factor_in']
-    towDF['Mass Density [kg/m]'] = outfitting * getter.get_tower_rho(prob)[0] * A
-    towDF['Fore-aft inertia [kg.m]'] = towDF['Side-side inertia [kg.m]'] = towDF['Mass Density [kg/m]'] * I/A
-    towDF['Fore-aft stiffness [N.m^2]'] = towDF['Side-side stiffness [N.m^2]'] = getter.get_tower_E(prob)[0] * I
-    towDF['Torsional stiffness [N.m^2]'] = getter.get_tower_G(prob)[0] * 2*I
-    towDF['Axial stiffness [N]'] = getter.get_tower_E(prob)[0] * A
-    #with open('tow.tbl','w') as f:
-    #    towDF.to_latex(f, index=False)
+    towDF = getter.get_tower_table(prob)
 
     # Frequency plot
     fn_tower = getter.get_tower_freqs(prob)
@@ -294,7 +153,10 @@ def run_15mw(fname_wt_input):
     overview['Nacelle mass [t]'] = 1e-3*prob['drivese.nacelle_mass']
     overview['RNA mass [t]'] = 1e-3*prob['drivese.rna_mass']
     overview['Tower mass [t]'] = 1e-3*getter.get_tower_mass(prob)
-    overview['Tower base diameter [m]'] = towdata[np.where(towdata[:,0] == h_trans)[0][0],1]
+    zvec = towDF['Height [m]'].to_numpy()
+    h_trans = getter.get_transition_height(prob)
+    idx = np.where(zvec == h_trans)[0][0]
+    overview['Tower base diameter [m]'] = towDF['OD [m]'].iloc[idx]
     overview['Transition piece height [m]'] = h_trans
     if not float_flag:
         overview['Monopile embedment depth [m]'] = prob['fixedse.suctionpile_depth']
@@ -315,6 +177,6 @@ def run_15mw(fname_wt_input):
 
 
 if __name__ == '__main__':
-    run_15mw( os.path.join(ontology_dir, "IEA-22-280-RWT.yaml") )
-    run_15mw( os.path.join(ontology_dir, "IEA-22-280-RWT_Floater.yaml") )
+    run_22mw( os.path.join(ontology_dir, "IEA-22-280-RWT.yaml") )
+    run_22mw( os.path.join(ontology_dir, "IEA-22-280-RWT_Floater.yaml") )