From 47715c9cf00b5f12c1acf431c384c153280e1f5d Mon Sep 17 00:00:00 2001 From: Romain-Gauthier Date: Wed, 23 Oct 2024 14:02:12 +0200 Subject: [PATCH] update for D150 --- ceasiompy/AeroFrame_new/aeroframe_run.py | 82 ++++++++-------- .../AeroFrame_new/func/aeroframe_config.py | 93 +++++++------------ .../AeroFrame_new/func/aeroframe_debbug.py | 4 +- .../AeroFrame_new/func/aeroframe_results.py | 3 +- .../AeroFrame_new/func/aeroframe_utils.py | 17 +++- 5 files changed, 95 insertions(+), 104 deletions(-) diff --git a/ceasiompy/AeroFrame_new/aeroframe_run.py b/ceasiompy/AeroFrame_new/aeroframe_run.py index 702c3ce3..afe44fbe 100644 --- a/ceasiompy/AeroFrame_new/aeroframe_run.py +++ b/ceasiompy/AeroFrame_new/aeroframe_run.py @@ -137,9 +137,37 @@ def aeroelastic_loop(cpacs_path, CASE_PATH, q, xyz, fxyz): wg_center_z_list[0] + wing_origin[2]]) fxyz_root = np.zeros(3) - xyz_tip = np.array([wg_center_x_list[-1] + wing_origin[0] + wg_chord_list[-1] / 2, - wg_center_y_list[-1] + wing_origin[1], - wg_center_z_list[-1] + wing_origin[2]]) + # xyz_tip = np.array([wg_center_x_list[-1] + wing_origin[0] + wg_chord_list[-1] / 2, + # wg_center_y_list[-1] + wing_origin[1], + # wg_center_z_list[-1] + wing_origin[2]]) + + Xle_list = [] + Yle_list = [] + Zle_list = [] + surface_count = 0 + with open(AVL_UNDEFORMED_PATH, "r") as f: + lines = f.readlines() + for i, line in enumerate(lines): + if "SURFACE" in line: + surface_count += 1 + if surface_count > 1: + break + + if "Xle" in line: + next_line = lines[i + 1].strip() + parts = next_line.split() + Xle_list.append(float(parts[0]) + wing_origin[0]) + Yle_list.append(float(parts[1]) + wing_origin[1]) + Zle_list.append(float(parts[2]) + wing_origin[2]) + + Xle_array = np.array(Xle_list) + Yle_array = np.array(Yle_list) + Zle_array = np.array(Zle_list) + + xyz_tip = np.array([Xle_array[-1] + wg_chord_list[-1] / 2, + Yle_array[-1], + Zle_array[-1]]) + # Get cross-section properties from CPACS file (if constants) area_const, Ix_const, Iy_const = get_section_properties(cpacs_path) @@ -245,8 +273,8 @@ def aeroelastic_loop(cpacs_path, CASE_PATH, q, xyz, fxyz): log.info(f"Iteration {n_iter} done!") log.info( - f"Wing tip deflection: {deflection:.3e} m ({percentage:.2%} of the semi-span length).") - log.info(f"Residual: {res[-1]:.3e}") + f"Wing tip deflection: {deflection:.3e} m ({percentage:.2%} of the semi-span length).") + log.info(f"Residual: {res[-1]:.3e}") # Run AVL with the new deformed geometry write_deformed_geometry(AVL_UNDEFORMED_PATH, AVL_DEFORMED_PATH, centerline_df, deformed_df) @@ -279,33 +307,6 @@ def compute_aero_work(row): wing_df['aero_work'] = wing_df.apply(compute_aero_work, axis=1) total_aero_work = wing_df['aero_work'].sum() - # E = 325e9 # Young's modulus in Pascals (example value) - # G = 125e9 # Shear modulus in Pascals (example value) - # A = 1.435e-3 # Cross-sectional area in square meters (example value) - # I_x = 2.01e-9 # Second moment of area about y-axis in meters^4 (example value) - # I_z = 1.465e-5 # Second moment of area about z-axis in meters^4 (example value) - # J = I_x + I_z # Polar moment of inertia in meters^4 (example value) - - # # Compute axial strain energy - # centerline_df['axial_strain_energy'] = 0.5 * (centerline_df['Fy']**2) / (E * A) - - # # Compute bending strain energy - # centerline_df['bending_strain_energy'] = 0.5 * \ - # ((centerline_df['Mx']**2 / (E * I_x)) + (centerline_df['Mz']**2 / (E * I_z))) - - # # Compute torsional strain energy - # centerline_df['torsional_strain_energy'] = 0.5 * (centerline_df['My']**2) / (G * J) - - # # Sum the strain energies to get the total strain energy for each node - # centerline_df['total_strain_energy'] = ( - # centerline_df['axial_strain_energy'] - # + centerline_df['bending_strain_energy'] - # + centerline_df['torsional_strain_energy'] - # ) - - # # Sum the total strain energy over all nodes to get the total strain energy of the beam - # total_strain_energy = centerline_df['total_strain_energy'].sum() - def compute_structural_work(row): force = np.array([row['Fx'], row['Fy'], row['Fz']]) moment = np.array([row['Mx'], row['My'], row['Mz']]) @@ -317,10 +318,10 @@ def compute_structural_work(row): centerline_df['structural_work'] = centerline_df.apply(compute_structural_work, axis=1) total_structural_work = centerline_df['structural_work'].sum() - log.info(f"Total aerodynamic work: {total_aero_work:.3e} J.") - log.info(f"Total structural work: {total_structural_work:.3e} J.") + log.info(f"Total aerodynamic work: {total_aero_work:.3e} J.") + log.info(f"Total structural work: {total_structural_work:.3e} J.") log.info( - f"Work variation: {((total_aero_work-total_structural_work)/total_aero_work):.2%}.") + f"Work variation: {((total_aero_work-total_structural_work)/total_aero_work):.2%}.") log.info("") log.info("----- Final results -----") @@ -333,13 +334,12 @@ def compute_structural_work(row): percentage = deflection / semi_span log.info( - f"Wing tip deflection: {deflection:.3e} m ({percentage:.2%} of the semi-span length).") - log.info(f"Wing tip twist: {tip_twist:.3e} degrees.") - log.info(f"Total aerodynamic work: {total_aero_work:.3e} J.") - log.info(f"Total structural work: {total_structural_work:.3e} J.") + f"Wing tip deflection: {deflection:.3e} m ({percentage:.2%} of the semi-span length).") + log.info(f"Wing tip twist: {tip_twist:.3e} degrees.") + log.info(f"Total aerodynamic work: {total_aero_work:.3e} J.") + log.info(f"Total structural work: {total_structural_work:.3e} J.") log.info( - f"Work variation: {((total_aero_work-total_structural_work)/total_aero_work):.2%}.") - # log.info(f"Total Strain Energy of the Beam: {total_strain_energy:.3e} J.") + f"Work variation: {((total_aero_work-total_structural_work)/total_aero_work):.2%}.") return delta_tip, res diff --git a/ceasiompy/AeroFrame_new/func/aeroframe_config.py b/ceasiompy/AeroFrame_new/func/aeroframe_config.py index 1f6d5965..1eb9e126 100644 --- a/ceasiompy/AeroFrame_new/func/aeroframe_config.py +++ b/ceasiompy/AeroFrame_new/func/aeroframe_config.py @@ -340,8 +340,14 @@ def create_framat_model(young_modulus, shear_modulus, material_density, beam.add('point_load', {'at': node_uid, 'load': load}) # ===== BOUNDARY CONDITIONS ===== - idx_to_fix = centerline_df["y"].idxmin() - bc.add('fix', {'node': "wing1_node" + str(idx_to_fix + 1), 'fix': ['all']}) + #idx_to_fix = centerline_df["y"].idxmin() + #bc.add('fix', {'node': "wing1_node" + str(idx_to_fix + 1), 'fix': ['all']}) + + # /!\ FOR D150 ONLY + idx_to_fix = centerline_df[centerline_df["y"] < 1.8].index + for idx in idx_to_fix: + bc.add('fix', {'node': "wing1_node" + str(idx + 1), 'fix': ['all']}) + # ===== POST-PROCESSING ===== pp = model.set_feature('post_proc') @@ -480,10 +486,11 @@ def create_wing_centerline(wing_df, centerline_df, N_beam, wg_origin, xyz_tot, f wing_df = pd.concat([wing_df, tip_row], ignore_index=True) - Xle, Yle, Zle = interpolate_leading_edge(AVL_UNDEFORMED_PATH, - CASE_PATH, n_iter, - wg_origin, - y_queries=wing_df["y"].unique()) + _, _, _, Xle, Yle, Zle = interpolate_leading_edge(AVL_UNDEFORMED_PATH, + CASE_PATH, + wg_origin, + y_queries=wing_df["y"].unique(), + n_iter=n_iter) leading_edge = [] trailing_edge = [] @@ -603,7 +610,7 @@ def compute_distance_and_moment(row): for force in ['Fx', 'Fy', 'Fz', 'Mx', 'My', 'Mz']: centerline_df.at[centerline_index, force] += wing_df.at[i, force] - log.info(f"Total aerodynamic force: {centerline_df['Fz'].sum()} N.") + log.info(f"Total aerodynamic force: {centerline_df['Fz'].sum():.2f} N.") return wing_df, centerline_df, internal_load_df @@ -653,10 +660,6 @@ def compute_cross_section(cpacs_path): Ix_list = [] Iy_list = [] - # x_beam_range = [[22.92, 31], - # [29.7, 34.6], - # [46.46, 48.3]] - for i_wing in range(1): wing_xpath = WINGS_XPATH + "/wing[" + str(i_wing + 1) + "]" wing_transf = Transformation() @@ -826,16 +829,6 @@ def compute_cross_section(cpacs_path): + pos_z_list[i_sec] ) * wing_transf.scaling.z - # beam_vect_x = [] - # beam_vect_z = [] - - # for x, z in zip(prof_vect_x, prof_vect_z): - # log.info(f"{x} ; {z}") - # if x > x_beam_range[i_sec][0] - x_le[i_sec] - # and x < x_beam_range[i_sec][1] - x_le[i_sec]: - # beam_vect_x.append(x) - # beam_vect_z.append(z) - wg_twist_list.append(wg_sec_twist) area_list.append(PolyArea(prof_vect_x, prof_vect_z)) Ix, Iy = second_moments_of_area(x=prof_vect_x, y=prof_vect_z) @@ -848,9 +841,10 @@ def compute_cross_section(cpacs_path): wg_chord_list.append(wg_sec_chord) log.info(f"Section number: {i_sec}") - log.info(f"Area of the sections: {area_list}") - log.info(f"Ix of the sections: {Ix_list}") - log.info(f"Iy of the sections: {Iy_list}") + log.info(f"Area of the sections: [{', '.join([f'{area:.2e}' for area in area_list])}] m^2.") + log.info(f"Ix of the sections: [{', '.join([f'{Ix:.2e}' for Ix in Ix_list])}] m^4.") + log.info(f"Iy of the sections: [{', '.join([f'{Iy:.2e}' for Iy in Iy_list])}] m^4.") + return ( wg_origin, wg_twist_list, area_list, Ix_list, Iy_list, @@ -887,42 +881,23 @@ def write_deformed_geometry(UNDEFORMED_PATH, DEFORMED_PATH, centerline_df, defor "0.0\t0.0\t0.0\n\n", "#---------------\n"]) - # target_pairs = [(20.870, 0), (28.050, 10.074), (45.930, 35.134)] - - # airfoils = ["/users/disk8/cfse3/Stage_Romain/CEASIOMpy/WKDIR/ - # Workflow_353/Results/PyAVL/Airfoil1.dat", - # "/users/disk8/cfse3/Stage_Romain/CEASIOMpy/WKDIR/ - # Workflow_353/Results/PyAVL/Airfoil2.dat", - # "/users/disk8/cfse3/Stage_Romain/CEASIOMpy/WKDIR/ - # Workflow_353/Results/PyAVL/Airfoil3.dat"] - - # closest_indices = [ - # deformed_df.apply(lambda row: np.linalg.norm( - # [row['x_leading'] - target[0], - # row['y_leading'] - target[1]]), axis=1).idxmin() - # for target in target_pairs] - - # for i_sec, airfoil in zip(closest_indices, airfoils): - # x_new = deformed_df.iloc[i_sec]["x_leading"] - # y_new = deformed_df.iloc[i_sec]["y_leading"] - # z_new = deformed_df.iloc[i_sec]["z_leading"] - # chord = deformed_df.iloc[i_sec]["chord"] - # AoA = twist_profile(y_new) # deformed_df.iloc[i_node]["AoA"] - # file_deformed.writelines( - # ["SECTION\n", - # "#Xle Yle Zle Chord Ainc\n", - # f"{x_new:.3f} {y_new:.3f} {z_new:.3e} {chord:.3f} {AoA:.3e}\n\n", - # "AFILE\n", - # airfoil + "\n\n", - # "#---------------\n"]) - - step = 1 + step = 7 + root_sec_added = False for i_node in range(0, len(deformed_df), step): x_new = deformed_df.iloc[i_node]["x_leading"] y_new = deformed_df.iloc[i_node]["y_leading"] z_new = deformed_df.iloc[i_node]["z_leading"] chord = deformed_df.iloc[i_node]["chord"] - AoA = twist_profile(y_new) # deformed_df.iloc[i_node]["AoA"] + AoA = 2#twist_profile(y_new) + + if y_new > 1.856 and root_sec_added == False: + file_deformed.writelines( + ["SECTION\n", + "#Xle Yle Zle Chord Ainc\n", + f"{12.746} {1.856} {-1.136} {6.076} {2}\n\n", + "#---------------\n"]) + root_sec_added = True + file_deformed.writelines( ["SECTION\n", "#Xle Yle Zle Chord Ainc\n", @@ -933,7 +908,7 @@ def write_deformed_geometry(UNDEFORMED_PATH, DEFORMED_PATH, centerline_df, defor y_new = deformed_df.iloc[-1]["y_leading"] z_new = deformed_df.iloc[-1]["z_leading"] chord = deformed_df.iloc[-1]["chord"] - AoA = centerline_df.iloc[-1]["AoA_new"] + AoA = 2#centerline_df.iloc[-1]["AoA_new"] file_deformed.writelines( ["SECTION\n", "#Xle Yle Zle Chord Ainc\n", @@ -957,7 +932,7 @@ def write_deformed_command(UNDEFORMED_COMMAND, DEFORMED_COMMAND): deformed.write(line) -def interpolate_leading_edge(AVL_UNDEFORMED_PATH, CASE_PATH, n_iter, wg_origin, y_queries): +def interpolate_leading_edge(AVL_UNDEFORMED_PATH, CASE_PATH, wg_origin, y_queries, n_iter): """Function to get the coordinates of the leading-edge points. Function 'interpolate_leading_edge' computes the coordinates of the leading-edge points @@ -966,9 +941,9 @@ def interpolate_leading_edge(AVL_UNDEFORMED_PATH, CASE_PATH, n_iter, wg_origin, Args: AVL_UNDEFORMED_PATH (Path): path to the undeformed AVL geometry. CASE_PATH (Path): path to the flight case directory. - n_iter (int): number of the current iteration. wg_origin (list): list of the coordinates of the origin of the wing geometry [m]. y_queries (list): list of unique spanwise location of the VLM panels [m]. + n_iter (int): number of the current iteration. Returns: interpolated_Xle (numpy array): array of the x-coordinates of the interpolated LE points. @@ -1060,7 +1035,7 @@ def interpolate_leading_edge_points(Xle_array, Yle_array, Zle_array, y_queries): interpolated_Yle = interpolated_points[:, 1] interpolated_Zle = interpolated_points[:, 2] - return interpolated_Xle, interpolated_Yle, interpolated_Zle + return Xle_array, Yle_array, Zle_array, interpolated_Xle, interpolated_Yle, interpolated_Zle # ================================================================================================= # MAIN diff --git a/ceasiompy/AeroFrame_new/func/aeroframe_debbug.py b/ceasiompy/AeroFrame_new/func/aeroframe_debbug.py index 965595b1..4c41f82b 100644 --- a/ceasiompy/AeroFrame_new/func/aeroframe_debbug.py +++ b/ceasiompy/AeroFrame_new/func/aeroframe_debbug.py @@ -3,7 +3,9 @@ Developed by CFS ENGINEERING, 1015 Lausanne, Switzerland -Script to... +Script to plot the FEM mesh and VLM panels used in the aeroelastic computations, +and plot the shape of the deformed wing. It helps to see if the geometry was +accurately captured and if the meshes are fine. Python version: >=3.8 diff --git a/ceasiompy/AeroFrame_new/func/aeroframe_results.py b/ceasiompy/AeroFrame_new/func/aeroframe_results.py index e9d115a4..64ea7e2e 100644 --- a/ceasiompy/AeroFrame_new/func/aeroframe_results.py +++ b/ceasiompy/AeroFrame_new/func/aeroframe_results.py @@ -3,7 +3,8 @@ Developed by CFS ENGINEERING, 1015 Lausanne, Switzerland -Script to... +Script to compute the wing deformation, plot the displacements and rotations, +and plot the convergence. Python version: >=3.8 diff --git a/ceasiompy/AeroFrame_new/func/aeroframe_utils.py b/ceasiompy/AeroFrame_new/func/aeroframe_utils.py index 0d833857..b0e2ffff 100644 --- a/ceasiompy/AeroFrame_new/func/aeroframe_utils.py +++ b/ceasiompy/AeroFrame_new/func/aeroframe_utils.py @@ -22,6 +22,9 @@ import numpy as np import math from shapely.geometry import Polygon +from ceasiompy.utils.ceasiomlogger import get_logger + +log = get_logger() # ================================================================================================= @@ -42,8 +45,9 @@ def second_moments_of_area(x, y): Ix = 0 Iy = 0 x_centr, y_centr = compute_centroid(x, y) - x -= x_centr - y -= y_centr + x = [xi - x_centr for xi in x] + y = [yi - y_centr for yi in y] + n = len(x) for i in range(n): j = (i + 1) % n @@ -104,3 +108,12 @@ def rotate_3D_points(x, y, z, angle_x, angle_y, angle_z): rotation_matrix[2, 1] + z * rotation_matrix[2, 2] return x_rot, y_rot, z_rot + + +# ================================================================================================= +# MAIN +# ================================================================================================= + +if __name__ == "__main__": + + log.info("Nothing to execute!") \ No newline at end of file