From 98853a130295e281821015fad2dd79dc7ff05763 Mon Sep 17 00:00:00 2001 From: Ben Preston <144227999+ben-l-p@users.noreply.github.com> Date: Thu, 26 Sep 2024 09:32:25 +0100 Subject: [PATCH 1/3] Add rewritten aeroforcescalculator.py using flow angle rotation --- sharpy/postproc/aeroforcescalculator.py | 259 +++++++++++------------- 1 file changed, 114 insertions(+), 145 deletions(-) diff --git a/sharpy/postproc/aeroforcescalculator.py b/sharpy/postproc/aeroforcescalculator.py index cb6833ff3..10e18c5db 100644 --- a/sharpy/postproc/aeroforcescalculator.py +++ b/sharpy/postproc/aeroforcescalculator.py @@ -5,15 +5,13 @@ from sharpy.utils.solver_interface import solver, BaseSolver import sharpy.utils.settings as settings_utils import sharpy.utils.algebra as algebra -import sharpy.aero.utils.mapping as mapping -import warnings @solver class AeroForcesCalculator(BaseSolver): """AeroForcesCalculator - Calculates the total aerodynamic forces and moments on the frame of reference ``A``. + Calculates the total aerodynamic forces and moments on the body ``A`` and inertial ``G`` frame of reference . """ solver_id = 'AeroForcesCalculator' @@ -63,6 +61,10 @@ class AeroForcesCalculator(BaseSolver): settings_default['c_ref'] = 1 settings_description['c_ref'] = 'Reference chord' + settings_types['u_inf_dir'] = 'list(float)' + settings_default['u_inf_dir'] = [1., 0., 0.] + settings_description['u_inf_dir'] = 'Flow direction' + settings_table = settings_utils.SettingsTable() __doc__ += settings_table.generate(settings_types, settings_default, settings_description) @@ -76,7 +78,6 @@ def __init__(self): self.caller = None self.table = None - self.rot = None self.moment_reference_location = np.array([0., 0., 0.]) def initialise(self, data, custom_settings=None, caller=None, restart=False): @@ -99,10 +100,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False): self.table.print_header(['tstep', 'fx_g', 'fy_g', 'fz_g', 'mx_g', 'my_g', 'mz_g']) def run(self, **kwargs): - - online = settings_utils.set_value_or_default(kwargs, 'online', False) - - if online: + if kwargs.get('online', False): self.ts_max = len(self.data.structure.timestep_info) self.calculate_forces(-1) @@ -118,174 +116,145 @@ def run(self, **kwargs): if self.settings['write_text_file']: self.file_output(self.settings['text_file_name']) return self.data - + + # generate forcing data for a given timestep def calculate_forces(self, ts): - # Forces per surface in G frame - self.rot = algebra.quat2rotation(self.data.structure.timestep_info[ts].quat) - force = self.data.aero.timestep_info[ts].forces - unsteady_force = self.data.aero.timestep_info[ts].dynamic_forces - for i_surf in range(self.data.aero.n_surf): - ( - self.data.aero.timestep_info[ts].inertial_steady_forces[i_surf, 0:3], - self.data.aero.timestep_info[ts].inertial_unsteady_forces[i_surf, 0:3], - self.data.aero.timestep_info[ts].body_steady_forces[i_surf, 0:3], - self.data.aero.timestep_info[ts].body_unsteady_forces[i_surf, 0:3] - ) = self.calculate_forces_for_isurf_in_g_frame(force[i_surf], unsteady_force=unsteady_force[i_surf]) - - if self.settings["nonlifting_body"]: - for i_surf in range(self.data.nonlifting_body.n_surf): - ( - self.data.nonlifting_body.timestep_info[ts].inertial_steady_forces[i_surf, 0:3], - self.data.nonlifting_body.timestep_info[ts].body_steady_forces[i_surf, 0:3], - ) = self.calculate_forces_for_isurf_in_g_frame(self.data.nonlifting_body.timestep_info[ts].forces[i_surf], nonlifting=True) - - # Convert to forces in B frame - try: - steady_forces_b = self.data.structure.timestep_info[ts].postproc_node['aero_steady_forces'] - except KeyError: - if self.settings["nonlifting_body"]: - warnings.warn('Nonlifting forces are not considered in aero forces calculation since forces cannot not be retrieved from postproc node.') - steady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, force) - - try: - unsteady_forces_b = self.data.structure.timestep_info[ts].postproc_node['aero_unsteady_forces'] - except KeyError: - unsteady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, unsteady_force) - # Convert to forces in A frame - steady_forces_a = self.data.structure.timestep_info[ts].nodal_b_for_2_a_for(steady_forces_b, - self.data.structure) - - unsteady_forces_a = self.data.structure.timestep_info[ts].nodal_b_for_2_a_for(unsteady_forces_b, - self.data.structure) - - # Express total forces in A frame - self.data.aero.timestep_info[ts].total_steady_body_forces = \ - mapping.total_forces_moments(steady_forces_a, - self.data.structure.timestep_info[ts].pos, - ref_pos=self.moment_reference_location) - self.data.aero.timestep_info[ts].total_unsteady_body_forces = \ - mapping.total_forces_moments(unsteady_forces_a, - self.data.structure.timestep_info[ts].pos, - ref_pos=self.moment_reference_location) - - # Express total forces in G frame - self.data.aero.timestep_info[ts].total_steady_inertial_forces = \ - np.block([[self.rot, np.zeros((3, 3))], - [np.zeros((3, 3)), self.rot]]).dot( - self.data.aero.timestep_info[ts].total_steady_body_forces) - - self.data.aero.timestep_info[ts].total_unsteady_inertial_forces = \ - np.block([[self.rot, np.zeros((3, 3))], - [np.zeros((3, 3)), self.rot]]).dot( - self.data.aero.timestep_info[ts].total_unsteady_body_forces) - - def calculate_forces_for_isurf_in_g_frame(self, force, unsteady_force = None, nonlifting = False): - """ - Forces for a surface in G frame - """ - # Forces per surface in G frame - total_steady_force = np.zeros((3,)) - total_unsteady_force = np.zeros((3,)) - _, n_rows, n_cols = force.shape - for i_m in range(n_rows): - for i_n in range(n_cols): - total_steady_force += force[0:3, i_m, i_n] - if not nonlifting: - total_unsteady_force += unsteady_force[0:3, i_m, i_n] - if not nonlifting: - return total_steady_force, total_unsteady_force, np.dot(self.rot.T, total_steady_force), np.dot(self.rot.T, total_unsteady_force) - else: - return total_steady_force, np.dot(self.rot.T, total_steady_force) + # body orientation rotation + rot_ag = algebra.quat2rotation(self.data.structure.timestep_info[ts].quat).T # R_AG + + # flow rotation angle from x and z components of flow direction + # WARNING: this will give incorrect results when there is sideslip (u_inf_dir[1] != 0.) + rot_xy = np.arctan2(-self.settings['u_inf_dir'][2], self.settings['u_inf_dir'][0]) + rmat_xy = algebra.euler2rot((0., rot_xy, 0.)) + + # number of aero surfaces + n_surf = self.data.aero.n_surf + # forces per node per wing surface + # naming convention: s=steady, u=unsteady, a=body for, g=inertial for + force_s_a = [np.squeeze(rmat_xy.T @ np.expand_dims(np.moveaxis( + self.data.aero.timestep_info[ts].forces[i], 0, -1)[..., :3], -1)) for i in range(n_surf)] + force_u_a = [np.squeeze(rmat_xy.T @ np.expand_dims(np.moveaxis( + self.data.aero.timestep_info[ts].dynamic_forces[i], 0, -1)[..., :3], -1)) for i in range(n_surf)] - def map_forces_beam_dof(self, aero_data, ts, force): - struct_tstep = self.data.structure.timestep_info[ts] - aero_forces_beam_dof = mapping.aero2struct_force_mapping(force, - aero_data.struct2aero_mapping, - aero_data.timestep_info[ts].zeta, - struct_tstep.pos, - struct_tstep.psi, - None, - self.data.structure.connectivities, - struct_tstep.cag()) - return aero_forces_beam_dof + # add nonlifting: TODO: check FoR of nonlifting forces + if self.settings["nonlifting_body"]: + n_surf_nonlift = self.data.nonlifting_body.n_surf + force_s_a.extend([np.squeeze(rmat_xy.T @ np.expand_dims(np.moveaxis( + self.data.nonlifting_body.timestep_info[ts].inertial_steady_forces[i], 0, -1)[..., :3], -1)) + for i in range(n_surf_nonlift)]) + + # list of steady and unsteady body and inertial forces per surface + f_s_a_surf = [np.sum(force_s_a[i], axis=(0, 1)) for i in range(n_surf)] + f_s_g_surf = [rot_ag.T @ f_s_a_surf[i] for i in range(n_surf)] + f_u_a_surf = [np.sum(force_u_a[i], axis=(0, 1)) for i in range(n_surf)] + f_u_g_surf = [rot_ag.T @ f_u_a_surf[i] for i in range(n_surf)] + + # calculate moments + # surface dimensions + m = [surf.shape[0] for surf in force_s_a] + n = [surf.shape[1] for surf in force_s_a] + + # aerogrid node body position + zeta_a = [np.squeeze(rot_ag @ np.expand_dims(np.moveaxis( + self.data.aero.timestep_info[ts].zeta[i], 0, -1), axis=-1)) for i in range(n_surf)] + + # add nonlifting body node positions + if self.settings["nonlifting_body"]: + zeta_a.extend([np.squeeze(rot_ag @ np.expand_dims(np.moveaxis( + self.data.nonlifting_body.timestep_info[ts].zeta[i], 0, -1), axis=-1)) for i in range(n_surf)]) + + # node body position relative to reference position and its skew + r_a = [zeta_a[i] - np.tile(self.moment_reference_location.reshape(1, 1, 3), (m[i], n[i], 1)) for i in + range(n_surf)] + r_a_tilde = [np.apply_along_axis(algebra.skew, 2, r_a[i]) for i in range(n_surf)] + + # steady moments per node per surface + mom_s_a = [np.squeeze(r_a_tilde[i] @ np.expand_dims(force_s_a[i], axis=-1)) for i in range(n_surf)] + + # list of steady body and inertial forces per surface + mom_s_a_surf = [np.sum(mom_s_a[i], axis=(0, 1)) for i in range(n_surf)] + mom_s_g_surf = [rot_ag.T @ mom_s_a_surf[i] for i in range(n_surf)] + + # unsteady moments per node per surface + mom_u_a = [np.squeeze(r_a_tilde[i] @ np.expand_dims(force_u_a[i], axis=-1)) for i in range(n_surf)] + + # list of unsteady body and inertial forces per surface + mom_u_a_surf = [np.sum(mom_u_a[i], axis=(0, 1)) for i in range(n_surf)] + mom_u_g_surf = [rot_ag.T @ mom_u_a_surf[i] for i in range(n_surf)] + + # save steady forcing to timestep data + self.data.aero.timestep_info[ts].body_steady_forces = np.hstack((f_s_a_surf, mom_s_a_surf)) + self.data.aero.timestep_info[ts].inertial_steady_forces = np.hstack((f_s_g_surf, mom_s_g_surf)) + self.data.aero.timestep_info[ts].total_steady_body_forces \ + = np.sum(np.hstack((f_s_a_surf, mom_s_a_surf)), axis=0) + self.data.aero.timestep_info[ts].total_steady_inertial_forces \ + = np.sum(np.hstack((f_s_g_surf, mom_s_g_surf)), axis=0) + + # save unsteady forcing to timestep data + self.data.aero.timestep_info[ts].body_unsteady_forces = np.hstack((f_u_a_surf, mom_u_a_surf)) + self.data.aero.timestep_info[ts].inertial_unsteady_forces = np.hstack((f_u_g_surf, mom_u_g_surf)) + self.data.aero.timestep_info[ts].total_unsteady_body_forces \ + = np.sum(np.hstack((f_u_a_surf, mom_u_a_surf)), axis=0) + self.data.aero.timestep_info[ts].total_unsteady_inertial_forces \ + = np.sum(np.hstack((f_u_g_surf, mom_u_g_surf)), axis=0) def calculate_coefficients(self, fx, fy, fz, mx, my, mz): - qS = self.settings['q_ref'] * self.settings['S_ref'] - return fx/qS, fy/qS, fz/qS, mx/qS/self.settings['b_ref'], my/qS/self.settings['c_ref'], \ - mz/qS/self.settings['b_ref'] + q_s = self.settings['q_ref'] * self.settings['S_ref'] + return (fx / q_s, + fy / q_s, + fz / q_s, + mx / (q_s * self.settings['b_ref']), + my / (q_s * self.settings['c_ref']), + mz / (q_s * self.settings['b_ref'])) def screen_output(self, ts): # print time step total aero forces - forces = np.zeros(3) - moments = np.zeros(3) - aero_tstep = self.data.aero.timestep_info[ts] - forces += aero_tstep.total_steady_inertial_forces[:3] + aero_tstep.total_unsteady_inertial_forces[:3] - moments += aero_tstep.total_steady_inertial_forces[3:] + aero_tstep.total_unsteady_inertial_forces[3:] + forces = aero_tstep.total_steady_inertial_forces[:3] + aero_tstep.total_unsteady_inertial_forces[:3] + moments = aero_tstep.total_steady_inertial_forces[3:] + aero_tstep.total_unsteady_inertial_forces[3:] - if self.settings['coefficients']: # TODO: Check if coefficients have to be computed differently for fuselages - Cfx, Cfy, Cfz, Cmx, Cmy, Cmz = self.calculate_coefficients(*forces, *moments) - self.table.print_line([ts, Cfx, Cfy, Cfz, Cmx, Cmy, Cmz]) + if self.settings['coefficients']: # TODO: Check if coefficients have to be computed differently for fuselages + coeffs = self.calculate_coefficients(*forces, *moments) + self.table.print_line([ts, *coeffs]) else: self.table.print_line([ts, *forces, *moments]) def file_output(self, filename): - # assemble forces/moments matrix - # (1 timestep) + (3+3 inertial steady+unsteady) + (3+3 body steady+unsteady) - force_matrix = np.zeros((self.ts_max, 1 + 3 + 3 + 3 + 3 )) - moment_matrix = np.zeros((self.ts_max, 1 + 3 + 3 + 3 + 3)) + # assemble forces/moments matrix: (1 timestep) + (3+3 inertial steady+unsteady) + (3+3 body steady+unsteady) + force_matrix = np.zeros((self.ts_max, 13)) + moment_matrix = np.zeros((self.ts_max, 13)) for ts in range(self.ts_max): aero_tstep = self.data.aero.timestep_info[ts] - i = 0 - force_matrix[ts, i] = ts - moment_matrix[ts, i] = ts - i += 1 - + # timestep column + force_matrix[ts, 0] = moment_matrix[ts, 0] = ts # Steady forces/moments G - force_matrix[ts, i:i+3] = aero_tstep.total_steady_inertial_forces[:3] - moment_matrix[ts, i:i+3] = aero_tstep.total_steady_inertial_forces[3:] - i += 3 - + force_matrix[ts, 1:4], moment_matrix[ts, 1:4] = np.split(aero_tstep.total_steady_inertial_forces, 2) # Unsteady forces/moments G - force_matrix[ts, i:i+3] = aero_tstep.total_unsteady_inertial_forces[:3] - moment_matrix[ts, i:i+3] = aero_tstep.total_unsteady_inertial_forces[3:] - i += 3 - + force_matrix[ts, 4:7], moment_matrix[ts, 4:7] = np.split(aero_tstep.total_unsteady_inertial_forces, 2) # Steady forces/moments A - force_matrix[ts, i:i+3] = aero_tstep.total_steady_body_forces[:3] - moment_matrix[ts, i:i+3] = aero_tstep.total_steady_body_forces[3:] - i += 3 - + force_matrix[ts, 7:10], moment_matrix[ts, 7:10] = np.split(aero_tstep.total_steady_body_forces, 2) # Unsteady forces/moments A - force_matrix[ts, i:i+3] = aero_tstep.total_unsteady_body_forces[:3] - moment_matrix[ts, i:i+3] = aero_tstep.total_unsteady_body_forces[3:] - + force_matrix[ts, 10:13], moment_matrix[ts, 10:13] = np.split(aero_tstep.total_unsteady_body_forces, 2) - header = '' - header += 'tstep, ' - header += 'fx_steady_G, fy_steady_G, fz_steady_G, ' - header += 'fx_unsteady_G, fy_unsteady_G, fz_unsteady_G, ' - header += 'fx_steady_a, fy_steady_a, fz_steady_a, ' - header += 'fx_unsteady_a, fy_unsteady_a, fz_unsteady_a' + # save linear force data + header_f = ('tstep, fx_steady_G, fy_steady_G, fz_steady_G, fx_unsteady_G, fy_unsteady_G, fz_unsteady_G, ' + 'fx_steady_a, fy_steady_a, fz_steady_a, fx_unsteady_a, fy_unsteady_a, fz_unsteady_a') np.savetxt(self.folder + 'forces_' + filename, force_matrix, fmt='%i' + ', %10e' * (np.shape(force_matrix)[1] - 1), delimiter=',', - header=header, + header=header_f, comments='#') - header = '' - header += 'tstep, ' - header += 'mx_steady_G, my_steady_G, mz_steady_G, ' - header += 'mx_unsteady_G, my_unsteady_G, mz_unsteady_G, ' - header += 'mx_steady_a, my_steady_a, mz_steady_a, ' - header += 'mx_unsteady_a, my_unsteady_a, mz_unsteady_a' + # save moment data + header_m = ('tstep, mx_steady_G, my_steady_G, mz_steady_G, mx_unsteady_G, my_unsteady_G, mz_unsteady_G, ' + 'mx_steady_a, my_steady_a, mz_steady_a, mx_unsteady_a, rot_agmy_unsteady_a, mz_unsteady_a') np.savetxt(self.folder + 'moments_' + filename, moment_matrix, fmt='%i' + ', %10e' * (np.shape(moment_matrix)[1] - 1), delimiter=',', - header=header, + header=header_m, comments='#') From faf51261d64ad2d74efe75ec50ab08241db5fc9f Mon Sep 17 00:00:00 2001 From: Ben Preston Date: Fri, 27 Sep 2024 14:41:36 +0100 Subject: [PATCH 2/3] Fixed FoR issue with aero forces --- sharpy/postproc/aeroforcescalculator.py | 264 +++++++++++++++--------- sharpy/utils/datastructures.py | 113 +--------- tests/uvlm/static/polars/test_polars.py | 17 +- 3 files changed, 176 insertions(+), 218 deletions(-) diff --git a/sharpy/postproc/aeroforcescalculator.py b/sharpy/postproc/aeroforcescalculator.py index 10e18c5db..92af0daaf 100644 --- a/sharpy/postproc/aeroforcescalculator.py +++ b/sharpy/postproc/aeroforcescalculator.py @@ -5,13 +5,15 @@ from sharpy.utils.solver_interface import solver, BaseSolver import sharpy.utils.settings as settings_utils import sharpy.utils.algebra as algebra +import sharpy.aero.utils.mapping as mapping +import warnings @solver class AeroForcesCalculator(BaseSolver): """AeroForcesCalculator - Calculates the total aerodynamic forces and moments on the body ``A`` and inertial ``G`` frame of reference . + Calculates the total aerodynamic forces and moments on the frame of reference ``A``. """ solver_id = 'AeroForcesCalculator' @@ -78,6 +80,7 @@ def __init__(self): self.caller = None self.table = None + self.rot = None self.moment_reference_location = np.array([0., 0., 0.]) def initialise(self, data, custom_settings=None, caller=None, restart=False): @@ -100,7 +103,10 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False): self.table.print_header(['tstep', 'fx_g', 'fy_g', 'fz_g', 'mx_g', 'my_g', 'mz_g']) def run(self, **kwargs): - if kwargs.get('online', False): + + online = settings_utils.set_value_or_default(kwargs, 'online', False) + + if online: self.ts_max = len(self.data.structure.timestep_info) self.calculate_forces(-1) @@ -117,144 +123,200 @@ def run(self, **kwargs): self.file_output(self.settings['text_file_name']) return self.data - # generate forcing data for a given timestep def calculate_forces(self, ts): - # body orientation rotation - rot_ag = algebra.quat2rotation(self.data.structure.timestep_info[ts].quat).T # R_AG + self.rot = algebra.quat2rotation(self.data.structure.timestep_info[ts].quat) # R_GA # flow rotation angle from x and z components of flow direction - # WARNING: this will give incorrect results when there is sideslip (u_inf_dir[1] != 0.) + # WARNING: this will give incorrect results when there is sideslip rot_xy = np.arctan2(-self.settings['u_inf_dir'][2], self.settings['u_inf_dir'][0]) rmat_xy = algebra.euler2rot((0., rot_xy, 0.)) - # number of aero surfaces + # Forces per surface in A frame relative to the flow + force = self.data.aero.timestep_info[ts].forces + unsteady_force = self.data.aero.timestep_info[ts].dynamic_forces + n_surf = self.data.aero.n_surf - # forces per node per wing surface - # naming convention: s=steady, u=unsteady, a=body for, g=inertial for - force_s_a = [np.squeeze(rmat_xy.T @ np.expand_dims(np.moveaxis( - self.data.aero.timestep_info[ts].forces[i], 0, -1)[..., :3], -1)) for i in range(n_surf)] - force_u_a = [np.squeeze(rmat_xy.T @ np.expand_dims(np.moveaxis( - self.data.aero.timestep_info[ts].dynamic_forces[i], 0, -1)[..., :3], -1)) for i in range(n_surf)] + rmat = np.block([[self.rot @ rmat_xy.T, np.zeros((3, 3))], [np.zeros((3, 3)), self.rot @ rmat_xy.T]]) + force_g = [np.moveaxis(np.squeeze(rmat @ np.expand_dims(np.moveaxis( + force[i], (0, 1, 2), (2, 0, 1)), -1)), source=(0, 1, 2), destination=(1, 2, 0)) + for i in range(n_surf)] + unsteady_force_g = [np.moveaxis(np.squeeze(rmat @ np.expand_dims(np.moveaxis( + unsteady_force[i], (0, 1, 2), (2, 0, 1)), -1)), source=(0, 1, 2), destination=(1, 2, 0)) + for i in range(n_surf)] + + for i_surf in range(n_surf): + ( + self.data.aero.timestep_info[ts].inertial_steady_forces[i_surf, 0:3], + self.data.aero.timestep_info[ts].inertial_unsteady_forces[i_surf, 0:3], + self.data.aero.timestep_info[ts].body_steady_forces[i_surf, 0:3], + self.data.aero.timestep_info[ts].body_unsteady_forces[i_surf, 0:3] + ) = self.calculate_forces_for_isurf_in_g_frame(force_g[i_surf], unsteady_force=unsteady_force_g[i_surf]) - # add nonlifting: TODO: check FoR of nonlifting forces if self.settings["nonlifting_body"]: - n_surf_nonlift = self.data.nonlifting_body.n_surf - force_s_a.extend([np.squeeze(rmat_xy.T @ np.expand_dims(np.moveaxis( - self.data.nonlifting_body.timestep_info[ts].inertial_steady_forces[i], 0, -1)[..., :3], -1)) - for i in range(n_surf_nonlift)]) - - # list of steady and unsteady body and inertial forces per surface - f_s_a_surf = [np.sum(force_s_a[i], axis=(0, 1)) for i in range(n_surf)] - f_s_g_surf = [rot_ag.T @ f_s_a_surf[i] for i in range(n_surf)] - f_u_a_surf = [np.sum(force_u_a[i], axis=(0, 1)) for i in range(n_surf)] - f_u_g_surf = [rot_ag.T @ f_u_a_surf[i] for i in range(n_surf)] - - # calculate moments - # surface dimensions - m = [surf.shape[0] for surf in force_s_a] - n = [surf.shape[1] for surf in force_s_a] - - # aerogrid node body position - zeta_a = [np.squeeze(rot_ag @ np.expand_dims(np.moveaxis( - self.data.aero.timestep_info[ts].zeta[i], 0, -1), axis=-1)) for i in range(n_surf)] - - # add nonlifting body node positions - if self.settings["nonlifting_body"]: - zeta_a.extend([np.squeeze(rot_ag @ np.expand_dims(np.moveaxis( - self.data.nonlifting_body.timestep_info[ts].zeta[i], 0, -1), axis=-1)) for i in range(n_surf)]) - - # node body position relative to reference position and its skew - r_a = [zeta_a[i] - np.tile(self.moment_reference_location.reshape(1, 1, 3), (m[i], n[i], 1)) for i in - range(n_surf)] - r_a_tilde = [np.apply_along_axis(algebra.skew, 2, r_a[i]) for i in range(n_surf)] - - # steady moments per node per surface - mom_s_a = [np.squeeze(r_a_tilde[i] @ np.expand_dims(force_s_a[i], axis=-1)) for i in range(n_surf)] - - # list of steady body and inertial forces per surface - mom_s_a_surf = [np.sum(mom_s_a[i], axis=(0, 1)) for i in range(n_surf)] - mom_s_g_surf = [rot_ag.T @ mom_s_a_surf[i] for i in range(n_surf)] - - # unsteady moments per node per surface - mom_u_a = [np.squeeze(r_a_tilde[i] @ np.expand_dims(force_u_a[i], axis=-1)) for i in range(n_surf)] - - # list of unsteady body and inertial forces per surface - mom_u_a_surf = [np.sum(mom_u_a[i], axis=(0, 1)) for i in range(n_surf)] - mom_u_g_surf = [rot_ag.T @ mom_u_a_surf[i] for i in range(n_surf)] - - # save steady forcing to timestep data - self.data.aero.timestep_info[ts].body_steady_forces = np.hstack((f_s_a_surf, mom_s_a_surf)) - self.data.aero.timestep_info[ts].inertial_steady_forces = np.hstack((f_s_g_surf, mom_s_g_surf)) - self.data.aero.timestep_info[ts].total_steady_body_forces \ - = np.sum(np.hstack((f_s_a_surf, mom_s_a_surf)), axis=0) - self.data.aero.timestep_info[ts].total_steady_inertial_forces \ - = np.sum(np.hstack((f_s_g_surf, mom_s_g_surf)), axis=0) - - # save unsteady forcing to timestep data - self.data.aero.timestep_info[ts].body_unsteady_forces = np.hstack((f_u_a_surf, mom_u_a_surf)) - self.data.aero.timestep_info[ts].inertial_unsteady_forces = np.hstack((f_u_g_surf, mom_u_g_surf)) - self.data.aero.timestep_info[ts].total_unsteady_body_forces \ - = np.sum(np.hstack((f_u_a_surf, mom_u_a_surf)), axis=0) - self.data.aero.timestep_info[ts].total_unsteady_inertial_forces \ - = np.sum(np.hstack((f_u_g_surf, mom_u_g_surf)), axis=0) + for i_surf in range(self.data.nonlifting_body.n_surf): + ( + self.data.nonlifting_body.timestep_info[ts].inertial_steady_forces[i_surf, 0:3], + self.data.nonlifting_body.timestep_info[ts].body_steady_forces[i_surf, 0:3], + ) = self.calculate_forces_for_isurf_in_g_frame( + self.data.nonlifting_body.timestep_info[ts].forces[i_surf], nonlifting=True) + + # Convert to forces in B frame + try: + steady_forces_b = self.data.structure.timestep_info[ts].postproc_node['aero_steady_forces'] + except KeyError: + if self.settings["nonlifting_body"]: + warnings.warn( + 'Nonlifting forces are not considered in aero forces calculation since forces cannot not be retrieved from postproc node.') + steady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, force_g) + + try: + unsteady_forces_b = self.data.structure.timestep_info[ts].postproc_node['aero_unsteady_forces'] + except KeyError: + unsteady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, unsteady_force_g) + # Convert to forces in A frame + steady_forces_a = self.data.structure.timestep_info[ts].nodal_b_for_2_a_for(steady_forces_b, + self.data.structure) + + unsteady_forces_a = self.data.structure.timestep_info[ts].nodal_b_for_2_a_for(unsteady_forces_b, + self.data.structure) + + # Express total forces in A frame + # self.data.aero.timestep_info[ts].total_steady_body_forces = \ + # mapping.total_forces_moments(steady_forces_a, + # self.data.structure.timestep_info[ts].pos, + # ref_pos=self.moment_reference_location) + # self.data.aero.timestep_info[ts].total_unsteady_body_forces = \ + # mapping.total_forces_moments(unsteady_forces_a, + # self.data.structure.timestep_info[ts].pos, + # ref_pos=self.moment_reference_location) + + self.data.aero.timestep_info[ts].total_steady_body_forces = \ + rmat @ mapping.total_forces_moments(steady_forces_a, + self.data.structure.timestep_info[ts].pos, + ref_pos=self.moment_reference_location) + self.data.aero.timestep_info[ts].total_unsteady_body_forces = \ + rmat @ mapping.total_forces_moments(unsteady_forces_a, + self.data.structure.timestep_info[ts].pos, + ref_pos=self.moment_reference_location) + + # Express total forces in G frame + self.data.aero.timestep_info[ts].total_steady_inertial_forces = \ + np.block([[self.rot, np.zeros((3, 3))], + [np.zeros((3, 3)), self.rot]]).dot( + self.data.aero.timestep_info[ts].total_steady_body_forces) + + self.data.aero.timestep_info[ts].total_unsteady_inertial_forces = \ + np.block([[self.rot, np.zeros((3, 3))], + [np.zeros((3, 3)), self.rot]]).dot( + self.data.aero.timestep_info[ts].total_unsteady_body_forces) + + def calculate_forces_for_isurf_in_g_frame(self, force, unsteady_force=None, nonlifting=False): + """ + Forces for a surface in G frame + """ + # Forces per surface in G frame + total_steady_force = np.zeros((3,)) + total_unsteady_force = np.zeros((3,)) + _, n_rows, n_cols = force.shape + for i_m in range(n_rows): + for i_n in range(n_cols): + total_steady_force += force[0:3, i_m, i_n] + if not nonlifting: + total_unsteady_force += unsteady_force[0:3, i_m, i_n] + if not nonlifting: + return total_steady_force, total_unsteady_force, np.dot(self.rot.T, total_steady_force), np.dot(self.rot.T, + total_unsteady_force) + else: + return total_steady_force, np.dot(self.rot.T, total_steady_force) + + def map_forces_beam_dof(self, aero_data, ts, force): + struct_tstep = self.data.structure.timestep_info[ts] + aero_forces_beam_dof = mapping.aero2struct_force_mapping(force, + aero_data.struct2aero_mapping, + aero_data.timestep_info[ts].zeta, + struct_tstep.pos, + struct_tstep.psi, + None, + self.data.structure.connectivities, + struct_tstep.cag()) + return aero_forces_beam_dof def calculate_coefficients(self, fx, fy, fz, mx, my, mz): - q_s = self.settings['q_ref'] * self.settings['S_ref'] - return (fx / q_s, - fy / q_s, - fz / q_s, - mx / (q_s * self.settings['b_ref']), - my / (q_s * self.settings['c_ref']), - mz / (q_s * self.settings['b_ref'])) + qS = self.settings['q_ref'] * self.settings['S_ref'] + return fx / qS, fy / qS, fz / qS, mx / qS / self.settings['b_ref'], my / qS / self.settings['c_ref'], \ + mz / qS / self.settings['b_ref'] def screen_output(self, ts): # print time step total aero forces + forces = np.zeros(3) + moments = np.zeros(3) + aero_tstep = self.data.aero.timestep_info[ts] - forces = aero_tstep.total_steady_inertial_forces[:3] + aero_tstep.total_unsteady_inertial_forces[:3] - moments = aero_tstep.total_steady_inertial_forces[3:] + aero_tstep.total_unsteady_inertial_forces[3:] + forces += aero_tstep.total_steady_inertial_forces[:3] + aero_tstep.total_unsteady_inertial_forces[:3] + moments += aero_tstep.total_steady_inertial_forces[3:] + aero_tstep.total_unsteady_inertial_forces[3:] if self.settings['coefficients']: # TODO: Check if coefficients have to be computed differently for fuselages - coeffs = self.calculate_coefficients(*forces, *moments) - self.table.print_line([ts, *coeffs]) + Cfx, Cfy, Cfz, Cmx, Cmy, Cmz = self.calculate_coefficients(*forces, *moments) + self.table.print_line([ts, Cfx, Cfy, Cfz, Cmx, Cmy, Cmz]) else: self.table.print_line([ts, *forces, *moments]) def file_output(self, filename): - # assemble forces/moments matrix: (1 timestep) + (3+3 inertial steady+unsteady) + (3+3 body steady+unsteady) - force_matrix = np.zeros((self.ts_max, 13)) - moment_matrix = np.zeros((self.ts_max, 13)) + # assemble forces/moments matrix + # (1 timestep) + (3+3 inertial steady+unsteady) + (3+3 body steady+unsteady) + force_matrix = np.zeros((self.ts_max, 1 + 3 + 3 + 3 + 3)) + moment_matrix = np.zeros((self.ts_max, 1 + 3 + 3 + 3 + 3)) for ts in range(self.ts_max): aero_tstep = self.data.aero.timestep_info[ts] - # timestep column - force_matrix[ts, 0] = moment_matrix[ts, 0] = ts + i = 0 + force_matrix[ts, i] = ts + moment_matrix[ts, i] = ts + i += 1 + # Steady forces/moments G - force_matrix[ts, 1:4], moment_matrix[ts, 1:4] = np.split(aero_tstep.total_steady_inertial_forces, 2) + force_matrix[ts, i:i + 3] = aero_tstep.total_steady_inertial_forces[:3] + moment_matrix[ts, i:i + 3] = aero_tstep.total_steady_inertial_forces[3:] + i += 3 + # Unsteady forces/moments G - force_matrix[ts, 4:7], moment_matrix[ts, 4:7] = np.split(aero_tstep.total_unsteady_inertial_forces, 2) + force_matrix[ts, i:i + 3] = aero_tstep.total_unsteady_inertial_forces[:3] + moment_matrix[ts, i:i + 3] = aero_tstep.total_unsteady_inertial_forces[3:] + i += 3 + # Steady forces/moments A - force_matrix[ts, 7:10], moment_matrix[ts, 7:10] = np.split(aero_tstep.total_steady_body_forces, 2) + force_matrix[ts, i:i + 3] = aero_tstep.total_steady_body_forces[:3] + moment_matrix[ts, i:i + 3] = aero_tstep.total_steady_body_forces[3:] + i += 3 + # Unsteady forces/moments A - force_matrix[ts, 10:13], moment_matrix[ts, 10:13] = np.split(aero_tstep.total_unsteady_body_forces, 2) + force_matrix[ts, i:i + 3] = aero_tstep.total_unsteady_body_forces[:3] + moment_matrix[ts, i:i + 3] = aero_tstep.total_unsteady_body_forces[3:] - # save linear force data - header_f = ('tstep, fx_steady_G, fy_steady_G, fz_steady_G, fx_unsteady_G, fy_unsteady_G, fz_unsteady_G, ' - 'fx_steady_a, fy_steady_a, fz_steady_a, fx_unsteady_a, fy_unsteady_a, fz_unsteady_a') + header = '' + header += 'tstep, ' + header += 'fx_steady_G, fy_steady_G, fz_steady_G, ' + header += 'fx_unsteady_G, fy_unsteady_G, fz_unsteady_G, ' + header += 'fx_steady_a, fy_steady_a, fz_steady_a, ' + header += 'fx_unsteady_a, fy_unsteady_a, fz_unsteady_a' np.savetxt(self.folder + 'forces_' + filename, force_matrix, fmt='%i' + ', %10e' * (np.shape(force_matrix)[1] - 1), delimiter=',', - header=header_f, + header=header, comments='#') - # save moment data - header_m = ('tstep, mx_steady_G, my_steady_G, mz_steady_G, mx_unsteady_G, my_unsteady_G, mz_unsteady_G, ' - 'mx_steady_a, my_steady_a, mz_steady_a, mx_unsteady_a, rot_agmy_unsteady_a, mz_unsteady_a') + header = '' + header += 'tstep, ' + header += 'mx_steady_G, my_steady_G, mz_steady_G, ' + header += 'mx_unsteady_G, my_unsteady_G, mz_unsteady_G, ' + header += 'mx_steady_a, my_steady_a, mz_steady_a, ' + header += 'mx_unsteady_a, my_unsteady_a, mz_unsteady_a' np.savetxt(self.folder + 'moments_' + filename, moment_matrix, fmt='%i' + ', %10e' * (np.shape(moment_matrix)[1] - 1), delimiter=',', - header=header_m, + header=header, comments='#') diff --git a/sharpy/utils/datastructures.py b/sharpy/utils/datastructures.py index 29222ad3d..2f90153df 100644 --- a/sharpy/utils/datastructures.py +++ b/sharpy/utils/datastructures.py @@ -431,6 +431,10 @@ class AeroTimeStepInfo(TimeStepInfo): ``[num_surf x chordwise panels x spanwise panels]`` dimensions_star (np.ndarray): Matrix defining the dimensions of the vortex grid on wakes ``[num_surf x streamwise panels x spanwise panels]`` + + Note, Ben Preston 27/09/24: forces and dynamic forces are stated to be in ``G`` however were actually + determined to be in ``A``. This issue may apply to other parameters, however errors due to this were only apparent + in the AeroForcesCalculator postprocessor. """ def __init__(self, dimensions, dimensions_star): super().__init__(dimensions) @@ -948,22 +952,6 @@ def change_to_local_AFoR(self, for0_pos, for0_vel, quat0): # Modify local rotations for ielem in range(self.psi.shape[0]): for inode in range(3): - # psi_master = self.psi[ielem, inode, :] + np.zeros((3,),) - # self.psi[ielem, inode, :] = np.dot(Csm, self.psi[ielem, inode, :]) - # self.psi_dot[ielem, inode, :] = (np.dot(Csm, self.psi_dot[ielem, inode, :] + algebra.cross3(for0_vel[3:6], psi_master)) - - # algebra.multiply_matrices(CAslaveG, algebra.skew(self.for_vel[3:6]), CGAmaster, psi_master)) - - # psi_master = self.psi[ielem,inode,:] + np.zeros((3,),) - # self.psi[ielem, inode, :] = algebra.rotation2crv(np.dot(Csm,algebra.crv2rotation(self.psi[ielem,inode,:]))) - # psi_slave = self.psi[ielem, inode, :] + np.zeros((3,),) - # cbam = algebra.crv2rotation(psi_master).T - # cbas = algebra.crv2rotation(psi_slave).T - # tm = algebra.crv2tan(psi_master) - # inv_ts = np.linalg.inv(algebra.crv2tan(psi_slave)) - - # self.psi_dot[ielem, inode, :] = np.dot(inv_ts, (np.dot(tm, self.psi_dot[ielem, inode, :]) + - # np.dot(cbam, for0_vel[3:6]) - - # np.dot(cbas, self.for_vel[3:6]))) self.psi[ielem, inode, :] = self.psi_local[ielem,inode,:].copy() self.psi_dot[ielem, inode, :] = self.psi_dot_local[ielem, inode, :].copy() @@ -997,12 +985,6 @@ def change_to_global_AFoR(self, for0_pos, for0_vel, quat0): for ielem in range(self.psi.shape[0]): for inode in range(3): - # psi_slave = self.psi[ielem,inode,:] + np.zeros((3,),) - # self.psi[ielem, inode, :] = np.dot(Cms, self.psi[ielem, inode, :]) - # self.psi_dot[ielem, inode, :] = (np.dot(Cms, self.psi_dot[ielem, inode, :] + algebra.cross3(self.for_vel[3:6], psi_slave)) - - # algebra.multiply_matrices(CAmasterG, algebra.skew(self.for0_vel[3:6]), CGAslave, psi_slave)) - - self.psi_local[ielem, inode, :] = self.psi[ielem, inode, :].copy() # Copy here the result from the structural computation self.psi_dot_local[ielem, inode, :] = self.psi_dot[ielem, inode, :].copy() # Copy here the result from the structural computation @@ -1020,93 +1002,6 @@ def change_to_global_AFoR(self, for0_pos, for0_vel, quat0): np.dot(cbas, self.for_vel[3:6]) - np.dot(cbam, for0_vel[3:6]))) - # def whole_structure_to_local_AFoR(self, beam, compute_psi_local=False): - # """ - # Same as change_to_local_AFoR but for a multibody structure - # - # Args: - # beam(sharpy.structure.models.beam.Beam): Beam structure of ``PreSharpy`` - # """ - # if not self.in_global_AFoR: - # raise NotImplementedError("Wrong managing of FoR") - # - # self.in_global_AFoR = False - # quat0 = self.quat.astype(dtype=ct.c_double, order='F', copy=True) - # for0_pos = self.for_pos.astype(dtype=ct.c_double, order='F', copy=True) - # for0_vel = self.for_vel.astype(dtype=ct.c_double, order='F', copy=True) - # - # MB_beam = [None]*beam.num_bodies - # MB_tstep = [None]*beam.num_bodies - # - # for ibody in range(beam.num_bodies): - # MB_beam[ibody] = beam.get_body(ibody = ibody) - # MB_tstep[ibody] = self.get_body(beam, MB_beam[ibody].num_dof, ibody = ibody) - # if compute_psi_local: - # MB_tstep[ibody].compute_psi_local_AFoR(for0_pos, for0_vel, quat0) - # MB_tstep[ibody].change_to_local_AFoR(for0_pos, for0_vel, quat0) - # - # first_dof = 0 - # for ibody in range(beam.num_bodies): - # # Renaming for clarity - # ibody_elems = MB_beam[ibody].global_elems_num - # ibody_nodes = MB_beam[ibody].global_nodes_num - # - # # Merge tstep - # self.pos[ibody_nodes,:] = MB_tstep[ibody].pos.astype(dtype=ct.c_double, order='F', copy=True) - # self.psi[ibody_elems,:,:] = MB_tstep[ibody].psi.astype(dtype=ct.c_double, order='F', copy=True) - # self.psi_local[ibody_elems,:,:] = MB_tstep[ibody].psi_local.astype(dtype=ct.c_double, order='F', copy=True) - # self.gravity_forces[ibody_nodes,:] = MB_tstep[ibody].gravity_forces.astype(dtype=ct.c_double, order='F', copy=True) - # - # self.pos_dot[ibody_nodes,:] = MB_tstep[ibody].pos_dot.astype(dtype=ct.c_double, order='F', copy=True) - # self.psi_dot[ibody_elems,:,:] = MB_tstep[ibody].psi_dot.astype(dtype=ct.c_double, order='F', copy=True) - # self.psi_dot_local[ibody_elems,:,:] = MB_tstep[ibody].psi_dot_local.astype(dtype=ct.c_double, order='F', copy=True) - # - # # TODO: Do I need a change in FoR for the following variables? Maybe for the FoR ones. - # # tstep.forces_constraints_nodes[ibody_nodes,:] = MB_tstep[ibody].forces_constraints_nodes.astype(dtype=ct.c_double, order='F', copy=True) - # # tstep.forces_constraints_FoR[ibody, :] = MB_tstep[ibody].forces_constraints_FoR[ibody, :].astype(dtype=ct.c_double, order='F', copy=True) - - # def whole_structure_to_global_AFoR(self, beam): - # """ - # Same as change_to_global_AFoR but for a multibody structure - # - # Args: - # beam(sharpy.structure.models.beam.Beam): Beam structure of ``PreSharpy`` - # """ - # if self.in_global_AFoR: - # raise NotImplementedError("Wrong managing of FoR") - # - # self.in_global_AFoR = True - # - # MB_beam = [None]*beam.num_bodies - # MB_tstep = [None]*beam.num_bodies - # quat0 = self.quat.astype(dtype=ct.c_double, order='F', copy=True) - # for0_pos = self.for_pos.astype(dtype=ct.c_double, order='F', copy=True) - # for0_vel = self.for_vel.astype(dtype=ct.c_double, order='F', copy=True) - # - # for ibody in range(beam.num_bodies): - # MB_beam[ibody] = beam.get_body(ibody = ibody) - # MB_tstep[ibody] = self.get_body(beam, MB_beam[ibody].num_dof, ibody = ibody) - # MB_tstep[ibody].change_to_global_AFoR(for0_pos, for0_vel, quat0) - # - # - # first_dof = 0 - # for ibody in range(beam.num_bodies): - # # Renaming for clarity - # ibody_elems = MB_beam[ibody].global_elems_num - # ibody_nodes = MB_beam[ibody].global_nodes_num - # - # # Merge tstep - # self.pos[ibody_nodes,:] = MB_tstep[ibody].pos.astype(dtype=ct.c_double, order='F', copy=True) - # self.psi[ibody_elems,:,:] = MB_tstep[ibody].psi.astype(dtype=ct.c_double, order='F', copy=True) - # self.psi_local[ibody_elems,:,:] = MB_tstep[ibody].psi_local.astype(dtype=ct.c_double, order='F', copy=True) - # self.gravity_forces[ibody_nodes,:] = MB_tstep[ibody].gravity_forces.astype(dtype=ct.c_double, order='F', - # copy=True) - # - # self.pos_dot[ibody_nodes,:] = MB_tstep[ibody].pos_dot.astype(dtype=ct.c_double, order='F', copy=True) - # self.psi_dot[ibody_elems,:,:] = MB_tstep[ibody].psi_dot.astype(dtype=ct.c_double, order='F', copy=True) - # self.psi_dot_local[ibody_elems,:,:] = MB_tstep[ibody].psi_dot_local.astype(dtype=ct.c_double, order='F', copy=True) - - def nodal_b_for_2_a_for(self, nodal, beam, filter=np.array([True]*6), ibody=None): """ Projects a nodal variable from the local, body-attached frame (B) to the reference A frame. diff --git a/tests/uvlm/static/polars/test_polars.py b/tests/uvlm/static/polars/test_polars.py index 31f50680d..43269a915 100644 --- a/tests/uvlm/static/polars/test_polars.py +++ b/tests/uvlm/static/polars/test_polars.py @@ -11,7 +11,7 @@ class InfiniteWing: area = 90000000.0 - chord = 3 + chord = 3. def force_coef(self, rho, uinf): return 0.5 * rho * uinf ** 2 * self.area @@ -53,8 +53,8 @@ def test_infinite_wing(self): results = postprocess(output_route + '/' + case_header + '/') - results[:, 1:3] /= wing.force_coef(1.225, 1) - results[:, -1] /= wing.moment_coef(1.225, 1) + results[:, 1:3] /= wing.force_coef(1.225, 1.) + results[:, -1] /= wing.moment_coef(1.225, 1.) with self.subTest('lift'): cl_polar = np.interp(results[:, 0], self.polar_data[:, 0], self.polar_data[:, 1]) @@ -95,12 +95,13 @@ def process_case(path_to_case): case_name = path_to_case.split('/')[-1] pmor = configobj.ConfigObj(path_to_case + f'/{case_name}.pmor.sharpy') alpha = pmor['parameters']['alpha'] - inertial_forces = np.loadtxt(f'{path_to_case}/forces/forces_aeroforces.txt', - skiprows=1, delimiter=',', dtype=float)[1:4] - inertial_moments = np.loadtxt(f'{path_to_case}/forces/moments_aeroforces.txt', - skiprows=1, delimiter=',', dtype=float)[1:4] - return alpha, inertial_forces[2], inertial_forces[0], inertial_moments[1] + body_forces = np.loadtxt(f'{path_to_case}/forces/forces_aeroforces.txt', + skiprows=1, delimiter=',', dtype=float)[7:10] + body_moments = np.loadtxt(f'{path_to_case}/forces/moments_aeroforces.txt', + skiprows=1, delimiter=',', dtype=float)[7:10] + + return alpha, body_forces[2], body_forces[0], body_moments[1] class TestStab(unittest.TestCase): From 86b8013c2460e37f9e281c64b3d7138132a63106 Mon Sep 17 00:00:00 2001 From: Ben Preston <144227999+ben-l-p@users.noreply.github.com> Date: Fri, 27 Sep 2024 15:47:41 +0100 Subject: [PATCH 3/3] Removed commented code in AeroForcesCalculator --- sharpy/postproc/aeroforcescalculator.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/sharpy/postproc/aeroforcescalculator.py b/sharpy/postproc/aeroforcescalculator.py index 92af0daaf..ea9511a7d 100644 --- a/sharpy/postproc/aeroforcescalculator.py +++ b/sharpy/postproc/aeroforcescalculator.py @@ -181,16 +181,6 @@ def calculate_forces(self, ts): unsteady_forces_a = self.data.structure.timestep_info[ts].nodal_b_for_2_a_for(unsteady_forces_b, self.data.structure) - # Express total forces in A frame - # self.data.aero.timestep_info[ts].total_steady_body_forces = \ - # mapping.total_forces_moments(steady_forces_a, - # self.data.structure.timestep_info[ts].pos, - # ref_pos=self.moment_reference_location) - # self.data.aero.timestep_info[ts].total_unsteady_body_forces = \ - # mapping.total_forces_moments(unsteady_forces_a, - # self.data.structure.timestep_info[ts].pos, - # ref_pos=self.moment_reference_location) - self.data.aero.timestep_info[ts].total_steady_body_forces = \ rmat @ mapping.total_forces_moments(steady_forces_a, self.data.structure.timestep_info[ts].pos,