diff --git a/sharpy/postproc/aeroforcescalculator.py b/sharpy/postproc/aeroforcescalculator.py index cb6833ff3..ea9511a7d 100644 --- a/sharpy/postproc/aeroforcescalculator.py +++ b/sharpy/postproc/aeroforcescalculator.py @@ -63,6 +63,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) @@ -118,39 +122,58 @@ def run(self, **kwargs): if self.settings['write_text_file']: self.file_output(self.settings['text_file_name']) return self.data - + def calculate_forces(self, ts): - # Forces per surface in G frame - self.rot = algebra.quat2rotation(self.data.structure.timestep_info[ts].quat) + 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 + rot_xy = np.arctan2(-self.settings['u_inf_dir'][2], self.settings['u_inf_dir'][0]) + rmat_xy = algebra.euler2rot((0., rot_xy, 0.)) + + # 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 - for i_surf in range(self.data.aero.n_surf): + + n_surf = self.data.aero.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[i_surf], unsteady_force=unsteady_force[i_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]) + 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) - + 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) + 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) + 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) @@ -158,15 +181,14 @@ 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) + 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 = \ - mapping.total_forces_moments(unsteady_forces_a, - self.data.structure.timestep_info[ts].pos, - ref_pos=self.moment_reference_location) + 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 = \ @@ -179,7 +201,7 @@ def calculate_forces(self, ts): [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): + def calculate_forces_for_isurf_in_g_frame(self, force, unsteady_force=None, nonlifting=False): """ Forces for a surface in G frame """ @@ -193,11 +215,11 @@ def calculate_forces_for_isurf_in_g_frame(self, force, unsteady_force = None, no 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) + 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, @@ -212,19 +234,19 @@ def map_forces_beam_dof(self, aero_data, ts, force): 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'] + 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) + 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:] - if self.settings['coefficients']: # TODO: Check if coefficients have to be computed differently for fuselages + 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]) else: @@ -233,7 +255,7 @@ def screen_output(self, ts): 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 )) + 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] @@ -243,24 +265,23 @@ def file_output(self, filename): i += 1 # 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:] + 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, i:i+3] = aero_tstep.total_unsteady_inertial_forces[:3] - moment_matrix[ts, i:i+3] = aero_tstep.total_unsteady_inertial_forces[3:] + 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, i:i+3] = aero_tstep.total_steady_body_forces[:3] - moment_matrix[ts, i:i+3] = aero_tstep.total_steady_body_forces[3:] + 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, 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, i:i + 3] = aero_tstep.total_unsteady_body_forces[:3] + moment_matrix[ts, i:i + 3] = aero_tstep.total_unsteady_body_forces[3:] header = '' header += 'tstep, ' 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):