From 3ca24b04ea473b2fb9f4140bc51674b34b076f72 Mon Sep 17 00:00:00 2001 From: yqliaohk Date: Wed, 4 Dec 2024 09:24:31 -0700 Subject: [PATCH] Updated to latest OWENS --- .../18_owens/analysis_options_owens_DVs.yaml | 65 ++- examples/18_owens/driver_weis_owens.py | 2 +- .../modeling_options_OWENS_windioExample.yaml | 108 +++-- examples/18_owens/owens_land.yaml | 10 +- weis/glue_code/gc_PoseOptimization.py | 37 ++ weis/glue_code/glue_code.py | 2 +- weis/inputs/modeling_schema.yaml | 437 +++++++++++++++--- weis/owens/openmdao_owens.py | 234 +++++----- weis/owens/unsteady_example.py | 9 +- 9 files changed, 643 insertions(+), 261 deletions(-) diff --git a/examples/18_owens/analysis_options_owens_DVs.yaml b/examples/18_owens/analysis_options_owens_DVs.yaml index 4316aab6f..5f653c3e5 100644 --- a/examples/18_owens/analysis_options_owens_DVs.yaml +++ b/examples/18_owens/analysis_options_owens_DVs.yaml @@ -3,54 +3,43 @@ general: fname_output: refturb_output design_variables: - tower: - outer_diameter: + blade: + aero_shape: + chord: + flag: True + # max: 20.0 + rotor_radius_vawt: + flag: True + n_opt: 7 + index_start: 1 + index_end: 6 + # max: 20.0 + + +merit_figure: vawt_pseudolcoe + +constraints: + owens: + SF: flag: True - lower_bound: 3.87 - upper_bound: 8.0 - layer_thickness: + lower_bound: 2.0 + fatigue: flag: True - lower_bound: 4e-3 - upper_bound: 2e-1 - -merit_figure: tower_mass + power: + flag: True # Enforce positive power -# constraints: -# tower: -# height_constraint: -# flag: False -# lower_bound: 1.e-2 -# upper_bound: 1.e-2 -# stress: -# flag: True -# global_buckling: -# flag: True -# shell_buckling: -# flag: True -# d_to_t: -# flag: True -# lower_bound: 120.0 -# upper_bound: 500.0 -# taper: -# flag: True -# lower_bound: 0.2 -# slope: -# flag: True -# frequency_1: -# flag: True -# lower_bound: 0.13 -# upper_bound: 0.40 driver: optimization: - flag: False # Flag to enable optimization + flag: True # Flag to enable optimization tol: 1.e-6 # Optimality tolerance # max_major_iter: 10 # Maximum number of major design iterations (SNOPT) # max_minor_iter: 100 # Maximum number of minor design iterations (SNOPT) - max_iter: 100 # Maximum number of iterations (SLSQP) - solver: SNOPT # Optimization solver. Other options are 'SLSQP' - 'CONMIN' - step_size: 1.e-6 # Step size for finite differencing + max_iter: 200 # Maximum number of iterations (SLSQP) + solver: SLSQP # Optimization solver. Other options are 'SLSQP' - 'CONMIN' + step_size: 1.e-4 # Step size for finite differencing form: forward # Finite differencing mode, either forward or central + setp_calc: rel_avg # design_of_experiments: # flag: True # Flag to enable design of experiments # run_parallel: False # Flag to run using parallel processing diff --git a/examples/18_owens/driver_weis_owens.py b/examples/18_owens/driver_weis_owens.py index af8ff58e4..6e6d14114 100644 --- a/examples/18_owens/driver_weis_owens.py +++ b/examples/18_owens/driver_weis_owens.py @@ -19,7 +19,7 @@ analysis_override = {} analysis_override['general'] = {} -analysis_override['general']['folder_output'] = os.path.join('outputs/18_OWENS_OptStudies/1_change_opt/',"slsqp") +analysis_override['general']['folder_output'] = os.path.join('outputs/18_OWENS_OptStudies/1_change_opt_chord_radius/',"slsqp") analysis_override['driver'] = {} analysis_override['driver']['optimization'] = {} analysis_override['driver']['optimization']['solver'] = "SLSQP" diff --git a/examples/18_owens/modeling_options_OWENS_windioExample.yaml b/examples/18_owens/modeling_options_OWENS_windioExample.yaml index 64881973b..4ffb9e146 100644 --- a/examples/18_owens/modeling_options_OWENS_windioExample.yaml +++ b/examples/18_owens/modeling_options_OWENS_windioExample.yaml @@ -21,52 +21,70 @@ WISDEM: # flag: True OWENS: - flag: True - OWENS_project_path: /Users/yliao/Documents/Projects/CT-OPT/OWENS - master_input: /Users/yliao/repos/OWENS.jl/docs/src/literate/sampleOWENS.yml - analysisType: unsteady # unsteady, steady, modal - turbineType: Darrieus #Darrieus, H-VAWT, ARCUS - controlStrategy: constantRPM # TODO: incorporate the others - RPM: 17.2 #RPM - numTS: 20 # - delta_t: 0.01 # s - Nslices: 30 # number of OWENSAero discritizations #TODO: AD parameters - ntheta: 30 # number of OWENSAero azimuthal discretizations - AModel: DMS # AD, DMS, AC - structuralModel: GX #GX, TNB, ROM - structuralNonlinear: false #TODO: propogate - ntelem: 10 #tower elements in each - nbelem: 60 #blade elements in each - ncelem: 10 #central cable elements in each if turbineType is ARCUS - nselem: 5 #strut elements in each if turbineType has struts - run_path: /Users/yliao/repos/OWENS.jl/docs/src/literate/ + flag: True + general: + OWENS_project_path: /Users/yliao/Documents/Projects/CT-OPT/OWENS + master_input: /Users/yliao/repos/OWENS.jl/examples/literate/sampleOWENS.yml + run_path: /Users/yliao/repos/WEIS/examples/18_owens # points to the data folders + analysisType: unsteady # unsteady, steady, modal + AeroModel: "DMS" # OWENSAero model "DMS" for double multiple streamtube or "AC" for actuator cylinder, or "AD" for aerodyn + structuralModel: "TNB" #Structural models available: TNB full timoshenko beam elements with time newmark beta time stepping, ROM reduced order modal model of the timoshenko elements, GX with GXBeam's methods for geometrically exact beam theory and more efficient methods and time stepping + controlStrategy: "prescribedRPM" # should be in WindIO?- yes, + numTS: 20 # number of time steps TODO: change to sim time and make this derived + delta_t: 0.01 # time step in seconds + dataOutputFilename: "./InitialDataOutputs.out" # data output filename with path, set to nothing or don't specify to not output anything + TOL: 1e-4 # gauss-seidel iteration tolerance - i.e. the two-way iteration tolerance + MAXITER: 300 # gauss-seidel max iterations - i.e. the two-way iterations + verbosity: 2 # verbosity where 0 is nothing, 1 is warnings, 2 is summary outputs, 3 is detailed outputs, and 4 is everything + VTKsaveName: "./vtk/windio" # Path and name of the VTK outputs, recommended to put it in its own folder (which it will automatically create if needed) + aeroLoadsOn: 2 # Level of aero coupling 0 structures only, 1 no deformation passed to the aero, 2 two-way coupling, 1.5 last time step's deformations passed to this timesteps aero and no internal iteration. + Prescribed_RPM_time_controlpoints: [0.0,100000.1] + Prescribed_RPM_RPM_controlpoints: [17.2,17.2] + Prescribed_Vinf_time_controlpoints: [0.0,100000.1] + Prescribed_Vinf_Vinf_controlpoints: [17.2,17.2] + + DLC_Options: + DLCs: ["1_1","2_1"] # name of DLC + Vinf_range: [17.2] # inflow Cutin to cutout and discretization + IEC_std: "\"1-ED3\"" # turbsim input file IEC standard + WindChar: "\"A\"" # turbsim wind charasteric + WindClass: 1 # turbsim wind class + turbsimsavepath: "./turbsimfiles" # path where the turbsim files are saved + pathtoturbsim: nothing # path to the turbsim executable + NumGrid_Z: 38 # turbsim vertical discretizations + NumGrid_Y: 26 # turbsim horizontal discretizations + Vref: 17.2 # reference/nominal wind speed m/s for turbsim or other inflow wind input file (depending on which DLC is selected) + Vdesign: 11.0 # Design or rated speed of turbine, used for certain DLC cases + grid_oversize: 1.1 # amount that the turbsim inflow is oversized compared to the turbine to allow for deflection + regenWindFiles: false #, force regeneration of turbsim files even if they already exist + delta_t_turbsim: 0.05 # turbsim timestep + simtime_turbsim: 600.0 # turbsim total time, which loops if simtime exceeds turbsim time -# designParameters: - eta: 0.5 # blade mount point ratio, 0.5 is the blade half chord is perpendicular with the axis of rotation, 0.25 is the quarter chord, etc - # move to geometry? - # Nbld: 3 # number of blades - # Blade_Radius: 54.01123056 # blade height m - # Blade_Height: 110.1829092 # blade radius m - # towerHeight: 3.0 # m tower extension height below blades + OWENSAero_Options: + Nslices: 30 # number of 3-D slices for the strip method to go from 2D to 3D considering curved deforming blades + ntheta: 30 # number of azimuthal discretizations + ifw: false # use the OpenFASTWrappers inflow wind coupling to get inflow velocities + DynamicStallModel: "BV" # dynamic stall model, should be under an OWENSAero options + RPI: true # rotating point iterative method (i.e. it just calculates at the blade positions and is much faster) + Aero_Buoyancy_Active: false # flag to turn buoyancy on for the blades. This is likely to be replaced by a different model -# operationParameters: -# rho: 1.225 # air density -# Vinf: 17.2 # m/s #optional, supersceeded if ifw=true + OWENSFEA_Options: + nlOn: true #nonlinear effects + RayleighAlpha: 0.05 #damping coefficient scalar on the stiffness matrix + RayleighBeta: 0.05 #damping coefficient scalar on the mass matrix + iterationType: "DI" #internal iteration type DI direct iteration, NR newton rhapson (which is less stable than DI) + numModes: 20 #if ROM model, number of modes used in the analysis type. Less is faster but less accurate + tolerance: 1.0e-06 #total mesh unsteady analysis convergence tolerance for a timestep within the structural model + maxIterations: 50 #total mesh unsteady analysis convergence max iterations for a timestep + AddedMass_Coeff_Ca: 0.0 #added mass coefficient, scaling factor (typically 0-1) on the cones of water mass applied to each structural element in the 22 and 33 diagonal terms. 0 turns this off -# turbulentInflow: - ifw: false - windType: 3 - windINPfilename: /Users/yliao/repos/OWENS.jl/examples/Optimization/data/turbsim/115mx115m_30x30_20.0msETM.bts - ifw_libfiles: None + Mesh_Options: + ntelem: 10 # number of tower elements in each blade, plus nodes wherever there is a component overlap + nbelem: 60 # number of blade elements in each blade, plus nodes wherever there is a component overlap + nselem: 5 # number of elements in each strut + angularOffset: -1.5707963267948966 + AD15hubR: 0.1 # parameter, used in aerodyn coupling for the hub radius so that the vortex sheets don't go within the hub + turbineType: Darrieus -# AeroParameters: - adi_lib: None - adi_rootname: "/ExampleB" - -# structuralParameters: -# NuMad_geom_xlscsv_file_twr: /data/NuMAD_Geom_SNL_5MW_D_TaperedTower.csv -# NuMad_mat_xlscsv_file_twr: /data/NuMAD_Materials_SNL_5MW.csv -# NuMad_geom_xlscsv_file_bld: /data/NuMAD_Geom_SNL_5MW_D_Carbon_LCDT_ThickFoils_ThinSkin.csv -# NuMad_mat_xlscsv_file_bld: /data/NuMAD_Materials_SNL_5MW.csv -# NuMad_geom_xlscsv_file_strut: /data/NuMAD_Geom_SNL_5MW_Struts.csv -# NuMad_mat_xlscsv_file_strut: /data/NuMAD_Materials_SNL_5MW.csv + OWENSOpenFASTWrappers_Options: + windINPfilename: /data/turbsim/115mx115m_30x30_20.0msETM.bts #OWENSOpenFASTWrappers If ifw or AeroDyn is being used, gets overwritten if using the DLC analysis type, the moordyn file location, like in the unit test \ No newline at end of file diff --git a/examples/18_owens/owens_land.yaml b/examples/18_owens/owens_land.yaml index 583a82488..a5bbd76d2 100644 --- a/examples/18_owens/owens_land.yaml +++ b/examples/18_owens/owens_land.yaml @@ -12,7 +12,7 @@ components: values: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] z: grid: *id000 - values: [0, 6.607421057, 13.21484211, 19.82226317, 26.42968423, 33.03710529, 39.64452634, 46.2519474, 52.85936846, 59.46678951, 66.07421057, 72.68163163, 79.28905268, 85.89647374, 92.5038948, 99.11131586, 105.7187369, 112.326158, 118.933579, 125.5410001, 132.1484211] + values: [0, 6.607421057, 13.21484211, 19.82226317, 26.42968423, 33.03710529, 39.64452634, 46.2519474, 52.85936846, 59.46678951, 66.07421057, 72.68163163, 79.28905268, 85.89647374, 92.5038948, 99.11131586, 105.7187369, 112.326158, 118.933579, 125.5410001, 132.1484211] # The reference axis z of tower cannot be zeros, otherwise the towerse in wisdem will give errors outer_diameter: grid: *id000 values: [5.888066637, 5.788249614, 5.684865834, 5.577578504, 5.465996871, 5.349663569, 5.22803793, 5.100473576, 4.966187737, 4.824218239, 4.673361587, 4.512081003, 4.338364592, 4.149496455, 3.941665819, 3.70924955, 3.443362258, 3.128502234, 2.732998436, 2.169182296, 2.0] @@ -51,7 +51,7 @@ components: material: Generic_Skin thickness: grid: *id000 - values: [0.089213, 0.087701, 0.086134, 0.084509, 0.082818, 0.081056, 0.079213, 0.07728 , 0.075245, 0.073094, 0.070809, 0.068365, 0.065733, 0.062871, 0.059722, 0.056201, 0.052172, 0.047402, 0.041409, 0.032866, 0.030303] + values: [0.008921, 0.00877 , 0.008613, 0.008451, 0.008282, 0.008106,0.007921, 0.007728, 0.007525, 0.007309, 0.007081, 0.006836, 0.006573, 0.006287, 0.005972, 0.00562 , 0.005217, 0.00474 , 0.004141, 0.003287, 0.00303 ] # n_plies: # grid: *id000 # values: [8.921, 8.770, 8.613, 8.451, 8.282, 8.106, 7.921, 7.728, 7.525, 7.309, 7.081, 6.836, 6.573, 6.287, 5.972, 5.620, 5.217, 4.740, 4.141, 3.287, 3.030] @@ -72,7 +72,7 @@ components: pitch_axis: grid: *id001 values: [0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3] - + blade_mountpoint: 0.5 # OWENS only uses the reference axis of internal structure 2d fem but in windIO this is usually defined equal to the refenrece axis in the outer shape bem, so move the values in internal sturcture 2d fem here reference_axis: x: &xid001 @@ -326,6 +326,8 @@ components: joint: {position: 0.0, mass: 0.0, cost: 0.0, bolt: M30, nonmaterial_cost: 0.0, reinforcement_layer_ss: joint_reinf_ss, reinforcement_layer_ps: joint_reinf_ps} struts: # - name: strut1 + mountfraction_blade: [0.11, 0.89] + mountfraction_tower: [0.11, 0.89] outer_shape_bem: airfoil_position: grid: &id002 [0, 0.33331652, 0.66668348, 1.0] # Changed to nondimensionalized @@ -750,7 +752,7 @@ materials: # m: [0.0,1.0,2.0,4.0,6.0,20.0] - {name: resin, description: epoxy, source: Described in https://doi.org/10.5194/wes-7-19-2022, E: 1000000.0, nu: 0.3, G: 312500.0, GIc: 0.0, GIIc: 0.0, alp0: 0.0, Xt: 0.0, Xc: 0.0, S: 0.0, rho: 1150.0, alpha: 0.0, orth: 0.0, unit_cost: 3.63} # No resin in the original OWENS yaml, but WEIS requires resin for composites - {name: Adhesive, orth: 0.0, rho: 1100.0, E: 4560000.0, G: 1450000.0, nu: 0.3, alpha: 0.0, Xt: 61510000.0, Xc: 65360000.0, S: 36610000.0, GIc: 0.0, GIIc: 0.0, alp0: 0.0, ply_t: 0.001, waste: 0.25, unit_cost: 7.23} # No adhesive in the original OWENS yaml, but WEIS requires adhesive for cost models -control: +control: # This has to be here for connection to work. Copied from IEA22 supervisory: {Vin: 0.5, Vout: 4.0, maxTS: 60} pitch: {PC_zeta: !!null '', PC_omega: !!null '', ps_percent: !!null '', max_pitch: !!null '', max_pitch_rate: 0.1745, min_pitch: 0.00088} torque: {control_type: !!null '', tsr: 7.64, VS_zeta: !!null '', VS_omega: !!null '', max_torque_rate: 1500000.0, VS_minspd: 0.0, VS_maxspd: 1.26711} diff --git a/weis/glue_code/gc_PoseOptimization.py b/weis/glue_code/gc_PoseOptimization.py index 53ac0daf0..247f41d8a 100644 --- a/weis/glue_code/gc_PoseOptimization.py +++ b/weis/glue_code/gc_PoseOptimization.py @@ -1,5 +1,6 @@ from wisdem.glue_code.gc_PoseOptimization import PoseOptimization import numpy as np +from scipy.interpolate import PchipInterpolator class PoseOptimizationWEIS(PoseOptimization): @@ -113,6 +114,9 @@ def set_objective(self, wt_opt): elif self.opt['merit_figure'] == 'OL2CL_pitch': wt_opt.model.add_objective('aeroelastic.OL2CL_pitch') + + elif self.opt['merit_figure'] == 'vawt_pseudolcoe': + wt_opt.model.add_objective('owens.lcoe') else: super(PoseOptimizationWEIS, self).set_objective(wt_opt) @@ -216,6 +220,30 @@ def set_design_variables(self, wt_opt, wt_init): upper=tmd_group['damping_ratio']['upper_bound'] ) + # OWENS radius + rotor_radius_vawt_opt = self.opt["design_variables"]["blade"]["aero_shape"]["rotor_radius_vawt"] + if rotor_radius_vawt_opt["flag"]: + if rotor_radius_vawt_opt["index_end"] > rotor_radius_vawt_opt["n_opt"]: + raise Exception( + "Check the analysis options yaml, index_end of the blade radius is higher than the number of DVs n_opt" + ) + elif rotor_radius_vawt_opt["index_end"] == 0: + rotor_radius_vawt_opt["index_end"] = rotor_radius_vawt_opt["n_opt"] + indices_radius = range(rotor_radius_vawt_opt["index_start"], rotor_radius_vawt_opt["index_end"]) + s_opt_radius = np.linspace(0.0, 1.0, rotor_radius_vawt_opt["n_opt"]) + radius_interpolator = PchipInterpolator( + wt_init["components"]["blade"]["outer_shape_bem"]["reference_axis"]["x"]["grid"], + wt_init["components"]["blade"]["outer_shape_bem"]["reference_axis"]["x"]["values"], + ) + init_radius_opt = radius_interpolator(s_opt_radius) + wt_opt.model.add_design_var( + "blade.opt_var.rotor_radius_vawt", + indices=indices_radius, + lower=init_radius_opt[indices_radius] * rotor_radius_vawt_opt["max_decrease"], + upper=init_radius_opt[indices_radius] * rotor_radius_vawt_opt["max_increase"], + ) + + return wt_opt @@ -396,6 +424,15 @@ def set_constraints(self, wt_opt): wt_opt.model.add_constraint('aeroelastic.damage_tower_base',upper = tower_base_damage_max) + # vawt owens constraints + owens_constr = self.opt["constraints"]["owens"] + if owens_constr["SF"]["flag"]: + wt_opt.model.add_constraint("owens.SF", lower=self.opt["constraints"]["owens"]["SF"]["lower_bound"]) + if owens_constr["fatigue"]["flag"]: + wt_opt.model.add_constraint("owens.fatigue_damage", upper=1.0) + if owens_constr["power"]["flag"]: + wt_opt.model.add_constraint("owens.power", lower=0.0) + return wt_opt diff --git a/weis/glue_code/glue_code.py b/weis/glue_code/glue_code.py index 690d24710..51a895b26 100644 --- a/weis/glue_code/glue_code.py +++ b/weis/glue_code/glue_code.py @@ -960,7 +960,7 @@ def setup(self): self.connect("blade.pa.chord_param", "owens.blade_chord_values") self.connect("blade.pa.twist_param", "owens.blade_twist_values") self.connect("blade.outer_shape_bem.pitch_axis", "owens.blade_pitch_axis_values") - self.connect("blade.outer_shape_bem.ref_axis", "owens.blade_ref_axis") + self.connect("blade.high_level_blade_props.blade_ref_axis", "owens.blade_ref_axis") self.connect("blade.internal_structure_2d_fem.web_start_nd", "owens.blade_web_start_nd_arc") self.connect("blade.internal_structure_2d_fem.web_end_nd", "owens.blade_web_end_nd_arc") diff --git a/weis/inputs/modeling_schema.yaml b/weis/inputs/modeling_schema.yaml index c3d58e50c..c4701705d 100644 --- a/weis/inputs/modeling_schema.yaml +++ b/weis/inputs/modeling_schema.yaml @@ -2879,64 +2879,385 @@ properties: flag: type: boolean default: False - analysisType: - type: string - default: unsteady # unsteady, steady, modal - turbineType: - type: string - default: Darrieus #Darrieus, H-VAWT, ARCUS - controlStrategy: - type: string - default: constantRPM - description: currently constant RPM (in work is to expose the scripting options of the discon and custom control options at the super simple file level) - RPM: - type: number - default: 0 - description: Rotational rate if constant RPM control, is also the initialized RPM - numTS: - type: integer - default: 1 - description: Total number of time steps - delta_t: - type: number - default: 0.01 - description: Time step, in seconds - Nslices: - type: integer - default: 1 - descrption: Number of vertical slices in the strip theory to get the simplified aero model to be 3D - ntheta: - type: integer - default: 1 - description: Number of azimuthal discretizations of the simplified aero model, if active, must by wholly divisible by NBlades - AModel: - type: string - default: DMS - description: Aero model, "AC" actuator cylinder (slower but slightly more accurate with array capability at the planar aero only level) "DMS" (faster and nearly as accurate) or "AD" (aerodyn interface, specifically OLAF, which is a relatively slow but higher fidelity vortex line wake model. This includes strut drag and tip losses.) - structuralModel: - type: string - default: GX - description: TNB time newmark beta time stepping linear or nonlinear, ROM reduced order modal linear or nonlinear, GX geometrically exact beam theory - structural_nonlinear: - type: boolean - default: false - description: Flag for nonlinear structural analysis - ntelem: - type: integer - default: 1 - description: Number of tower elements, in the automatically generated mesh - nbelem: - type: integer - default: 1 - description: Number of blade elements, in the automatically generated mesh - ncelem: - type: integer - default: 1 - description: Number of cable elements, in the automatically generated mesh - nselem: - type: integer - default: 1 - description: Number of strut elements, in the automatically generated mesh + general: + type: object + default: {} + properties: + OWENS_project_path: + type: string + default: none + desciption: Path to OWENS project + master_input: + type: string + default: nothing + desciption: A OWENS mater input file for initialization + run_path: + type: string + default: nothing + desciption: Points to the data folders + analysisType: + type: string + default: unsteady # unsteady, steady, modal + AeroModel: + type: string + default: DMS + description: Aero model, "AC" actuator cylinder (slower but slightly more accurate with array capability at the planar aero only level) "DMS" (faster and nearly as accurate) or "AD" (aerodyn interface, specifically OLAF, which is a relatively slow but higher fidelity vortex line wake model. This includes strut drag and tip losses.) + structuralModel: + type: string + default: GX + description: TNB time newmark beta time stepping linear or nonlinear, ROM reduced order modal linear or nonlinear, GX geometrically exact beam theory + controlStrategy: + type: string + default: constantRPM + description: currently constant RPM (in work is to expose the scripting options of the discon and custom control options at the super simple file level) + platformActive: + type: boolean + default: false + description: Flag to indicate if the floating platform model is active. + topsideOn: + type: boolean + default: true + desciption: Flag to be able to turn off the rotor and just run the floating portions + InterpOrder: + type: integer + default: 2 + description: if platformActive, order used for extrapolating inputs and states, 0 flat, 1 linear, 2 quadratic + defaultOutputFilename: + type: string + default: nothing + desciption: Data output filename with path, set to nothing or don't specify to not output anything + rigid: + type: boolean + default: false + description: This bypasses the structural solve and just mapps the applied loads as the reaction loads, and the deflections remain 0 + TOL: + type: number + default: 1e-4 + desciption: Gauss-seidel iteration tolerance - i.e. the two-way iteration tolerance + MAXITER: + type: integer + default: 300 + desciption: Gauss-seidel max iterations - i.e. the two-way iterations + verbosity: + type: integer + default: 2 + description: Verbosity where 0 is nothing, 1 is warnings, 2 is summary outputs, 3 is detailed outputs, and 4 is everything + VTKsaveName: + type: string + default: ./vtk/windio + desciption: Path and name of the VTK outputs, recommended to put it in its own folder (which it will automatically create if needed) + aeroLoadsOn: + type: integer + default: 2 + description: Level of aero coupling 0 structures only, 1 no deformation passed to the aero, 2 two-way coupling, 1.5 last time step's deformations passed to this timesteps aero and no internal iteration + Prescribed_RPM_time_controlpoints: + type: array + default: [0.0,100000.1] + description: If controlStrategy is "fixedRPM", array of time control points for the internal spline + Prescribed_RPM_RPM_controlpoints: + type: array + default: [17.2,17.2] + desciption: If controlStrategy is "fixedRPM", array of RPM control points for the internal spline + Prescribed_Vinf_time_controlpoints: + type: array + default: [0.0,100000.1] + desciption: If AeroModel is "DMS" or "AC, and ifw is false, array of time control points for the internal spline + Prescribed_Vinf_Vinf_controlpoints: + type: array + default: [17.2,17.2] + desciption: If AeroModel is "DMS" or "AC, and ifw is false, array of Vinf control points for the internal spline + numTS: + type: integer + default: 1 + description: Total number of time steps + delta_t: + type: number + default: 0.01 + description: Time step, in seconds + # turbineType: + # type: string + # default: Darrieus #Darrieus, H-VAWT, ARCUS + + DLC: + type: object + default: {} + properties: + DLCs: + type: array + default: ["none"] + description: name of DLC + Vinf_range: + type: array + default: [17.2] + description: inflow Cutin to cutout and discretization + IEC_std: + type: string + default: 1-ED3 + description: turbsim input file IEC standard + WindChar: + type: string + default: A + description: turbsim wind charasteric + WindClass: + type: integer + default: 1 + description: DLC turbsim wind class + turbsimsavepath: + type: string + default: ./turbsimfiles + description: path where the turbsim files are saved + pathtoturbsim: + type: string + default: nothing + description: path to the turbsim executable + NumGrid_Z: + type: integer + default: 8 + description: turbsim vertical discretizations + NumGrid_Y: + type: integer + default: 6 + description: turbsim horizontal discretizations + Vref: + type: number + default: 0.0 + description: reference/nominal wind speed m/s for turbsim or other inflow wind input file (depending on which DLC is selected) + Vdesign: + type: number + default: 1.0 + description: Design or rated speed of turbine, used for certain DLC cases + grid_oversize: + type: number + default: 0.1 + description: amount that the turbsim inflow is oversized compared to the turbine to allow for deflection + regenWindFiles: + type: boolean + default: false + desciption: force regeneration of turbsim files even if they already exist + delta_t_turbsim: + type: number + default: 0.05 + description: turbsim timestep + simtime_turbsim: + type: number + default: 0.0 + description: turbsim total time, which loops if simtime exceeds turbsim time + RandSeed1: + type: integer + default: 40071 + description: turbsim random seed number + AeroOptions: + type: object + default: {} + properties: + Nslices: + type: integer + default: 1 + descrption: Number of 3-D slices for the strip method to go from 2D to 3D considering curved deforming blades + ntheta: + type: integer + default: 1 + description: Number of azimuthal discretizations of the simplified aero model, if active, must by wholly divisible by NBlades + ifw: + type: boolean + default: false + description: Use the inflow wind coupling to get inflow velocities + DynamicStallModel: + type: string + default: "BV" + description: Dynamic stall model, should be under an OWENSAero options + RPI: + type: boolean + default: true + description: rotating point iterative method (i.e. it just calculates at the blade positions and is much faster) + Aero_Buoyancy_Active: + type: boolean + default: false + description: flag to turn buoyancy on for the blades. This is likely to be replaced by a different model + Aero_AddedMass_Active: + type: boolean + default: false + description: flag to turn added mass forces on, don't turn on if the added mass in the structures are on + Aero_RotAccel_Active: + type: boolean + default: false + desciption: Flag to turn added mass forces on, don't turn on if the added mass in the structures are on + FEAOptions: + nlOn: + type: boolean + default: true + desciption: Nonlinear effects + RayleighAlpha: + type: number + default: 0.05 + desciption: damping coefficient scalar on the stiffness matrix + RayleighBeta: + type: number + default: 0.05 + description: damping coefficient scalar on the mass matrix + iterationType: + type: string + default: DI + description: internal iteration type DI direct iteration, NR newton rhapson (which is less stable than DI) + guessFreq: + type: number + default: 0.0 + description: for the built in flutter model frequency guessed for the flutter frequency + numModes: + type: integer + default: 20 + desciption: ROM model, number of modes used in the analysis type. Less is faster but less accurate + adaptiveLoadSteppingFlag: + type: boolean + default: true + desciption: For steady analysis if convergence fails, it will reduce the load and retry then increase the load + minLoadStepDelta: + type: number + default: 0.05 + description: minimum change in load step + minLoadStep: + type: number + default: 0.05 + desciption: minimum value of reduced load + prescribedLoadStep: + type: number + default: 0.0 + description: optional prescribed load fraction + maxNumLoadSteps: + type: interger + default: 20 + desciption: used in static (steady state) analysis, max load steps for adaptive load stepping + tolerance: + type: number + defaul: 1.0000e-06 + desciption: total mesh unsteady analysis convergence tolerance for a timestep within the structural model + maxIterations: + type: integer + default: 50 + desciption: total mesh unsteady analysis convergence max iterations for a timestep + elementOrder: + type: integer + default: 1 + desciption: Element order, 1st order, 2nd order etc; determines the number of nodes per element (order +1). Orders above 1 have not been tested in a long time + alpha: + type: number + default: 0.5 + desciption: newmark time integration alpha parameter + gamma: + type: number + default: 0.5 + decription: newmark time integration gamma parameter + AddedMass_Coeff_Ca: + type: number + default: 0.0 + desciption: added mass coefficient, scaling factor (typically 0-1) on the cones of water mass applied to each structural element in the 22 and 33 diagonal terms. 0 turns this off + platformTurbineConnectionNodeNumber: + type: integer + default: 1 + aeroElasticOn: + type: boolean + default: false + description: OWENSFEA for the built in flutter model + spinUpOn: + type: boolean + default: true + desciption: OWENSFEA modal analysis, calculates steady centrifugal strain stiffening and then passes that model to the modal analysis + predef: + type: boolean + default: false + desciption: Predeformation flag for two state analysis where a portion of the blade is deformed and the nonlinear strain stiffening terms are "update"-d, then "use"-d in two different analysis + Mesh_Options: + type: object + default: {} + properties: + ntelem: + type: integer + default: 20 + description: number of tower elements in each blade, plus nodes wherever there is a component overlap + nbelem: + type: integer + default: 30 + desciption: number of blade elements in each blade, plus nodes wherever there is a component overlap + ncelem: + type: integer + default: 30 + desciption: number of cable elements in each cable if ARCUS + nselem: + type: integer + default: 10 + description: number of elements in each strut + angularOffset: + type: number + default: 0.0 + decription: moves the structure to align with the aero model + joint_type: + type: integer + default: 0 + desciption: optionally can specify the strut to blade joints to be pinned about different axes, or 0 for welded + c_mount_ratio: + type: number + default: 0.05 + description: for ARCUS, where the cable mounts on the lower side of the blade + AD15hubR: + type: number + default: 0.1 + description: parameter, used in aerodyn coupling for the hub radius so that the vortex sheets don't go within the hub + cables_connected_to_blade_base: + type: boolean + default: true + desciption: for ARCUS, for the two part simulation of the blade bending + turbineType: + type: string + default: Darrieus" + desciption: mesh Darrieus, H-VAWT, controls if the tips of the blades are joined to the tower in the mesh or not. + OWENSOpenFASTWrappers: + type: object + default: {} + properties: + windINPfilename: + type: string + default: nothing + decription: If ifw or AeroDyn is being used, gets overwritten if using the DLC analysis type, the moordyn file location, like in the unit test + ifw_libfile: + type: string + default: nothing + desciption: location of the respective OpenFAST library, if nothing it will use the internal OWENS installation + hd_lib: + type: string + default: nothing + description: location of the respective OpenFAST library, if nothing it will use the internal OWENS installation + md_lib: + type: string + default: nothing + desciption: location of the respective OpenFAST library, if nothing it will use the internal OWENS installation + adi_lib: + type: string + default: nothing + desciption: location of the respective OpenFAST library, if nothing it will use the internal OWENS installation + adi_rootname: + type: string + default: /aerodyn + desciption: location of the respective OpenFAST library, if nothing it will use the internal OWENS installation + hd_input_file: + type: string + default: none + decription: If platformActive, the hydrodyn file location, like in the unit test + ss_input_file: + type: string + default: none + desciption: If platformActive, the sea state file location, like in the unit test + md_input_file: + type: string + default: none + decription: If platformActive, the moordyn file location, like in the unit test + potflowfile: + type: string + default: nothing + desciption: If platformActive, the potential flow files location, like in the unit test + WindType: + type: integer + default: 3 + decription: Derived parameter, inflowwind wind file type when DLC generator is active, matches inflowwind WindType HydroDyn: *ofhydrodyn MoorDyn: *ofmoordyn diff --git a/weis/owens/openmdao_owens.py b/weis/owens/openmdao_owens.py index 0a336608a..738401e24 100644 --- a/weis/owens/openmdao_owens.py +++ b/weis/owens/openmdao_owens.py @@ -68,7 +68,7 @@ def setup(self): rotorse_options = self.options["rotorse_options"] towerse_options = self.options["towerse_options"] strut_options = self.options["strut_options"] - OWENS_path = modopt["OWENS"]["OWENS_project_path"] + OWENS_path = modopt["OWENS"]["general"]["OWENS_project_path"] jlPkg.activate(OWENS_path) @@ -78,11 +78,6 @@ def setup(self): jl.seval("using OrderedCollections") - master_file = modopt["OWENS"]["master_input"] - - # Initialize an OWENS model - self.model = jl.OWENS.MasterInput(master_file) - self.n_span = rotorse_options["n_span"] self.n_layers = rotorse_options["n_layers"] self.n_webs = rotorse_options["n_webs"] @@ -98,32 +93,63 @@ def setup(self): n_webs_strut = strut_options["n_webs"] - # analysisType: unsteady # unsteady, steady, modal - # turbineType: Darrieus #Darrieus, H-VAWT, ARCUS - # These should be in modeling options? - # Currently, they are just normal openmdao component options, - # Not considering any WEIS data/input structure - # Probably needs to be like modopt["crossflow"]["analysis_type"] - self.analysis_type = modopt["OWENS"]["analysisType"] - self.turbine_type = modopt["OWENS"]["turbineType"] - self.eta = modopt["OWENS"]["eta"] # Now eta is a modeling option - self.control_strategy = modopt["OWENS"]["controlStrategy"] - self.aeroModel = modopt["OWENS"]["AModel"] - self.adi_lib = modopt["OWENS"]["adi_lib"] - self.adi_rootname = modopt["OWENS"]["adi_rootname"] - # self.NumadSpec = modopt["OWENS"]["NumadSpec"] # All files go here - # self.Turbulence = modopt["OWENS"]["TurbulenceSpec"] # All files go here - self.ifw = modopt["OWENS"]["ifw"] - self.windType = modopt["OWENS"]["windType"] - self.windINPfilename = modopt["OWENS"]["windINPfilename"] - self.ifw_libfiles = modopt["OWENS"]["ifw_libfiles"] - self.n_blades = modopt["assembly"]["number_of_blades"] - self.structuralModel = modopt["OWENS"]["structuralModel"] - self.structuralNonlinear = modopt["OWENS"]["structuralNonlinear"] - self.run_path = modopt["OWENS"]["run_path"] # What exactly is this path? - - # Reinitialize the model with the inputs from modeling options - self.initialize_model() + # Initialize OWENS modeling options + self.owens_modeling_options = {} + self.owens_modeling_options["OWENS_Options"] = {} + self.owens_modeling_options["OWENS_Options"]["analysisType"] = modopt["OWENS"]["general"]["analysisType"] + self.owens_modeling_options["OWENS_Options"]["AeroModel"] = modopt["OWENS"]["general"]["AeroModel"] + self.owens_modeling_options["OWENS_Options"]["structuralModel"] = modopt["OWENS"]["general"]["structuralModel"] + self.owens_modeling_options["OWENS_Options"]["controlStrategy"] = modopt["OWENS"]["general"]["controlStrategy"] + self.owens_modeling_options["OWENS_Options"]["numTS"] = modopt["OWENS"]["general"]["numTS"] + self.owens_modeling_options["OWENS_Options"]["delta_t"] = modopt["OWENS"]["general"]["delta_t"] + self.owens_modeling_options["OWENS_Options"]["dataOutputFilename"] = modopt["OWENS"]["general"]["dataOutputFilename"] + self.owens_modeling_options["OWENS_Options"]["MAXITER"] = modopt["OWENS"]["general"]["MAXITER"] + self.owens_modeling_options["OWENS_Options"]["TOL"] = modopt["OWENS"]["general"]["TOL"] + self.owens_modeling_options["OWENS_Options"]["verbosity"] = modopt["OWENS"]["general"]["verbosity"] + self.owens_modeling_options["OWENS_Options"]["VTKsaveName"] = modopt["OWENS"]["general"]["VTKsaveName"] + self.owens_modeling_options["OWENS_Options"]["aeroLoadsOn"] = modopt["OWENS"]["general"]["aeroLoadsOn"] + self.owens_modeling_options["OWENS_Options"]["Prescribed_RPM_time_controlpoints"] = modopt["OWENS"]["general"]["Prescribed_RPM_time_controlpoints"] + self.owens_modeling_options["OWENS_Options"]["Prescribed_RPM_RPM_controlpoints"] = modopt["OWENS"]["general"]["Prescribed_RPM_RPM_controlpoints"] + self.owens_modeling_options["OWENS_Options"]["Prescribed_Vinf_time_controlpoints"] = modopt["OWENS"]["general"]["Prescribed_Vinf_time_controlpoints"] + self.owens_modeling_options["OWENS_Options"]["Prescribed_Vinf_Vinf_controlpoints"] = modopt["OWENS"]["general"]["Prescribed_Vinf_Vinf_controlpoints"] + + + # Write out modeling option file + self.owens_modeling_options["DLC_Options"] = modopt["OWENS"]["DLC_Options"] + self.owens_modeling_options["OWENSAero_Options"] = modopt["OWENS"]["OWENSAero_Options"] + self.owens_modeling_options["OWENSFEA_Options"] = modopt["OWENS"]["OWENSFEA_Options"] + self.owens_modeling_options["Mesh_Options"] = modopt["OWENS"]["Mesh_Options"] + self.owens_modeling_options["OWENSOpenFASTWrappers_Options"] = modopt["OWENS"]["OWENSOpenFASTWrappers"] + FileTools.save_yaml(outdir="./", fname="OWENS_Opt.yml", data_out=self.owens_modeling_options) + + + # self.analysis_type = modopt["OWENS"]["general"]["analysisType"] + # self.turbine_type = modopt["OWENS"]["general"]["turbineType"] + # self.control_strategy = modopt["OWENS"]["controlParameters"]["controlStrategy"] + # self.aeroModel = modopt["OWENS"]["AeroParameters"]["AeroModel"] + # self.adi_lib = modopt["OWENS"]["AeroParameters"]["adi_lib"] + # self.adi_rootname = modopt["OWENS"]["AeroParameters"]["adi_rootname"] + # self.Nslices = modopt["OWENS"]["AeroParameters"]["Nslices"] + # self.ntheta = modopt["OWENS"]["AeroParameters"]["ntheta"] + # # self.NumadSpec = modopt["OWENS"]["NumadSpec"] # All files go here + # # self.Turbulence = modopt["OWENS"]["TurbulenceSpec"] # All files go here + # self.ifw = modopt["OWENS"]["Inflow"]["ifw"] + # self.windType = modopt["OWENS"]["Inflow"]["windType"] + # self.windINPfilename = modopt["OWENS"]["Inflow"]["windINPfilename"] + # self.ifw_libfiles = modopt["OWENS"]["Inflow"]["ifw_libfiles"] + # self.n_blades = modopt["assembly"]["number_of_blades"] + # self.structuralModel = modopt["OWENS"]["structuralParameters"]["structuralModel"] + # self.structuralNonlinear = modopt["OWENS"]["structuralParameters"]["nonlinear"] + # self.run_path = modopt["OWENS"]["general"]["run_path"] # What exactly is this path? + + # # set other modeling options + # self.numTS = modopt["OWENS"]["general"]["numTS"] + # self.delta_t = modopt["OWENS"]["general"]["delta_t"] + # self.RPM = modopt["OWENS"]["controlParameters"]["RPM"] + + + # # Reinitialize the model with the inputs from modeling options + # self.initialize_model() # number_of_grid_pts = modopt["number_of_grid_pts"] @@ -142,14 +168,8 @@ def setup(self): # Solver options (come from modeling options) # Control inputs - if self.control_strategy == "constantRPM": - self.add_input("RPM", val=17.2, desc="RPM") - self.add_input("numTS", val=20) - self.add_input("delta_t", val=0.01, units="s") - - # Aero parameters - self.add_input("NSlices", val=30, desc="number of VAWTAero discritizations") - self.add_input("ntheta", val=30, desc="number of VAWTAero azimuthal discretizations") + if self.owens_modeling_options["OWENS_Options"]["controlStrategy"] == "prescribedRPM": + self.add_input("RPM", val=15.0, desc="RPM") # Environmental conditions # Maybe this fits better into load cases? @@ -303,10 +323,13 @@ def initialize_model(self): self.model.analysisType = self.analysis_type self.model.turbineType = self.turbine_type self.model.controlStrategy = self.control_strategy - self.model.AModel = self.aeroModel + self.model.RPM = self.RPM + self.model.AeroModel = self.aeroModel self.model.adi_lib = self.adi_lib self.model.adi_rootname = self.adi_rootname self.model.Nbld = int(self.n_blades) + self.model.Nslices = self.Nslices + self.model.ntheta = self.ntheta # This live in initialize model # TODO: If at any point we want to change the material parameters, we need to take this out to compute @@ -327,6 +350,10 @@ def initialize_model(self): self.model.structuralModel = self.structuralModel # self.model.nonlinear = self.structuralNonlinear + # Simulation time + self.model.numTS = self.numTS + self.model.delta_t = self.delta_t + # def update_model_with_yaml(self, inputs): # TODO: this is for updating the model using the yaml files, will be called in compute # First take the inputs and write out the yaml file @@ -372,37 +399,19 @@ def compute(self, inputs, outputs): material_name = modopt["materials"]["mat_name"] - path = self.run_path - owensopt = modopt["OWENS"] - eta = owensopt["eta"] + path = modopt["OWENS"]["general"]["run_path"] number_of_blades = int(inputs["Nbld"][0]) + + # Below numbers should not be set separately from blade shape and tower shape (ref_axis) Blade_Radius = inputs["Blade_Radius"][0] Blade_Height = inputs["Blade_Height"][0] towerHeight = inputs["towerHeight"][0] - # Blade dicretization - blade_x = inputs["blade_ref_axis"][:,0] - - # The reference axis is usually the same between outer shape bem and internal sturctgure fem - blade_y = inputs["blade_ref_axis"][:,1] - - # The reference axis is usually the same between outer shape bem and internal sturctgure fem - blade_z = inputs["blade_ref_axis"][:,2] rho = inputs["rho"][0] Vinf = inputs["Vinf"][0] RPM = inputs["RPM"][0] - numTS = int(inputs["numTS"][0]) - delta_t = inputs["delta_t"][0] - - Nslices = int(inputs["NSlices"][0]) - ntheta = int(inputs["ntheta"][0]) - - ntelem = int(owensopt["ntelem"]) - nbelem = int(owensopt["nbelem"]) - ncelem = int(owensopt["ncelem"]) - nselem = int(owensopt["nselem"]) # TODO: depending on the owens_yaml option, we can either update the model options directly, or write the intermediate yaml # and then update the model inputs @@ -410,14 +419,12 @@ def compute(self, inputs, outputs): # self.update_model_with_yaml(inputs): # else update the model directly as follows: - AModel = self.aeroModel - mesh_type = self.turbine_type # Assemble the ordered dictionaries for the readNumad # NuMad_geom_xlscsv_file_twr should contain everything under tower from the yaml inputs # ---------- blade structure -------------- - nd_blade_grid = inputs["tower_grid"] # WEIS gives non-dimensional grid + nd_blade_grid = inputs["blade_structure_grid"] # WEIS gives non-dimensional grid blade_accum_distances = arc_length(inputs["blade_ref_axis"]) # OWENS wants dimensional grid blade_grid = blade_accum_distances # it should be the same too nd_blade_grid* np.maximum(blade_accum_distances) @@ -425,17 +432,18 @@ def compute(self, inputs, outputs): blade_geo_dict = {} blade_geo_dict["outer_shape_bem"] = {} + blade_geo_dict["outer_shape_bem"]["blade_mountpoint"] = modopt["OWENS"]["blade"]["blade_mountpoint"] blade_geo_dict["outer_shape_bem"]["airfoil_position"] = {} - blade_geo_dict["outer_shape_bem"]["airfoil_position"]["grid"] = blade_grid_rescale # not sure how those numbers from + blade_geo_dict["outer_shape_bem"]["airfoil_position"]["grid"] = nd_blade_grid # not sure how those numbers from blade_geo_dict["outer_shape_bem"]["airfoil_position"]["labels"] = rotorse_options["airfoil_labels"] blade_geo_dict["outer_shape_bem"]["chord"] = {} - blade_geo_dict["outer_shape_bem"]["chord"]["grid"] = blade_grid_rescale + blade_geo_dict["outer_shape_bem"]["chord"]["grid"] = nd_blade_grid blade_geo_dict["outer_shape_bem"]["chord"]["values"] = inputs["blade_chord_values"] blade_geo_dict["outer_shape_bem"]["twist"] = {} - blade_geo_dict["outer_shape_bem"]["twist"]["grid"] = blade_grid_rescale + blade_geo_dict["outer_shape_bem"]["twist"]["grid"] = nd_blade_grid blade_geo_dict["outer_shape_bem"]["twist"]["values"] = inputs["blade_twist_values"] blade_geo_dict["outer_shape_bem"]["pitch_axis"] = {} - blade_geo_dict["outer_shape_bem"]["pitch_axis"]["grid"] = blade_grid_rescale + blade_geo_dict["outer_shape_bem"]["pitch_axis"]["grid"] = nd_blade_grid blade_geo_dict["outer_shape_bem"]["pitch_axis"]["values"] = inputs["blade_pitch_axis_values"] # put the outer shape bem here but OWENS does not use outer shape bem reference axis # OWENS uses that from the internal structure 2d fem @@ -443,14 +451,14 @@ def compute(self, inputs, outputs): blade_geo_dict["outer_shape_bem"]["reference_axis"]["x"] = {} blade_geo_dict["outer_shape_bem"]["reference_axis"]["y"] = {} blade_geo_dict["outer_shape_bem"]["reference_axis"]["z"] = {} - blade_geo_dict["outer_shape_bem"]["reference_axis"]["x"]["grid"] = blade_grid_rescale - blade_geo_dict["outer_shape_bem"]["reference_axis"]["y"]["grid"] = blade_grid_rescale - blade_geo_dict["outer_shape_bem"]["reference_axis"]["z"]["grid"] = blade_grid_rescale + blade_geo_dict["outer_shape_bem"]["reference_axis"]["x"]["grid"] = nd_blade_grid + blade_geo_dict["outer_shape_bem"]["reference_axis"]["y"]["grid"] = nd_blade_grid + blade_geo_dict["outer_shape_bem"]["reference_axis"]["z"]["grid"] = nd_blade_grid blade_geo_dict["outer_shape_bem"]["reference_axis"]["x"]["values"] = inputs["blade_ref_axis"][:,0] blade_geo_dict["outer_shape_bem"]["reference_axis"]["y"]["values"] = inputs["blade_ref_axis"][:,1] blade_geo_dict["outer_shape_bem"]["reference_axis"]["z"]["values"] = inputs["blade_ref_axis"][:,2] - # The reference axis is usually the same between outer shape bem and internal sturctgure fem blade_geo_dict["outer_shape_bem"]["reference_axis"] = inputs["blade_ref_axis"] # Use the same as structural reference axis + # The reference axis is usually the same between outer shape bem and internal sturctgure fem blade_geo_dict["outer_shape_bem"]["reference_axis"] = inputs["blade_ref_axis"] # Use the same as structural reference axis blade_geo_dict["internal_structure_2d_fem"] = {} blade_geo_dict["internal_structure_2d_fem"]["reference_axis"] = {} blade_geo_dict["internal_structure_2d_fem"]["reference_axis"]["x"] = {} @@ -471,10 +479,10 @@ def compute(self, inputs, outputs): blade_i = {} blade_i["name"] = blade_web_names[i] blade_i["start_nd_arc"] = {} - blade_i["start_nd_arc"]["grid"] = blade_grid_rescale + blade_i["start_nd_arc"]["grid"] = nd_blade_grid blade_i["start_nd_arc"]["values"] = inputs["blade_web_start_nd_arc"][i,:] blade_i["end_nd_arc"] = {} - blade_i["end_nd_arc"]["grid"] = blade_grid_rescale + blade_i["end_nd_arc"]["grid"] = nd_blade_grid blade_i["end_nd_arc"]["values"] = inputs["blade_web_end_nd_arc"][i,:] blade_geo_dict["internal_structure_2d_fem"]["webs"].append(blade_i) @@ -488,10 +496,10 @@ def compute(self, inputs, outputs): else: blade_i["start_nd_arc"] = {} blade_i["start_nd_arc"]["values"] = inputs["blade_layer_start_nd_arc"][i,:] - blade_i["start_nd_arc"]["grid"] = blade_grid_rescale + blade_i["start_nd_arc"]["grid"] = nd_blade_grid blade_i["end_nd_arc"] = {} blade_i["end_nd_arc"]["values"] = inputs["blade_layer_end_nd_arc"][i,:] - blade_i["end_nd_arc"]["grid"] = blade_grid_rescale + blade_i["end_nd_arc"]["grid"] = nd_blade_grid # loop to find the material ply thickness and compute the n_plies for OWENS for m, name in enumerate(material_name): @@ -501,10 +509,10 @@ def compute(self, inputs, outputs): blade_i["n_plies"] = {} blade_i["n_plies"]["values"] = n_plies - blade_i["n_plies"]["grid"] = blade_grid_rescale + blade_i["n_plies"]["grid"] = nd_blade_grid blade_i["fiber_orientation"] = {} blade_i["fiber_orientation"]["values"] = inputs["blade_layer_fiber_orientation"][i,:] - blade_i["fiber_orientation"]["grid"] = blade_grid_rescale + blade_i["fiber_orientation"]["grid"] = nd_blade_grid blade_geo_dict["internal_structure_2d_fem"]["layers"].append(blade_i) @@ -520,24 +528,24 @@ def compute(self, inputs, outputs): tower_geo_dict = {} tower_geo_dict["outer_shape_bem"] = {} tower_geo_dict["outer_shape_bem"]["airfoil_position"] = {} - tower_geo_dict["outer_shape_bem"]["airfoil_position"]["grid"] = tower_grid + tower_geo_dict["outer_shape_bem"]["airfoil_position"]["grid"] = nd_tower_grid tower_geo_dict["outer_shape_bem"]["airfoil_position"]["labels"] = n_height_tower *["Circular"] tower_geo_dict["outer_shape_bem"]["chord"] = {} - tower_geo_dict["outer_shape_bem"]["chord"]["grid"] = tower_grid + tower_geo_dict["outer_shape_bem"]["chord"]["grid"] = nd_tower_grid tower_geo_dict["outer_shape_bem"]["chord"]["values"] = tower_diameter tower_geo_dict["outer_shape_bem"]["twist"] = {} - tower_geo_dict["outer_shape_bem"]["twist"]["grid"] = tower_grid + tower_geo_dict["outer_shape_bem"]["twist"]["grid"] = nd_tower_grid tower_geo_dict["outer_shape_bem"]["twist"]["values"] = np.zeros(len(nd_tower_grid)) tower_geo_dict["outer_shape_bem"]["pitch_axis"] = {} - tower_geo_dict["outer_shape_bem"]["pitch_axis"]["grid"] = tower_grid + tower_geo_dict["outer_shape_bem"]["pitch_axis"]["grid"] = nd_tower_grid tower_geo_dict["outer_shape_bem"]["pitch_axis"]["values"] = 0.5*np.ones(len(nd_tower_grid)) tower_geo_dict["outer_shape_bem"]["reference_axis"] = {} tower_geo_dict["outer_shape_bem"]["reference_axis"]["x"] = {} tower_geo_dict["outer_shape_bem"]["reference_axis"]["y"] = {} tower_geo_dict["outer_shape_bem"]["reference_axis"]["z"] = {} - tower_geo_dict["outer_shape_bem"]["reference_axis"]["x"]["grid"] = tower_grid - tower_geo_dict["outer_shape_bem"]["reference_axis"]["y"]["grid"] = tower_grid - tower_geo_dict["outer_shape_bem"]["reference_axis"]["z"]["grid"] = tower_grid + tower_geo_dict["outer_shape_bem"]["reference_axis"]["x"]["grid"] = nd_tower_grid + tower_geo_dict["outer_shape_bem"]["reference_axis"]["y"]["grid"] = nd_tower_grid + tower_geo_dict["outer_shape_bem"]["reference_axis"]["z"]["grid"] = nd_tower_grid tower_geo_dict["outer_shape_bem"]["reference_axis"]["x"]["values"] = inputs["tower_ref_axis"][:,0] tower_geo_dict["outer_shape_bem"]["reference_axis"]["y"]["values"] = inputs["tower_ref_axis"][:,1] tower_geo_dict["outer_shape_bem"]["reference_axis"]["z"]["values"] = inputs["tower_ref_axis"][:,2] @@ -547,9 +555,9 @@ def compute(self, inputs, outputs): tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["x"] = {} tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["y"] = {} tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["z"] = {} - tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["x"]["grid"] = tower_grid - tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["y"]["grid"] = tower_grid - tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["z"]["grid"] = tower_grid + tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["x"]["grid"] = nd_tower_grid + tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["y"]["grid"] = nd_tower_grid + tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["z"]["grid"] = nd_tower_grid tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["x"]["values"] = inputs["tower_ref_axis"][:,0] tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["y"]["values"] = inputs["tower_ref_axis"][:,1] tower_geo_dict["internal_structure_2d_fem"]["reference_axis"]["z"]["values"] = inputs["tower_ref_axis"][:,2]# Not needed in readNuMad @@ -570,10 +578,10 @@ def compute(self, inputs, outputs): tower_i["material"] = tower_layer_materials[i] tower_i["start_nd_arc"] = {} tower_i["start_nd_arc"]["values"] = np.zeros(n_height_tower) - tower_i["start_nd_arc"]["grid"] = tower_grid + tower_i["start_nd_arc"]["grid"] = nd_tower_grid tower_i["end_nd_arc"] = {} tower_i["end_nd_arc"]["values"] = np.ones(n_height_tower) - tower_i["end_nd_arc"]["grid"] = tower_grid + tower_i["end_nd_arc"]["grid"] = nd_tower_grid # loop to find the material ply thickness and compute the n_plies for OWENS for m, name in enumerate(material_name): @@ -583,10 +591,10 @@ def compute(self, inputs, outputs): tower_i["n_plies"] = {} tower_i["n_plies"]["values"] = n_plies - tower_i["n_plies"]["grid"] = tower_grid + tower_i["n_plies"]["grid"] = nd_tower_grid tower_i["fiber_orientation"] = {} tower_i["fiber_orientation"]["values"] = np.zeros(n_height_tower) - tower_i["fiber_orientation"]["grid"] = tower_grid + tower_i["fiber_orientation"]["grid"] = nd_tower_grid tower_geo_dict["internal_structure_2d_fem"]["layers"].append(tower_i) @@ -600,28 +608,31 @@ def compute(self, inputs, outputs): strut_geo_dict = {} + strut_geo_dict["mountfraction_blade"] = {} + strut_geo_dict["mountfraction_blade"] = strut_options["mountfraction_blade"] + strut_geo_dict["mountfraction_tower"] = strut_options["mountfraction_tower"] strut_geo_dict["outer_shape_bem"] = {} strut_geo_dict["outer_shape_bem"]["airfoil_position"] = {} # strut grid in OWENS yaml is dimensional but it doesn't matter # the strut dimension is determined by the strut mountpoint to tower and blade - strut_geo_dict["outer_shape_bem"]["airfoil_position"]["grid"] = strut_grid_rescale + strut_geo_dict["outer_shape_bem"]["airfoil_position"]["grid"] = nd_strut_grid strut_geo_dict["outer_shape_bem"]["airfoil_position"]["labels"] = strut_options["airfoils"] strut_geo_dict["outer_shape_bem"]["chord"] = {} - strut_geo_dict["outer_shape_bem"]["chord"]["grid"] = strut_grid_rescale + strut_geo_dict["outer_shape_bem"]["chord"]["grid"] = nd_strut_grid strut_geo_dict["outer_shape_bem"]["chord"]["values"] = inputs["strut_chord"] strut_geo_dict["outer_shape_bem"]["twist"] = {} - strut_geo_dict["outer_shape_bem"]["twist"]["grid"] = strut_grid_rescale + strut_geo_dict["outer_shape_bem"]["twist"]["grid"] = nd_strut_grid strut_geo_dict["outer_shape_bem"]["twist"]["values"] = inputs["strut_twist"] strut_geo_dict["outer_shape_bem"]["pitch_axis"] = {} - strut_geo_dict["outer_shape_bem"]["pitch_axis"]["grid"] = strut_grid_rescale + strut_geo_dict["outer_shape_bem"]["pitch_axis"]["grid"] = nd_strut_grid strut_geo_dict["outer_shape_bem"]["pitch_axis"]["values"] = inputs["strut_pitch_axis"] strut_geo_dict["outer_shape_bem"]["reference_axis"] = {} strut_geo_dict["outer_shape_bem"]["reference_axis"]["x"] = {} strut_geo_dict["outer_shape_bem"]["reference_axis"]["y"] = {} strut_geo_dict["outer_shape_bem"]["reference_axis"]["z"] = {} - strut_geo_dict["outer_shape_bem"]["reference_axis"]["x"]["grid"] = strut_grid_rescale - strut_geo_dict["outer_shape_bem"]["reference_axis"]["y"]["grid"] = strut_grid_rescale - strut_geo_dict["outer_shape_bem"]["reference_axis"]["z"]["grid"] = strut_grid_rescale + strut_geo_dict["outer_shape_bem"]["reference_axis"]["x"]["grid"] = nd_strut_grid + strut_geo_dict["outer_shape_bem"]["reference_axis"]["y"]["grid"] = nd_strut_grid + strut_geo_dict["outer_shape_bem"]["reference_axis"]["z"]["grid"] = nd_strut_grid strut_geo_dict["outer_shape_bem"]["reference_axis"]["x"]["values"] = inputs["strut_ref_axis"][:,0] strut_geo_dict["outer_shape_bem"]["reference_axis"]["y"]["values"] = inputs["strut_ref_axis"][:,1] strut_geo_dict["outer_shape_bem"]["reference_axis"]["z"]["values"] = inputs["strut_ref_axis"][:,2] @@ -633,9 +644,9 @@ def compute(self, inputs, outputs): strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["x"] = {} strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["y"] = {} strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["z"] = {} - strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["x"]["grid"] = strut_grid_rescale - strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["y"]["grid"] = strut_grid_rescale - strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["z"]["grid"] = strut_grid_rescale + strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["x"]["grid"] = nd_strut_grid + strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["y"]["grid"] = nd_strut_grid + strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["z"]["grid"] = nd_strut_grid strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["x"]["values"] = inputs["strut_ref_axis"][:,0] strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["y"]["values"] = inputs["strut_ref_axis"][:,1] strut_geo_dict["internal_structure_2d_fem"]["reference_axis"]["z"]["values"] = inputs["strut_ref_axis"][:,2] @@ -657,10 +668,10 @@ def compute(self, inputs, outputs): strut_i["name"] = strut_web_names[i] strut_i["start_nd_arc"] = {} - strut_i["start_nd_arc"]["grid"] = strut_grid_rescale + strut_i["start_nd_arc"]["grid"] = nd_strut_grid strut_i["start_nd_arc"]["values"] = inputs["strut_web_start_nd_arc"][i,:] strut_i["end_nd_arc"] = {} - strut_i["end_nd_arc"]["grid"] = strut_grid_rescale + strut_i["end_nd_arc"]["grid"] = nd_strut_grid strut_i["end_nd_arc"]["values"] = inputs["strut_web_end_nd_arc"][i,:] strut_geo_dict["internal_structure_2d_fem"]["webs"].append(strut_i) @@ -673,10 +684,10 @@ def compute(self, inputs, outputs): strut_i["material"] = strut_layer_materials[i] strut_i["start_nd_arc"] = {} strut_i["start_nd_arc"]["values"] = inputs["strut_layer_start_nd_arc"][i,:] - strut_i["start_nd_arc"]["grid"] = strut_grid_rescale + strut_i["start_nd_arc"]["grid"] = nd_strut_grid strut_i["end_nd_arc"] = {} strut_i["end_nd_arc"]["values"] = inputs["strut_layer_end_nd_arc"][i,:] - strut_i["end_nd_arc"]["grid"] = strut_grid_rescale + strut_i["end_nd_arc"]["grid"] = nd_strut_grid # strut_i["n_plies"]["values"] = inputs["strut_layer_nplies"] # loop to find the material ply thickness and compute the n_plies for OWENS for m, name in enumerate(material_name): @@ -684,10 +695,10 @@ def compute(self, inputs, outputs): n_plies = inputs["strut_layer_thickness"][i,:]/inputs["ply_t"][m] strut_i["n_plies"] = {} strut_i["n_plies"]["values"] = n_plies - strut_i["n_plies"]["grid"] = strut_grid_rescale + strut_i["n_plies"]["grid"] = nd_strut_grid strut_i["fiber_orientation"] = {} strut_i["fiber_orientation"]["values"] = inputs["strut_layer_fiber_orientation"][i,:] - strut_i["fiber_orientation"]["grid"] = strut_grid_rescale + strut_i["fiber_orientation"]["grid"] = nd_strut_grid strut_geo_dict["internal_structure_2d_fem"]["layers"].append(strut_i) # materials dict @@ -730,7 +741,8 @@ def compute(self, inputs, outputs): # jl_ordered_dict = convert(jl.OrderedDict, yaml_dict) FileTools.save_yaml(outdir="/Users/yliao/repos/WEIS/examples/18_owens", fname="data.yml", data_out=yaml_dict) - jl.OWENS.runOWENSWINDIO("/Users/yliao/repos/WEIS/examples/18_owens/data.yml",self.model,path,verbosity=2) + + jl.OWENS.runOWENSWINDIO("./OWENS_Opt.yml", "/Users/yliao/repos/WEIS/examples/18_owens/data.yml",path) # setup_outputs = jl.OWENS.setupOWENS(jl.OWENSAero, path, # rho=rho, @@ -945,6 +957,7 @@ def compute(self, inputs, outputs): output_h5 = OWENSOutput(output_path, output_channels=["t", "FReactionHist", "massOwens", "topDamage_blade_U", "SF_ult_U"]) massOwens = output_h5["massOwens"] FReactionHist = output_h5["FReactionHist"] + # print("shape of FReactionHist: ", FReactionHist.shape) topDamage_blade_U = output_h5["topDamage_blade_U"] SF_ult_U= output_h5["SF_ult_U"] t = output_h5["t"] @@ -965,6 +978,7 @@ def compute(self, inputs, outputs): # # Other outputs for constraints outputs["SF"] = minSF outputs["fatigue_damage"] = maxFatiguePer20yr # - 1.0 + print("power: ", outputs["power"]) print("fatigue damage: ", outputs["fatigue_damage"]) # # power constraint can be imposed elsewhere # # since it is already an output diff --git a/weis/owens/unsteady_example.py b/weis/owens/unsteady_example.py index 876754520..90b1fcc2c 100644 --- a/weis/owens/unsteady_example.py +++ b/weis/owens/unsteady_example.py @@ -1,6 +1,6 @@ import numpy as np import openmdao.api as om -from openmdao_owens import * +from openmdao_owens_standalone import * from scipy.interpolate import PchipInterpolator, Akima1DInterpolator class discretization_x(om.ExplicitComponent): @@ -28,7 +28,8 @@ def compute(self, inputs, outputs): max_radius = inputs["max_radius"][0] x_control_values = np.array([0, quarter_span, max_radius, quarter_span, 0]) - blade_x = Akima1DInterpolator(x_control_pts_grid, x_control_values)(x_grid) + # blade_x = Akima1DInterpolator(x_control_pts_grid, x_control_values)(x_grid) + blade_x = np.array([0.0, 6.961447494399997, 13.442795161599998, 19.444043001599994, 24.965191014399995, 30.006239199999996, 34.5671875584, 38.6480360896, 42.248784793599995, 45.3694336704, 48.00998272, 50.1704319424, 51.8507813376, 53.0510309056, 53.7711806464, 54.01123056, 53.7711806464, 53.0510309056, 51.8507813376, 50.1704319424, 48.009982720000004, 45.36943367040001, 42.24878479360001, 38.6480360896, 34.56718755839999, 30.006239199999996, 24.965191014399995, 19.444043001599994, 13.442795161599998, 6.961447494399997, 0.0]) outputs["blade_x"] = blade_x @@ -69,7 +70,7 @@ def compute(self, inputs, outputs): } model.add_subsystem("discretization", discretization_x(number_of_control_pts=5, number_of_grid_pts=len(x)), promotes=["*"]) -model.add_subsystem("crossflow", OWENSUnsteadySetup(modeling_options=options), promotes=["*"]) +model.add_subsystem("crossflow", OWENSUnsteadySetup2(modeling_options=options), promotes=["*"]) # Note: Change the return line in topRunDLC in OWENS to "return mass_breakout_twr, genPower, massOwens" for this example to work @@ -91,7 +92,7 @@ def compute(self, inputs, outputs): print("quarter_span of 24 m and max radius of 42 m") print("Mean power is : "+str(prob.get_val("power"))) print("lcoe is : "+str(prob.get_val("lcoe"))) -# exit() +exit() # run opt prob.driver = om.pyOptSparseDriver() prob.driver.options["optimizer"] = "SLSQP"