Source code for sharpy.postproc.liftdistribution

import os

import numpy as np

from sharpy.utils.solver_interface import solver, BaseSolver
import sharpy.utils.settings as settings_utils
import sharpy.aero.utils.mapping as mapping
import sharpy.utils.algebra as algebra
import sharpy.aero.utils.utils as aeroutils


[docs]@solver class LiftDistribution(BaseSolver): """LiftDistribution Calculates and exports the lift distribution on lifting surfaces """ solver_id = 'LiftDistribution' solver_classification = 'post-processor' settings_types = dict() settings_default = dict() settings_description = dict() settings_types['text_file_name'] = 'str' settings_default['text_file_name'] = 'liftdistribution' settings_description['text_file_name'] = 'Text file name' settings_default['coefficients'] = True settings_types['coefficients'] = 'bool' settings_description['coefficients'] = 'Calculate aerodynamic lift coefficients' settings_types['rho'] = 'float' settings_default['rho'] = 1.225 settings_description['rho'] = 'Reference freestream density [kg/m³]' settings_table = settings_utils.SettingsTable() __doc__ += settings_table.generate(settings_types, settings_default, settings_description) def __init__(self): self.settings = None self.data = None self.folder = None self.caller = None def initialise(self, data, custom_settings=None, restart=False, caller=None): self.data = data self.settings = data.settings[self.solver_id] settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default) self.caller = caller self.folder = data.output_folder + '/liftdistribution/' if not os.path.exists(self.folder): os.makedirs(self.folder) def run(self, **kwargs): self.lift_distribution(self.data.structure.timestep_info[self.data.ts], self.data.aero.timestep_info[self.data.ts]) return self.data def lift_distribution(self, struct_tstep, aero_tstep): # Force mapping forces = mapping.aero2struct_force_mapping( aero_tstep.forces + aero_tstep.dynamic_forces, self.data.aero.struct2aero_mapping, aero_tstep.zeta, struct_tstep.pos, struct_tstep.psi, self.data.structure.node_master_elem, self.data.structure.connectivities, struct_tstep.cag(), self.data.aero.data_dict) # Prepare output matrix and file N_nodes = self.data.structure.num_node numb_col = 6 header = "x,y,z,fx,fy,fz" # get aero forces # get rotation matrix cga = algebra.quat2rotation(struct_tstep.quat) if self.settings["coefficients"]: # TODO: add nondimensional spanwise column y/s header += ", cfx, cfy, cfz" numb_col += 3 lift_distribution = np.zeros((N_nodes, numb_col)) for inode in range(N_nodes): if self.data.aero.data_dict['aero_node'][inode]: local_node = self.data.aero.struct2aero_mapping[inode][0]["i_n"] ielem, inode_in_elem = self.data.structure.node_master_elem[inode] i_surf = int(self.data.aero.surface_distribution[ielem]) # get c_gb cab = algebra.crv2rotation(struct_tstep.psi[ielem, inode_in_elem, :]) cgb = np.dot(cga, cab) # Get c_bs urel, dir_urel = aeroutils.magnitude_and_direction_of_relative_velocity(struct_tstep.pos[inode, :], struct_tstep.pos_dot[inode, :], struct_tstep.for_vel[:], cga, aero_tstep.u_ext[i_surf][:, :, local_node]) dir_span, span, dir_chord, chord = aeroutils.span_chord(local_node, aero_tstep.zeta[i_surf]) # Stability axes - projects forces in B onto S c_bs = aeroutils.local_stability_axes(cgb.T.dot(dir_urel), cgb.T.dot(dir_chord)) aero_forces = c_bs.T.dot(forces[inode, :3]) # Store data in export matrix lift_distribution[inode, 3:6] = aero_forces lift_distribution[inode, 2] = struct_tstep.pos[inode, 2] # z lift_distribution[inode, 1] = struct_tstep.pos[inode, 1] # y lift_distribution[inode, 0] = struct_tstep.pos[inode, 0] # x if self.settings["coefficients"]: # Get lift coefficient for idim in range(3): lift_distribution[inode, 6+idim] = np.sign(aero_forces[idim]) * np.linalg.norm(aero_forces[idim]) \ / (0.5 * self.settings['rho'] \ * np.linalg.norm(urel) ** 2 * span * chord) # Check if shared nodes from different surfaces exist (e.g. two wings joining at symmetry plane) # Leads to error since panel area just donates for half the panel size while lift forces is summed up lift_distribution[inode, 6+idim] /= len(self.data.aero.struct2aero_mapping[inode]) # Export lift distribution data np.savetxt(os.path.join(self.folder, self.settings['text_file_name'] + '_ts{}'.format(str(self.data.ts)) + '.txt'), lift_distribution, fmt='%10e,' * (numb_col - 1) + '%10e', delimiter=", ", header=header)