From 7a88e6f804c6e7093d7e7285ae84d88ead8c3fcd Mon Sep 17 00:00:00 2001 From: mlauer154 Date: Thu, 14 Sep 2023 10:08:22 -0500 Subject: [PATCH] Updates to shock dissipation calculation --- pymead/analysis/calc_aero_data.py | 108 +++++++++++++++++++++------ pymead/analysis/compressible_flow.py | 8 ++ pymead/analysis/read_aero_data.py | 8 +- pymead/post/post_process.py | 3 +- pymead/tests/aero_tests/Edot_test.py | 10 +-- 5 files changed, 106 insertions(+), 31 deletions(-) create mode 100644 pymead/analysis/compressible_flow.py diff --git a/pymead/analysis/calc_aero_data.py b/pymead/analysis/calc_aero_data.py index f7da112a..36403594 100644 --- a/pymead/analysis/calc_aero_data.py +++ b/pymead/analysis/calc_aero_data.py @@ -13,6 +13,7 @@ from pymead.analysis.read_aero_data import read_aero_data_from_xfoil, read_Cp_from_file_xfoil, read_bl_data_from_mses, \ read_forces_from_mses, read_grid_stats_from_mses, read_field_from_mses, read_streamline_grid_from_mses, \ flow_var_idx, convert_blade_file_to_3d_array, read_actuator_disk_data_mses, read_Mach_from_mses_file +from pymead.analysis.compressible_flow import calculate_normal_shock_total_pressure_ratio from pymead.utils.file_conversion import convert_ps_to_svg from pymead.utils.geometry import check_airfoil_self_intersection, convert_numpy_array_to_shapely_LineString, \ calculate_area_triangle_heron @@ -252,7 +253,7 @@ def calculate_aero_data(airfoil_coord_dir: str, airfoil_name: str, coords: typin outputs_CPK = calculate_CPK_power_consumption(os.path.join(airfoil_coord_dir, airfoil_name)) except Exception as e: print(f"{e = }") - outputs_CPK = {"CPK": 1e9, "diss_shock": 1e9, "diss_surf": 1e9, "Edot": 1e9} + outputs_CPK = {"CPK": 1e9, "diss_shock": 1e9, "diss_surf": 1e9, "Edot": 1e9, "Cd": 1e9} aero_data = {**aero_data, **outputs_CPK} t2 = time.time() @@ -1117,7 +1118,7 @@ def surface_dissipation(last_BL_point: typing.Dict[str, np.ndarray]): def shock_dissipation(x_edge: typing.List[np.ndarray], y_edge: typing.List[np.ndarray], u: np.ndarray, - v: np.ndarray, dCpt: np.ndarray): + v: np.ndarray, M: np.ndarray, Cpt: np.ndarray, dCpt: np.ndarray, dCp: np.ndarray, gam: float): Shock = namedtuple("Shock", ("j", "i_start", "i_end")) shock_indices = np.argwhere(dCpt < -0.002) shock_indices = np.fliplr(shock_indices) @@ -1145,10 +1146,44 @@ def shock_dissipation(x_edge: typing.List[np.ndarray], y_edge: typing.List[np.nd edge_split_i.append([i1, i2]) i1 = i2 + 1 + # print(f"{edge_split_i = }") + def integral(s: Shock): - total_dCpt = np.sum(dCpt[s.i_start:s.i_end, s.j]) - _u = u[s.i_start - 1, s.j] - _v = v[s.i_start - 1, s.j] + + _dCp = dCp[s.i_start:s.i_end, s.j] + + if np.max(_dCp) < 0.1: # If the static pressure rise is too small, do not consider this to be a shock wave + return 0.0, 0.0 + + # print(f"{_dCp = }, {s.i_start = }, {s.i_end = }") + # + # M_temp = M[s.i_start:s.i_end, s.j] + # print(f"{M_temp = }") + + upstream_shock_cell_index = np.arange(s.i_start, s.i_end)[np.argwhere(_dCp > 0.0).flatten()[0]] + # print(f"{upstream_shock_cell_index = }, {s.j = }") + + # max_Mach_number_search_iter = 7 + # iter_count = 0 + # old_M = M[upstream_shock_cell_index, s.j] + # while True: + # print(f"{iter_count = }") + # if iter_count > max_Mach_number_search_iter: + # break + # upstream_shock_cell_index -= 1 + # current_M = M[upstream_shock_cell_index, s.j] + # print(f"{old_M = }, {current_M = }") + # if current_M < old_M: + # upstream_shock_cell_index += 1 + # break + # + # old_M = current_M + # iter_count += 1 + + _u = u[upstream_shock_cell_index, s.j] + _v = v[upstream_shock_cell_index, s.j] + _M = M[upstream_shock_cell_index, s.j] + _Cpt = Cpt[upstream_shock_cell_index, s.j] j_flow = 0 j_grid = 0 @@ -1158,21 +1193,43 @@ def integral(s: Shock): j_grid = s.j - split_idx[0] break - x1 = x_edge[j_flow][s.i_start, j_grid] - x2 = x_edge[j_flow][s.i_start, j_grid + 1] - y1 = y_edge[j_flow][s.i_start, j_grid] - y2 = y_edge[j_flow][s.i_start, j_grid + 1] + x1 = x_edge[j_flow][upstream_shock_cell_index + 1, j_grid] + x2 = x_edge[j_flow][upstream_shock_cell_index + 1, j_grid + 1] + y1 = y_edge[j_flow][upstream_shock_cell_index + 1, j_grid] + y2 = y_edge[j_flow][upstream_shock_cell_index + 1, j_grid + 1] n_hat_angle = np.arctan2(y2 - y1, x2 - x1) - np.pi / 2 n_hat = np.array([np.cos(n_hat_angle), np.sin(n_hat_angle)]) + + V_dir = np.array([_u, _v]) / np.hypot(_u, _v) + M_vec = _M * V_dir + M_normal_up = np.dot(M_vec, n_hat) + + if isinstance(M_normal_up, np.ndarray): + M_normal_up = M_normal_up[0] + + if M_normal_up <= 1.0: + return 0.0, 0.0 + else: + # print(f"{M_normal_up = }, {x1 = }, {y1 = }, {s.j = }, {upstream_shock_cell_index = }") + Cpt2_Cpt1 = calculate_normal_shock_total_pressure_ratio(M_normal_up, gam) + _dCpt = _Cpt * Cpt2_Cpt1 - _Cpt + dot_product = np.dot(np.array([_u, _v]), n_hat) dl = np.hypot(x2 - x1, y2 - y1) - shock_int = np.abs(total_dCpt) * dot_product * dl - return shock_int + shock_int = np.abs(_dCpt) * dot_product * dl + return shock_int, dl + + # print(f"{shock_indices = }") shock_diss = 0.0 + dl_total = 0.0 for shock_list in shocks.values(): for shock in shock_list: - shock_diss += integral(shock) + s_int, incremental_dl = integral(shock) + shock_diss += s_int + dl_total += incremental_dl + + # print(f"{dl_total = }") return shock_diss @@ -1470,9 +1527,11 @@ def calculate_CPK_power_consumption(analysis_subdir: str): mses_file = os.path.join(analysis_subdir, f"mses.{airfoil_system_name}") mses_log_file = os.path.join(analysis_subdir, "mses.log") coords = convert_blade_file_to_3d_array(blade_file) + mplot_log_file = os.path.join(analysis_subdir, "mplot.log") M_inf = read_Mach_from_mses_file(mses_file) gam = 1.4 # Hard-coded specific heat ratio + forces = read_forces_from_mses(mplot_log_file) field = read_field_from_mses(field_file, M_inf=M_inf, gam=gam) bl_data = read_bl_data_from_mses(bl_file) grid_stats = read_grid_stats_from_mses(grid_stats_file) @@ -1498,14 +1557,12 @@ def calculate_CPK_power_consumption(analysis_subdir: str): # Integrate over the propulsor outlet for the given flow section - inlet_integral = line_integral_Edot_inviscid(Cp[0, :], rho[0, :], u[0, :], v[0, :], V[0, :], x[0, :], - y[0, :], n_hat_dir="left") - Edot += inlet_integral - outlet_integral = line_integral_Edot_inviscid(Cp[-1, :], rho[-1, :], u[-1, :], v[-1, :], V[-1, :], x[-1, :], y[-1, :], n_hat_dir="right") Edot += outlet_integral + # print(f"{outlet_integral = }") + # The side cylinder contributions are negligible # if flow_section_idx == len(x_grid) - 1: # @@ -1553,26 +1610,29 @@ def calculate_CPK_power_consumption(analysis_subdir: str): last_BL_point[k] = np.array(v) break - print(f"{Edot = }") + # print(f"{Edot = }") Edot_visc = viscous_Edot(last_BL_point=last_BL_point) Edot += Edot_visc - print(f"{Edot_visc = }") + # print(f"{Edot_visc = }") surface_diss = surface_dissipation(last_BL_point=last_BL_point) - print(f"{surface_diss = }") + # print(f"{surface_diss = }") shock_diss = shock_dissipation(x_edge=x_grid, y_edge=y_grid, u=field[flow_var_idx["u"]][:, :], - v=field[flow_var_idx["v"]][:, :], dCpt=field[flow_var_idx["dCpt"]][:, :]) + v=field[flow_var_idx["v"]][:, :], M=field[flow_var_idx["M"]][:, :], + Cpt=field[flow_var_idx["Cpt"]][:, :], dCpt=field[flow_var_idx["dCpt"]][:, :], + dCp=field[flow_var_idx["dCp"]][:, :], gam=gam) - print(f"{shock_diss = }") + # print(f"{shock_diss = }") + # print(f"{forces['Cdw'] = }, {forces['Cd'] = }") - CPK = Edot + surface_diss + shock_diss + CPK = Edot + surface_diss + shock_diss - forces["Cd"] - print(f"{CPK = }") + # print(f"{CPK = }") if np.isnan(CPK): CPK = 1e9 @@ -1584,7 +1644,7 @@ def calculate_CPK_power_consumption(analysis_subdir: str): # else: # return CPK - return {"CPK": CPK, "Edot": Edot, "diss_surf": surface_diss, "diss_shock": shock_diss} + return {"CPK": CPK, "Edot": Edot, "diss_surf": surface_diss, "diss_shock": shock_diss, "Cd": forces["Cd"]} def calculate_CPK_mses(analysis_subdir: str, configuration: str = "underwing_te"): diff --git a/pymead/analysis/compressible_flow.py b/pymead/analysis/compressible_flow.py new file mode 100644 index 00000000..9d4cdf84 --- /dev/null +++ b/pymead/analysis/compressible_flow.py @@ -0,0 +1,8 @@ + + +def calculate_normal_shock_total_pressure_ratio(M_up: float, gam: float): + A = (gam + 1) / 2 * M_up**2 + B = 1 + (gam - 1) / 2 * M_up**2 + C = 2 * gam * M_up**2 / (gam + 1) + D = (gam - 1) / (gam + 1) + return (A / B) ** (gam / (gam - 1)) * (C - D) ** (1 / (1 - gam)) diff --git a/pymead/analysis/read_aero_data.py b/pymead/analysis/read_aero_data.py index 31c2a3da..4224b610 100644 --- a/pymead/analysis/read_aero_data.py +++ b/pymead/analysis/read_aero_data.py @@ -6,7 +6,7 @@ import numpy as np -flow_var_idx = {'M': 7, 'Cp': 8, 'p': 3, 'rho': 2, 'u': 4, 'v': 5, 'q': 6, "Cpt": 9, "dCpt": 10} +flow_var_idx = {'M': 7, 'Cp': 8, 'p': 3, 'rho': 2, 'u': 4, 'v': 5, 'q': 6, "Cpt": 9, "dCpt": 10, "dCp": 11} def read_Cp_from_file_xfoil(fname: str): @@ -350,13 +350,19 @@ def read_field_from_mses(src_file: str, M_inf: float = None, gam: float = None): if M_inf is not None and gam is not None: Cpt = field_array[3][:, :] * (1 + (gam - 1) / 2 * field_array[7][:, :] ** 2) ** (gam / (gam - 1)) * (2 / gam / M_inf**2) + Cp = field_array[flow_var_idx["Cp"]][:, :] Cpt = np.array([Cpt]) + Cp = np.array([Cp]) field_array = np.concatenate((field_array, Cpt)) dCpt = np.zeros(Cpt[0].shape) + dCp = np.zeros(Cpt[0].shape) for idx in range(1, Cpt[0].shape[0]): dCpt[idx, :] = Cpt[0][idx, :] - Cpt[0][idx - 1, :] + dCp[idx, :] = Cp[0][idx, :] - Cp[0][idx - 1, :] dCpt = np.array([dCpt]) + dCp = np.array([dCp]) field_array = np.concatenate((field_array, dCpt)) + field_array = np.concatenate((field_array, dCp)) return field_array diff --git a/pymead/post/post_process.py b/pymead/post/post_process.py index 1e69ac77..cdf241da 100644 --- a/pymead/post/post_process.py +++ b/pymead/post/post_process.py @@ -812,7 +812,8 @@ def generate_single_field(self, var: str, index: int, index_list, cmap_field: mp 'v': r'Velocity-y ($v/V_\infty$)', 'q': r'Speed of Sound ($q/V_\infty$)', "Cpt": r"Total Pressure Over P_inf", - "dCpt": r"Delta Total Pressure"} + "dCpt": r"Delta Total Pressure", + "dCp": r"Delta Pressure Coefficient"} post_process_forces = load_data( os.path.splitext( self.post_process_force_file)[0] + weight_str + os.path.splitext(self.post_process_force_file)[1]) diff --git a/pymead/tests/aero_tests/Edot_test.py b/pymead/tests/aero_tests/Edot_test.py index 3d41c90e..00581f4f 100644 --- a/pymead/tests/aero_tests/Edot_test.py +++ b/pymead/tests/aero_tests/Edot_test.py @@ -33,11 +33,11 @@ def test_Edot_opt_oblique_shock_4(self): analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_82\analysis\analysis_127_w70-30" Edot = calculate_CPK_power_consumption(analysis_dir) print(f"Opt: {Edot = }") - - def test_Edot_opt_oblique_shock_5(self): - analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_82\analysis\analysis_127_w100-0" - Edot = calculate_CPK_power_consumption(analysis_dir) - print(f"Opt: {Edot = }") + # + # def test_Edot_opt_oblique_shock_5(self): + # analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_82\analysis\analysis_127_w100-0" + # Edot = calculate_CPK_power_consumption(analysis_dir) + # print(f"Opt: {Edot = }") def test_Edot_opt_6(self): analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\analysis_temp\ga_airfoil_7"