Skip to content

Commit

Permalink
update for D150
Browse files Browse the repository at this point in the history
  • Loading branch information
@Romain-Gauthier
Romain-Gauthier committed Oct 23, 2024
1 parent d5298be commit 47715c9
Show file tree
Hide file tree
Showing 5 changed files with 95 additions and 104 deletions.
82 changes: 41 additions & 41 deletions ceasiompy/AeroFrame_new/aeroframe_run.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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']])
Expand All @@ -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 -----")
Expand All @@ -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

Expand Down
93 changes: 34 additions & 59 deletions ceasiompy/AeroFrame_new/func/aeroframe_config.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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')
Expand Down Expand Up @@ -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 = []

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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)
Expand All @@ -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,
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand All @@ -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
Expand All @@ -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.
Expand Down Expand Up @@ -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
Expand Down
4 changes: 3 additions & 1 deletion ceasiompy/AeroFrame_new/func/aeroframe_debbug.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion ceasiompy/AeroFrame_new/func/aeroframe_results.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
17 changes: 15 additions & 2 deletions ceasiompy/AeroFrame_new/func/aeroframe_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()


# =================================================================================================
Expand All @@ -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
Expand Down Expand Up @@ -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!")

0 comments on commit 47715c9

Please sign in to comment.