From 32bcdc6543860eac70a8298960c479486bad1935 Mon Sep 17 00:00:00 2001 From: Ernesto Camarena Date: Thu, 19 Oct 2023 09:17:41 -0600 Subject: [PATCH 01/11] fixing ansys example --- src/pynumad/analysis/ansys/run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pynumad/analysis/ansys/run.py b/src/pynumad/analysis/ansys/run.py index 43396eb..2b7faa1 100644 --- a/src/pynumad/analysis/ansys/run.py +++ b/src/pynumad/analysis/ansys/run.py @@ -5,7 +5,7 @@ def call_ansys(script_name,log=False,script_out='ansys.out',ncpus=1): for f in glob.glob("*.lock"): os.remove(f) - ansys_path=SOFTWARE_PATHS['ansys_path'] + ansys_path=SOFTWARE_PATHS['ansys'] MAXnLicenceTries=100 try: From 628fb7fe7167701b92334801cbccabb2757e36ad Mon Sep 17 00:00:00 2001 From: Ernesto Camarena Date: Thu, 19 Oct 2023 10:30:22 -0600 Subject: [PATCH 02/11] making sure ansys and cubit example work. adding to gitignore. adding to software_paths --- .gitignore | 19 +++++++++++++++++++ docs/user-guide/meshing.rst | 2 +- examples/cubit_beam.py | 19 +++++++++++++++---- examples/cubit_solid.py | 4 ++-- src/pynumad/analysis/make_models.py | 12 ++++++------ src/pynumad/software_paths.json | 3 ++- 6 files changed, 45 insertions(+), 14 deletions(-) diff --git a/.gitignore b/.gitignore index 82fa191..034f892 100644 --- a/.gitignore +++ b/.gitignore @@ -183,3 +183,22 @@ cython_debug/ # pynumad src/pynumad/software_paths.json + +#VABS +*.in +*.in.* + +#Cubit +*.cub +*.g +*.jou + +#BeamDyn +*BeamDyn*dat + +#Sierra +sd.i +sm.i +euler +mat_ori.py + diff --git a/docs/user-guide/meshing.rst b/docs/user-guide/meshing.rst index c42d9eb..b387f62 100644 --- a/docs/user-guide/meshing.rst +++ b/docs/user-guide/meshing.rst @@ -76,7 +76,7 @@ Then add Cubit to the path import sys sys.path.append(pynumad.SOFTWARE_PATHS['cubit']) - sys.path.append(pynumad.SOFTWARE_PATHS['cubitEnhancements']) + sys.path.append(pynumad.SOFTWARE_PATHS['cubit_enhancements']) import cubit diff --git a/examples/cubit_beam.py b/examples/cubit_beam.py index bcd27c6..cd24876 100644 --- a/examples/cubit_beam.py +++ b/examples/cubit_beam.py @@ -2,14 +2,14 @@ import pynumad sys.path.append(pynumad.SOFTWARE_PATHS['cubit']) -sys.path.append(pynumad.SOFTWARE_PATHS['cubitEnhancements']) +sys.path.append(pynumad.SOFTWARE_PATHS['cubit_enhancements']) import cubit from pynumad.analysis.cubit.make_blade import * import numpy as np from pynumad.analysis.make_models import write_beam_model - +import logging def get_cs_params(): @@ -75,12 +75,23 @@ def get_cs_params(): cs_params=get_cs_params() settings={} -settings['make_input_for']='VABS' #SM, VABS, ANBA, or None +settings['make_input_for']='VABSbeamdyn' #SM, VABS, ANBA, or None settings['export']='cubg' #cub, g, or None cubit_make_cross_sections(blade,wt_name,settings,cs_params,'2D',stationList=[],directory=dirName) #Note an empty list for stationList will make all cross sections. - +#Proportional damping values for BeamDyn file. mu=[0.00257593, 0.0017469, 0.0017469, 0.0017469, 0.00257593, 0.0017469] + + +#Set up logging +log =logging.getLogger(__name__) +log.setLevel(logging.DEBUG) +fh=logging.FileHandler(wt_name+'driver.log',mode='w') +log.addHandler(fh) + +#Read in a fresh new blade +blade=pynumad.Blade() +blade.read_yaml('example_data/'+yamlName+'.yaml') file_names=write_beam_model(wt_name,settings,blade,mu,log,directory=dirName) \ No newline at end of file diff --git a/examples/cubit_solid.py b/examples/cubit_solid.py index e387a3c..ff2d482 100644 --- a/examples/cubit_solid.py +++ b/examples/cubit_solid.py @@ -2,7 +2,7 @@ import pynumad sys.path.append(pynumad.SOFTWARE_PATHS['cubit']) -sys.path.append(pynumad.SOFTWARE_PATHS['cubitEnhancements']) +sys.path.append(pynumad.SOFTWARE_PATHS['cubit_enhancements']) import cubit from pynumad.analysis.cubit.make_blade import * @@ -75,7 +75,7 @@ def get_cs_params(): cs_params=get_cs_params() settings={} -settings['make_input_for']='Sd' #SM, VABS, ANBA, or None +settings['make_input_for']='SmSd' #SM, VABS, ANBA, or None settings['export']='cubg' #cub, g, or None materials_used=cubit_make_solid_blade(blade, wt_name, settings, cs_params, stationList=[2,3]) diff --git a/src/pynumad/analysis/make_models.py b/src/pynumad/analysis/make_models.py index 1b27ccb..b250a78 100644 --- a/src/pynumad/analysis/make_models.py +++ b/src/pynumad/analysis/make_models.py @@ -4,7 +4,7 @@ import glob import numpy as np from pynumad.utils.misc_utils import copy_and_replace - +from pynumad.paths import SOFTWARE_PATHS def write_beam_model(wt_name,settings,blade,mu,log,directory='.'): import pynumad.analysis.beam_utils as beam_utils @@ -31,7 +31,7 @@ def write_beam_model(wt_name,settings,blade,mu,log,directory='.'): for filePath in glob.glob(directory+'/'+wt_name+'*.in'): fileCount+=1 try: - this_cmd = 'VABS ' +filePath + this_cmd = SOFTWARE_PATHS['vabs']+' ' +filePath log.info(f' running: {this_cmd}') licenseAvailable=False @@ -100,7 +100,7 @@ def write_sierra_model(wt_name,settings,blade,materials_used,directory='.'): # #Runs VABS or OpenSG to homogenize # #Makes beamDyn or GEBT files - + template_path=SOFTWARE_PATHS['pynumad']+'src/data/templates/' # radial_stations=blade.ispan/blade.ispan[-1] # nStations=len(radial_stations) @@ -109,7 +109,7 @@ def write_sierra_model(wt_name,settings,blade,materials_used,directory='.'): materials = blade.definition.materials solver_string=settings['make_input_for'].lower() if 'sm' in solver_string or 'sd' in solver_string: - templateFileName='mat_ori.py.template' + templateFileName=template_path+'mat_ori.py.template' mat_ori_file_name='mat_ori.py' copy_and_replace(templateFileName, mat_ori_file_name, @@ -118,7 +118,7 @@ def write_sierra_model(wt_name,settings,blade,materials_used,directory='.'): }) if 'sm' in settings['make_input_for'].lower(): - templateFileName='sm.i.template' + templateFileName=template_path+'sm.i.template' adagioFileName='sm.i' materialLines=f'' @@ -167,7 +167,7 @@ def write_sierra_model(wt_name,settings,blade,materials_used,directory='.'): 'BLADE_BLOCKS': blockLines, }) if 'sd' in settings['make_input_for'].lower(): - templateFileName='sd.i.template' + templateFileName=template_path+'sd.i.template' adagioFileName='sd.i' materialLines=f'' diff --git a/src/pynumad/software_paths.json b/src/pynumad/software_paths.json index 755693e..8a18e88 100644 --- a/src/pynumad/software_paths.json +++ b/src/pynumad/software_paths.json @@ -1,6 +1,7 @@ { - "numad": "", + "pynumad": "", "ansys": "", + "vabs": "", "cubit": "", "cubit_enhancements": "" } From 415089617e22b7e21b04d4d1226a29ac73316afc Mon Sep 17 00:00:00 2001 From: Camarena Date: Tue, 10 Sep 2024 08:54:32 -0600 Subject: [PATCH 03/11] updated buckling to create .g file for sd and updated sd input file for new mat ori implementation in sd --- src/pynumad/analysis/cubit/make_blade.py | 99 ++++++++++++++++++++---- src/pynumad/analysis/make_models.py | 2 + 2 files changed, 88 insertions(+), 13 deletions(-) diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index 286fe38..a8954d3 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -185,6 +185,42 @@ def get_hex_orientations_two_points(volume_id): return global_el_ids_in_vol,spanwise_directions_in_vol,perimeter_directions_in_vol +def get_hex_orientations_vectors(volume_id): + global_el_ids_in_vol=[] + + + + spanwise_directions_in_vol = [] + perimeter_directions_in_vol = [] + surface_normal_directions_in_vol = [] + + surf_id_for_mat_ori,sign = get_mat_ori_surface(volume_id) + #volume_name = cubit.get_entity_name("volume", volume_id) + #t0 = time.time() + + for el_id in get_volume_hexes(volume_id): + coords = cubit.get_center_point("hex", el_id) + + surface_normal = vectNorm( + list(sign*np.array(get_surface_normal_at_coord(surf_id_for_mat_ori, coords)))) + + ref_line_direction = [0,0,1] + #https://www.maplesoft.com/support/help/maple/view.aspx?path=MathApps%2FProjectionOfVectorOntoPlane + spanwise_direction = vectNorm(np.array(ref_line_direction)-np.dot(ref_line_direction,surface_normal)*np.array(surface_normal)) + + perimeter_direction = vectNorm(np.cross(surface_normal, spanwise_direction)) + + global_id=get_global_element_id('hex',el_id) + + global_el_ids_in_vol.append(global_id) + + + spanwise_directions_in_vol.append(spanwise_direction) + perimeter_directions_in_vol.append(perimeter_direction) + surface_normal_directions_in_vol.append(surface_normal) + + return global_el_ids_in_vol,spanwise_directions_in_vol,perimeter_directions_in_vol,surface_normal_directions_in_vol + def get_tet_orientations(volume_id): global_el_ids_in_vol=[] theta1s_in_vol=[] @@ -231,21 +267,40 @@ def get_tet_orientations(volume_id): return global_el_ids_in_vol,theta1s_in_vol,theta2s_in_vol,theta3s_in_vol,spanwise_directions_in_vol,perimeter_directions_in_vol -def assign_material_orientations(orientation_data): +def assign_material_orientations(orientation_data,output_format = 'euler'): #Apply Material Orientation global_ids=orientation_data[0] n_el = len(global_ids) - theta1s=orientation_data[1] - theta2s=orientation_data[2] - theta3s=orientation_data[3] - cubit.set_element_variable(global_ids, 'rotation_angle_one', theta1s) - cubit.set_element_variable(global_ids, 'rotation_angle_two', theta2s) - cubit.set_element_variable(global_ids, 'rotation_angle_three', theta3s) + if output_format =='euler': + theta1s=orientation_data[1] + theta2s=orientation_data[2] + theta3s=orientation_data[3] + + cubit.set_element_variable(global_ids, 'rotation_angle_one', theta1s) + cubit.set_element_variable(global_ids, 'rotation_angle_two', theta2s) + cubit.set_element_variable(global_ids, 'rotation_angle_three', theta3s) + + cubit.set_element_variable(global_ids, 'rotation_axis_one', 1*np.ones(n_el)) + cubit.set_element_variable(global_ids, 'rotation_axis_two', 2*np.ones(n_el)) + cubit.set_element_variable(global_ids, 'rotation_axis_three', 3*np.ones(n_el)) + elif output_format =='vectors': + + one_axis=np.array(orientation_data[1]) + two_axis=np.array(orientation_data[2]) + three_axis=np.array(orientation_data[3]) + + cubit.set_element_variable(global_ids, 'matCoord_1_x', one_axis[:,0]) + cubit.set_element_variable(global_ids, 'matCoord_1_y', one_axis[:,1]) + cubit.set_element_variable(global_ids, 'matCoord_1_z', one_axis[:,2]) - cubit.set_element_variable(global_ids, 'rotation_axis_one', 1*np.ones(n_el)) - cubit.set_element_variable(global_ids, 'rotation_axis_two', 2*np.ones(n_el)) - cubit.set_element_variable(global_ids, 'rotation_axis_three', 3*np.ones(n_el)) + cubit.set_element_variable(global_ids, 'matCoord_2_x', two_axis[:,0]) + cubit.set_element_variable(global_ids, 'matCoord_2_y', two_axis[:,1]) + cubit.set_element_variable(global_ids, 'matCoord_2_z', two_axis[:,2]) + + cubit.set_element_variable(global_ids, 'matCoord_3_x', three_axis[:,0]) + cubit.set_element_variable(global_ids, 'matCoord_3_y', three_axis[:,1]) + cubit.set_element_variable(global_ids, 'matCoord_3_z', three_axis[:,2]) return @@ -269,6 +324,9 @@ def compute_material_orientations(element_shape,output_format = 'euler',ncpus = elif 'two_points' in output_format: for vol_id in all_volume_ids: ans.append(get_hex_orientations_two_points(vol_id)) + elif 'vectors' in output_format: + for vol_id in all_volume_ids: + ans.append(get_hex_orientations_vectors(vol_id)) else: raise NameError(f'Material Orientation output format: {output_format} is not supported') else: @@ -277,6 +335,8 @@ def compute_material_orientations(element_shape,output_format = 'euler',ncpus = ans = pool_obj.map(get_hex_orientations_euler,all_volume_ids) elif 'two_points' in output_format: ans = pool_obj.map(get_hex_orientations_two_points,all_volume_ids) + elif 'vectors' in output_format: + ans = pool_obj.map(get_hex_orientations_vectors,all_volume_ids) else: raise NameError(f'Material Orientation output format: {output_format} is not supported') @@ -313,6 +373,19 @@ def compute_material_orientations(element_shape,output_format = 'euler',ncpus = return [global_ids,spanwise_directions,perimiter_directions] + elif 'vectors' in output_format: + spanwise_directions = [] + perimiter_directions = [] + normal_directions = [] + + for i in range(len(all_volume_ids)): + global_ids+=list(ans[i][0]) + spanwise_directions+=list(ans[i][1]) + perimiter_directions+=list(ans[i][2]) + normal_directions+=list(ans[i][3]) + + return [global_ids,spanwise_directions,perimiter_directions,normal_directions] + def order_path_points(points, ind): points_new = [ points.pop(ind) ] # initialize a new list of points with the known first point @@ -1373,9 +1446,9 @@ def cubit_make_solid_blade( cubit.cmd(f"nodeset {nodeset_id} add surface {l2s(surface_ids)} ") cubit.cmd(f'nodeset {nodeset_id} name "{node_set_name}_ns"') - sideset_id=cubit.get_next_sideset_id() - cubit.cmd(f"sideset {sideset_id} add surface {l2s(surface_ids)} ") - cubit.cmd(f'sideset {sideset_id} name "{node_set_name}_ss"') + # sideset_id=cubit.get_next_sideset_id() + # cubit.cmd(f"sideset {sideset_id} add surface {l2s(surface_ids)} ") + # cubit.cmd(f'sideset {sideset_id} name "{node_set_name}_ss"') # Outer mold-line sideset diff --git a/src/pynumad/analysis/make_models.py b/src/pynumad/analysis/make_models.py index 1fa54dc..1632c43 100644 --- a/src/pynumad/analysis/make_models.py +++ b/src/pynumad/analysis/make_models.py @@ -222,11 +222,13 @@ def write_sierra_sd_model(template_file,wt_name,station_list,blade,materials_use blockLines+=f'block {material.name}\n' blockLines+=f' material {material.name}\n' + blockLines+=f' coordinate from_transfer\n' blockLines+=f'end \n' blockLines+='\n\n' copy_and_replace(template_file, sierra_file_name, { + 'ROOT_STATION': 'station'+str(station_list[0]).zfill(3), 'BLADE_MATERIALS': materialLines, 'IN_MESH':wt_name+'.g', 'BLADE_BLOCKS': blockLines, From b9258c8c821b86470288a165a93042d499183d78 Mon Sep 17 00:00:00 2001 From: Camarena Date: Wed, 16 Oct 2024 11:37:46 -0600 Subject: [PATCH 04/11] adding precomp capabiliy to cubit workflow --- src/pynumad/analysis/cubit/make_blade.py | 445 +++++++++--------- .../analysis/cubit/make_cross_sections.py | 356 ++++++++++++++ 2 files changed, 581 insertions(+), 220 deletions(-) diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index a8954d3..fac32d6 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -760,253 +760,258 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati else: make_webs = hasWebs[i_station] - # Only save birds_mouth_verts for the right cross-section - if i_station == i_station_first_web: - birds_mouth_verts = make_a_cross_section(wt_name, - surface_dict, - i_station, - i_station_geometry, - blade, - make_webs, - aft_web_stack, - fore_web_stack, - iLE, - cs_params, - geometry_scaling, - thickness_scaling, - last_round_station, - last_flat_station, - materials_used, - cs_normal, - ) + if 'precomp' in settings['make_input_for']: + make_a_precomp_cross_section(wt_name, + surface_dict,i_station,i_station_geometry,blade,make_webs,aft_web_stack,fore_web_stack,iLE,cs_params, + geometry_scaling,thickness_scaling,last_round_station,last_flat_station,materials_used,cs_normal) else: - make_a_cross_section(wt_name, - surface_dict, - i_station, - i_station_geometry, - blade, - make_webs, - aft_web_stack, - fore_web_stack, - iLE, - cs_params, - geometry_scaling, - thickness_scaling, - last_round_station, - last_flat_station, - materials_used, - cs_normal, - ) - - - cubit.cmd(f"delete curve all with Is_Free except {spanwise_mat_ori_curve}") - - # Chord line for rotation of cross-section for homogenization - if model2Dor3D.lower() == "2d": - # #Blocks - if 'd_tube' in cs_params.keys() and cs_params['d_tube']: - keep_list=[] - - cubit.cmd(f'delete surface with x_coord < 0"') - cubit.cmd(f'delete surface with name "*layer9*"') - cubit.cmd(f'delete surface with name "*layer10*"') - cubit.cmd(f'delete surface with name "*layer11*"') - - delete_list=[] - parse_string = f'with name "*layer3*"' - delete_list += list(parse_cubit_list("surface", parse_string)) - parse_string = f'with name "*layer4*"' - delete_list += list(parse_cubit_list("surface", parse_string)) - - keep_list=[] - #LE - for i in [121,122,123]: - parse_string = f'with name "shell*Station*surface{i}"' - keep_list += list(parse_cubit_list("surface", parse_string)) - - #Web - for i in [1,2,3,4,5,6,17,18,19,20,21,22]: - parse_string = f'with name "web_web*surface{i}"' - keep_list += list(parse_cubit_list("surface", parse_string)) - - vol_ids=set(delete_list).difference(set(keep_list)) - - cubit.cmd(f'delete vol {l2s(vol_ids)}') - cubit.cmd(f'delete vol with name "*Station005*"') - - for imat, material_name in enumerate(materials_used): - cubit.cmd(f'block {imat+1} add surface with name "*{material_name}*"') - cubit.cmd(f'block {imat+1} name "{material_name}"') - - addColor(blade, "surface") - - # create_vertex(blade.geometry[0,0,i_station]*geometry_scaling,blade.geometry[0,1,i_station]*geometry_scaling,blade.geometry[0,2,i_station]*geometry_scaling) - # TEvert=get_last_id("vertex") - # create_vertex(blade.geometry[iLE-1,0,i_station]*geometry_scaling,blade.geometry[iLE-1,1,i_station]*geometry_scaling,blade.geometry[iLE-1,2,i_station]*geometry_scaling) - # LEvert=get_last_id("vertex") - - # cubit.cmd(f'create curve vertex {TEvert} {LEvert}') - # coords=cubit.vertex(TEvert).coordinates() - # tangent=cubit.curve(get_last_id("curve")).tangent(coords) - # tangent_direction=vectNorm(list(tangent)) #Unit vector of tangent. - # crossSectionRotationAngle=math.atan2(tangent_direction[1],tangent_direction[0])*180/pi - - parse_string = f'with name "*Station{str(i_station).zfill(3)}*"' - volume_ids = parse_cubit_list("surface", parse_string) - - # Undo initial twist - cubit.cmd( - f"rotate Surface {l2s(volume_ids)} angle {definition.degreestwist[i_station]} about Z include_merged " - ) - - # Undo prebend - if definition.prebend[i_station] != 0: - cubit.cmd(f"move surface {l2s(volume_ids)} y {-1*definition.prebend[i_station]} include_merged") - - # Undo sweep - if definition.sweep[i_station] != 0: - raise ValueError("Presweep is untested for cross-sectional meshing") - - if 'ref_line_type' in cs_params and 'centroid' in cs_params['ref_line_type'].lower(): - centroidal_vert_id, centroidal_ref_line_coords=get_locus_of_cross_sectional_centroids([i_station]) - cubit.cmd(f"move surface {l2s(volume_ids)} x {-1*centroidal_ref_line_coords[0][0]} include_merged") - cubit.cmd(f"move surface {l2s(volume_ids)} y {-1*centroidal_ref_line_coords[0][1]} include_merged") - #cubit.cmd(f"move surface {l2s(volume_ids)} z {-1*centroidal_ref_line_coords[0][2]} include_merged") + # Only save birds_mouth_verts for the right cross-section + if i_station == i_station_first_web: + birds_mouth_verts = make_a_cross_section(wt_name, + surface_dict, + i_station, + i_station_geometry, + blade, + make_webs, + aft_web_stack, + fore_web_stack, + iLE, + cs_params, + geometry_scaling, + thickness_scaling, + last_round_station, + last_flat_station, + materials_used, + cs_normal, + ) + else: + make_a_cross_section(wt_name, + surface_dict, + i_station, + i_station_geometry, + blade, + make_webs, + aft_web_stack, + fore_web_stack, + iLE, + cs_params, + geometry_scaling, + thickness_scaling, + last_round_station, + last_flat_station, + materials_used, + cs_normal, + ) - centroidal_vert_id, centroidal_ref_line_coords=get_locus_of_cross_sectional_centroids([i_station]) + cubit.cmd(f"delete curve all with Is_Free except {spanwise_mat_ori_curve}") + + # Chord line for rotation of cross-section for homogenization + if model2Dor3D.lower() == "2d": + # #Blocks + if 'd_tube' in cs_params.keys() and cs_params['d_tube']: + keep_list=[] + + cubit.cmd(f'delete surface with x_coord < 0"') + cubit.cmd(f'delete surface with name "*layer9*"') + cubit.cmd(f'delete surface with name "*layer10*"') + cubit.cmd(f'delete surface with name "*layer11*"') + + delete_list=[] + parse_string = f'with name "*layer3*"' + delete_list += list(parse_cubit_list("surface", parse_string)) + parse_string = f'with name "*layer4*"' + delete_list += list(parse_cubit_list("surface", parse_string)) + + keep_list=[] + #LE + for i in [121,122,123]: + parse_string = f'with name "shell*Station*surface{i}"' + keep_list += list(parse_cubit_list("surface", parse_string)) + + #Web + for i in [1,2,3,4,5,6,17,18,19,20,21,22]: + parse_string = f'with name "web_web*surface{i}"' + keep_list += list(parse_cubit_list("surface", parse_string)) + + vol_ids=set(delete_list).difference(set(keep_list)) + + cubit.cmd(f'delete vol {l2s(vol_ids)}') + cubit.cmd(f'delete vol with name "*Station005*"') + + for imat, material_name in enumerate(materials_used): + cubit.cmd(f'block {imat+1} add surface with name "*{material_name}*"') + cubit.cmd(f'block {imat+1} name "{material_name}"') + + addColor(blade, "surface") + + # create_vertex(blade.geometry[0,0,i_station]*geometry_scaling,blade.geometry[0,1,i_station]*geometry_scaling,blade.geometry[0,2,i_station]*geometry_scaling) + # TEvert=get_last_id("vertex") + # create_vertex(blade.geometry[iLE-1,0,i_station]*geometry_scaling,blade.geometry[iLE-1,1,i_station]*geometry_scaling,blade.geometry[iLE-1,2,i_station]*geometry_scaling) + # LEvert=get_last_id("vertex") + + # cubit.cmd(f'create curve vertex {TEvert} {LEvert}') + # coords=cubit.vertex(TEvert).coordinates() + # tangent=cubit.curve(get_last_id("curve")).tangent(coords) + # tangent_direction=vectNorm(list(tangent)) #Unit vector of tangent. + # crossSectionRotationAngle=math.atan2(tangent_direction[1],tangent_direction[0])*180/pi + + parse_string = f'with name "*Station{str(i_station).zfill(3)}*"' + volume_ids = parse_cubit_list("surface", parse_string) + + # Undo initial twist + cubit.cmd( + f"rotate Surface {l2s(volume_ids)} angle {definition.degreestwist[i_station]} about Z include_merged " + ) + + # Undo prebend + if definition.prebend[i_station] != 0: + cubit.cmd(f"move surface {l2s(volume_ids)} y {-1*definition.prebend[i_station]} include_merged") + + # Undo sweep + if definition.sweep[i_station] != 0: + raise ValueError("Presweep is untested for cross-sectional meshing") + + if 'ref_line_type' in cs_params and 'centroid' in cs_params['ref_line_type'].lower(): + centroidal_vert_id, centroidal_ref_line_coords=get_locus_of_cross_sectional_centroids([i_station]) + cubit.cmd(f"move surface {l2s(volume_ids)} x {-1*centroidal_ref_line_coords[0][0]} include_merged") + cubit.cmd(f"move surface {l2s(volume_ids)} y {-1*centroidal_ref_line_coords[0][1]} include_merged") + #cubit.cmd(f"move surface {l2s(volume_ids)} z {-1*centroidal_ref_line_coords[0][2]} include_merged") + + centroidal_vert_id, centroidal_ref_line_coords=get_locus_of_cross_sectional_centroids([i_station]) - # Mesh the cross-section - cubit.cmd(f'curve with name "face_thickness*" interval {cs_params["nel_per_layer"]}') - cubit.cmd(f'curve with name "*face_web_thickness*" interval {cs_params["nel_per_layer"]}') - cubit.cmd(f'curve with name "core_thickness*" interval {cs_params["nel_per_core_layer"]}') - cubit.cmd(f'curve with name "*core_web_thickness*" interval {cs_params["nel_per_core_layer"]}') + # Mesh the cross-section + cubit.cmd(f'curve with name "face_thickness*" interval {cs_params["nel_per_layer"]}') + cubit.cmd(f'curve with name "*face_web_thickness*" interval {cs_params["nel_per_layer"]}') + cubit.cmd(f'curve with name "core_thickness*" interval {cs_params["nel_per_core_layer"]}') + cubit.cmd(f'curve with name "*core_web_thickness*" interval {cs_params["nel_per_core_layer"]}') - cubit.cmd(f'curve with name "*hoop*" in surface with name "roundTEadhesive*" interval {cs_params["nel_per_layer"]}') - # cubit.cmd(f'imprint volume {l2s(surface_ids)}') - cubit.cmd(f"merge volume {l2s(volume_ids)}") - cubit.cmd(f"set default autosize on") + cubit.cmd(f'curve with name "*hoop*" in surface with name "roundTEadhesive*" interval {cs_params["nel_per_layer"]}') - if cs_params["element_shape"].lower() == "tri": - cubit.cmd(f"surface {l2s(volume_ids)} scheme tri") - else: - cubit.cmd(f"surface {l2s(volume_ids)} scheme map") + # cubit.cmd(f'imprint volume {l2s(surface_ids)}') + cubit.cmd(f"merge volume {l2s(volume_ids)}") + cubit.cmd(f"set default autosize on") + if cs_params["element_shape"].lower() == "tri": + cubit.cmd(f"surface {l2s(volume_ids)} scheme tri") + else: + cubit.cmd(f"surface {l2s(volume_ids)} scheme map") - t_1=get_mean_layer_thickness_at_station(i_station) #Find mean layer thickness for first station - e_size_1=t_1/cs_params['nel_per_layer']*cs_params['element_ar'] #Get spanwise element size at station_id cross section - - cubit.cmd(f"surface all size {e_size_1}") - cubit.cmd(f"mesh surface {l2s(volume_ids)}") - - file_name = wt_name + "-" + str(i_station) + "-t-0.in" - - if not os.path.exists(directory): - os.makedirs(directory) - - if get_mesh_error_count() ==0: - if settings["make_input_for"] is not None: - if "vabs" in settings["make_input_for"].lower(): - write_vabs_input( - surface_dict, - blade, - cs_params, - directory, - file_name, - volume_ids, - materials_used, - cs_normal, - ) - elif "anba" in settings["make_input_for"].lower(): - raise ValueError("ANBA currently not supported") - else: - raise NameError( - f'Unknown beam cross-sectional solver: {settings["make_input_for"]}' - ) - else: - with open(f"{wt_name}.log", "a") as logFile: - logFile.write(f" Warning: {get_mesh_error_count()} cross section mesh errors exist in station {i_station}\n") - - if 'd_tube' in cs_params.keys() and cs_params['d_tube']: - - set_verts={} - set_verts[f'thickness_{str(i_station).zfill(3)}_s1']=[2804, 9817, 9823, 9831] - set_verts[f'thickness_{str(i_station).zfill(3)}_s2']=[9985, 10001, 10039, 10083] - set_verts[f'spanwise_s2']=[9985] - set_verts[f'circumferential_{str(i_station).zfill(3)}']=[9650,2801,2802,2804,2806,2808,2810,2812,2814,9979,9983,9985,10117,10114,10115,6396,6395,6398,6266,6264,6260,2764,2762,2760,2758,2756,2754,2752,2751,5931,5951,5989,6033,11530,11538,11714,11706,9752,9708,9670,9650] - set_verts[f'p1']=[9650] - set_verts[f'p2']=[5931] - - file_name=f'beam_{str(i_station).zfill(3)}.nodes' - write_path_node_ids_to_file(set_verts,file_name,directory) - #write_path_coords_to_file(set_verts,prepend,dir_name) - - file_name=f'{directory}/beam_{str(i_station).zfill(3)}.abscissa' - write_path_abscissas_to_file(set_verts,file_name) + t_1=get_mean_layer_thickness_at_station(i_station) #Find mean layer thickness for first station + e_size_1=t_1/cs_params['nel_per_layer']*cs_params['element_ar'] #Get spanwise element size at station_id cross section + cubit.cmd(f"surface all size {e_size_1}") + cubit.cmd(f"mesh surface {l2s(volume_ids)}") + + file_name = wt_name + "-" + str(i_station) + "-t-0.in" + + if not os.path.exists(directory): + os.makedirs(directory) + + if get_mesh_error_count() ==0: + if settings["make_input_for"] is not None: + if "vabs" in settings["make_input_for"].lower(): + write_vabs_input( + surface_dict, + blade, + cs_params, + directory, + file_name, + volume_ids, + materials_used, + cs_normal, + ) + + elif "anba" in settings["make_input_for"].lower(): + raise ValueError("ANBA currently not supported") + else: + raise NameError( + f'Unknown beam cross-sectional solver: {settings["make_input_for"]}' + ) + else: + with open(f"{wt_name}.log", "a") as logFile: + logFile.write(f" Warning: {get_mesh_error_count()} cross section mesh errors exist in station {i_station}\n") + + if 'd_tube' in cs_params.keys() and cs_params['d_tube']: + + set_verts={} + set_verts[f'thickness_{str(i_station).zfill(3)}_s1']=[2804, 9817, 9823, 9831] + set_verts[f'thickness_{str(i_station).zfill(3)}_s2']=[9985, 10001, 10039, 10083] + set_verts[f'spanwise_s2']=[9985] + set_verts[f'circumferential_{str(i_station).zfill(3)}']=[9650,2801,2802,2804,2806,2808,2810,2812,2814,9979,9983,9985,10117,10114,10115,6396,6395,6398,6266,6264,6260,2764,2762,2760,2758,2756,2754,2752,2751,5931,5951,5989,6033,11530,11538,11714,11706,9752,9708,9670,9650] + set_verts[f'p1']=[9650] + set_verts[f'p2']=[5931] + + file_name=f'beam_{str(i_station).zfill(3)}.nodes' + write_path_node_ids_to_file(set_verts,file_name,directory) + #write_path_coords_to_file(set_verts,prepend,dir_name) + + file_name=f'{directory}/beam_{str(i_station).zfill(3)}.abscissa' + write_path_abscissas_to_file(set_verts,file_name) + - # nodeset_id= cubit.get_next_nodeset_id() - # cubit.cmd(f'nodeset {nodeset_id} add curve 925 928 932') - # cubit.cmd(f'nodeset {nodeset_id} name "{node_set_name}"') + # nodeset_id= cubit.get_next_nodeset_id() + # cubit.cmd(f'nodeset {nodeset_id} add curve 925 928 932') + # cubit.cmd(f'nodeset {nodeset_id} name "{node_set_name}"') - # path_type='thickness' - #node_order=get_path_node_order(node_set_name,path_type) - #def get_path_node_order(node_set_name,path_type): + # path_type='thickness' + #node_order=get_path_node_order(node_set_name,path_type) + #def get_path_node_order(node_set_name,path_type): - # nodeset_nodes = get_nodeset_nodes_from_name(node_set_name) - # coords=get_nodal_coordinates_from_set_of_nodes(nodeset_nodes) + # nodeset_nodes = get_nodeset_nodes_from_name(node_set_name) + # coords=get_nodal_coordinates_from_set_of_nodes(nodeset_nodes) - # pointer=order_path_points(coords, ind) + # pointer=order_path_points(coords, ind) - # file = open(directory +'/'+ axisFileName, 'w') - # file.write('--------- BEAMDYN with OpenFAST INPUT FILE -------------------------------------------\n') - - # nodeset_id= cubit.get_next_nodeset_id() - # cubit.cmd(f'nodeset {nodeset_id} add curve 1040 1065 1091') - # cubit.cmd(f'nodeset {nodeset_id} name "s2_thickness_path"') + # file = open(directory +'/'+ axisFileName, 'w') + # file.write('--------- BEAMDYN with OpenFAST INPUT FILE -------------------------------------------\n') + + # nodeset_id= cubit.get_next_nodeset_id() + # cubit.cmd(f'nodeset {nodeset_id} add curve 1040 1065 1091') + # cubit.cmd(f'nodeset {nodeset_id} name "s2_thickness_path"') - # nodeset_id= cubit.get_next_nodeset_id() - # cubit.cmd(f'nodeset {nodeset_id} add vertex 9985') - # cubit.cmd(f'nodeset {nodeset_id} name "spanwise_path"') + # nodeset_id= cubit.get_next_nodeset_id() + # cubit.cmd(f'nodeset {nodeset_id} add vertex 9985') + # cubit.cmd(f'nodeset {nodeset_id} name "spanwise_path"') - # nodeset_id= cubit.get_next_nodeset_id() - # cubit.cmd(f'nodeset {nodeset_id} add curve 73 74 75 76 77 78 79 92 93 94 95 96 97 98 335 361 476 478 479 556 557 558 824 848 873 899 1014 1016 1017 1094 1095 1096 1196 1220 1224 1324 1328 1422') - # cubit.cmd(f'nodeset {nodeset_id} name "circumferential_path"') - + # nodeset_id= cubit.get_next_nodeset_id() + # cubit.cmd(f'nodeset {nodeset_id} add curve 73 74 75 76 77 78 79 92 93 94 95 96 97 98 335 361 476 478 479 556 557 558 824 848 873 899 1014 1016 1017 1094 1095 1096 1196 1220 1224 1324 1328 1422') + # cubit.cmd(f'nodeset {nodeset_id} name "circumferential_path"') + - # nodeset_id= cubit.get_next_nodeset_id() - # cubit.cmd(f'nodeset {nodeset_id} add curve with name "*oml*"') - # cubit.cmd(f'nodeset {nodeset_id} name "{node_set_name}"') - # oml_nodes = get_nodeset_nodes_from_name(node_set_name) - - - if settings["export"] is not None: - if ( - "g" in settings["export"].lower() - or "cub" in settings["export"].lower() - ): - if "g" in settings["export"].lower(): - cubit.cmd(f'export mesh "{path_name}-{str(i_station)}.g" overwrite') - if "cub" in settings["export"].lower(): - cubit.cmd(f"delete curve {spanwise_mat_ori_curve}") - cubit.cmd(f'save as "{path_name}-{str(i_station)}.cub" overwrite') - print('') - elif len(settings["export"]) == 0: - pass - else: - raise NameError( - f'Unknown model export format: {settings["export"]}' - ) + # nodeset_id= cubit.get_next_nodeset_id() + # cubit.cmd(f'nodeset {nodeset_id} add curve with name "*oml*"') + # cubit.cmd(f'nodeset {nodeset_id} name "{node_set_name}"') + # oml_nodes = get_nodeset_nodes_from_name(node_set_name) + + + if settings["export"] is not None: + if ( + "g" in settings["export"].lower() + or "cub" in settings["export"].lower() + ): + if "g" in settings["export"].lower(): + cubit.cmd(f'export mesh "{path_name}-{str(i_station)}.g" overwrite') + if "cub" in settings["export"].lower(): + cubit.cmd(f"delete curve {spanwise_mat_ori_curve}") + cubit.cmd(f'save as "{path_name}-{str(i_station)}.cub" overwrite') + print('') + elif len(settings["export"]) == 0: + pass + else: + raise NameError( + f'Unknown model export format: {settings["export"]}' + ) diff --git a/src/pynumad/analysis/cubit/make_cross_sections.py b/src/pynumad/analysis/cubit/make_cross_sections.py index a86cea9..70cd527 100644 --- a/src/pynumad/analysis/cubit/make_cross_sections.py +++ b/src/pynumad/analysis/cubit/make_cross_sections.py @@ -1640,6 +1640,362 @@ def make_cs_web_layer_areas( cubit.cmd(f'surface {surf_id} rename "{surf_name}"') return part_name_id, (vHP, vLP) +def make_a_precomp_cross_section(wt_name, + surface_dict, + i_station, + i_station_geometry, + blade, + hasWebs, + aft_web_stack, + fore_web_stack, + iLE, + cs_params, + geometry_scaling, + thickness_scaling, + last_round_station, + last_flat_station, + materials_used, + cs_normal, +): + geometry = blade.geometry + stackdb = blade.stackdb + keypoints = blade.keypoints + definition = blade.definition + + if i_station > last_round_station: + is_flatback=True + else: + is_flatback=False + + with open(f"{wt_name}.log", "a") as logFile: + logFile.write(f"Working on Station: {i_station}\n") + + part_name_id = 0 + + + #### Step one create outer mold line + xyz = get_blade_geometry_for_station(blade, i_station_geometry) * geometry_scaling + xyz[:,0] = xyz[:,0]/blade.geometry.ichord[i_station_geometry] + xyz[:,1] = xyz[:,1]/blade.geometry.ichord[i_station_geometry] + + # Start indexing from 1 (not 0) to ignore first point: because first point is not on the LP or HP surface but rather is the midpoint at the TE + #xyz = np.flip(xyz, 0) + xyz[:, 0]=-1*xyz[:, 0] + spline_points = xyz[iLE:len(xyz)-1, :] + + write_spline_from_coordinate_points(cubit, spline_points) + lp_key_curve = get_last_id("curve") + LE_vert, v2 = selCurveVerts(lp_key_curve) + + xyz = np.flip(xyz, 0) + spline_points = np.flip(xyz[iLE-1:len(xyz)-1, :],0) + write_spline_from_coordinate_points(cubit, spline_points) + hp_key_curve = get_last_id("curve") + v3, v4 = selCurveVerts(hp_key_curve) + + + # Undo initial twist + cubit.cmd(f"rotate curve {lp_key_curve} {hp_key_curve} angle {-1*geometry.idegreestwist[i_station_geometry]} about Z include_merged ") + x_move,y_move,_=cubit.vertex(LE_vert).coordinates() + cubit.cmd(f"move curve {lp_key_curve} {hp_key_curve} x {-x_move} y {-y_move} include_merged") + + #Check of HP LE is less than 0 + + x,_,_=cubit.vertex(v4).coordinates() + if x <= 0.0: + cubit.cmd(f"move curve {hp_key_curve} x {-x+0.0001} y {-y_move} include_merged") + + + cubit.cmd(f"delete vertex with is_free") + n_start = get_last_id("vertex")+1 + cubit.cmd(f"create vertex on curve {lp_key_curve} segment 50") + n_end = get_last_id("vertex") + + lp_vert_ids = [LE_vert]+list(range(n_start, n_end + 1))+[v2] + + n_start = get_last_id("vertex")+1 + cubit.cmd(f"create vertex on curve {hp_key_curve} segment 50") + n_end = get_last_id("vertex") + + hp_vert_ids = [v3]+list(range(n_start, n_end + 1))+[v4] + + vert_ids=lp_vert_ids+hp_vert_ids + cubit.cmd(f'save as "Debug.cub" overwrite') + # n_start = get_last_id("vertex")+1 + + # LE_vert = n_start + + # for coords in spline_points: + # create_vertex(coords[0], coords[1], coords[2]) + # flatback_vTop = get_last_id("vertex") + # flatback_vBot = flatback_vTop+1 + + + + # for coords in spline_points: + # create_vertex(coords[0], coords[1], coords[2]) + + # n_end = get_last_id("vertex") + # vert_ids = list(range(n_start, n_end + 1)) + + # # Undo initial twist + # cubit.cmd(f"rotate vertex {l2s(vert_ids)} angle {-1*geometry.idegreestwist[i_station_geometry]} about Z include_merged ") + + # # Move LE to origin + # x_move,y_move,_=cubit.vertex(LE_vert).coordinates() + # cubit.cmd(f"move vertex {l2s(vert_ids)} x {-x_move} y {-y_move} include_merged") + + # #### Step one create outer mold line + # xyz = get_blade_geometry_for_station(blade, i_station_geometry) * geometry_scaling + # xyz[:,0] = xyz[:,0]/blade.geometry.ichord[i_station_geometry] + # xyz[:,1] = xyz[:,1]/blade.geometry.ichord[i_station_geometry] + + # # Start indexing from 1 (not 0) to ignore first point: because first point is not on the LP or HP surface but rather is the midpoint at the TE + # #xyz = np.flip(xyz, 0) + # xyz[:, 0]=-1*xyz[:, 0] + # spline_points = xyz[iLE:len(xyz)-1, :] + + # n_start = get_last_id("vertex")+1 + + # LE_vert = n_start + + # for coords in spline_points: + # create_vertex(coords[0], coords[1], coords[2]) + # flatback_vTop = get_last_id("vertex") + # flatback_vBot = flatback_vTop+1 + + # xyz = np.flip(xyz, 0) + # spline_points = np.flip(xyz[iLE-1:len(xyz)-1, :],0) + + # for coords in spline_points: + # create_vertex(coords[0], coords[1], coords[2]) + + # n_end = get_last_id("vertex") + # vert_ids = list(range(n_start, n_end + 1)) + + # # Undo initial twist + # cubit.cmd(f"rotate vertex {l2s(vert_ids)} angle {-1*geometry.idegreestwist[i_station_geometry]} about Z include_merged ") + + # # Move LE to origin + # x_move,y_move,_=cubit.vertex(LE_vert).coordinates() + # cubit.cmd(f"move vertex {l2s(vert_ids)} x {-x_move} y {-y_move} include_merged") + + + ### Write Shape File + directory='.' + + if not os.path.exists(directory): + os.makedirs(directory) + + file = open(f'{directory}/shape_{i_station}.inp', 'w') + file.write(f'{len(vert_ids)} N_af_nodes :no of airfoil nodes, counted clockwise starting\n') + file.write(f' with leading edge (see users" manual, fig xx)\n\n') + file.write(f' Xnode Ynode !! chord-normalized coordinated of the airfoil nodes\n') + + for vertex in vert_ids: + x,y,_=cubit.vertex(vertex).coordinates() + file.write(f'{min(x,1.00)} {y}\n') + file.close() + + + + flatbackCurve = cubit.create_curve(cubit.vertex(v2), cubit.vertex(v3)) + flatback_curve_id = flatbackCurve.id() + + ## Create chord line + cubit.cmd(f'create vertex on curve {flatback_curve_id} fraction 0.5 from start') + chord_line = cubit.create_curve(cubit.vertex(LE_vert), cubit.vertex(get_last_id('vertex'))) + chord_line_id = chord_line.id() + + #HP Side + key_points = keypoints.key_points[1:5, :, i_station_geometry] + key_points[:, 0]=-1*key_points[:, 0] + key_points[:,:2] = key_points[:,:2]/blade.geometry.ichord[i_station_geometry] + key_points=np.flip(key_points,0) + + n_start = get_last_id("vertex")+1 + for coords in key_points: + create_vertex(coords[0], coords[1], coords[2]) + n_end = get_last_id("vertex") + vert_ids = list(range(n_start, n_end + 1)) + + # Undo initial twist + cubit.cmd(f"rotate vertex {l2s(vert_ids)} angle {-1*geometry.idegreestwist[i_station_geometry]} about Z include_merged " ) + cubit.cmd(f"move vertex {l2s(vert_ids)} x {-x_move} y {-y_move} include_merged") + + + hp_sector_boundaries = [0.0] + for vertex in vert_ids: + x,_,_=cubit.vertex(vertex).coordinates() + hp_sector_boundaries.append(x) + hp_sector_boundaries.append(1.0) + + ## LP Side + temp = np.flip(keypoints.key_points[:, :, i_station_geometry], 0) + + key_points =temp[1:5, :] + key_points[:, 0]=-1*key_points[:, 0] + key_points[:,:2] = key_points[:,:2]/geometry.ichord[i_station_geometry] + key_points=np.flip(key_points,0) + + n_start = get_last_id("vertex")+1 + for coords in key_points: + create_vertex(coords[0], coords[1], coords[2]) + n_end = get_last_id("vertex") + vert_ids = list(range(n_start, n_end + 1)) + + # Undo initial twist + cubit.cmd(f"rotate vertex {l2s(vert_ids)} angle {-1*geometry.idegreestwist[i_station_geometry]} about Z include_merged " ) + cubit.cmd(f"move vertex {l2s(vert_ids)} x {-x_move} y {-y_move} include_merged") + + + lp_sector_boundaries = [0.0] + for vertex in vert_ids: + x,_,_=cubit.vertex(vertex).coordinates() + lp_sector_boundaries.append(x) + lp_sector_boundaries.append(1.0) + + + + + file = open(f'{directory}/layup_{i_station}.inp', 'w') + file.write(f'Composite laminae lay-up inside the blade section\n\n') + + N_scts=len(lp_sector_boundaries)-1 + file.write(f'*************************** TOP SURFACE ****************************\n') + file.write(f'{N_scts} N_scts(1): no of sectors on top surface\n\n') + file.write(f'normalized chord location of nodes defining airfoil sectors boundaries (xsec_node)\n') + file.write(' '.join(str(e) for e in lp_sector_boundaries)+'\n') + #file.write('') + + blade_material_names = list(blade.definition.materials.keys()) + + temp = stackdb.stacks[:, i_station] + temp = np.flip(temp) + + #TOP SURFACE IS temp = stackdb.stacks[:, i_station] temp = np.flip(temp) + #BOTTOM SURFACE IS: stackdb.stacks[1:6, i_station] + for i_stack,stack in enumerate(np.flip(temp[1:6])): + file.write(f'...........................................................................\n') + file.write(f'Sect_num no of laminae (N_laminas)\n') + file.write(f'{i_stack+1} {len(stack.plygroups)}\n\n') + + file.write(f'lamina num of thickness fibers_direction composite_material ID\n') + file.write(f'number plies of ply (m) (deg) (-)\n') + file.write(f'lam_num N_plies Tply Tht_lam Mat_id\n') + + + for i_ply,ply in enumerate(stack.plygroups): + material_name =ply.materialid + mat_id = blade_material_names.index(material_name)+1 + file.write(f'{i_ply+1} {1} {ply.thickness*ply.nPlies/1000.0} {ply.angle} {mat_id} ({ply.materialid})\n') + + N_scts=len(lp_sector_boundaries)-1 + file.write(f'\n\n*************************** BOTTOM SURFACE ****************************\n') + file.write(f'{N_scts} N_scts(1): no of sectors on top surface\n\n') + file.write(f'normalized chord location of nodes defining airfoil sectors boundaries (xsec_node)\n') + file.write(' '.join(str(e) for e in hp_sector_boundaries)+'\n') + + # stackdb.stacks[1:6, i_station] + # temp = stackdb.stacks[:, i_station] + # temp = np.flip(temp) + for i_stack,stack in enumerate(np.flip(stackdb.stacks[1:6, i_station])): + file.write(f'...........................................................................\n') + file.write(f'Sect_num no of laminae (N_laminas)\n') + file.write(f'{i_stack+1} {len(stack.plygroups)}\n\n') + + file.write(f'lamina num of thickness fibers_direction composite_material ID\n') + file.write(f'number plies of ply (m) (deg) (-)\n') + file.write(f'lam_num N_plies Tply Tht_lam Mat_id\n') + + for i_ply,ply in enumerate(stack.plygroups): + material_name =ply.materialid + mat_id = blade_material_names.index(material_name)+1 + file.write(f'{i_ply+1} {1} {ply.thickness*ply.nPlies/1000} {ply.angle} {mat_id} ({ply.materialid})\n') + + + if hasWebs: + file.write(f'\n\n**********************************************************************\n') + file.write(f'Laminae schedule for webs (input required only if webs exist at this section):\n') + n_webs=len(stackdb.swstacks) + #Foreweb + for web_number in range(n_webs): + web_stack = stackdb.swstacks[web_number][i_station] + file.write(f'\nweb_num no of laminae (N_weblams)\n') + file.write(f'{web_number+1} {len(web_stack.plygroups)}\n\n') + + file.write(f'lamina num of thickness fibers_direction composite_material ID\n') + file.write(f'number plies of ply (m) (deg) (-)\n') + file.write(f'wlam_num N_Plies w_tply Tht_Wlam Wmat_Id\n') + + for i_ply,ply in enumerate(web_stack.plygroups): + material_name =ply.materialid + mat_id = blade_material_names.index(material_name)+1 + file.write(f'{i_ply+1} {1} {ply.thickness*ply.nPlies/1000} {ply.angle} {mat_id} ({ply.materialid})\n') + + + file.close() + + #### Materials ##### + from decimal import Decimal + file = open(f'{directory}/materials.inp', 'w') + + file.write(f'\nMat_Id E1 E2 G12 Nu12 Density Mat_Name\n') + file.write(f' (-) (Pa) (Pa) (Pa) (-) (Kg/m^3) (-)\n') + + for i_mat,material_name in enumerate(blade_material_names): + E1=blade.definition.materials[material_name].ex + E2=blade.definition.materials[material_name].ey + G12=blade.definition.materials[material_name].gxy + nu12=blade.definition.materials[material_name].prxy + rho=blade.definition.materials[material_name].density + file.write(f'{i_mat+1} {Decimal(str(E1)):10.6E} {Decimal(str(E2)):10.6E} {Decimal(str(G12)):10.6E} {Decimal(str(nu12)):10.6E} {Decimal(str(rho)):10.6E} ({material_name})\n') + + file.close() + + file = open(f'{directory}/input_{i_station}.pci', 'w') + + file.write(f'***************** main input file for PreComp *****************************\n') + file.write(f'Sample Composite Blade Section Properties\n') + file.write(f'\n') + file.write(f'General information -----------------------------------------------\n') + file.write(f'1 Bl_length : blade length (m)\n') + file.write(f'2 N_sections : no of blade sections (-)\n') + file.write(f'{len(blade_material_names)} N_materials : no of materials listed in the materials table (material.inp)\n') + file.write(f'3 Out_format : output file (1: general format, 2: BModes-format, 3: both)\n') + file.write(f'f TabDelim (true: tab-delimited table; false: space-delimited table)\n') + file.write(f'\n') + file.write(f'Blade-sections-specific data --------------------------------------\n') + file.write(f'Sec span l.e. chord aerodynamic af_shape int str layup\n') + file.write(f'location position length twist file file\n') + file.write(f'Span_loc Le_loc Chord Tw_aero Af_shape_file Int_str_file\n') + file.write(f' (-) (-) (m) (degrees) (-) (-)\n\n') + + file.write(f'0.00 {-1*x_move} {geometry.ichord[i_station_geometry]} {geometry.idegreestwist[i_station_geometry]} shape_{i_station}.inp layup_{i_station}.inp\n') + file.write(f'1.00 {-1*x_move} {geometry.ichord[i_station_geometry]} {geometry.idegreestwist[i_station_geometry]} shape_{i_station}.inp layup_{i_station}.inp\n') + + if hasWebs: + n_webs =2 + else: + n_webs = 0 + file.write(f'\nWebs (spars) data --------------------------------------------------\n\n') + file.write(f'{n_webs} Nweb : number of webs (-) ! enter 0 if the blade has no webs\n') + file.write(f'1 Ib_sp_stn : blade station number where inner-most end of webs is located (-)\n') + file.write(f'2 Ob_sp_stn : blade station number where outer-most end of webs is located (-)\n\n') + + + file.write(f'Web_num Inb_end_ch_loc Oub_end_ch_loc (fraction of chord length)\n') + if hasWebs: + i_loc = 2 + web_loc= mean([lp_sector_boundaries[i_loc],hp_sector_boundaries[i_loc]]) + file.write(f'1.0000 {web_loc} {web_loc}\n') + + i_loc = 3 + web_loc= mean([lp_sector_boundaries[i_loc],hp_sector_boundaries[i_loc]]) + file.write(f'2.0000 {web_loc} {web_loc}\n') + file.write(f'\n') + + file.close() def make_a_cross_section(wt_name, surface_dict, From 69408f78698687e1a689409ddb32bbafe6fee67d Mon Sep 17 00:00:00 2001 From: Camarena Date: Wed, 23 Oct 2024 09:15:59 -0600 Subject: [PATCH 05/11] computing last_flat_station bug fix --- src/pynumad/analysis/cubit/make_blade.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index fac32d6..b90a426 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -685,9 +685,10 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati athickness=cs_params["te_adhesive_thickness"][last_10deg_station+1+i_length] print(f'i {last_10deg_station+1+i_length},excess_length {excess_length*1000}, athickness{athickness*1000}') if (excess_length-athickness)/excess_length > 0.025: - last_flat_station=last_flat_station+1+i_length + last_flat_station+=1 else: break + last_flat_station+=1 ######################## From d51a8afe2d838c49d3ffed8ae7b29b6cafd2e43e Mon Sep 17 00:00:00 2001 From: Camarena Date: Fri, 15 Nov 2024 15:40:55 -0700 Subject: [PATCH 06/11] first pass at allwoing for more than one web --- .../analysis/cubit/make_cross_sections.py | 4 +- src/pynumad/io/yaml_to_blade.py | 349 ++++++++++++------ src/pynumad/objects/definition.py | 3 + src/pynumad/objects/geometry.py | 4 +- src/pynumad/objects/keypoints.py | 344 +++++++++-------- src/pynumad/objects/stackdb.py | 10 + src/pynumad/utils/misc_utils.py | 16 + 7 files changed, 474 insertions(+), 256 deletions(-) diff --git a/src/pynumad/analysis/cubit/make_cross_sections.py b/src/pynumad/analysis/cubit/make_cross_sections.py index 70cd527..0fa0790 100644 --- a/src/pynumad/analysis/cubit/make_cross_sections.py +++ b/src/pynumad/analysis/cubit/make_cross_sections.py @@ -879,9 +879,7 @@ def make_cs_perimeter_layer_areas(wt_name, bottom_left_vertex_curve_left, bottom_right_vertex_curve_left = selCurveVerts( left_bottom_curve ) - bottom_left_vertex_curve_right, bottom_right_vertex_curve_right = selCurveVerts( - right_bottom_curve - ) + bottom_left_vertex_curve_right, bottom_right_vertex_curve_right = selCurveVerts(right_bottom_curve) # This if statement prepares all layer curves such that they taper at the TE if i_perimeter == 0: diff --git a/src/pynumad/io/yaml_to_blade.py b/src/pynumad/io/yaml_to_blade.py index d017fcd..7848b55 100644 --- a/src/pynumad/io/yaml_to_blade.py +++ b/src/pynumad/io/yaml_to_blade.py @@ -2,12 +2,13 @@ import numpy as np import logging from scipy.stats import mode +import string from pynumad.utils.misc_utils import ( LARCetaT, LARCetaL, _parse_data, - full_keys_from_substrings, + full_keys_from_substrings, sort_based_on_root_values ) from pynumad.io.airfoil_to_xml import airfoil_to_xml from pynumad.utils.interpolation import interpolator_wrap @@ -59,6 +60,9 @@ def yaml_to_blade(blade, filename: str, write_airfoils: bool = False): # obtain blade internal structure blade_internal_structure = data["components"]["blade"]["internal_structure_2d_fem"] + #Get number of webs + blade.definition.n_webs = len(blade_internal_structure['webs']) + # obtain airfoil data af_data = data["airfoils"] @@ -82,9 +86,7 @@ def yaml_to_blade(blade, filename: str, write_airfoils: bool = False): ## Blade Components # Update "grid" and "values" keys to cover the whole span of the blade - blade_internal_structure = update_internal_structure( - blade_internal_structure, blade_outer_shape_bem - ) + blade_internal_structure = update_internal_structure(blade_internal_structure, blade_outer_shape_bem) blade_structure_dict = { blade_internal_structure["layers"][i]["name"].lower(): blade_internal_structure[ @@ -92,24 +94,79 @@ def yaml_to_blade(blade, filename: str, write_airfoils: bool = False): ][i] for i in range(len(blade_internal_structure["layers"])) } + + #Initialize + definition.key_arcs = [] + definition.key_arcs = [] + + # TE Bands + _add_te_le_bands(definition, blade_structure_dict,'te') + # Spar caps _add_spar_caps(definition, blade_structure_dict) - # TE Bands - _add_te_bands(definition, blade_structure_dict) + # Webs Bands + i_sorted = _add_webs(definition,blade_structure_dict,blade_internal_structure['webs']) #also returns i_sorted. The web sorting from te to le. # LE Bands - _add_le_bands(definition, blade_structure_dict) + _add_te_le_bands(definition, blade_structure_dict,'le') + + definition.key_arcs = np.transpose(np.array(definition.key_arcs)) + definition.key_arcs,_= sort_based_on_root_values(definition.key_arcs) + +########### TEMP + leReinfKeys = full_keys_from_substrings(blade_structure_dict.keys(), ["le", "reinf"]) + if len(leReinfKeys) == 1: + definition.leband = ( + blade_structure_dict[leReinfKeys[0]]["width"]["values"] * 1000 / 2 + ) + elif len(leReinfKeys) == 2: + definition.leband = ( + ( + blade_structure_dict[leReinfKeys[0]]["width"]["values"] + + blade_structure_dict[leReinfKeys[1]]["width"]["values"] + ) + * 1000 + / 2 + ) +########### TEMP + + # number of areas around airfoil profile; must be even (see calc of web areas) + _,num_defined_arcs = np.shape(definition.key_arcs) + + num_areas = num_defined_arcs+4 + + blade.keypoints.key_labels.append('te') + letter_list = list(string.ascii_lowercase[0:int(num_areas/2)-1]) + + for letter in np.flip(letter_list): + blade.keypoints.key_labels.append(letter) + blade.keypoints.key_labels.append('le') + for letter in letter_list: + blade.keypoints.key_labels.append(letter) + blade.keypoints.key_labels.append('te') + ### COMPONENTS - _add_components(definition, blade_internal_structure, blade_structure_dict) + _add_components(definition, blade_internal_structure,letter_list,blade.keypoints.key_labels,i_sorted,blade_internal_structure['webs']) + + + #Sort based on root values + + + + + + + + + # definition.key_arcs = np.transpose(np.array(definition.key_arcs)) blade.update_blade() - # save(blade_name) - # BladeDef_to_NuMADfile(obj,numad_name,matdb_name,numad_af_folder) - return blade + return blade + def _add_stations( definition, blade_outer_shape_bem, @@ -344,12 +401,14 @@ def _add_materials(definition, material_data): return -def _add_components(definition, blade_internal_structure, blade_structure_dict): +def _add_components(definition, blade_internal_structure, letter_list,key_labels,i_sorted,web_data_structure): N_layer_comp = len(blade_internal_structure["layers"]) component_list = list() + web_component_indecies=list() for i in range(N_layer_comp): i_component_data = blade_internal_structure["layers"][i] cur_comp = Component() + cur_comp._keylabels = key_labels cur_comp.group = 0 cur_comp.name = i_component_data["name"] # comp['material'] = blade_internal_structure['layers']{i}['material']; @@ -384,6 +443,9 @@ def _add_components(definition, blade_internal_structure, blade_structure_dict): cur_comp.pinnedends = 0 component_list.append(cur_comp) + if 'web' in i_component_data["name"].lower(): + web_component_indecies.append(i) + component_dict = dict() for comp in component_list: component_dict[comp.name] = comp @@ -421,18 +483,18 @@ def _add_components(definition, blade_internal_structure, blade_structure_dict): # TE Band(s) key_list = full_keys_from_substrings(component_dict.keys(), ["te", "reinf"]) if len(key_list) == 1: - component_dict[key_list[0]].hpextents = ["d", "te"] - component_dict[key_list[0]].lpextents = ["d", "te"] + component_dict[key_list[0]].hpextents = [letter_list[-2], "te"] + component_dict[key_list[0]].lpextents = [letter_list[-2], "te"] elif len(key_list) == 2: tempKeyList = full_keys_from_substrings(key_list, ["ss"]) if len(tempKeyList) == 1: - component_dict[tempKeyList[0]].lpextents = ["d", "te"] + component_dict[tempKeyList[0]].lpextents = [letter_list[-2], "te"] else: ValueError("Incorrect number of te reinf ss components") tempKeyList = full_keys_from_substrings(key_list, ["ps"]) if len(tempKeyList) == 1: - component_dict[tempKeyList[0]].hpextents = ["d", "te"] + component_dict[tempKeyList[0]].hpextents = [letter_list[-2], "te"] else: ValueError("Incorrect number of te reinf ps components") else: @@ -461,7 +523,7 @@ def _add_components(definition, blade_internal_structure, blade_structure_dict): # Trailing edge suction-side panel key_list = full_keys_from_substrings(component_dict.keys(), ["te_", "ss", "filler"]) if len(key_list) == 1: - component_dict[key_list[0]].lpextents = ["c", "d"] + component_dict[key_list[0]].lpextents = ["c", letter_list[-2]] else: raise ValueError("Invalid number of trailing edge suction-side panels") @@ -482,41 +544,56 @@ def _add_components(definition, blade_internal_structure, blade_structure_dict): # Leading edge suction-side panel key_list = full_keys_from_substrings(component_dict.keys(), ["te_", "ps", "filler"]) if len(key_list) == 1: - component_dict[key_list[0]].hpextents = ["c", "d"] + component_dict[key_list[0]].hpextents = ["c", letter_list[-2]] else: raise ValueError("Invalid number of trailing edge pressure-side panels") - # Web - - for comp in component_dict: - logging.debug(comp) - key_list = full_keys_from_substrings(component_dict.keys(), ["web", "fore"]) # Try 1 - if len(key_list) == 0: - key_list = full_keys_from_substrings(component_dict.keys(), ["web", "1"]) # Try 2 - - if len(key_list) > 0: - for key in key_list: - component_dict[key].hpextents = ["b"] - component_dict[key].lpextents = ["b"] - component_dict[key].group = 1 - elif len(key_list) == 0: - raise ValueError("No fore web layers found found") - - key_list = full_keys_from_substrings(component_dict.keys(), ["web", "aft"]) # Try 1 - if len(key_list) == 0: - key_list = full_keys_from_substrings(component_dict.keys(), ["web", "0"]) # Try 2 - if len(key_list) == 0: - key_list = full_keys_from_substrings( - component_dict.keys(), ["web", "rear"] - ) # Try 3 - - if len(key_list) > 0: - for key in key_list: - component_dict[key].hpextents = ["c"] - component_dict[key].lpextents = ["c"] - component_dict[key].group = 2 - elif len(key_list) == 0: - raise ValueError("No rear web layers found found") + # Webs + for i_web_index, i_web in enumerate(np.flip(i_sorted)): #Work from LE to TE + web_data= web_data_structure[i_web] + target_web_name = web_data['name'] + + target_web_components =[] + for web_component_index in web_component_indecies: + this_web_name = blade_internal_structure['layers'][web_component_index]['web'] + if this_web_name == target_web_name: + target_web_components.append(web_component_index) + + for web_component_index in target_web_components: + component_name = blade_internal_structure['layers'][web_component_index]['name'] + component_dict[component_name].hpextents = letter_list[i_web_index+1] + component_dict[component_name].lpextents = letter_list[i_web_index+1] + component_dict[component_name].group = i_web_index+1 + + # for comp in component_dict: + # logging.debug(comp) + # key_list = full_keys_from_substrings(component_dict.keys(), ["web", "fore"]) # Try 1 + # if len(key_list) == 0: + # key_list = full_keys_from_substrings(component_dict.keys(), ["web", "1"]) # Try 2 + + # if len(key_list) > 0: + # for key in key_list: + # component_dict[key].hpextents = ["b"] + # component_dict[key].lpextents = ["b"] + # component_dict[key].group = 1 + # elif len(key_list) == 0: + # raise ValueError("No fore web layers found found") + + # key_list = full_keys_from_substrings(component_dict.keys(), ["web", "aft"]) # Try 1 + # if len(key_list) == 0: + # key_list = full_keys_from_substrings(component_dict.keys(), ["web", "0"]) # Try 2 + # if len(key_list) == 0: + # key_list = full_keys_from_substrings( + # component_dict.keys(), ["web", "rear"] + # ) # Try 3 + + # if len(key_list) > 0: + # for key in key_list: + # component_dict[key].hpextents = ["c"] + # component_dict[key].lpextents = ["c"] + # component_dict[key].group = 2 + # elif len(key_list) == 0: + # raise ValueError("No rear web layers found found") ### add components to blade definition.components = component_dict @@ -572,87 +649,145 @@ def update_internal_structure(blade_internal_structure, blade_outer_shape_bem): ] = fullSpanValues return blade_internal_structure +def _add_webs(definition,blade_structure_dict,web_data_structure): + key_names = full_keys_from_substrings(blade_structure_dict.keys(), ["spar"]) + if len(key_names) != 2: + raise ValueError("Incorrect number of spar cap components") + + for i_side in range(2): + if "suc" in blade_structure_dict[key_names[i_side]]["side"].lower(): + lp_side_index = i_side + if "pres" in blade_structure_dict[key_names[i_side]]["side"].lower(): + hp_side_index = i_side + + + spar_cap_edges =np.transpose(np.array([blade_structure_dict[key_names[lp_side_index]]["start_nd_arc"]["values"], + blade_structure_dict[key_names[lp_side_index]]["end_nd_arc"]["values"]])) + + web_edges = [] + for i_web in range(definition.n_webs): + current_web = web_data_structure[i_web] + web_edges.append(current_web["start_nd_arc"]["values"][:]) + web_edges = np.transpose(np.array(web_edges)) + web_edges,i_sorted = sort_based_on_root_values(web_edges) + + if len(i_sorted)>2: #Only do the following if greater than two webs + mean_diff = [] + for i_web in range(definition.n_webs-1): + diff = abs(spar_cap_edges-web_edges[:,i_web:i_web+2]) + mean_diff.append(np.mean(diff)) + + matching_index = np.argmin(mean_diff) + keep_indecies = list(range(3)) + keep_indecies.pop(matching_index+1) + keep_indecies.pop(matching_index) #answer is zero + keep_webs = np.array(i_sorted)[keep_indecies] + + for i_web in keep_webs: + current_web = web_data_structure[i_web] + definition.key_arcs.append(current_web["start_nd_arc"]["values"]) + definition.key_arcs.append(current_web["end_nd_arc"]["values"]) + return i_sorted + def _add_spar_caps(definition, blade_structure_dict): - sparCapKeys = full_keys_from_substrings(blade_structure_dict.keys(), ["spar"]) - if len(sparCapKeys) != 2: + key_names = full_keys_from_substrings(blade_structure_dict.keys(), ["spar"]) + if len(key_names) != 2: raise ValueError("Incorrect number of spar cap components") - for iSparCap in range(2): - if "suc" in blade_structure_dict[sparCapKeys[iSparCap]]["side"].lower(): - lpSideIndex = iSparCap - if "pres" in blade_structure_dict[sparCapKeys[iSparCap]]["side"].lower(): - hpSideIndex = iSparCap + for i_side in range(2): + if "suc" in blade_structure_dict[key_names[i_side]]["side"].lower(): + lp_side_index = i_side + if "pres" in blade_structure_dict[key_names[i_side]]["side"].lower(): + hp_side_index = i_side - definition.sparcapwidth_lp = ( - blade_structure_dict[sparCapKeys[lpSideIndex]]["width"]["values"] * 1000 - ) + definition.sparcapwidth_lp = (blade_structure_dict[key_names[lp_side_index]]["width"]["values"] * 1000) + + # definition.sparcap_start_nd_arc = blade_structure_dict[key_names[lp_side_index]]["start_nd_arc"]["values"] + # definition.sparcap_end_nd_arc = blade_structure_dict[key_names[lp_side_index]]["end_nd_arc"]["values"] try: definition.sparcapoffset_lp = ( - blade_structure_dict[sparCapKeys[lpSideIndex]]["offset_y_pa"]["values"] + blade_structure_dict[key_names[lp_side_index]]["offset_y_pa"]["values"] * 1000 ) except KeyError: - definition.sparcap_start_nd_arc = blade_structure_dict[ - sparCapKeys[lpSideIndex] - ]["start_nd_arc"]["values"] - definition.sparcap_end_nd_arc = blade_structure_dict[sparCapKeys[lpSideIndex]][ - "end_nd_arc" - ]["values"] + pass definition.sparcapwidth_hp = ( - blade_structure_dict[sparCapKeys[hpSideIndex]]["width"]["values"] * 1000 + blade_structure_dict[key_names[hp_side_index]]["width"]["values"] * 1000 ) + + definition.key_arcs.append(blade_structure_dict[key_names[lp_side_index]]["start_nd_arc"]["values"]) + definition.key_arcs.append(blade_structure_dict[key_names[lp_side_index]]["end_nd_arc"]["values"]) + + definition.key_arcs.append(blade_structure_dict[key_names[hp_side_index]]["start_nd_arc"]["values"]) + definition.key_arcs.append(blade_structure_dict[key_names[hp_side_index]]["end_nd_arc"]["values"]) + try: definition.sparcapoffset_hp = ( - blade_structure_dict[sparCapKeys[hpSideIndex]]["offset_y_pa"]["values"] + blade_structure_dict[key_names[hp_side_index]]["offset_y_pa"]["values"] * 1000 ) except KeyError: - definition.sparcap_start_nd_arc = blade_structure_dict[ - sparCapKeys[hpSideIndex] - ]["start_nd_arc"]["values"] - definition.sparcap_end_nd_arc = blade_structure_dict[sparCapKeys[hpSideIndex]][ - "end_nd_arc" - ]["values"] + pass return definition -def _add_te_bands(definition, blade_structure_dict): - teReinfKeys = full_keys_from_substrings(blade_structure_dict.keys(), ["te", "reinf"]) - if len(teReinfKeys) == 1: - definition.teband = ( - blade_structure_dict[teReinfKeys[0]]["width"]["values"] * 1000 / 2 - ) - elif len(teReinfKeys) == 2: +def _add_te_le_bands(definition, blade_structure_dict,te_or_le): + key_names = full_keys_from_substrings(blade_structure_dict.keys(), [te_or_le, "reinf"]) + + if len(key_names) == 1: + definition.teband = (blade_structure_dict[key_names[0]]["width"]["values"] * 1000 / 2) + + # le_array = np.ones(np.shape(definition.key_arcs)[1],)*0.5 #Approximate location of LE + + definition.key_arcs.append(blade_structure_dict[key_names[0]]["start_nd_arc"]["values"]) + # definition.key_arcs.append(le_array) + # definition.key_arcs.append(le_array) + definition.key_arcs.append(blade_structure_dict[key_names[0]]["end_nd_arc"]["values"]) + + elif len(key_names) == 2: definition.teband = ( - ( - blade_structure_dict[teReinfKeys[0]]["width"]["values"] - + blade_structure_dict[teReinfKeys[1]]["width"]["values"] - ) - * 1000 - / 2 - ) - else: - raise ValueError("Unknown number of TE reinforcements") - return definition + (blade_structure_dict[key_names[0]]["width"]["values"] + + blade_structure_dict[key_names[1]]["width"]["values"])* 1000/ 2) + + for i_side in range(2): + if "ss" in blade_structure_dict[key_names[i_side]]["name"].lower(): + lp_side_index = i_side + if "ps" in blade_structure_dict[key_names[i_side]]["name"].lower(): + hp_side_index = i_side -def _add_le_bands(definition, blade_structure_dict): - leReinfKeys = full_keys_from_substrings(blade_structure_dict.keys(), ["le", "reinf"]) - if len(leReinfKeys) == 1: - definition.leband = ( - blade_structure_dict[leReinfKeys[0]]["width"]["values"] * 1000 / 2 - ) - elif len(leReinfKeys) == 2: - definition.leband = ( - ( - blade_structure_dict[leReinfKeys[0]]["width"]["values"] - + blade_structure_dict[leReinfKeys[1]]["width"]["values"] - ) - * 1000 - / 2 - ) + definition.key_arcs.append(blade_structure_dict[key_names[lp_side_index]]["end_nd_arc"]["values"]) + definition.key_arcs.append(blade_structure_dict[key_names[hp_side_index]]["start_nd_arc"]["values"]) + + + # definition.key_arcs.append(blade_structure_dict[key_names[1]]["start_nd_arc"]["values"])) + # definition.key_arcs.append(blade_structure_dict[key_names[1]]["end_nd_arc"]["values"])) else: - raise ValueError("Invalid number of LE reinforcements") + raise ValueError(f"Unknown number of {te_or_le} reinforcements") return definition + + +# def _add_le_bands(definition, blade_structure_dict): +# leReinfKeys = full_keys_from_substrings(blade_structure_dict.keys(), ["le", "reinf"]) + + + +# if len(leReinfKeys) == 1: +# definition.leband = (blade_structure_dict[leReinfKeys[0]]["width"]["values"] * 1000 / 2) + +# definition.key_arcs.append(blade_structure_dict[leReinfKeys[0]]["start_nd_arc"]["values"])) +# definition.key_arcs.append(blade_structure_dict[leReinfKeys[0]]["end_nd_arc"]["values"])) + +# elif len(leReinfKeys) == 2: +# definition.leband = ( +# (blade_structure_dict[leReinfKeys[0]]["width"]["values"] +# + blade_structure_dict[leReinfKeys[1]]["width"]["values"])* 1000/ 2) + + +# definition.key_arcs.append(blade_structure_dict[leReinfKeys[0]]["start_nd_arc"]["values"])) +# definition.key_arcs.append(blade_structure_dict[leReinfKeys[1]]["end_nd_arc"]["values"])) +# else: +# raise ValueError("Invalid number of LE reinforcements") +# return definition diff --git a/src/pynumad/objects/definition.py b/src/pynumad/objects/definition.py index 16e1606..d581175 100644 --- a/src/pynumad/objects/definition.py +++ b/src/pynumad/objects/definition.py @@ -77,6 +77,9 @@ def __init__(self): self.sweep: ndarray = None self.teband: float = None self.leband: float = None + self.n_webs: float = None + self.key_arcs: ndarray = None + # init properties self._natural_offset: int = 1 diff --git a/src/pynumad/objects/geometry.py b/src/pynumad/objects/geometry.py index 4867d36..70f6a1b 100644 --- a/src/pynumad/objects/geometry.py +++ b/src/pynumad/objects/geometry.py @@ -292,10 +292,10 @@ def generate(self, definition): self.arclength[:, k] = self.arclength[:, k] - LEarcsum # find where x=0 intersects the surface - self.HParcx0[0, k] = ( + self.LParcx0[0, k] = ( interpolator_wrap(xx[1 : le_index + 1], arclen[1 : le_index + 1], 0) - LEarcsum ) - self.LParcx0[0, k] = ( + self.HParcx0[0, k] = ( interpolator_wrap(xx[-2 : le_index - 1 : -1], arclen[-2 : le_index - 1 : -1], 0) - LEarcsum ) diff --git a/src/pynumad/objects/keypoints.py b/src/pynumad/objects/keypoints.py index 338e5c7..54efb47 100644 --- a/src/pynumad/objects/keypoints.py +++ b/src/pynumad/objects/keypoints.py @@ -6,6 +6,7 @@ from pynumad.utils.interpolation import interpolator_wrap from pynumad.objects.definition import Definition from pynumad.objects.geometry import Geometry +from pynumad.utils.misc_utils import find_frist_fully_populated_row @@ -38,21 +39,7 @@ class KeyPoints: """ def __init__(self): - self.key_labels = [ - "te", - "e", - "d", - "c", - "b", - "a", - "le", - "a", - "b", - "c", - "d", - "e", - "te", - ] + self.key_labels = [] self.key_points: ndarray = None self.key_arcs: ndarray = None self.key_cpos: ndarray = None @@ -116,14 +103,54 @@ def generate(self, definition: Definition, geometry: Geometry): num_istations = geometry.ispan.size # number of areas around airfoil profile; must be even (see calc of web areas) - num_areas = 12 + _,num_defined_arcs = np.shape(definition.key_arcs) + + num_areas = num_defined_arcs+4 + num_kp = num_defined_arcs+2 + self.initialize(num_areas, num_istations) # initialize keypoints + # start and finish indices in geometry/arcs ns = 1 nf = geometry.coordinates.shape[0] - 2 + + + i_start = find_frist_fully_populated_row(definition.key_arcs) + arc_copy = definition.key_arcs[i_start,:] + + for k_station in range(i_start): + + for i_arc in range(num_defined_arcs): + if definition.key_arcs[k_station,i_arc] == 0: + definition.key_arcs[k_station,i_arc] = arc_copy[i_arc]#*total_arc_len/arc_copy_arclen + + + i_start = find_frist_fully_populated_row(np.flip(definition.key_arcs)) + arc_copy = definition.key_arcs[i_start-1,:] + + for k_station in range(i_start,num_istations): + for i_arc in range(num_defined_arcs): + if definition.key_arcs[k_station,i_arc] == 0: + definition.key_arcs[k_station,i_arc] = arc_copy[i_arc]#*total_arc_len/arc_copy_arclen + + if num_areas>12: #Make third web equidistant (target_percent = 0.5)from spar cap and te reinf. + target_percent =0.5 + for k_station in range(i_start,num_istations): + #LP + prior_percent = (definition.key_arcs[k_station-1,1]-definition.key_arcs[k_station-1,0])/(definition.key_arcs[k_station-1,2]-definition.key_arcs[k_station-1,0]) + next_percent = min(prior_percent+0.2*(target_percent-prior_percent),0.5) + definition.key_arcs[k_station,1] = (definition.key_arcs[k_station,2]-definition.key_arcs[k_station,0])*next_percent+definition.key_arcs[k_station,0] + + #HP + prior_percent = (definition.key_arcs[k_station-1,-2]-definition.key_arcs[k_station-1,-1])/(definition.key_arcs[k_station-1,-3]-definition.key_arcs[k_station-1,-1]) + next_percent = min(prior_percent+0.2*(target_percent-prior_percent),0.5) + definition.key_arcs[k_station,-2] = (definition.key_arcs[k_station,-3]-definition.key_arcs[k_station,-1])*next_percent+definition.key_arcs[k_station,-1] #LP + # definition.key_arcs[k_station,-2] = (definition.key_arcs[k_station,-1]+definition.key_arcs[k_station,-3])/2.0 #HP + + # keypoints, keyarcs, keycpos te_types = [] # reset te_type i_leband_start=np.min(np.nonzero(definition.leband)) @@ -132,11 +159,11 @@ def generate(self, definition: Definition, geometry: Geometry): i_leband_end=np.max(np.nonzero(definition.leband)) i_teband_end=np.max(np.nonzero(definition.teband)) - for k in range(num_istations): + for k_station in range(num_istations): # allow for separate definitions of HP and LP spar cap # width and offset [HP LP] - n1 = mm_to_m * definition.leband[k] # no foam width - n2 = mm_to_m * definition.teband[k] # no foam width + n1 = mm_to_m * definition.leband[k_station] # no foam width + n2 = mm_to_m * definition.teband[k_station] # no foam width ### #In order to avoid abrupt changes in geometry when the le/te bands @@ -145,10 +172,10 @@ def generate(self, definition: Definition, geometry: Geometry): #le/te band widths since it is based on nonzero values. - if k < i_leband_start : + if k_station < i_leband_start : n1= mm_to_m * definition.leband[i_leband_start] - if k < i_teband_start : + if k_station < i_teband_start : n2= mm_to_m * definition.teband[i_teband_start] #In order to avoid abrupt changes in geometry when the le/te bands @@ -156,22 +183,22 @@ def generate(self, definition: Definition, geometry: Geometry): #value. This algorithm will not work as well if small numbers exist in the #le/te band widths since it is based on nonzero values. - if k > i_leband_end : - n1= mm_to_m * definition.leband[k-1]*0.75 - definition.leband[k]=n1/mm_to_m + if k_station > i_leband_end : + n1= mm_to_m * definition.leband[k_station-1]*0.75 + definition.leband[k_station]=n1/mm_to_m - if k > i_teband_end : - n2= mm_to_m * definition.teband[k-1]*0.75 - definition.teband[k]=n2/mm_to_m + if k_station > i_teband_end : + n2= mm_to_m * definition.teband[k_station-1]*0.75 + definition.teband[k_station]=n2/mm_to_m ### - scwidth_hp = mm_to_m * definition.sparcapwidth_hp[k] # type: float - scwidth_lp = mm_to_m * definition.sparcapwidth_lp[k] # type: float + # scwidth_hp = mm_to_m * definition.sparcapwidth_hp[k_station] # type: float + # scwidth_lp = mm_to_m * definition.sparcapwidth_lp[k_station] # type: float - scoffset_hp = mm_to_m * definition.sparcapoffset_hp[k] # type: float - scoffset_lp = mm_to_m * definition.sparcapoffset_lp[k] # type: float + # scoffset_hp = mm_to_m * definition.sparcapoffset_hp[k_station] # type: float + # scoffset_lp = mm_to_m * definition.sparcapoffset_lp[k_station] # type: float - tempTE = geometry.get_profile_te_type(k) + tempTE = geometry.get_profile_te_type(k_station) if te_types: te_types.append(tempTE) else: @@ -181,118 +208,147 @@ def generate(self, definition: Definition, geometry: Geometry): # get angle of each xy pair w.r.t. pitch axis (0,0) xyangle = np.zeros(geometry.coordinates.shape[0]) for j in range(len(xyangle)): - xy = geometry.coordinates[j, 0:2, k] + xy = geometry.coordinates[j, 0:2, k_station] xyangle[j] = np.arctan2(definition.rotorspin * xy[1], xy[0]) # unwrap and center around 0 xyangle = np.unwrap(xyangle) xyangle = xyangle - np.pi * np.round(xyangle[self.LEindex] / np.pi) - k_arclen = geometry.arclength[ns : nf + 1, k] - k_geom = geometry.coordinates[ns : nf + 1, :, k] - k_cpos = geometry.cpos[ns : nf + 1, k] + k_arclen = geometry.arclength[ns : nf + 1, k_station] + k_geom = geometry.coordinates[ns : nf + 1, :, k_station] + k_cpos = geometry.cpos[ns : nf + 1, k_station] # ==================== HP surface ==================== if definition.swtwisted: # find arclength where xyangle equals normal to chord # angle normal to chord line - twistnorm = np.pi / 180 * (-self.idegreestwist[k] - 90) + twistnorm = np.pi / 180 * (-self.idegreestwist[k_station] - 90) z = interpolator_wrap(xyangle[ns : nf + 1], k_arclen, twistnorm) else: - z = geometry.HParcx0[0, k] - z0 = z - z = z - scoffset_hp - a = np.amax(((0 - n1), 0.1 * geometry.arclength[ns, k])) # type: float - a = np.amin((a, 0.01 * geometry.arclength[ns, k])) - b = np.amin(((z + 0.5 * scwidth_hp), 0.15 * geometry.arclength[ns, k])) - c = np.amax(((z - 0.5 * scwidth_hp), 0.8 * geometry.arclength[ns, k])) - d = np.amin( - ((geometry.arclength[0, k] + n2), 0.85 * geometry.arclength[ns, k]) - ) - d = np.amax((d, 0.98 * geometry.arclength[ns, k])) - if str(te_types[k]) == "flat": - e = geometry.arclength[ns, k] - self.key_points[0, :, k] = geometry.coordinates[ns, :, k] - self.key_cpos[1, k] = -1 + z = geometry.HParcx0[0, k_station] + # z0 = z + # z = z - scoffset_hp + # a = np.amax(((0 - n1), 0.1 * geometry.arclength[ns, k_station])) # type: float + # a = np.amin((a, 0.01 * geometry.arclength[ns, k_station])) + # b = np.amin(((z + 0.5 * scwidth_hp), 0.15 * geometry.arclength[ns, k_station])) + # c = np.amax(((z - 0.5 * scwidth_hp), 0.8 * geometry.arclength[ns, k_station])) + # d = np.amin( + # ((geometry.arclength[0, k_station] + n2), 0.85 * geometry.arclength[ns, k_station]) + # ) + # d = np.amax((d, 0.98 * geometry.arclength[ns, k_station])) + if str(te_types[k_station]) == "flat": + e = geometry.arclength[ns, k_station] + self.key_points[0, :, k_station] = geometry.coordinates[ns, :, k_station] + self.key_cpos[1, k_station] = -1 else: - # e = 0.5 * (d + geometry.arclength(ns,k)); - e = 0.99 * geometry.arclength[ns, k] - self.key_points[0, :, k] = interpolator_wrap(k_arclen, k_geom, e) - self.key_cpos[1, k] = interpolator_wrap(k_arclen, k_cpos, e) + # e = 0.5 * (d + geometry.arclength(ns,k_station)); + e = 0.99 * geometry.arclength[ns, k_station] + self.key_points[0, :, k_station] = interpolator_wrap(k_arclen, k_geom, e) + self.key_cpos[1, k_station] = interpolator_wrap(k_arclen, k_cpos, e) # 1 -> e - self.key_points[1, :, k] = interpolator_wrap(k_arclen, k_geom, d) - self.key_points[2, :, k] = interpolator_wrap(k_arclen, k_geom, c) - # self.key_points( ,:,k) = interpolator_wrap(geometry.arclength(ns:nf,k),self.geometry(ns:nf,:,k),z); - self.key_points[3, :, k] = interpolator_wrap(k_arclen, k_geom, b) - self.key_points[4, :, k] = interpolator_wrap(k_arclen, k_geom, a) - self.key_arcs[0, k] = geometry.arclength[ns, k] - self.key_arcs[1, k] = e - self.key_arcs[2, k] = d - self.key_arcs[3, k] = c - # self.key_arcs( ,k) = z; - self.key_arcs[4, k] = b - self.key_arcs[5, k] = a - self.key_arcs[6, k] = 0 # le - self.key_cpos[0, k] = geometry.cpos[ns, k] # te, hp surface + arclen = k_arclen-k_arclen[0] + total_arc_len = (arclen)[-1] + + LEarcsum = arclen[geometry.LEindex] + reversed_scaled_yaml_key_arcs= np.flip((np.ones(np.shape(definition.key_arcs[k_station,:]))-definition.key_arcs[k_station,:])*total_arc_len) - LEarcsum + + #Make sure LE band is centered about leading edge + le_band = reversed_scaled_yaml_key_arcs[int(num_areas/2)-2]-reversed_scaled_yaml_key_arcs[int(num_areas/2)-3] + reversed_scaled_yaml_key_arcs[int(num_areas/2)-3]=-le_band/2 #HP side + reversed_scaled_yaml_key_arcs[int(num_areas/2)-2]=le_band/2 #LP side + + self.key_arcs[0, k_station] = geometry.arclength[ns, k_station] + self.key_arcs[1, k_station] = e + self.key_arcs[-1, k_station] = geometry.arclength[nf, k_station] + + + self.key_cpos[0, k_station] = geometry.cpos[ns, k_station] # te, hp surface + self.key_cpos[-1, k_station] = geometry.cpos[nf, k_station] # te, lp surface + + for i_kp in range(1,len(reversed_scaled_yaml_key_arcs) +1): + + self.key_points[i_kp, :, k_station] = interpolator_wrap(k_arclen, k_geom, reversed_scaled_yaml_key_arcs[i_kp-1]) + self.key_arcs[i_kp+1, k_station] = reversed_scaled_yaml_key_arcs[i_kp-1] + self.key_cpos[i_kp+1, k_station] = interpolator_wrap(k_arclen, k_cpos, reversed_scaled_yaml_key_arcs[i_kp-1]) + + self.key_cpos[int(num_areas/2)+1:num_areas-2+1, k_station] = self.key_cpos[int(num_areas/2):num_areas-2, k_station] #Shift second half by one + self.key_cpos[int(num_areas/2), k_station] = 0.0 + + self.key_arcs[int(num_areas/2)+1:num_areas-2+1, k_station] = self.key_arcs[int(num_areas/2):num_areas-2, k_station] #Shift second half by one + self.key_arcs[int(num_areas/2), k_station] = 0.0 + # self.key_points[1, :, k_station] = interpolator_wrap(k_arclen, k_geom, d) + # self.key_points[2, :, k_station] = interpolator_wrap(k_arclen, k_geom, c) + # # self.key_points( ,:,k_station) = interpolator_wrap(geometry.arclength(ns:nf,k_station),self.geometry(ns:nf,:,k_station),z); + # self.key_points[3, :, k_station] = interpolator_wrap(k_arclen, k_geom, b) + # self.key_points[4, :, k_station] = interpolator_wrap(k_arclen, k_geom, a) + # self.key_arcs[0, k_station] = geometry.arclength[ns, k_station] + # self.key_arcs[1, k_station] = e + # self.key_arcs[2, k_station] = d + # self.key_arcs[3, k_station] = c + # # self.key_arcs( ,k_station) = z; + # self.key_arcs[4, k_station] = b + # self.key_arcs[5, k_station] = a + # self.key_arcs[6, k_station] = 0 # le + # self.key_cpos[0, k_station] = geometry.cpos[ns, k_station] # te, hp surface # 2 -> e - self.key_cpos[2, k] = interpolator_wrap(k_arclen, k_cpos, d) - self.key_cpos[3, k] = interpolator_wrap(k_arclen, k_cpos, c) - # self.key_cpos( ,k) = interpolator_wrap(geometry.arclength(ns:nf,k),self.cpos(ns:nf,k),z); - self.key_cpos[4, k] = interpolator_wrap(k_arclen, k_cpos, b) - self.key_cpos[5, k] = interpolator_wrap(k_arclen, k_cpos, a) - self.key_cpos[6, k] = interpolator_wrap(k_arclen, k_cpos, 0) + # self.key_cpos[2, k_station] = interpolator_wrap(k_arclen, k_cpos, d) + # self.key_cpos[3, k_station] = interpolator_wrap(k_arclen, k_cpos, c) + # # self.key_cpos( ,k_station) = interpolator_wrap(geometry.arclength(ns:nf,k_station),self.cpos(ns:nf,k_station),z); + # self.key_cpos[4, k_station] = interpolator_wrap(k_arclen, k_cpos, b) + # self.key_cpos[5, k_station] = interpolator_wrap(k_arclen, k_cpos, a) + # self.key_cpos[6, k_station] = interpolator_wrap(k_arclen, k_cpos, 0) # ==================== LP surface ==================== - if definition.swtwisted: - # angle normal to chord line - twistnorm = np.pi / 180 * (-self.idegreestwist[k] + 90) - z = interpolator_wrap(xyangle[ns : nf + 1], k_arclen, twistnorm) + # if definition.swtwisted: + # # angle normal to chord line + # twistnorm = np.pi / 180 * (-self.idegreestwist[k_station] + 90) + # z = interpolator_wrap(xyangle[ns : nf + 1], k_arclen, twistnorm) + # else: + # z = geometry.LParcx0[0, k_station] + # z0 = z # ble: location where airfoil surface crosses Xglobal=0 + # z = z + scoffset_lp # positive scoffset moves z toward t.e. + # a = np.amin(((0 + n1), 0.1 * geometry.arclength[nf, k_station])) + # a = np.amax((a, 0.01 * geometry.arclength[nf, k_station])) + # b = np.amax(((z - 0.5 * scwidth_lp), 0.15 * geometry.arclength[nf, k_station])) + # c = np.amin((z + 0.5 * scwidth_lp, 0.8 * geometry.arclength[nf, k_station])) + # d = np.amax( + # (geometry.arclength[-1, k_station] - n2, 0.85 * geometry.arclength[nf, k_station]) + # ) + # d = np.amin((d, 0.96 * geometry.arclength[nf, k_station])) + if str(te_types[k_station]) == str("flat"): + e = geometry.arclength[nf, k_station] + self.key_points[-1, :, k_station] = geometry.coordinates[nf, :, k_station] + self.key_cpos[-2, k_station] = 1 else: - z = geometry.LParcx0[0, k] - z0 = z # ble: location where airfoil surface crosses Xglobal=0 - z = z + scoffset_lp # positive scoffset moves z toward t.e. - a = np.amin(((0 + n1), 0.1 * geometry.arclength[nf, k])) - a = np.amax((a, 0.01 * geometry.arclength[nf, k])) - b = np.amax(((z - 0.5 * scwidth_lp), 0.15 * geometry.arclength[nf, k])) - c = np.amin((z + 0.5 * scwidth_lp, 0.8 * geometry.arclength[nf, k])) - d = np.amax( - (geometry.arclength[-1, k] - n2, 0.85 * geometry.arclength[nf, k]) - ) - d = np.amin((d, 0.96 * geometry.arclength[nf, k])) - if str(te_types[k]) == str("flat"): - e = geometry.arclength[nf, k] - self.key_points[9, :, k] = geometry.coordinates[nf, :, k] - self.key_cpos[11, k] = 1 - else: - # e = 0.5 * (d + geometry.arclength(nf,k)); - e = 0.98 * geometry.arclength[nf, k] - self.key_points[9, :, k] = interpolator_wrap(k_arclen, k_geom, e) - self.key_cpos[11, k] = interpolator_wrap(k_arclen, k_cpos, e) - self.key_points[5, :, k] = interpolator_wrap(k_arclen, k_geom, a) - self.key_points[6, :, k] = interpolator_wrap(k_arclen, k_geom, b) - # self.key_points( ,:,k) = interpolator_wrap(geometry.arclength(ns:nf,k),self.geometry(ns:nf,:,k),z); - self.key_points[7, :, k] = interpolator_wrap(k_arclen, k_geom, c) - self.key_points[8, :, k] = interpolator_wrap(k_arclen, k_geom, d) + # e = 0.5 * (d + geometry.arclength(nf,k_station)); + e = 0.98 * geometry.arclength[nf, k_station] + self.key_points[-1, :, k_station] = interpolator_wrap(k_arclen, k_geom, e) + self.key_cpos[-2, k_station] = interpolator_wrap(k_arclen, k_cpos, e) + # self.key_points[5, :, k_station] = interpolator_wrap(k_arclen, k_geom, a) + # self.key_points[6, :, k_station] = interpolator_wrap(k_arclen, k_geom, b) + # # self.key_points( ,:,k_station) = interpolator_wrap(geometry.arclength(ns:nf,k_station),self.geometry(ns:nf,:,k_station),z); + # self.key_points[7, :, k_station] = interpolator_wrap(k_arclen, k_geom, c) + # self.key_points[8, :, k_station] = interpolator_wrap(k_arclen, k_geom, d) # 10 -> e - self.key_arcs[7, k] = a - self.key_arcs[8, k] = b - # self.key_arcs( ,k) = z; - self.key_arcs[9, k] = c - self.key_arcs[10, k] = d - self.key_arcs[11, k] = e - self.key_arcs[12, k] = geometry.arclength[nf, k] - self.key_cpos[7, k] = interpolator_wrap(k_arclen, k_cpos, a) - self.key_cpos[8, k] = interpolator_wrap(k_arclen, k_cpos, b) - # self.key_cpos( ,k) = interpolator_wrap(geometry.arclength(ns:nf,k),self.cpos(ns:nf,k),z); - self.key_cpos[9, k] = interpolator_wrap(k_arclen, k_cpos, c) - self.key_cpos[10, k] = interpolator_wrap(k_arclen, k_cpos, d) - # 12 -> e - self.key_cpos[12, k] = geometry.cpos[nf, k] # te, lp surface + # self.key_arcs[7, k_station] = a + # self.key_arcs[8, k_station] = b + # # self.key_arcs( ,k_station) = z; + # self.key_arcs[9, k_station] = c + # self.key_arcs[10, k_station] = d + self.key_arcs[-2, k_station] = e + # self.key_arcs[11, k_station] = e + # self.key_arcs[12, k_station] = geometry.arclength[nf, k_station] + # self.key_cpos[7, k_station] = interpolator_wrap(k_arclen, k_cpos, a) + # self.key_cpos[8, k_station] = interpolator_wrap(k_arclen, k_cpos, b) + # # self.key_cpos( ,k_station) = interpolator_wrap(geometry.arclength(ns:nf,k_station),self.cpos(ns:nf,k_station),z); + # self.key_cpos[9, k_station] = interpolator_wrap(k_arclen, k_cpos, c) + # self.key_cpos[10, k_station] = interpolator_wrap(k_arclen, k_cpos, d) + # # 12 -> e + # self.key_cpos[12, k_station] = geometry.cpos[nf, k_station] # te, lp surface # find the points used by each shear web - component_groups = [ - definition.components[name].group for name in definition.components - ] + component_groups = [definition.components[name].group for name in definition.components] self.web_indices = [] self.web_arcs = [] self.web_cpos = [] @@ -369,12 +425,12 @@ def generate(self, definition: Definition, geometry: Geometry): p2 = self.key_arcs[n2, :] p = (1 - f) * p1 + f * p2 self.web_arcs[ksw][0, :] = p - for k in range(num_istations): - self.web_cpos[ksw][0, k] = interpolator_wrap( - k_arclen, k_cpos, p[k] + for k_station in range(num_istations): + self.web_cpos[ksw][0, k_station] = interpolator_wrap( + k_arclen, k_cpos, p[k_station] ) - self.web_points[ksw][0, :, k] = interpolator_wrap( - k_arclen, k_geom, p[k] + self.web_points[ksw][0, :, k_station] = interpolator_wrap( + k_arclen, k_geom, p[k_station] ) elif hp["pt3"]: try: @@ -397,14 +453,14 @@ def generate(self, definition: Definition, geometry: Geometry): ) p[np.abs(p) < np.abs(pMin)] = pMin[np.abs(p) < np.abs(pMin)] self.web_cpos[ksw][0, :] = p - for k in range(num_istations): - self.web_arcs[ksw][0, k] = interpolator_wrap( - self.cpos[ns : nf + 1, :, k], - geometry.arclength[ns : nf + 1, :, k], - p[k], + for k_station in range(num_istations): + self.web_arcs[ksw][0, k_station] = interpolator_wrap( + self.cpos[ns : nf + 1, :, k_station], + geometry.arclength[ns : nf + 1, :, k_station], + p[k_station], ) - self.web_points[ksw][0, :, k] = interpolator_wrap( - k_cpos, k_geom, p[k] + self.web_points[ksw][0, :, k_station] = interpolator_wrap( + k_cpos, k_geom, p[k_station] ) else: raise Exception( @@ -440,12 +496,12 @@ def generate(self, definition: Definition, geometry: Geometry): p2 = self.key_arcs[n2, :] p = (1 - f) * p1 + f * p2 self.web_arcs[ksw][1, :] = p - for k in range(num_istations): - self.web_cpos[ksw][1, k] = interpolator_wrap( - k_arclen, k_cpos, p[k] + for k_station in range(num_istations): + self.web_cpos[ksw][1, k_station] = interpolator_wrap( + k_arclen, k_cpos, p[k_station] ) - self.web_points[ksw][1, :, k] = interpolator_wrap( - k_arclen, k_geom, p[k] + self.web_points[ksw][1, :, k_station] = interpolator_wrap( + k_arclen, k_geom, p[k_station] ) elif lp["pt3"]: try: @@ -466,12 +522,12 @@ def generate(self, definition: Definition, geometry: Geometry): ) p[np.abs(p) < np.abs(pMin)] = pMin[np.abs(p) < np.abs(pMin)] self.web_cpos[ksw][1, :] = p - for k in range(num_istations): - self.web_arcs[ksw][1, k] = interpolator_wrap( - k_cpos, k_arclen, p[k] + for k_station in range(num_istations): + self.web_arcs[ksw][1, k_station] = interpolator_wrap( + k_cpos, k_arclen, p[k_station] ) - self.web_points[ksw][1, :, k] = interpolator_wrap( - k_cpos, k_geom, p[k] + self.web_points[ksw][1, :, k_station] = interpolator_wrap( + k_cpos, k_geom, p[k_station] ) else: raise Exception( diff --git a/src/pynumad/objects/stackdb.py b/src/pynumad/objects/stackdb.py index a2a9f3d..d11a7d0 100644 --- a/src/pynumad/objects/stackdb.py +++ b/src/pynumad/objects/stackdb.py @@ -50,6 +50,16 @@ def generate(self, keypoints: KeyPoints, bom: BillOfMaterials): "LP_TE_REINF", "LP_TE_FLAT", ] + + if n_webs > 2: + temp_segment_labels = segment_labels.copy() + ct = 0 + for i_segment,segment_label in enumerate(temp_segment_labels): + prepend = segment_label.split('_')[0] + if 'TE_PANEL' in segment_label: + for i in range(n_webs-2): + segment_labels.insert(i_segment+1+ct, f'{prepend}_TE_PANEL') + ct+=1 # initialize stack array self.stacks = np.empty(shape=(n_segments, n_stations), dtype=object) diff --git a/src/pynumad/utils/misc_utils.py b/src/pynumad/utils/misc_utils.py index a50cd74..695e0ea 100644 --- a/src/pynumad/utils/misc_utils.py +++ b/src/pynumad/utils/misc_utils.py @@ -46,6 +46,22 @@ def full_keys_from_substrings(key_list, subtring_list,ignore_case=True): return matched_keys +def find_frist_fully_populated_row(matrix): + n_rows,_ = np.shape(matrix) + for i_row in range(n_rows): + n_zeros = len(np.where(matrix[i_row,:] == 0)[0]) + if n_zeros == 0: + break + if i_row == n_rows: + raise ValueError('Was not able to sort matrix. All rows had at least one zero.') + return i_row +def sort_based_on_root_values(matrix): + + i_index = find_frist_fully_populated_row(matrix) + values = matrix[i_index,:] + i_sorted = sorted(range(len(values)), key=lambda i: values[i]) + return matrix[:,i_sorted],i_sorted + # SED-like substitution def copy_and_replace(fin, fout, replacements): inf = open(fin, "r") From e32b2be93f769e2cf5ae9d798f984f179f64c312 Mon Sep 17 00:00:00 2001 From: Camarena Date: Tue, 17 Dec 2024 10:49:54 -0700 Subject: [PATCH 07/11] multi web capability from cubit --- src/pynumad/analysis/cubit/make_blade.py | 97 +++-- .../analysis/cubit/make_cross_sections.py | 405 +++++++++--------- 2 files changed, 261 insertions(+), 241 deletions(-) diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index b90a426..ec270c7 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -565,12 +565,13 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati all_layer_thicknesses = [[] for _ in range(3)] #three modeled layers - _,n_stations = np.shape(stackdb.stacks) + n_areas,n_stations = np.shape(stackdb.stacks) for i_station in range(n_stations): temp = stackdb.stacks[:, i_station] - temp = np.flip(temp) + # temp = np.flip(temp) - stacks = list(stackdb.stacks[1:6, i_station]) + list(temp[1:6]) + # stacks = list(stackdb.stacks[1:int(n_areas/2), i_station]) + list(temp[1:int(n_areas/2)]) + stacks = list(temp[1:int(n_areas/2)]) + list(temp[-2:int(n_areas/2)-1:-1]) for stack in stacks: for i_layer,layer_thicknesses in enumerate(all_layer_thicknesses): @@ -593,16 +594,20 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati - hasWebs = [] - webNumber = 1 - for i_station in range(len(stackdb.swstacks[webNumber])): - if not len(stackdb.swstacks[webNumber][i_station].plygroups) == 0: - hasWebs.append(True) - else: - hasWebs.append(False) + has_webs = [] + n_webs, n_stations = np.shape(stackdb.swstacks) + for i_station in range(n_stations-1): + temp=[] + for i_web in range(n_webs): + if not len(stackdb.swstacks[i_web][i_station].plygroups) == 0: + temp.append(True) + else: + temp.append(False) + has_webs.append(temp) # WARNING - Last station never has webs. Fix later - hasWebs.append(False) + has_webs.append([False for x in range(n_webs)]) + has_webs=np.array(has_webs) # WARNING - Last station never has webs. Fix later # Create Referece line as a spline @@ -723,29 +728,22 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati #is_flatback=is_station_flatback[i_station_geometry] + web_stacks = [] + for i_web in range(n_webs-1,-1,-1): #stackdb.swstacks is arranged from LE to TE. Need TE to LE order - i_station_first_web = np.argwhere(hasWebs)[0][0] - i_station_last_web = np.argwhere(hasWebs)[-1][0] + i_station_first_web = np.argwhere(has_webs[:,i_web])[0][0] + i_station_last_web = np.argwhere(has_webs[:,i_web])[-1][0] - if hasWebs[i_station] == True: - webNumber = 1 - aft_web_stack = stackdb.swstacks[webNumber][i_station] - webNumber = 0 - fore_web_stack = stackdb.swstacks[webNumber][i_station] - else: - if i_station < i_station_first_web: - iWebStation = i_station_first_web + if has_webs[i_station,i_web] == True: - # elif i_station_last_web == len(blade.ispan) - 1-1: + web_stacks.append(stackdb.swstacks[i_web,i_station]) else: - iWebStation = i_station_last_web - # else: - # raise Exception('assuming web ends at last station for now. ') + if i_station < i_station_first_web: + iWebStation = i_station_first_web + else: + iWebStation = i_station_last_web - webNumber = 1 - aft_web_stack = stackdb.swstacks[webNumber][iWebStation] - webNumber = 0 - fore_web_stack = stackdb.swstacks[webNumber][iWebStation] + web_stacks.append(stackdb.swstacks[i_web,iWebStation]) cs_normal = get_cs_normal_vector( np.array( @@ -757,9 +755,9 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati ) ) if model2Dor3D.lower() == "3d": - make_webs = True + make_webs = [True for x in range(n_webs)] else: - make_webs = hasWebs[i_station] + make_webs = np.flip(has_webs[i_station,:]) #Order from TE web to LE Web if 'precomp' in settings['make_input_for']: make_a_precomp_cross_section(wt_name, @@ -774,8 +772,7 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati i_station_geometry, blade, make_webs, - aft_web_stack, - fore_web_stack, + web_stacks, iLE, cs_params, geometry_scaling, @@ -792,8 +789,7 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati i_station_geometry, blade, make_webs, - aft_web_stack, - fore_web_stack, + web_stacks, iLE, cs_params, geometry_scaling, @@ -1034,7 +1030,7 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati # Remove unnecessary files to save space # for filePath in glob.glob(f"{path_name}-*.cub"): # os.remove(filePath) - return (cubit,blade,surface_dict,birds_mouth_verts,i_station_first_web,i_station_last_web,materials_used,spanwise_mat_ori_curve,hasWebs) + return (cubit,blade,surface_dict,birds_mouth_verts,i_station_first_web,i_station_last_web,materials_used,spanwise_mat_ori_curve,has_webs) def cubit_make_solid_blade( @@ -1071,7 +1067,7 @@ def cubit_make_solid_blade( raise ValueError("Need more than one cross section to make a solid model") (cubit,blade,surface_dict,birds_mouth_verts,i_station_first_web, - i_station_last_web,materials_used,spanwise_mat_ori_curve,hasWebs) = cubit_make_cross_sections( + i_station_last_web,materials_used,spanwise_mat_ori_curve,has_webs) = cubit_make_cross_sections( blade, wt_name, settings, cs_params, "3D", stationList) @@ -1090,13 +1086,15 @@ def cubit_make_solid_blade( else: shell_vol_list=[] - part_name = "web" - ordered_list = get_ordered_list(part_name) - ordered_list_web = ordered_list.copy() - if ordered_list and len(ordered_list[0]) > 1: - web_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, ordered_list, i_station_end,spanwise_splines) - else: - web_vol_list=[] + n_webs=len(blade.stackdb.swstacks) + for i_web in range(n_webs): + part_name = f'web{i_web}' + ordered_list = get_ordered_list(part_name) + ordered_list_web = ordered_list.copy() + if ordered_list and len(ordered_list[0]) > 1: + web_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, ordered_list, i_station_end,spanwise_splines) + else: + web_vol_list=[] part_name = "roundTEadhesive" ordered_list = get_ordered_list(part_name) @@ -1133,9 +1131,10 @@ def cubit_make_solid_blade( cubit.cmd(f"delete surface with Is_Free") for i_station in stationList[0:-1]: - #for i_keep,keep_web in enumerate(hasWebs): - if not hasWebs[i_station]: - cubit.cmd(f"delete volume with name 'web*tation{str(i_station).zfill(3)}*'") + temp = np.flip(has_webs[i_station]) + for i_web in range(n_webs): + if not temp[i_web]: + cubit.cmd(f"delete volume with name 'web{i_web}*tation{str(i_station).zfill(3)}*'") @@ -1155,7 +1154,7 @@ def cubit_make_solid_blade( elif float(cs_params['element_ar']) != 0.0: cubit.cmd("set default autosize on") omit_surf_mesh=[] - hplp_max_stack_ct = 14 #The number of stacks in HP surface. Hard code for now. + hplp_max_stack_ct = 14 + (n_webs - 2)#The number of stacks in HP surface. Hard code for now. #cubit.cmd(f'surface with name "*shell*Station*layer0_bottomFace*" scheme map') cubit.cmd(f"surf with name '*Station*surface*' scheme map") #All cross sections are map @@ -1269,7 +1268,7 @@ def cubit_make_solid_blade( parse_string = f'with name "hoop_direction{str(station_id).zfill(3)}_stack{str(i_stack).zfill(3)}*"' curve_ids = parse_cubit_list("curve", parse_string) - if i_stack == 12 or i_stack == 26: + if i_stack == 12+(n_webs-2) or i_stack == 26+2*(n_webs-2): cubit.cmd(f"curve {l2s(curve_ids)} interval {cs_params['nel_per_layer']}") else: cubit.cmd(f"curve {curve_ids[0]} scheme bias fine size {e_size_1} coarse size {e_size_1} ") ### SAME FOR NOW diff --git a/src/pynumad/analysis/cubit/make_cross_sections.py b/src/pynumad/analysis/cubit/make_cross_sections.py index 0fa0790..89e9e34 100644 --- a/src/pynumad/analysis/cubit/make_cross_sections.py +++ b/src/pynumad/analysis/cubit/make_cross_sections.py @@ -626,74 +626,95 @@ def split_curve_at_coordinte_points(coordinates_to_split_curve, curve_idToSplit) return key_curves -def split_key_curves(key_curves, aft_web_stack, fore_web_stack, web_adhesive_width): +def split_key_curves(key_curves, web_stacks, web_adhesive_width): ###Do not split TE reinf temp_base_curve_ids = [key_curves[0]] ###split TE panel in half - cubit.cmd(f"split curve {key_curves[1]} fraction 0.5") - - temp_base_curve_ids.append(get_last_id("curve") - 1) - temp_base_curve_ids.append(get_last_id("curve")) + if len(web_stacks) == 2: + cubit.cmd(f"split curve {key_curves[1]} fraction 0.5") + temp_base_curve_ids.append(get_last_id("curve") - 1) + temp_base_curve_ids.append(get_last_id("curve")) + else: + temp_base_curve_ids.append(key_curves[1]) - ###Partition sparcap curve - vertex_list = [] - web_layer_thickness = 0 - n_start = get_last_id("vertex") + 1 - for iLayer in reversed(range(len(aft_web_stack.plygroups))): - web_layer_thickness += ( - aft_web_stack.plygroups[iLayer].thickness - * aft_web_stack.plygroups[iLayer].nPlies - / 1000 - ) + ###Partition curves with webs for all but the fore web + overwrap_base_curves=[] + for i_web,web_stack in enumerate(web_stacks[:-1]): + vertex_list = [] + web_layer_thickness = 0 + n_start = get_last_id("vertex") + 1 + for i_layer in reversed(range(len(web_stack.plygroups))): + web_layer_thickness += ( + web_stack.plygroups[i_layer].thickness* web_stack.plygroups[i_layer].nPlies/ 1000) + cubit.cmd( + f"create vertex on curve {key_curves[i_web+2]} distance {web_layer_thickness} from start") cubit.cmd( - f"create vertex on curve {key_curves[2]} distance {web_layer_thickness} from start" - ) - cubit.cmd( - f"create vertex on curve {key_curves[2]} distance {web_layer_thickness+web_adhesive_width} from start" - ) + f"create vertex on curve {key_curves[i_web+2]} distance {web_layer_thickness+web_adhesive_width} from start") + + n_end = get_last_id("vertex") + + vertex_list += list(range(n_start, n_end + 1)) + + n_start = get_last_id("curve") + 1 + cubit.cmd(f"split curve {key_curves[i_web+2]} at vertex {l2s(vertex_list)}") + n_end = get_last_id("curve") + + cubit.cmd(f'curve {n_start} {n_start+2} rename "face_web_thickness"') + # cubit.cmd(f'curve {n_end-2} {n_end} rename "face_web_thickness"') + + cubit.cmd(f'curve {n_start+1} rename "core_web_thickness"') + # cubit.cmd(f'curve {n_end-1} rename "core_web_thickness"') + + # overwrap_base_curves += list(range(n_start + 1, n_end)) + # temp_base_curve_ids.append(n_start) + if i_web < len(web_stacks) -2: #Webs that sit on panels + overwrap_base_curves += list(range(n_start, n_end)) + temp_base_curve_ids.append(n_end) + else: + overwrap_base_curves += list(range(n_start + 1, n_end)) + temp_base_curve_ids.append(n_start) + key_curves[i_web+2]=n_end #aft panel sitting on sparcap + # get total foreweb thickness + vertex_list = [] + fore_web_stack = web_stacks[-1] web_layer_thickness = sum(fore_web_stack.layer_thicknesses()) / 1000 + n_start = get_last_id("vertex") + 1 cubit.cmd( - f"create vertex on curve {key_curves[2]} distance {web_layer_thickness+web_adhesive_width} from end" - ) + f"create vertex on curve {key_curves[i_web+2]} distance {web_layer_thickness+web_adhesive_width} from end") for iLayer in reversed(range(len(fore_web_stack.plygroups))): cubit.cmd( - f"create vertex on curve {key_curves[2]} distance {web_layer_thickness} from end" - ) + f"create vertex on curve {key_curves[i_web+2]} distance {web_layer_thickness} from end") web_layer_thickness -= ( - fore_web_stack.plygroups[iLayer].thickness - * fore_web_stack.plygroups[iLayer].nPlies - / 1000 - ) - + fore_web_stack.plygroups[iLayer].thickness* fore_web_stack.plygroups[iLayer].nPlies/ 1000) n_end = get_last_id("vertex") vertex_list += list(range(n_start, n_end + 1)) n_start = get_last_id("curve") + 1 - cubit.cmd(f"split curve {key_curves[2]} at vertex {l2s(vertex_list)}") + cubit.cmd(f"split curve {key_curves[i_web+2]} at vertex {l2s(vertex_list)}") n_end = get_last_id("curve") - cubit.cmd(f'curve {n_start} {n_start+2} rename "face_web_thickness"') + # cubit.cmd(f'curve {n_start} {n_start+2} rename "face_web_thickness"') cubit.cmd(f'curve {n_end-2} {n_end} rename "face_web_thickness"') - cubit.cmd(f'curve {n_start+1} rename "core_web_thickness"') + # cubit.cmd(f'curve {n_start+1} rename "core_web_thickness"') cubit.cmd(f'curve {n_end-1} rename "core_web_thickness"') - temp_base_curve_ids.append(n_start) + sparcap_base_curve = n_start temp_base_curve_ids.append(n_end) - spar_cap_base_curves = list(range(n_start + 1, n_end)) + overwrap_base_curves += list(range(n_start + 1, n_end)) ###split LE panel in half - cubit.cmd(f"split curve {key_curves[3]} fraction 0.5") + cubit.cmd(f"split curve {key_curves[-2]} fraction 0.5") temp_base_curve_ids.append(get_last_id("curve") - 1) temp_base_curve_ids.append(get_last_id("curve")) ###Do not split LE reinf temp_base_curve_ids.append(key_curves[-1]) - return temp_base_curve_ids, spar_cap_base_curves + return temp_base_curve_ids, overwrap_base_curves,sparcap_base_curve def get_mid_line(blade, iLE, i_station, geometry_scaling): @@ -751,6 +772,7 @@ def make_cs_perimeter_layer_areas(wt_name, i_station, station_stacks, cs_params, + has_webs, thickness_scaling, lp_hp_side, last_round_station, @@ -783,13 +805,18 @@ def make_cs_perimeter_layer_areas(wt_name, ) base_curve_index_ct = 0 - station_stacks.shape - nStationLayups = len(station_stacks) - station_stacks.shape + n_perimeter_layups = int(len(lp_hp_dict["base_curve_ids"][0])/2) - last_perimeter = nStationLayups - 2 + last_perimeter = n_perimeter_layups - 1 - for i_perimeter in range(nStationLayups - 1): # Skip the last stack since the current and the next stack are generated at the same time. + if len(station_stacks) > 5: + if len(station_stacks) ==6: + station_stacks=list(station_stacks) + station_stacks.pop(2) #Remove reduntant panel stack. + else: + raise ValueError('More than three webs is unsupported.') + + for i_perimeter in range(n_perimeter_layups): # Skip the last stack since the current and the next stack are generated at the same time. with open(f"{wt_name}.log", "a") as logFile: logFile.write(f"\tlp_hp_side {lp_hp_side}, i_perimeter={i_perimeter}\n") @@ -800,13 +827,12 @@ def make_cs_perimeter_layer_areas(wt_name, next_stack_layer_thicknesses = np.array(next_stack.layer_thicknesses()) / 1000 cubit.cmd( - f'curve {lp_hp_dict["base_curve_ids"][lp_hp_side_index][base_curve_index_ct]} copy' - ) + f'curve {lp_hp_dict["base_curve_ids"][lp_hp_side_index][base_curve_index_ct]} copy') current_base_curve_id = get_last_id("curve") base_curve_index_ct += 1 cubit.cmd( - f'curve {lp_hp_dict["base_curve_ids"][lp_hp_side_index][base_curve_index_ct]} copy' - ) + f'curve {lp_hp_dict["base_curve_ids"][lp_hp_side_index][base_curve_index_ct]} copy') + next_base_curve_id = get_last_id("curve") base_curve_index_ct += 1 @@ -876,9 +902,7 @@ def make_cs_perimeter_layer_areas(wt_name, else: raise ValueError(f"i_perimeter {i_perimeter} not recognized") - bottom_left_vertex_curve_left, bottom_right_vertex_curve_left = selCurveVerts( - left_bottom_curve - ) + bottom_left_vertex_curve_left, bottom_right_vertex_curve_left = selCurveVerts(left_bottom_curve) bottom_left_vertex_curve_right, bottom_right_vertex_curve_right = selCurveVerts(right_bottom_curve) # This if statement prepares all layer curves such that they taper at the TE @@ -1380,52 +1404,47 @@ def make_cs_perimeter_layer_areas(wt_name, bottom_left_vertex_curve_right = top_left_vertex_curve_right bottom_right_vertex_curve_right = top_right_vertex_curve_right + stack_ct+=3 # Build spar caps - if i_perimeter == 1: - lp_hp_dict["web_interface_curves"][lp_hp_side_index] = [right_top_curve] - temp_ct=0 - for ic, current_curveID in enumerate( - lp_hp_dict["spar_cap_base_curves"][lp_hp_side_index] - ): - bottom_curve = current_curveID - offSetSign = get_curve_offset_direction( - bottom_curve, lp_hp_side, cs_normal - ) - temp_ct+=1 - for it, thickness in enumerate(next_stack_layer_thicknesses): - cubit.cmd( - f"create curve offset curve {bottom_curve} distance {offSetSign*thickness} extended" - ) - top_curve = get_last_id("curve") + if i_perimeter==0 and len(lp_hp_dict["overwrap_base_curves"][lp_hp_side_index])==10: #Build web base on TE panel + stack_ct+=1 + curve_ids = lp_hp_dict["overwrap_base_curves"][lp_hp_side_index][:4] + top_curves, part_name_id = make_stack_of_surfaces_from_curve_list(curve_ids,next_stack,lp_hp_side,cs_normal,surface_dict,i_station,part_name,stack_ct,part_name_id,lp_hp_dict,materials_used) + i_web = 0 + if has_webs[i_web]: + lp_hp_dict["web_interface_curves"][lp_hp_side][i_web]+=top_curves + + if i_perimeter==1: - material_name = next_stack.plygroups[it].materialid - ply_angle = next_stack.plygroups[it].angle - part_name_id, materials_used = make_cross_section_surface(lp_hp_side, - surface_dict, - i_station, - part_name, - top_curve, - bottom_curve, - material_name, - ply_angle, - part_name_id, - it, - materials_used, - stack_ct+temp_ct - ) - next_stack_surface_list.append(get_last_id("surface")) + i_web =-2 + #Build web base on aft sparcap edge + lp_hp_dict["web_interface_curves"][lp_hp_side][i_web].append(right_top_curve) + curve_ids = lp_hp_dict["overwrap_base_curves"][lp_hp_side_index][-6:-3] + top_curves, part_name_id = make_stack_of_surfaces_from_curve_list(curve_ids,next_stack,lp_hp_side,cs_normal,surface_dict,i_station,part_name,stack_ct,part_name_id,lp_hp_dict,materials_used) + + if has_webs[i_web]: + lp_hp_dict["web_interface_curves"][lp_hp_side][i_web]+=top_curves + stack_ct+=1 - if it == 2 and ic != 3: - lp_hp_dict["web_interface_curves"][lp_hp_side_index].append( - top_curve - ) - bottom_curve = top_curve - + #Build center region of sparcap + curve_ids = [lp_hp_dict["sparcap_base_curve"][lp_hp_side_index]] + top_curves, part_name_id = make_stack_of_surfaces_from_curve_list(curve_ids,next_stack,lp_hp_side,cs_normal,surface_dict,i_station,part_name,stack_ct,part_name_id,lp_hp_dict,materials_used) stack_ct+=1 + + #Build web base on fore sparcap edge + i_web =-1 + curve_ids = lp_hp_dict["overwrap_base_curves"][lp_hp_side_index][-3:] + top_curves, part_name_id = make_stack_of_surfaces_from_curve_list(curve_ids,next_stack,lp_hp_side,cs_normal,surface_dict,i_station,part_name,stack_ct,part_name_id,lp_hp_dict,materials_used) + if has_webs[i_web]: + top_curves.reverse() + lp_hp_dict["web_interface_curves"][lp_hp_side][i_web]+=top_curves + stack_ct-=1 + elif i_perimeter == 2: - lp_hp_dict["web_interface_curves"][lp_hp_side_index].append(leftTopCurve) + i_web =-1 + lp_hp_dict["web_interface_curves"][lp_hp_side][i_web].insert(0,leftTopCurve) - stack_ct+=3 + return part_name_id, lp_hp_dict, stack_ct @@ -1435,7 +1454,39 @@ def make_cs_perimeter_layer_areas(wt_name, #################################################### #################################################### +def make_stack_of_surfaces_from_curve_list(curve_ids,stack_data,lp_hp_side,cs_normal,surface_dict,i_station,part_name,stack_ct,part_name_id,lp_hp_dict,materials_used): + top_curves=[] + layer_thicknesses = np.array(stack_data.layer_thicknesses()) / 1000 + # temp_ct=0 + for ic, current_curveID in enumerate(curve_ids): + bottom_curve = current_curveID + offSetSign = get_curve_offset_direction(bottom_curve, lp_hp_side, cs_normal) + # temp_ct+=1 + for it, thickness in enumerate(layer_thicknesses): + cubit.cmd( + f"create curve offset curve {bottom_curve} distance {offSetSign*thickness} extended") + top_curve = get_last_id("curve") + + material_name = stack_data.plygroups[it].materialid + ply_angle = stack_data.plygroups[it].angle + part_name_id, materials_used = make_cross_section_surface(lp_hp_side, + surface_dict, + i_station, + part_name, + top_curve, + bottom_curve, + material_name, + ply_angle, + part_name_id, + it, + materials_used, + stack_ct) + + bottom_curve = top_curve + + top_curves.append(top_curve) + return top_curves,part_name_id def create_simplist_surface_for_TE_or_LE_adhesive( i_station, surface_dict, @@ -1510,53 +1561,37 @@ def print_sine_curve_between_two_verts(vBot, vTop, amplitude, direction): return get_last_id("curve"),vertex_list -def make_cs_web_layer_areas( +def make_cs_web_layer_areas(part_name, surface_dict, i_station, - aft_web_stack, - fore_web_stack, + web_stack, web_interface_curves, cs_params, part_name_id, cs_normal, n_modeled_layers, - materials_used, -): + materials_used): - aft_web_overwrap_thickness = ( - aft_web_stack.layer_thicknesses()[0] + aft_web_stack.layer_thicknesses()[-1] - ) / 1000 - fore_web_overwrap_thickness = ( - fore_web_stack.layer_thicknesses()[0] + fore_web_stack.layer_thicknesses()[-1] - ) / 1000 - part_name = "web" + overwrap_thickness = ( + web_stack.layer_thicknesses()[0] + web_stack.layer_thicknesses()[-1]) / 1000 + + lp_hp_side_dict = {0:'HP',1:'LP'} ### First create the first two layers. The first layer is the adhesive. The second layer is the web overwrap layer for i_curveList, curveList in enumerate(web_interface_curves): n_base_curves_web = len(curveList) - if i_curveList == 0: - lp_hp_side = "HP" - else: - lp_hp_side = "LP" + lp_hp_side = lp_hp_side_dict[i_curveList] + for i_curve, bottom_curve in enumerate(curveList): - offSetSign = get_curve_offset_direction( - bottom_curve, lp_hp_side, cs_normal - ) + + offSetSign = get_curve_offset_direction(bottom_curve, lp_hp_side, cs_normal) - if i_curve < n_base_curves_web / 2: - layer_thicknesses = [ - cs_params["web_aft_adhesive_thickness"][i_station], - aft_web_overwrap_thickness, - ] - else: - layer_thicknesses = [ - cs_params["web_fore_adhesive_thickness"][i_station], - fore_web_overwrap_thickness, - ] + layer_thicknesses = [ + cs_params["web_adhesive_thickness"][i_station],overwrap_thickness] for it, thickness in enumerate(layer_thicknesses): cubit.cmd( - f"create curve offset curve {bottom_curve} distance {offSetSign*thickness} extended" - ) + f"create curve offset curve {bottom_curve} distance {offSetSign*thickness} extended") + top_curve = get_last_id("curve") if it == 0: @@ -1564,12 +1599,9 @@ def make_cs_web_layer_areas( ply_angle = 0 else: - if i_curve < n_base_curves_web / 2: - material_name = aft_web_stack.plygroups[0].materialid - ply_angle = aft_web_stack.plygroups[0].angle - else: - material_name = fore_web_stack.plygroups[0].materialid - ply_angle = fore_web_stack.plygroups[0].angle + material_name = web_stack.plygroups[0].materialid + ply_angle = web_stack.plygroups[0].angle + part_name_id, materials_used = make_cross_section_surface(lp_hp_side, surface_dict, @@ -1582,8 +1614,7 @@ def make_cs_web_layer_areas( part_name_id, n_modeled_layers + it, materials_used, - -1 - ) + -1) bottom_curve = top_curve @@ -1592,32 +1623,21 @@ def make_cs_web_layer_areas( ### Create vertical web regions lp_hp_side='' - # remove curves that are not going to be part of the vertical web - for i_curveList, curveList in enumerate(web_interface_curves): - curveList.pop(3) - curveList.pop(3) n_base_curves_web = len(web_interface_curves[0]) - for i_curve in range(n_base_curves_web): + for i_curve in range(n_base_curves_web-1): #Minus 1 since the overwrap is not part of vertical web vHP, _ = selCurveVerts(web_interface_curves[0][i_curve]) vLP, _ = selCurveVerts(web_interface_curves[1][i_curve]) top_curve,_ = print_sine_curve_between_two_verts( - vHP, vLP, cs_params["max_web_imperfection_distance"][i_station], "x" - ) + vHP, vLP, cs_params["max_web_imperfection_distance"][i_station], "x") _, vHP = selCurveVerts(web_interface_curves[0][i_curve]) _, vLP = selCurveVerts(web_interface_curves[1][i_curve]) bottom_curve,_ = print_sine_curve_between_two_verts( - vHP, vLP, cs_params["max_web_imperfection_distance"][i_station], "x" - ) + vHP, vLP, cs_params["max_web_imperfection_distance"][i_station], "x") + + material_name = web_stack.plygroups[i_curve].materialid + ply_angle = web_stack.plygroups[i_curve].angle - if i_curve < n_base_curves_web / 2: - material_name = aft_web_stack.plygroups[i_curve].materialid - ply_angle = aft_web_stack.plygroups[i_curve].angle - else: - material_name = fore_web_stack.plygroups[ - i_curve - int(n_base_curves_web / 2) - ].materialid - ply_angle = fore_web_stack.plygroups[i_curve - int(n_base_curves_web / 2)].angle part_name_id, materials_used = make_cross_section_surface(lp_hp_side, surface_dict, i_station, @@ -1629,13 +1649,14 @@ def make_cs_web_layer_areas( part_name_id, n_modeled_layers + it + 2 + i_curve, materials_used, - -1 - ) + -1) + surf_id = get_last_id("surface") surf_name = cubit.get_entity_name("surface", surf_id).split('_') surf_name.insert(-1,'vertical') surf_name = '_'.join(surf_name) cubit.cmd(f'surface {surf_id} rename "{surf_name}"') + return part_name_id, (vHP, vLP) def make_a_precomp_cross_section(wt_name, @@ -1643,7 +1664,7 @@ def make_a_precomp_cross_section(wt_name, i_station, i_station_geometry, blade, - hasWebs, + has_webs, aft_web_stack, fore_web_stack, iLE, @@ -1912,7 +1933,7 @@ def make_a_precomp_cross_section(wt_name, file.write(f'{i_ply+1} {1} {ply.thickness*ply.nPlies/1000} {ply.angle} {mat_id} ({ply.materialid})\n') - if hasWebs: + if has_webs: file.write(f'\n\n**********************************************************************\n') file.write(f'Laminae schedule for webs (input required only if webs exist at this section):\n') n_webs=len(stackdb.swstacks) @@ -1972,7 +1993,7 @@ def make_a_precomp_cross_section(wt_name, file.write(f'0.00 {-1*x_move} {geometry.ichord[i_station_geometry]} {geometry.idegreestwist[i_station_geometry]} shape_{i_station}.inp layup_{i_station}.inp\n') file.write(f'1.00 {-1*x_move} {geometry.ichord[i_station_geometry]} {geometry.idegreestwist[i_station_geometry]} shape_{i_station}.inp layup_{i_station}.inp\n') - if hasWebs: + if has_webs: n_webs =2 else: n_webs = 0 @@ -1983,7 +2004,7 @@ def make_a_precomp_cross_section(wt_name, file.write(f'Web_num Inb_end_ch_loc Oub_end_ch_loc (fraction of chord length)\n') - if hasWebs: + if has_webs: i_loc = 2 web_loc= mean([lp_sector_boundaries[i_loc],hp_sector_boundaries[i_loc]]) file.write(f'1.0000 {web_loc} {web_loc}\n') @@ -2000,9 +2021,8 @@ def make_a_cross_section(wt_name, i_station, i_station_geometry, blade, - hasWebs, - aft_web_stack, - fore_web_stack, + has_webs, + web_stacks, iLE, cs_params, geometry_scaling, @@ -2044,9 +2064,7 @@ def make_a_cross_section(wt_name, flatback_vTop, _ = selCurveVerts(lp_key_curve) - flatbackCurve = cubit.create_curve( - cubit.vertex(flatback_vBot), cubit.vertex(flatback_vTop) - ) + flatbackCurve = cubit.create_curve(cubit.vertex(flatback_vBot), cubit.vertex(flatback_vTop)) flatback_curve_id = flatbackCurve.id() #### Extend flatback ### @@ -2167,20 +2185,19 @@ def make_a_cross_section(wt_name, flatback_curve_id = get_last_id("curve") - n_stacks = len(stackdb.stacks) + n_areas = len(stackdb.stacks) le_hp_stack_thickness = ( - sum(stackdb.stacks[int(n_stacks / 2.0) - 1, i_station].layer_thicknesses()) / 1000 - ) + sum(stackdb.stacks[int(n_areas / 2.0) - 1, i_station].layer_thicknesses()) / 1000) le_lp_stack_thickness = ( - sum(stackdb.stacks[int(n_stacks / 2.0), i_station].layer_thicknesses()) / 1000 - ) + sum(stackdb.stacks[int(n_areas / 2.0), i_station].layer_thicknesses()) / 1000) # Define variables with HP side in index 0, LP side in index 1 - + n_webs = len(web_stacks) lp_hp_dict = {} - lp_hp_dict["spar_cap_base_curves"] = [[], []] - lp_hp_dict["web_interface_curves"] = [[], []] + lp_hp_dict["overwrap_base_curves"] = [[], []] + lp_hp_dict["sparcap_base_curve"] = [[], []] + lp_hp_dict["web_interface_curves"] = {'LP':[[] for x in range(n_webs)],'HP':[[] for x in range(n_webs)]} lp_hp_dict["base_curve_ids"] = [[], []] lp_hp_dict["round_te_adhesive_curve_list"] = [[], []] lp_hp_dict["flat_te_adhesive_curve_list"] = [[], []] @@ -2196,21 +2213,23 @@ def make_a_cross_section(wt_name, cs_normal, ) - key_curves = split_curve_at_coordinte_points( - keypoints.key_points[1:5, :, i_station_geometry], hp_key_curve - ) + key_points = keypoints.key_points[1:int(n_areas/2), :, i_station_geometry] + key_curves = split_curve_at_coordinte_points(key_points, hp_key_curve) web_adhesive_width = cs_params["web_adhesive_width"][i_station] ( lp_hp_dict["base_curve_ids"][0], - lp_hp_dict["spar_cap_base_curves"][0], - ) = split_key_curves(key_curves, aft_web_stack, fore_web_stack, web_adhesive_width) + lp_hp_dict["overwrap_base_curves"][0], + lp_hp_dict["sparcap_base_curve"][0] + ) = split_key_curves(key_curves, web_stacks, web_adhesive_width) - temp = np.flip(keypoints.key_points[:, :, i_station_geometry], 0) - key_curves = split_curve_at_coordinte_points(temp[1:5, :], lp_key_curve) + # temp = np.flip(keypoints.key_points[:, :, i_station_geometry], 0) + key_points = keypoints.key_points[-2:int(n_areas/2)-2:-1, :, i_station_geometry] + key_curves = split_curve_at_coordinte_points(key_points, lp_key_curve) ( lp_hp_dict["base_curve_ids"][1], - lp_hp_dict["spar_cap_base_curves"][1], - ) = split_key_curves(key_curves, aft_web_stack, fore_web_stack, web_adhesive_width) + lp_hp_dict["overwrap_base_curves"][1], + lp_hp_dict["sparcap_base_curve"][1] + ) = split_key_curves(key_curves, web_stacks, web_adhesive_width) # Make sure that the adhesive width is the same on HP and LP sides. Also @@ -2320,8 +2339,9 @@ def make_a_cross_section(wt_name, part_name_id, lp_hp_dict,stack_ct = make_cs_perimeter_layer_areas(wt_name, surface_dict, i_station, - stackdb.stacks[1:6, i_station], + stackdb.stacks[1:int(n_areas/2), i_station], cs_params, + has_webs, thickness_scaling, lp_hp_side, last_round_station, @@ -2336,13 +2356,12 @@ def make_a_cross_section(wt_name, stack_ct+=1 lp_hp_side = "LP" - temp = stackdb.stacks[:, i_station] - temp = np.flip(temp) part_name_id, lp_hp_dict,stack_ct = make_cs_perimeter_layer_areas(wt_name, surface_dict, i_station, - temp[1:6], + stackdb.stacks[-2:int(n_areas/2)-1:-1, i_station], cs_params, + has_webs, thickness_scaling, lp_hp_side, last_round_station, @@ -2425,21 +2444,23 @@ def make_a_cross_section(wt_name, cubit.cmd(f'curve {curve_id} rename "{curve_name}"') birds_mouth_verts = [] - if hasWebs: - part_name_id = 0 # Reset since outer areoshell is complete (LE adhesive is accouted for as aeroshell) - - part_name_id, birds_mouth_verts = make_cs_web_layer_areas( - surface_dict, - i_station, - aft_web_stack, - fore_web_stack, - lp_hp_dict["web_interface_curves"], - cs_params, - part_name_id, - cs_normal, - n_modeled_layers, - materials_used, - ) + for i_web in range(n_webs): + + if has_webs[i_web]: + web_stack = web_stacks[i_web] + web_interface_curves=[lp_hp_dict["web_interface_curves"]['HP'][i_web],lp_hp_dict["web_interface_curves"]['LP'][i_web]] + part_name_id = 0 # Reset since outer areoshell is complete (LE adhesive is accouted for as aeroshell) + part_name = f'web{i_web}' + part_name_id, birds_mouth_verts = make_cs_web_layer_areas(part_name, + surface_dict, + i_station, + web_stack, + web_interface_curves, + cs_params, + part_name_id, + cs_normal, + n_modeled_layers, + materials_used) parse_string = f'with name "*station{str(i_station).zfill(3)}*"' cs_surfaces = parse_cubit_list("surface", parse_string) From cb415f480ba29c2a4c7a66f588d011e0ebfead53 Mon Sep 17 00:00:00 2001 From: Camarena Date: Tue, 7 Jan 2025 10:59:34 -0700 Subject: [PATCH 08/11] bug fix to iea22 failed flatback TE --- src/pynumad/analysis/cubit/make_blade.py | 3 +- .../analysis/cubit/make_cross_sections.py | 32 ++++++++++++------- 2 files changed, 22 insertions(+), 13 deletions(-) diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index ec270c7..b0c54c9 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -694,7 +694,8 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati else: break last_flat_station+=1 - + cubit.cmd(f'delete curve all') + cubit.cmd(f'delete vertex all') ######################## diff --git a/src/pynumad/analysis/cubit/make_cross_sections.py b/src/pynumad/analysis/cubit/make_cross_sections.py index 89e9e34..e303c7c 100644 --- a/src/pynumad/analysis/cubit/make_cross_sections.py +++ b/src/pynumad/analysis/cubit/make_cross_sections.py @@ -1085,19 +1085,12 @@ def make_cs_perimeter_layer_areas(wt_name, last_layer_offset, curve_start_or_end, lp_hp_dict["flatback_curve_id"] ) - if ( - is_flatback - and abs( - min(current_stack.layer_thicknesses()) - - max(current_stack.layer_thicknesses()) - ) - > 0.0001 - ): + if (is_flatback and abs(min(current_stack.layer_thicknesses()) + - max(current_stack.layer_thicknesses()))> 0.0001): layer_offset_dist = current_stack_layer_thicknesses[-1] curve_start_or_end = "end" last_layer_offset = get_adjustment_curve( - curve_ids, layer_offset_dist, curve_start_or_end, end_layer_taper_curve - ) + curve_ids, layer_offset_dist, curve_start_or_end, end_layer_taper_curve) # cubit.cmd(f'split curve {first_layer_offset} at vertex {lp_hp_dict["perimeter_verts_for_te_adhesive"][lp_hp_side]}') @@ -1125,14 +1118,29 @@ def make_cs_perimeter_layer_areas(wt_name, for vertex_id in lp_hp_dict["perimeter_verts_for_te_adhesive"][lp_hp_side]: temp_list_left=[] temp_list_right=[] + break_flag = False for curve_id in current_stack_right_curves: + cubit.cmd(f'curve {curve_id} copy') n_start = get_last_id("curve") - cubit.cmd(f'split curve {curve_id} at vertex {vertex_id}') + cubit.cmd(f'split curve {n_start} at vertex {vertex_id}') n_end = get_last_id("curve") if n_end - n_start< 2: - print('') + break_flag = True + break + temp_list_left.append(get_last_id("curve") - 1) temp_list_right.append(get_last_id("curve")) + + if break_flag: + temp_list_left=[] + temp_list_right=[] + for curve_id in current_stack_right_curves: + cubit.cmd(f'curve {curve_id} copy') + cubit.cmd(f'split curve {get_last_id("curve")} distance {cs_params["te_adhesive_width"][i_station]} from start') + temp_list_left.append(get_last_id("curve") - 1) + temp_list_right.append(get_last_id("curve")) + + current_stack_right_curves=temp_list_right current_stack_left_curves_splited.append(temp_list_left) #current_stack_left_curves_splited.append(temp_list_right) From 4a8c7dd6955f77d945ca7d65c39672010b2aa777 Mon Sep 17 00:00:00 2001 From: Camarena Date: Fri, 10 Jan 2025 09:39:45 -0700 Subject: [PATCH 09/11] bug fix to beam utils when trapizodal integration is used. --- src/pynumad/analysis/beam_utils.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/pynumad/analysis/beam_utils.py b/src/pynumad/analysis/beam_utils.py index 95029cf..c71d3fa 100644 --- a/src/pynumad/analysis/beam_utils.py +++ b/src/pynumad/analysis/beam_utils.py @@ -358,8 +358,11 @@ def get_bd_source_radial_stations(bd_sum_file_name): with open(bd_sum_file_name) as blade_yaml: data = yaml.load(blade_yaml,Loader=yaml.FullLoader) - node_location=list(np.array(data['Init_Nodes_E1'])[:,2]/np.array(data['Init_Nodes_E1'])[-1,2]) - return node_location + if 'trap' in data['Quadrature_method'].lower(): + return None #The the source locations are target locations + else: + node_location=list(np.array(data['Init_Nodes_E1'])[:,2]/np.array(data['Init_Nodes_E1'])[-1,2]) + return node_location # def get_bd_source_radial_stations(qp_location,blade_length,nqp): #Adds root and top to qp data From 30e19277c994138cea91e6aba7178db418992653 Mon Sep 17 00:00:00 2001 From: Camarena Date: Mon, 13 Jan 2025 12:24:17 -0700 Subject: [PATCH 10/11] adding stack for tip station --- src/pynumad/objects/stackdb.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/pynumad/objects/stackdb.py b/src/pynumad/objects/stackdb.py index d11a7d0..ad840fe 100644 --- a/src/pynumad/objects/stackdb.py +++ b/src/pynumad/objects/stackdb.py @@ -34,7 +34,7 @@ def __eq__(self, other): def generate(self, keypoints: KeyPoints, bom: BillOfMaterials): # build the material stack for each area n_segments = keypoints.key_areas.shape[0] - n_stations = keypoints.key_areas.shape[1] + n_stations = keypoints.key_areas.shape[1] +1 n_webs = len(keypoints.web_points) segment_labels = [ "HP_TE_FLAT", @@ -84,14 +84,14 @@ def generate(self, keypoints: KeyPoints, bom: BillOfMaterials): cur_ply.thickness = bom["hp"][k].thickness cur_ply.angle = 0 # TODO, set to 0 for now, bom['lp'](k, ); cur_ply.nPlies = 1 # default to 1, modified in addply() if necessary - # ... and add the ply to every area that is part of the region ind = bom.indices["hp"][k] for k_seg in range(ind[2], ind[3]): for k_stat in range(ind[0], ind[1]): # deepcopy is important to keep make ply object unique in each stack self.stacks[k_seg, k_stat].addply(deepcopy(cur_ply)) - + if ind[1] == keypoints.key_areas.shape[1]: + self.stacks[k_seg, ind[1]].addply(deepcopy(cur_ply)) for k in range(len(bom["lp"])): # for each row in the BOM, get the ply definition ... cur_ply = Ply() @@ -106,7 +106,8 @@ def generate(self, keypoints: KeyPoints, bom: BillOfMaterials): for k_seg in range(ind[2], ind[3]): for k_stat in range(ind[0], ind[1]): self.stacks[k_seg, k_stat].addply(deepcopy(cur_ply)) - + if ind[1] == keypoints.key_areas.shape[1]: + self.stacks[k_seg, ind[1]].addply(deepcopy(cur_ply)) self.swstacks = np.empty(shape=(n_webs, n_stations), dtype=object) for k_web in range(n_webs): for k_stat in range(n_stations): From bd52487bdc8ff99ccd08b53d3852c86da4d13cf2 Mon Sep 17 00:00:00 2001 From: Camarena Date: Tue, 14 Jan 2025 13:21:42 -0700 Subject: [PATCH 11/11] Deleting i_station_geometry and last station hack. Bug fix on the deletion of last webs. --- src/pynumad/analysis/cubit/make_blade.py | 44 +++++++------------ .../analysis/cubit/make_cross_sections.py | 33 ++++++++------ 2 files changed, 37 insertions(+), 40 deletions(-) diff --git a/src/pynumad/analysis/cubit/make_blade.py b/src/pynumad/analysis/cubit/make_blade.py index b0c54c9..8df787f 100644 --- a/src/pynumad/analysis/cubit/make_blade.py +++ b/src/pynumad/analysis/cubit/make_blade.py @@ -640,8 +640,8 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati #### Step one create outer mold line excess_lengths=[] te_angles=[] - for i_station_geometry in range(len(blade.geometry.ispan)-1): #-1 b/c fewer stacks than stations - xyz = get_blade_geometry_for_station(blade, i_station_geometry) * geometry_scaling + for i_station in range(len(blade.geometry.ispan)): + xyz = get_blade_geometry_for_station(blade, i_station) * geometry_scaling npts=5 # Start indexing from 1 (not 0) to ignore first point: because first point is not on the LP or HP surface but rather is the midpoint at the TE @@ -660,9 +660,9 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati flatback_length=np.linalg.norm(second_point - first_point) - athickness=cs_params["te_adhesive_thickness"][i_station_geometry] - stack_thicknesses_hp=sum(stackdb.stacks[1, i_station_geometry].layer_thicknesses())/1000 - stack_thicknesses_lp=sum(stackdb.stacks[-2, i_station_geometry].layer_thicknesses())/1000 + athickness=cs_params["te_adhesive_thickness"][i_station] + stack_thicknesses_hp=sum(stackdb.stacks[1, i_station].layer_thicknesses())/1000 + stack_thicknesses_lp=sum(stackdb.stacks[-2, i_station].layer_thicknesses())/1000 excess_lengths.append(flatback_length-(stack_thicknesses_lp+stack_thicknesses_hp+athickness)) @@ -712,22 +712,9 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati for i_station in stationList: if model2Dor3D.lower() == "2d": - cubit.cmd( - "reset " - ) # This is needed to restart node numbering for VABS. VABS neeeds every element and node starting from 1 to nelem/nnode should be present + cubit.cmd("reset ") # This is needed to restart node numbering for VABS. + # VABS neeeds every element and node starting from 1 to nelem/nnode should be present write_spline_from_coordinate_points(cubit, ref_line_coords) - i_station_geometry = i_station - if i_station == len(geometry.ispan) - 1: # Only do this for the last station - blade.add_interpolated_station(geometry.ispan[-1] * 0.999) - stackdb.edit_stacks_for_solid_mesh() - expandTEthicknesses.append(expandTEthicknesses[-1]) - blade.expand_blade_geometry_te(expandTEthicknesses) - - # adjustLastStackAfterNewTipStation(i_station) - - i_station_geometry = i_station + 1 - - #is_flatback=is_station_flatback[i_station_geometry] web_stacks = [] for i_web in range(n_webs-1,-1,-1): #stackdb.swstacks is arranged from LE to TE. Need TE to LE order @@ -749,9 +736,9 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati cs_normal = get_cs_normal_vector( np.array( [ - keypoints.key_points[2, :, i_station_geometry], - keypoints.key_points[3, :, i_station_geometry], - keypoints.key_points[7, :, i_station_geometry], + keypoints.key_points[2, :, i_station], + keypoints.key_points[3, :, i_station], + keypoints.key_points[7, :, i_station], ] ) ) @@ -770,7 +757,6 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati birds_mouth_verts = make_a_cross_section(wt_name, surface_dict, i_station, - i_station_geometry, blade, make_webs, web_stacks, @@ -787,7 +773,6 @@ def cubit_make_cross_sections(blade,wt_name,settings,cs_params,model2Dor3D,stati make_a_cross_section(wt_name, surface_dict, i_station, - i_station_geometry, blade, make_webs, web_stacks, @@ -1086,7 +1071,6 @@ def cubit_make_solid_blade( shell_vol_list,spanwise_splines = make_all_volumes_for_a_part(surface_dict, ordered_list, i_station_end,spanwise_splines) else: shell_vol_list=[] - n_webs=len(blade.stackdb.swstacks) for i_web in range(n_webs): part_name = f'web{i_web}' @@ -1133,10 +1117,16 @@ def cubit_make_solid_blade( for i_station in stationList[0:-1]: temp = np.flip(has_webs[i_station]) + for i_web in range(n_webs): if not temp[i_web]: cubit.cmd(f"delete volume with name 'web{i_web}*tation{str(i_station).zfill(3)}*'") - + #If the last cross section does not have webs, delete last volumes + temp = np.flip(has_webs[-1]) + i_station = stationList[-2] + for i_web in range(n_webs): + if not temp[i_web]: + cubit.cmd(f"delete volume with name 'web{i_web}*tation{str(i_station).zfill(3)}*'") cubit.cmd(f"merge volume all") diff --git a/src/pynumad/analysis/cubit/make_cross_sections.py b/src/pynumad/analysis/cubit/make_cross_sections.py index e303c7c..02ee64f 100644 --- a/src/pynumad/analysis/cubit/make_cross_sections.py +++ b/src/pynumad/analysis/cubit/make_cross_sections.py @@ -1136,9 +1136,17 @@ def make_cs_perimeter_layer_areas(wt_name, temp_list_right=[] for curve_id in current_stack_right_curves: cubit.cmd(f'curve {curve_id} copy') - cubit.cmd(f'split curve {get_last_id("curve")} distance {cs_params["te_adhesive_width"][i_station]} from start') - temp_list_left.append(get_last_id("curve") - 1) - temp_list_right.append(get_last_id("curve")) + curve_id = get_last_id("curve") + if cs_params["te_adhesive_width"][i_station] < cubit.curve(get_last_id("curve")).length(): + cubit.cmd(f'split curve {curve_id} distance {cs_params["te_adhesive_width"][i_station]} from start') + temp_list_left.append(get_last_id("curve") - 1) + temp_list_right.append(get_last_id("curve")) + else: + #TE adhesive width is too wide at station {i_station}. Making sure adhesive width fits in width of the reinf. + cubit.cmd(f'split curve {curve_id} fraction 0.1 from end') + temp_list_left.append(get_last_id("curve") - 1) + temp_list_right.append(get_last_id("curve")) + # raise ValueError(f'TE adhesive width is too wide at station {i_station}. Widen TE reif. or make adhesive width more narrow. ') current_stack_right_curves=temp_list_right @@ -2027,7 +2035,6 @@ def make_a_precomp_cross_section(wt_name, def make_a_cross_section(wt_name, surface_dict, i_station, - i_station_geometry, blade, has_webs, web_stacks, @@ -2043,7 +2050,7 @@ def make_a_cross_section(wt_name, geometry = blade.geometry stackdb = blade.stackdb keypoints = blade.keypoints - + if i_station > last_round_station: is_flatback=True else: @@ -2055,7 +2062,7 @@ def make_a_cross_section(wt_name, part_name_id = 0 #### Step one create outer mold line - xyz = get_blade_geometry_for_station(blade, i_station_geometry) * geometry_scaling + xyz = get_blade_geometry_for_station(blade, i_station) * geometry_scaling # Start indexing from 1 (not 0) to ignore first point: because first point is not on the LP or HP surface but rather is the midpoint at the TE splinePoints = xyz[1:iLE, :] @@ -2078,7 +2085,7 @@ def make_a_cross_section(wt_name, #### Extend flatback ### curve_start_or_end = "start" extension_length = ( - 100 * geometry.ichord[i_station_geometry] * cubit.curve(flatback_curve_id).length() + 100 * geometry.ichord[i_station] * cubit.curve(flatback_curve_id).length() ) flatback_curve_id = extend_curve_at_vertex_to_length( flatback_curve_id, extension_length, curve_start_or_end @@ -2094,7 +2101,7 @@ def make_a_cross_section(wt_name, # Crate camber line offset_distance = 0 npts = 100 - spacing = geometry.ichord[i_station_geometry] * 0.5 / npts + spacing = geometry.ichord[i_station] * 0.5 / npts cubit.cmd(f"curve {flatback_curve_id} copy") flatback_offset_curve_id = get_last_id("curve") @@ -2152,7 +2159,7 @@ def make_a_cross_section(wt_name, # # cubit.cmd(f'save as "Debug.cub" overwrite') # # foo else: - xyz = get_mid_line(blade, iLE, i_station_geometry, geometry_scaling) + xyz = get_mid_line(blade, iLE, i_station, geometry_scaling) npts, _ = xyz.shape npts = round( npts * 0.75 @@ -2174,7 +2181,7 @@ def make_a_cross_section(wt_name, #### Extend flatback ### curve_start_or_end = "start" extension_length = ( - 100 * geometry.ichord[i_station_geometry] * cubit.curve(flatback_curve_id).length() + 100 * geometry.ichord[i_station] * cubit.curve(flatback_curve_id).length() ) cubit.cmd(f"curve {flatback_curve_id} copy") curve_id = extend_curve_at_vertex_to_length( @@ -2221,7 +2228,7 @@ def make_a_cross_section(wt_name, cs_normal, ) - key_points = keypoints.key_points[1:int(n_areas/2), :, i_station_geometry] + key_points = keypoints.key_points[1:int(n_areas/2), :, i_station] key_curves = split_curve_at_coordinte_points(key_points, hp_key_curve) web_adhesive_width = cs_params["web_adhesive_width"][i_station] ( @@ -2230,8 +2237,8 @@ def make_a_cross_section(wt_name, lp_hp_dict["sparcap_base_curve"][0] ) = split_key_curves(key_curves, web_stacks, web_adhesive_width) - # temp = np.flip(keypoints.key_points[:, :, i_station_geometry], 0) - key_points = keypoints.key_points[-2:int(n_areas/2)-2:-1, :, i_station_geometry] + # temp = np.flip(keypoints.key_points[:, :, i_station], 0) + key_points = keypoints.key_points[-2:int(n_areas/2)-2:-1, :, i_station] key_curves = split_curve_at_coordinte_points(key_points, lp_key_curve) ( lp_hp_dict["base_curve_ids"][1],