From 8d2166253cf4cce2d0c863918e0db1611f46132e Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Wed, 28 Apr 2021 19:23:37 +0200 Subject: [PATCH 01/19] wip_lin_polar --- sharpy/generators/polaraeroforces.py | 18 +++ sharpy/linear/assembler/linearaeroelastic.py | 134 +++++++++++++++++++ sharpy/linear/utils/derivatives.py | 2 + sharpy/postproc/stabilityderivatives.py | 2 + 4 files changed, 156 insertions(+) diff --git a/sharpy/generators/polaraeroforces.py b/sharpy/generators/polaraeroforces.py index cff836f65..0d948b7a6 100644 --- a/sharpy/generators/polaraeroforces.py +++ b/sharpy/generators/polaraeroforces.py @@ -154,6 +154,9 @@ def generate(self, **params): cga = algebra.quat2rotation(structural_kstep.quat) pos_g = np.array([cga.dot(structural_kstep.pos[inode]) + np.array([0, 0, 0]) for inode in range(nnode)]) + aero_kstep.polar_coeff = np.zeros((nnode, 6)) + aero_kstep.stability_transform = np.zeros((nnode, 3, 3)) + for inode in range(nnode): new_struct_forces[inode, :] = struct_forces[inode, :].copy() if aero_dict['aero_node'][inode]: @@ -266,6 +269,21 @@ def generate(self, **params): new_struct_forces[inode, 3:6] = c_bs.dot(moment_s) + linear_force_coeffs = forces_s / c_bs.T.dot(struct_forces[inode, :3]) + linear_moment_coeffs = moment_s / c_bs.T.dot(struct_forces[inode, 3:]) + + # print('New forces = ', forces_s / coef) + # print('Old forces = ', c_bs.T.dot(struct_forces[inode, :3]) / coef) + + print('Forces Coef', linear_force_coeffs) + print('Moments coef', linear_moment_coeffs) + aero_kstep.polar_coeff[inode] = np.concatenate((linear_force_coeffs, linear_moment_coeffs)) + orig_forces = np.concatenate((c_bs.T.dot(struct_forces[inode, :3]), c_bs.T.dot(struct_forces[inode, 3:]))) + for ind in range(6): + if np.abs(orig_forces[ind]) < 1e-4: + aero_kstep.polar_coeff[inode, ind] = 0 + aero_kstep.stability_transform[inode] = c_bs + return new_struct_forces diff --git a/sharpy/linear/assembler/linearaeroelastic.py b/sharpy/linear/assembler/linearaeroelastic.py index bd451dd5e..75359b434 100644 --- a/sharpy/linear/assembler/linearaeroelastic.py +++ b/sharpy/linear/assembler/linearaeroelastic.py @@ -131,6 +131,7 @@ def initialise(self, data): self.linearisation_vectors[k] = v self.get_gebm2uvlm_gains(data) + self.data = data def assemble(self): r""" @@ -227,6 +228,12 @@ def assemble(self): self.couplings['Ksa'] = gain_ksa self.couplings['Kas'] = gain_kas + polar_matrix_disp, polar_matrix_vel = self.polar_gains(self.data, uvlm.ss) + uvlm.ss.D *= 0 + uvlm.ss.C *= 0 + uvlm.ss.D[:, :polar_matrix_disp.shape[1]] = polar_matrix_disp + uvlm.ss.D[:, -polar_matrix_vel.shape[1]:] = polar_matrix_vel + if self.settings['beam_settings']['modal_projection'] is True and \ self.settings['beam_settings']['inout_coords'] == 'modes': # Project UVLM onto modal space @@ -996,3 +1003,130 @@ def load_uvlm(filename): uvlm_ss_read = read_data return libss.StateSpace(uvlm_ss_read.A, uvlm_ss_read.B, uvlm_ss_read.C, uvlm_ss_read.D, dt=uvlm_ss_read.dt) + + def polar_gains(self, data, ss): + + import sharpy.aero.utils.utils as aeroutils + aero = data.aero + structure = data.structure + tsaero = self.uvlm.tsaero0 + tsstr = self.beam.tsstruct0 + + Kzeta = self.uvlm.sys.Kzeta + num_dof_str = self.beam.sys.num_dof_str + num_dof_rig = self.beam.sys.num_dof_rig + num_dof_flex = self.beam.sys.num_dof_flex + use_euler = self.beam.sys.use_euler + + # allocate output + Kdisp = np.zeros((3 * Kzeta, num_dof_str)) + Kdisp_vel = np.zeros((3 * Kzeta, num_dof_str)) # Orientation is in velocity DOFs + Kvel_disp = np.zeros((3 * Kzeta, num_dof_str)) + Kvel_vel = np.zeros((3 * Kzeta, num_dof_str)) + Kforces = np.zeros((num_dof_str, 3 * Kzeta)) + + Kss = np.zeros((num_dof_flex, num_dof_flex)) + Csr = np.zeros((num_dof_flex, num_dof_rig)) + Crs = np.zeros((num_dof_rig, num_dof_flex)) + Crr = np.zeros((num_dof_rig, num_dof_rig)) + Krs = np.zeros((num_dof_rig, num_dof_flex)) + + # get projection matrix A->G + # (and other quantities indep. from nodal position) + Cga = algebra.quat2rotation(tsstr.quat) # NG 6-8-19 removing .T + Cag = Cga.T + + # for_pos=tsstr.for_pos + for_vel = tsstr.for_vel[:3] + for_rot = tsstr.for_vel[3:] + skew_for_rot = algebra.skew(for_rot) + Der_vel_Ra = np.dot(Cga, skew_for_rot) + + Faero = np.zeros(3) + FaeroA = np.zeros(3) + + # GEBM degrees of freedom + jj_for_tra = range(num_dof_str - num_dof_rig, + num_dof_str - num_dof_rig + 3) + jj_for_rot = range(num_dof_str - num_dof_rig + 3, + num_dof_str - num_dof_rig + 6) + + if use_euler: + jj_euler = range(num_dof_str - 3, num_dof_str) + euler = algebra.quat2euler(tsstr.quat) + tsstr.euler = euler + else: + jj_quat = range(num_dof_str - 4, num_dof_str) + + polar_matrix_vel = np.zeros((ss.outputs, num_dof_str)) + polar_matrix_disp = np.zeros((ss.outputs, num_dof_str)) + + jj = 0 # nodal dof index + for node_glob in range(structure.num_node): + + ### detect bc at node (and no. of dofs) + bc_here = structure.boundary_conditions[node_glob] + + if bc_here == 1: # clamp (only rigid-body) + dofs_here = 0 + jj_tra, jj_rot = [], [] + # continue + + elif bc_here == -1 or bc_here == 0: # (rigid+flex body) + dofs_here = 6 + jj_tra = 6 * structure.vdof[node_glob] + np.array([0, 1, 2], dtype=int) + jj_rot = 6 * structure.vdof[node_glob] + np.array([3, 4, 5], dtype=int) + else: + raise NameError('Invalid boundary condition (%d) at node %d!' \ + % (bc_here, node_glob)) + + jj += dofs_here + + # retrieve element and local index + ee, node_loc = structure.node_master_elem[node_glob, :] + + # get position, crv and rotation matrix + Ra = tsstr.pos[node_glob, :] # in A FoR, w.r.t. origin A-G + Rg = np.dot(Cag.T, Ra) # in G FoR, w.r.t. origin A-G + psi = tsstr.psi[ee, node_loc, :] + psi_dot = tsstr.psi_dot[ee, node_loc, :] + Cab = algebra.crv2rotation(psi) + Cba = Cab.T + Cbg = np.dot(Cab.T, Cag) + Tan = algebra.crv2tan(psi) + + track_body = self.settings['track_body'] + isurf = data.aero.struct2aero_mapping[node_glob][0]['i_surf'] + i_n = data.aero.struct2aero_mapping[node_glob][0]['i_n'] + + cas = Cab.dot(tsaero.stability_transform[node_glob]) + + dir_span, span, dir_chord, chord = aeroutils.span_chord(i_n, tsaero.zeta[isurf]) + + # Define the relative velocity and its direction + urel, dir_urel = aeroutils.magnitude_and_direction_of_relative_velocity(tsstr.pos[node_glob, :], + tsstr.pos_dot[node_glob, :], + tsstr.for_vel[:], + Cga, + tsaero.u_ext[isurf][:, :, i_n]) + coef = 0.5 * 1.225 * np.linalg.norm(urel) + cla = np.zeros((3, 3)) + cla[0, :] = 0.5 + cla[2, :] = 18 + u_rel_norm = np.linalg.norm(urel) + ldir = np.zeros((3, 3)) + ldir[:, 2] = 1 + if bc_here != 1: + polar_matrix_vel[np.ix_(jj_tra, jj_tra)] -= cas.dot(coef * cla * ldir).dot(cas.T) # eta_dot + polar_matrix_vel[np.ix_(jj_tra, jj_for_tra)] -= cas.dot(coef * cla * ldir).dot(cas.T) # va + polar_matrix_vel[np.ix_(jj_tra, jj_for_rot)] += cas.dot(coef * cla * ldir).dot(cas.T).dot(algebra.skew(Ra)) # eta * Ra + + polar_matrix_vel[np.ix_(jj_for_tra, jj_tra)] -= cas.dot(coef * cla * ldir).dot(cas.T) # eta_dot + polar_matrix_vel[np.ix_(jj_for_tra, jj_for_tra)] -= cas.dot(coef * cla * ldir).dot(cas.T) # va + polar_matrix_vel[np.ix_(jj_for_tra, jj_for_rot)] += cas.dot(coef * cla * ldir).dot(cas.T).dot(algebra.skew(Ra)) # eta * Ra + # polar_matrix_disp[np.ix_(jj_for_tra, jj_euler)] += cas.dot(coef * cla * ldir).dot(cas.T) + + + + import pdb; pdb.set_trace() + return polar_matrix_disp, polar_matrix_vel diff --git a/sharpy/linear/utils/derivatives.py b/sharpy/linear/utils/derivatives.py index 812835ffe..ceecc34dc 100644 --- a/sharpy/linear/utils/derivatives.py +++ b/sharpy/linear/utils/derivatives.py @@ -349,6 +349,8 @@ def angle_derivatives(self): stab_der_trans2 = cga.dot(delta_nodal_forces[:3, :]) stab_der_mom2 = cga.dot(delta_nodal_forces[3:, :]) + import pdb; pdb.set_trace() + self.matrix[:3, :] = stab_der_trans + stab_der_trans2 self.matrix[3:6, :] = stab_der_mom + stab_der_mom2 diff --git a/sharpy/postproc/stabilityderivatives.py b/sharpy/postproc/stabilityderivatives.py index 82942c7df..ca90a21f4 100644 --- a/sharpy/postproc/stabilityderivatives.py +++ b/sharpy/postproc/stabilityderivatives.py @@ -149,6 +149,7 @@ def run(self, online=False): self.data.linear.linear_system.update(self.settings['u_inf']) for target_system in ['aerodynamic', 'aeroelastic']: + print('-------- CURRENT SYSTEM ---------', target_system) state_space = self.get_state_space(target_system) current_derivative = derivatives[target_system] @@ -364,6 +365,7 @@ def angle_derivatives(self, H0): # second term in the stability derivative expression stab_der_trans2 = cga.dot(H0[:3, :3].real.dot(cga.T.dot(algebra.der_Peuler_by_v(euler0 * 0, v0)))) + import pdb; pdb.set_trace() stab_der_mom = algebra.der_Ceuler_by_v(euler0, m0a) stab_der_mom2 = cga.dot(H0[3:6, :3].real.dot(cga.T.dot(algebra.der_Peuler_by_v(euler0 * 0, v0)))) From ab55199acd39b9c7ee095dd75daa268e63967348 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Sun, 9 May 2021 12:28:07 +0200 Subject: [PATCH 02/19] Generate linear gain on airfoil polar coefficients --- sharpy/aero/utils/airfoilpolars.py | 20 +++ sharpy/generators/polaraeroforces.py | 141 ++++++++++++++++--- sharpy/linear/assembler/linearaeroelastic.py | 34 ++++- 3 files changed, 169 insertions(+), 26 deletions(-) diff --git a/sharpy/aero/utils/airfoilpolars.py b/sharpy/aero/utils/airfoilpolars.py index 64b6b1bb0..8b8227e08 100644 --- a/sharpy/aero/utils/airfoilpolars.py +++ b/sharpy/aero/utils/airfoilpolars.py @@ -50,6 +50,26 @@ def get_coefs(self, aoa_deg): return cl, cd, cm + def get_derivatives_at_aoa(self, aoa): + """ + Get Derivatives at angle of attack + + Args: + aoa (float): Angle of attack in radians + + Returns: + tuple: Containing ``cl_a``, ``cd_a``, ``cm_a`` + """ + + delta_aoa = 1e-5 + p1 = np.array(self.get_coefs(aoa + delta_aoa)) + p = np.array(self.get_coefs(aoa)) + m1 = np.array(self.get_coefs(aoa - delta_aoa)) + + der = (3 * p1 - 4 * p + m1) / 2 / delta_aoa + + return der + def get_aoa_deg_from_cl_2pi(self, cl): return cl/2/np.pi/deg2rad + self.aoa_cl0_deg diff --git a/sharpy/generators/polaraeroforces.py b/sharpy/generators/polaraeroforces.py index 1783d6267..00ffb261e 100644 --- a/sharpy/generators/polaraeroforces.py +++ b/sharpy/generators/polaraeroforces.py @@ -128,6 +128,8 @@ def generate(self, **params): aero_kstep.polar_coeff = np.zeros((nnode, 6)) aero_kstep.stability_transform = np.zeros((nnode, 3, 3)) + structural_kstep.postproc_node['aoa'] = np.zeros(nnode) + for inode in range(nnode): new_struct_forces[inode, :] = struct_forces[inode, :].copy() if aero_dict['aero_node'][inode]: @@ -165,18 +167,19 @@ def generate(self, **params): cl = np.sign(lift_force) * np.linalg.norm(lift_force) / coef cd_sharpy = np.linalg.norm(drag_force) / coef + # ii) Compute the effective angle of attack from potential flow theory. The local lift curve + # slope is 2pi and the zero-lift angle of attack is given by thin airfoil theory. From this, + # the effective angle of attack is computed for the section and includes 3D effects. + aoa_0cl = get_aoacl0_from_camber(airfoil_coords[:, 0], airfoil_coords[:, 1]) + aoa = cl / 2 / np.pi + aoa_0cl + structural_kstep.postproc_node['aoa'][inode] = aoa + if cd_from_cl: # Compute the drag from the UVLM computed lift cd, cm = polar.get_cdcm_from_cl(cl) else: - # Compute L, D, M from polar depending on: - # ii) Compute the effective angle of attack from potential flow theory. The local lift curve - # slope is 2pi and the zero-lift angle of attack is given by thin airfoil theory. From this, - # the effective angle of attack is computed for the section and includes 3D effects. - aoa_0cl = get_aoacl0_from_camber(airfoil_coords[:, 0], airfoil_coords[:, 1]) - aoa = cl / 2 / np.pi + aoa_0cl - # Compute the coefficients associated to that angle of attack + # Compute L, D, M from polar depending on the coefficients associated to that angle of attack cl_polar, cd, cm = polar.get_coefs(aoa) if correct_lift: @@ -213,22 +216,117 @@ def generate(self, **params): new_struct_forces[inode, 3:6] = c_bs.dot(moment_s) - linear_force_coeffs = forces_s / c_bs.T.dot(struct_forces[inode, :3]) - linear_moment_coeffs = moment_s / c_bs.T.dot(struct_forces[inode, 3:]) + return new_struct_forces - # print('New forces = ', forces_s / coef) - # print('Old forces = ', c_bs.T.dot(struct_forces[inode, :3]) / coef) + def generate_linear(self, **params): + """ + Generate Linear Gain matrix with correction factors (lift only at the moment) - print('Forces Coef', linear_force_coeffs) - print('Moments coef', linear_moment_coeffs) - aero_kstep.polar_coeff[inode] = np.concatenate((linear_force_coeffs, linear_moment_coeffs)) - orig_forces = np.concatenate((c_bs.T.dot(struct_forces[inode, :3]), c_bs.T.dot(struct_forces[inode, 3:]))) - for ind in range(6): - if np.abs(orig_forces[ind]) < 1e-4: - aero_kstep.polar_coeff[inode, ind] = 0 - aero_kstep.stability_transform[inode] = c_bs + Keyword Args: + beam (sharpy.linear.assembler.linearbeam.LinearBeam): Beam object + tsstruct0 (sharpy.utils.datastructures.StructTimestepInfo): Ref structural time step + tsaero0 (sharpy.utils.datastructures.AeroTimestepInfo): Ref aero time step info - return new_struct_forces + Returns: + np.array: Gain matrix + """ + beam = params['beam'] + tsstruct0 = params['tsstruct0'] + tsaero0 = params['tsaero0'] + + aerogrid = self.aero + structure = self.structure + + structural_kstep = tsstruct0 + aero_kstep = tsaero0 + nnode = tsstruct0.pos.shape[0] + + aero_dict = aerogrid.aero_dict + if aerogrid.polars is None: + return None + + cga = algebra.quat2rotation(structural_kstep.quat) + + aero_kstep.polar_coeff = np.zeros((nnode, 6)) + aero_kstep.stability_transform = np.zeros((nnode, 3, 3)) + + # Matrix allocation + num_dof_str = beam.sys.num_dof_str + num_dof_rig = beam.sys.num_dof_rig + use_euler = beam.sys.use_euler + + # GEBM degrees of freedom + jj_for_tra = range(num_dof_str - num_dof_rig, + num_dof_str - num_dof_rig + 3) + jj_for_rot = range(num_dof_str - num_dof_rig + 3, + num_dof_str - num_dof_rig + 6) + + if use_euler: + jj_euler = range(num_dof_str - 3, num_dof_str) + euler = algebra.quat2euler(tsstruct0.quat) + tsstruct0.euler = euler + else: + jj_quat = range(num_dof_str - 4, num_dof_str) + + polar_gain = np.zeros((num_dof_str, num_dof_str)) + + jj = 0 # Global DOF index + for inode in range(nnode): + if aero_dict['aero_node'][inode]: + + ### detect bc at node (and no. of dofs) + bc_here = structure.boundary_conditions[inode] + + if bc_here == 1: # clamp (only rigid-body) + dofs_here = 0 + jj_tra, jj_rot = [], [] + + elif bc_here == -1 or bc_here == 0: # (rigid+flex body) + dofs_here = 6 + jj_tra = 6 * structure.vdof[inode] + np.array([0, 1, 2], dtype=int) + jj_rot = 6 * structure.vdof[inode] + np.array([3, 4, 5], dtype=int) + else: + raise NameError('Invalid boundary condition (%d) at node %d!' \ + % (bc_here, inode)) + + jj += dofs_here + + ielem, inode_in_elem = structure.node_master_elem[inode] + iairfoil = aero_dict['airfoil_distribution'][ielem, inode_in_elem] + isurf = aerogrid.struct2aero_mapping[inode][0]['i_surf'] + i_n = aerogrid.struct2aero_mapping[inode][0]['i_n'] + N = aerogrid.aero_dimensions[isurf, 1] + polar = aerogrid.polars[iairfoil] + cab = algebra.crv2rotation(structural_kstep.psi[ielem, inode_in_elem, :]) + cgb = np.dot(cga, cab) + + dir_span, span, dir_chord, chord = span_chord(i_n, aero_kstep.zeta[isurf]) + + # Define the relative velocity and its direction + urel, dir_urel = magnitude_and_direction_of_relative_velocity(structural_kstep.pos[inode, :], + structural_kstep.pos_dot[inode, :], + structural_kstep.for_vel[:], + cga, + aero_kstep.u_ext[isurf][:, :, i_n]) + + # Stability axes - projects forces in B onto S + c_bs = local_stability_axes(cgb.T.dot(dir_urel), cgb.T.dot(dir_chord)) + cas = cab.dot(c_bs) + + local_correction = np.eye(6) + + cla, cda, cma = polar.get_derivatives_at_aoa(tsstruct0.postproc_node['aoa'][inode]) + local_correction[2, 2] = cla / 2 / np.pi + + weight = 1 / np.sum(aero_dict['aero_node']) + if bc_here != 1: + polar_gain[np.ix_(jj_tra, jj_tra)] += cas.dot(local_correction[:3, :3]).dot(cas.T) + polar_gain[np.ix_(jj_rot, jj_rot)] += cas.dot(local_correction[3:, 3:]).dot(cas.T) + + polar_gain[np.ix_(jj_for_tra, jj_for_tra)] += weight * cas.dot(local_correction[:3, :3]).dot(cas.T) + polar_gain[np.ix_(jj_for_rot, jj_for_rot)] += weight * cas.dot(local_correction[3:, 3:]).dot(cas.T) + + return polar_gain @generator_interface.generator @@ -301,3 +399,6 @@ def generate(self, **params): new_struct_forces[inode, 3:6] *= moment_efficiency[i_elem, i_local_node, 0, :] new_struct_forces[inode, 3:6] += moment_efficiency[i_elem, i_local_node, 1, :] return new_struct_forces + + def generate_linear(self, **params): + raise NotImplementedError('Efficiency factors not yet implemented') diff --git a/sharpy/linear/assembler/linearaeroelastic.py b/sharpy/linear/assembler/linearaeroelastic.py index 75359b434..bd2bf9d95 100644 --- a/sharpy/linear/assembler/linearaeroelastic.py +++ b/sharpy/linear/assembler/linearaeroelastic.py @@ -88,6 +88,9 @@ def __init__(self): self.Crs = None self.Crr = None + self.correct_forces = False + self.correct_forces_generator = None + def initialise(self, data): try: @@ -131,7 +134,22 @@ def initialise(self, data): self.linearisation_vectors[k] = v self.get_gebm2uvlm_gains(data) - self.data = data + + # correct forces generators + try: + data.settings['StaticCoupled']['correct_forces_method'] + except KeyError: + self.correct_forces = False + else: + if data.settings['StaticCoupled']['correct_forces_method'] is not '': + import sharpy.utils.generator_interface as gen_interface + self.correct_forces = True + self.correct_forces_generator = gen_interface.generator_from_string(data.settings['StaticCoupled']['correct_forces_method'])() + self.correct_forces_generator.initialise(in_dict=data.settings['StaticCoupled']['correct_forces_settings'], + aero=data.aero, + structure=data.structure, + rho=self.settings['aero_settings']['density'], + vortex_radius=self.settings['aero_settings']['vortex_radius']) def assemble(self): r""" @@ -228,11 +246,15 @@ def assemble(self): self.couplings['Ksa'] = gain_ksa self.couplings['Kas'] = gain_kas - polar_matrix_disp, polar_matrix_vel = self.polar_gains(self.data, uvlm.ss) - uvlm.ss.D *= 0 - uvlm.ss.C *= 0 - uvlm.ss.D[:, :polar_matrix_disp.shape[1]] = polar_matrix_disp - uvlm.ss.D[:, -polar_matrix_vel.shape[1]:] = polar_matrix_vel + if self.correct_forces: + polar_gain_value = self.correct_forces_generator.generate_linear(beam=beam, + tsstruct0=beam.sys.tsstruct0, + tsaero0=uvlm.tsaero0) + polar_gain = libss.Gain(polar_gain_value, + input_vars=LinearVector.transform(uvlm.ss.output_variables, + to_type=InputVariable), + output_vars=uvlm.ss.output_variables.copy()) + uvlm.ss.addGain(polar_gain, where='out') if self.settings['beam_settings']['modal_projection'] is True and \ self.settings['beam_settings']['inout_coords'] == 'modes': From 101b41bbefe6bfec751388278d8c6f1712843905 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Sun, 9 May 2021 12:28:44 +0200 Subject: [PATCH 03/19] Clean up and remove unused code --- sharpy/linear/utils/derivatives.py | 69 ---- sharpy/postproc/stabilityderivatives.py | 425 +----------------------- 2 files changed, 4 insertions(+), 490 deletions(-) diff --git a/sharpy/linear/utils/derivatives.py b/sharpy/linear/utils/derivatives.py index ceecc34dc..ee7ace461 100644 --- a/sharpy/linear/utils/derivatives.py +++ b/sharpy/linear/utils/derivatives.py @@ -228,70 +228,6 @@ def save(self, derivative_name, output_name): with h5py.File(output_name + '.stability.h5', 'w') as f: f.create_dataset(derivative_name, data=self.matrix) - # @classmethod - # def initialise_derivatives(cls, state_space, steady_forces, quat, v0, phi=None): - # """ - # Initialises the required class attributes for all derivative calculations/ - # - # Args: - # state_space (sharpy.linear.src.libss.StateSpace): State-space object for the target system - # steady_forces (np.array): Array of steady forces (at the linearisation) expressed in the beam degrees of - # freedom and with size equal to the number of structural degrees of freedom - # quat (np.array): Quaternion at the linearisation - # v0 (np.array): Free stream velocity vector at the linearisation condition - # phi (np.array (optional)): Mode shape matrix for modal systems - # - # """ - # cls.quat = quat - # cls.cga = algebra.quat2rotation(cls.quat) - # cls.v0 = v0 - # - # if phi is not None: - # cls.modal = True - # cls.phi = phi[-9:-3, :6] - # cls.inv_phi_forces = np.linalg.inv(phi[-9:-3, :6].T) - # cls.inv_phi_vel = np.linalg.inv(phi[-9:-3, :6]) - # else: - # cls.modal = False - # cls.steady_forces = steady_forces - # - # cls.transfer_function = cls.calculate_transfer_function(state_space) - # - # @classmethod - # def calculate_transfer_function(cls, ss): - # H0 = ss.freqresp(np.array([1e-3]))[:, :, 0] - # # A, B, C, D = ss.get_mats() - # # H0 = C.dot(np.linalg.inv(np.eye(ss.states) - A).dot(B)) + D - # # np.savetxt('./nodal_aeroelastic_static_manual.txt', H0.real) - # if cls.modal: - # vel_inputs_variables = ss.input_variables.get_variable_from_name('q_dot') - # output_indices = ss.output_variables.get_variable_from_name('Q').rows_loc[:6] - # cls.steady_forces = cls.inv_phi_forces.dot(cls.steady_forces[output_indices]) - # else: - # vel_inputs_variables = ss.input_variables.get_variable_from_name('beta') - # output_indices = ss.output_variables.get_variable_from_name('forces_n').rows_loc[-9:-3] - # cls.steady_forces = cls.steady_forces[output_indices] - # rbm_indices = vel_inputs_variables.cols_loc[:9] - # - # # look for control surfaces - # try: - # cs_input_variables = ss.input_variables.get_variable_from_name('control_surface_deflection') - # dot_cs_input_variables = ss.input_variables.get_variable_from_name('dot_control_surface_deflection') - # except ValueError: - # cs_indices = np.array([], dtype=int) - # dot_cs_indices = np.array([], dtype=int) - # cls.n_control_surfaces = 0 - # else: - # cs_indices = cs_input_variables.cols_loc - # dot_cs_indices = dot_cs_input_variables.cols_loc - # cls.n_control_surfaces = cs_input_variables.size - # finally: - # input_indices = np.concatenate((rbm_indices, cs_indices, dot_cs_indices)) - # - # H0 = H0[np.ix_(output_indices, input_indices)].real - # - # return H0 - def angle_derivatives(self): r""" Stability derivatives against aerodynamic angles (angle of attack and sideslip) expressed in stability axes, i.e @@ -349,8 +285,6 @@ def angle_derivatives(self): stab_der_trans2 = cga.dot(delta_nodal_forces[:3, :]) stab_der_mom2 = cga.dot(delta_nodal_forces[3:, :]) - import pdb; pdb.set_trace() - self.matrix[:3, :] = stab_der_trans + stab_der_trans2 self.matrix[3:6, :] = stab_der_mom + stab_der_mom2 @@ -368,9 +302,6 @@ def angle_derivatives_tb(self): f0a = self.steady_forces[:3] m0a = self.steady_forces[-3:] - print(f0a) - print(m0a) - euler0 = algebra.quat2euler(self.quat) cga = self.cga diff --git a/sharpy/postproc/stabilityderivatives.py b/sharpy/postproc/stabilityderivatives.py index ca90a21f4..11b15f594 100644 --- a/sharpy/postproc/stabilityderivatives.py +++ b/sharpy/postproc/stabilityderivatives.py @@ -115,26 +115,9 @@ def initialise(self, data, custom_settings=None, caller=None): def run(self, online=False): # TODO: consider running all required solvers inside this one to keep the correct settings - # i.e: run Modal, Linear Ass + # i.e: run Modal, Linear Ass... COMING SOON with SHARPy routines derivatives = self.data.linear.derivatives - # Y_freq = self.uvlm_steady_state_transfer_function() - # angle_ders = self.angle_derivatives(Y_freq) - # derivatives.dict_of_derivatives['force_angle_velocity'] = angle_ders[0] - # derivatives.dict_of_derivatives['force_angle'] = angle_ders[1] - # derivatives.dict_of_derivatives['force_angle_body'] = angle_ders[2] - # derivatives.dict_of_derivatives['force_velocity'] = self.body_derivatives(Y_freq) - # derivatives.dict_of_derivatives['force_cs_body'] = self.control_surface_derivatives(Y_freq) - # derivatives.save(self.settings['folder']) - - # derivatives_dimensional, derivatives_coeff = self.derivatives(Y_freq) - - # self.export_derivatives(np.hstack((derivatives_coeff[:, :6], derivatives_coeff[:, -2:]))) - # self.data.linear.derivatives.savetxt(self.folder) - # angle_derivative_set.print(derivative_filename=self.folder + '/stability_derivatives.txt') - - - # >>>>>>>> New methods if self.data.linear.linear_system.beam.sys.modal: phi = self.data.linear.linear_system.linearisation_vectors['mode_shapes'].real else: @@ -149,7 +132,7 @@ def run(self, online=False): self.data.linear.linear_system.update(self.settings['u_inf']) for target_system in ['aerodynamic', 'aeroelastic']: - print('-------- CURRENT SYSTEM ---------', target_system) + cout.cout_wrap('-------- {:s} SYSTEM DERIVATIVES ---------'.format(target_system.upper())) state_space = self.get_state_space(target_system) current_derivative = derivatives[target_system] @@ -277,409 +260,9 @@ def uvlm_steady_state_transfer_function(self): return H0 - def angle_derivatives(self, H0): - r""" - Stability derivatives against aerodynamic angles (angle of attack and sideslip) expressed in stability axes, i.e - forces are lift, drag... - - Linearised forces in stability axes are expressed as - - .. math:: - F^S = F_0^S + \frac{\partial}{\partial \alpha}\left(C^{GA}(\alpha)F_0^A\right)\delta\alpha + C_0^{GA}\delta F^A - - Therefore, the stability derivative becomes - - .. math:: \frac{\partial\F^S}{\partial\alpha} =\frac{\partial}{\partial \alpha}\left(C^{GA}(\alpha)F_0^A\right) + - C_0^{GA}\frac{\partial F^A}{\partial\alpha} - - where - - .. math:: \frac{\partial F^A}{\partial\alpha} = \frac{\partial F^A}{\partial v^A}\frac{\partial v^A}{\partial\alpha} - - and - - .. math:: \frac{\partial v^A}{\partial\alpha} = C^{AG}\frac{\partial}{\partial\alpha}\left(C(0)V_0^G\right). - - The term :math:`\frac{\partial F^A}{\partial v^A}` is obtained directly from the steady state transfer - function of the linear UVLM expressed in the beam degrees of freedoms. - - Args: - H0 (np.ndarray): Steady state gain transfer function of the linear UVLM expressed in the beam rigid - degrees of freedom - - Returns: - sharpy.linear.utils.derivatives.DerivativeSet: containing the derivatives. - """ - derivative_set = DerivativeSet('stability') - derivative_set.labels_in = ['phi', 'alpha', 'beta'] - derivative_set.labels_out = ['CD', 'CY', 'CL', 'Cl', 'Cm', 'Cn'] - derivative_set.matrix = np.zeros((6, 3)) - - modal = self.data.linear.linear_system.beam.sys.modal - phi = self.data.linear.linear_system.linearisation_vectors['mode_shapes'].real - - # Get free stream velocity direction - v0 = self.get_freestream_velocity() - - # Steady Forces - # in beam dof - f0a = self.data.linear.linear_system.linearisation_vectors['forces_aero_beam_dof'][:3].copy().real - m0a = self.data.linear.linear_system.linearisation_vectors['forces_aero_beam_dof'][3:6].copy().real - - print(f0a) - print(m0a) - - if self.ppal_axes: - cap = self.data.linear.tsstruct0.modal['t_pa'] # rotates A onto P, thus projects Vp to Va - rap = self.data.linear.tsstruct0.modal['r_pa'] - # f0a = cap.dot(f0a)/ - # m0a = cap.dot(m0a) - # mod_forces = rap.dot(np.concatenate((f0a, m0a))) - # f0a = mod_forces[:3] - # m0a = mod_forces[3:] - - print('After PA') - print(f0a) - print(m0a) - - if modal: - forces_vec = np.concatenate((f0a, m0a)) # modal forces Q - phirr = phi[-9:-3, :6] - inv_phirr = np.linalg.inv(phirr.T) - dim_forces = inv_phirr.dot(forces_vec) - f0a = dim_forces[:3] - m0a = dim_forces[-3:] - phi_vel = phi[-9:-6, :3] - inv_phi_vel = np.linalg.inv(phi_vel) - - print('After Modal') - print(f0a) - print(m0a) - # breakpoint() - - euler0 = self.data.linear.tsstruct0.euler_angles() - cga = self.data.linear.tsstruct0.cga() - - # first term in the stability derivative expression - stab_der_trans = algebra.der_Ceuler_by_v(euler0, f0a) - # second term in the stability derivative expression - stab_der_trans2 = cga.dot(H0[:3, :3].real.dot(cga.T.dot(algebra.der_Peuler_by_v(euler0 * 0, v0)))) - - import pdb; pdb.set_trace() - stab_der_mom = algebra.der_Ceuler_by_v(euler0, m0a) - stab_der_mom2 = cga.dot(H0[3:6, :3].real.dot(cga.T.dot(algebra.der_Peuler_by_v(euler0 * 0, v0)))) - - if modal: - if self.ppal_axes: - delta_nodal_vel = inv_phi_vel.dot(cga.T.dot(algebra.der_Peuler_by_v(euler0 * 0, v0))) - delta_nodal_forces = inv_phirr.dot(H0[:6, :3].real.dot(delta_nodal_vel)) # this is in A - - stab_der_trans2 = cga.dot(delta_nodal_forces[:3, :]) - stab_der_mom2 = cga.dot(delta_nodal_forces[3:, :]) - else: - delta_nodal_vel = inv_phi_vel.dot(cga.T.dot(algebra.der_Peuler_by_v(euler0 * 0, v0))) - delta_nodal_forces = inv_phirr.dot(H0[:6, :3].real.dot(delta_nodal_vel)) - - stab_der_trans2 = cga.dot(delta_nodal_forces[:3, :]) - stab_der_mom2 = cga.dot(delta_nodal_forces[3:, :]) - - derivative_set.matrix[:3, :] = (stab_der_trans + stab_der_trans2) / self.coefficients['force'] - derivative_set.matrix[3:6, :] = (stab_der_mom + stab_der_mom2) - derivative_set.matrix[np.ix_([3, 5]), :] /= self.coefficients['moment_lat'] - derivative_set.matrix[4, :] /= self.coefficients['moment_lon'] - - ######################################## - - # for debugging and checking purposes at the moment - # derivative_set.print() - derivative_set.name = 'Force derivatives to angles - velocity perturbation' - derivative_set.save('force_angle', self.settings['folder'] + '/force_angle') - - angle_derivative_set = DerivativeSet('stability') - angle_derivative_set.name = 'Force derivatives to angles - body perturbation' - angle_derivative_set.labels_in = ['phi', 'alpha', 'beta'] - angle_derivative_set.labels_out = ['CD', 'CY', 'CL', 'Cl', 'Cm', 'Cn'] - # These are onto the original stability axes at the linearisation - # The above take the stability axes to rotate with the perturbation!! - angles = stab_der_trans / self.coefficients['force'] - delta_nodal_forces_a = inv_phirr.dot(H0[:6, 6:9]) - # angles += cga.dot(inv_phi_vel.dot(H0[:3, 6:9])) / self.coefficients['force'] - angles += cga.dot(delta_nodal_forces_a[:3, :]) / self.coefficients['force'] - mom_angles = stab_der_mom - # mom_angles += cga.dot(inv_phirr.dot(H0[:6, 6:9])[3:]) - mom_angles += cga.dot(delta_nodal_forces_a[3:, :]) - mom_angles[np.ix_([0, 2]), :] /= self.coefficients['moment_lat'] - mom_angles[1, :] /= self.coefficients['moment_lon'] - angle_derivative_set.matrix = np.vstack((angles, mom_angles)) - - ############ - cout.cout_wrap('Body axes') - angle_derivative_body = DerivativeSet('body') - angle_derivative_set.name = 'Force derivatives to angles - body perturbation' - angle_derivative_body.labels_in = ['phi', 'alpha', 'beta'] - angle_derivative_body.labels_out = ['C_XA', 'C_YA', 'C_ZA', 'C_LA', 'C_MA', 'C_NA'] - # These are onto the original stability axes at the linearisation - # The above take the stability axes to rotate with the perturbation!! - # angles = H0[:3, 6:9] / phi[-9, 0] / self.coefficients['force'] - angles = inv_phirr.dot(H0[:6, 6:9])[:3, :] / self.coefficients['force'] - mom_angles = inv_phirr.dot(H0[:6, 6:9])[3:, :] - mom_angles[np.ix_([0, 2]), :] /= self.coefficients['moment_lat'] - mom_angles[1, :] /= self.coefficients['moment_lon'] - angle_derivative_body.matrix = np.vstack((angles, mom_angles)) - # angle_derivative_body.print() - - return derivative_set, angle_derivative_set, angle_derivative_body - - def body_derivatives(self, H0): - derivative_set = DerivativeSet('body') - derivative_set.name = 'Force derivatives to rigid body velocities - Body derivatives' - derivative_set.labels_in = ['uA', 'vA', 'wA', 'pA', 'qA', 'rA'] - derivative_set.labels_out = ['C_XA', 'C_YA', 'C_ZA', 'C_LA', 'C_MA', 'C_NA'] - derivative_set.matrix = np.zeros((6, 6)) - - modal = self.data.linear.linear_system.beam.sys.modal - - body_derivatives = H0[:6, :6] - - if modal: - phi = self.data.linear.linear_system.linearisation_vectors['mode_shapes'].real - phirr = phi[-9:-3, :6] - inv_phirr = np.linalg.inv(phirr.T) - - body_derivatives = inv_phirr.dot(body_derivatives).dot(np.linalg.inv(phirr)) - - derivative_set.matrix = body_derivatives - derivative_set.matrix[:3, :] /= self.coefficients['force'] - derivative_set.matrix[np.ix_([3, 5]), :] /= self.coefficients['moment_lat'] - derivative_set.matrix[4, :] /= self.coefficients['moment_lon'] - # derivative_set.print() - - return derivative_set - - def control_surface_derivatives(self, H0): - derivative_set = DerivativeSet('body') - derivative_set.name = 'Force derivatives wrt control surface inputs - Body axes' - derivative_set.labels_out = ['C_XA', 'C_YA', 'C_ZA', 'C_LA', 'C_MA', 'C_NA'] - n_control_surfaces = self.n_control_surfaces - if n_control_surfaces == 0: - return None - labels_in_deflection = [] - labels_in_rate = [] - for i in range(n_control_surfaces): - labels_in_deflection.append('delta_{:g}'.format(i)) - labels_in_rate.append('dot(delta)_{:g}'.format(i)) - derivative_set.labels_in = labels_in_deflection + labels_in_rate - - body_derivatives = H0[:6, 9:] - assert body_derivatives.shape == (6, 2 * self.n_control_surfaces), 'Incorrect TF shape' - modal = self.data.linear.linear_system.beam.sys.modal - - if modal: - phi = self.data.linear.linear_system.linearisation_vectors['mode_shapes'].real - phirr = phi[-9:-3, :6] - inv_phirr = np.linalg.inv(phirr.T) - - body_derivatives = inv_phirr.dot(body_derivatives) - - derivative_set.matrix = body_derivatives - derivative_set.matrix[:3, :] /= self.coefficients['force'] - derivative_set.matrix[np.ix_([3, 5]), :] /= self.coefficients['moment_lat'] - derivative_set.matrix[4, :] /= self.coefficients['moment_lon'] - # derivative_set.print() - - return derivative_set - - # # Get rigid body + control surface inputs - # try: - # n_ctrl_sfc = self.data.linear.linear_system.uvlm.control_surface.n_control_surfaces - # except AttributeError: - # n_ctrl_sfc = 0 - # - # self.inputs = rig_dof + n_ctrl_sfc - # - # in_matrix = np.zeros((ssuvlm.inputs, self.inputs)) - # out_matrix = np.zeros((nout, ssuvlm.outputs)) - # - # if modal: - # # Modal scaling - # raise NotImplementedError('Not yet implemented in modal space') - # else: - # in_matrix[-self.inputs:, :] = np.eye(self.inputs) - # out_matrix[:, -rig_dof:-rig_dof+6] = np.eye(nout) - # - # ssuvlm.addGain(in_matrix, where='in') - # ssuvlm.addGain(out_matrix, where='out') - # - # A, B, C, D = ssuvlm.get_mats() - # if type(A) == libsp.csc_matrix: - # Y_freq = C.dot(scsp.linalg.inv(scsp.eye(ssuvlm.states, format='csc') - A).dot(B)) + D - # else: - # Y_freq = C.dot(np.linalg.inv(np.eye(ssuvlm.states) - A).dot(B)) + D - # Yf = ssuvlm.freqresp(np.array([0])) - # - # return Y_freq - - def a_derivatives(self, Y_freq): - - Cng = np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]]) # Project SEU on NED - TODO implementation - u_inf = self.settings['u_inf'].value - s_ref = self.settings['S_ref'].value - b_ref = self.settings['b_ref'].value - c_ref = self.settings['c_ref'].value - rho = self.data.linear.tsaero0.rho - - # Inertial frame - try: - euler = self.data.linear.tsstruct0.euler - Pga = algebra.euler2rot(euler) - rig_dof = 9 - except AttributeError: - quat = self.data.linear.tsstruct0.quat - Pga = algebra.quat2rotation(quat) - rig_dof = 10 - - derivatives_g = np.zeros((6, Y_freq.shape[1] + 2)) - coefficients = {'force': 0.5*rho*u_inf**2*s_ref, - 'moment_lon': 0.5*rho*u_inf**2*s_ref*c_ref, - 'moment_lat': 0.5*rho*u_inf**2*s_ref*b_ref, - 'force_angular_vel': 0.5*rho*u_inf**2*s_ref*c_ref/u_inf, - 'moment_lon_angular_vel': 0.5*rho*u_inf**2*s_ref*c_ref*c_ref/u_inf} # missing rates - - for in_channel in range(Y_freq.shape[1]): - derivatives_g[:3, in_channel] = Pga.dot(Y_freq[:3, in_channel]) - derivatives_g[3:, in_channel] = Pga.dot(Y_freq[3:, in_channel]) - - derivatives_g[:3, :3] /= coefficients['force'] - derivatives_g[:3, 3:6] /= coefficients['force_angular_vel'] - derivatives_g[4, :3] /= coefficients['moment_lon'] - derivatives_g[4, 3:6] /= coefficients['moment_lon_angular_vel'] - derivatives_g[[3, 5], :] /= coefficients['moment_lat'] - - derivatives_g[:, -2] = derivatives_g[:, 2] * u_inf # ders wrt alpha - derivatives_g[:, -1] = derivatives_g[:, 1] * u_inf # ders wrt beta - - der_matrix = np.zeros((6, self.inputs - (rig_dof - 6))) - der_col = 0 - for i in list(range(6))+list(range(rig_dof, self.inputs)): - der_matrix[:3, der_col] = Y_freq[:3, i] - der_matrix[3:6, der_col] = Y_freq[3:6, i] - der_col += 1 - - labels_force = {0: 'X', - 1: 'Y', - 2: 'Z', - 3: 'L', - 4: 'M', - 5: 'N'} - - labels_velocity = {0: 'u', - 1: 'v', - 2: 'w', - 3: 'p', - 4: 'q', - 5: 'r', - 6: 'flap1', - 7: 'flap2', - 8: 'flap3'} - - table = cout.TablePrinter(n_fields=7, field_length=12, field_types=['s', 'f', 'f', 'f', 'f', 'f', 'f']) - table.print_header(['der'] + list(labels_force.values())) - for i in range(der_matrix.shape[1]): - table.print_line([labels_velocity[i]] + list(der_matrix[:, i])) - - table_coeff = cout.TablePrinter(n_fields=7, field_length=12, field_types=['s']+6*['f']) - labels_out = {0: 'C_D', - 1: 'C_Y', - 2: 'C_L', - 3: 'C_l', - 4: 'C_m', - 5: 'C_n'} - labels_der = {0: 'u', - 1: 'v', - 2: 'w', - 3: 'p', - 4: 'q', - 5: 'r', - 6: 'alpha', - 7: 'beta'} - table_coeff.print_header(['der'] + list(labels_out.values())) - for i in range(6): - table_coeff.print_line([labels_der[i]] + list(derivatives_g[:, i])) - table_coeff.print_line([labels_der[6]] + list(derivatives_g[:, -2])) - table_coeff.print_line([labels_der[7]] + list(derivatives_g[:, -1])) - - return der_matrix, derivatives_g - - def export_derivatives(self, der_matrix_g): - - folder = self.settings['folder'] + '/' + self.data.settings['SHARPy']['case'] + '/stability/' - if not os.path.exists(folder): - os.makedirs(folder) - filename = 'stability_derivatives.txt' - - u_inf = self.settings['u_inf'].value - s_ref = self.settings['S_ref'].value - b_ref = self.settings['b_ref'].value - c_ref = self.settings['c_ref'].value - rho = self.data.linear.tsaero0.rho - euler_orient = algebra.quat2euler(self.data.linear.tsstruct0.quat) * 180/np.pi - - labels_der = {0: 'u', - 1: 'v', - 2: 'w', - 3: 'p', - 4: 'q', - 5: 'r', - 6: 'alpha', - 7: 'beta'} - - labels_out = {0: 'C_D', - 1: 'C_Y', - 2: 'C_L', - 3: 'C_l', - 4: 'C_m', - 5: 'C_n'} - - separator = '\n' + 80*'#' + '\n' - - with open(folder + '/' + filename, mode='w') as outfile: - outfile.write('SHARPy Stability Derivatives Analysis\n') - - outfile.write('State:\n') - outfile.write('\t%.4f\t\t\t # Free stream velocity\n' % u_inf) - outfile.write('\t%.4f\t\t\t # Free stream density\n' % rho) - outfile.write('\t%.4f\t\t\t # Alpha [deg]\n' % euler_orient[1]) - outfile.write('\t%.4f\t\t\t # Beta [deg]\n' % euler_orient[2]) - - outfile.write(separator) - outfile.write('\nReference Dimensions:\n') - outfile.write('\t%.4f\t\t\t # Reference planform area\n' % s_ref) - outfile.write('\t%.4f\t\t\t # Reference chord\n' % c_ref) - outfile.write('\t%.4f\t\t\t # Reference span\n' % b_ref) - - outfile.write(separator) - outfile.write('\nCoefficients:\n') - coeffs = self.static_state() - for i in range(3): - outfile.write('\t%.4f\t\t\t # %s\n' % (coeffs[i], labels_out[i])) - - outfile.write(separator) - for k, v in labels_der.items(): - outfile.write('%s derivatives:\n' % v) - for i in range(6): - outfile.write('\t%.4f\t\t\t # %s_%s derivative\n' % (der_matrix_g[i, k], labels_out[i], labels_der[k])) - outfile.write(separator) - def steady_aero_forces(self): - fx = np.sum(self.data.aero.timestep_info[0].inertial_steady_forces[:, 0], 0) + \ - np.sum(self.data.aero.timestep_info[0].inertial_unsteady_forces[:, 0], 0) - - fy = np.sum(self.data.aero.timestep_info[0].inertial_steady_forces[:, 1], 0) + \ - np.sum(self.data.aero.timestep_info[0].inertial_unsteady_forces[:, 1], 0) - - fz = np.sum(self.data.aero.timestep_info[0].inertial_steady_forces[:, 2], 0) + \ - np.sum(self.data.aero.timestep_info[0].inertial_unsteady_forces[:, 2], 0) - + # Find ref forces in G + fx, fy, fz = self.data.linear.tsaero0.total_steady_inertial_forces[:3] return fx, fy, fz def static_state(self): From 969194ca584d58de85aef898662f013e98b13347 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Sun, 9 May 2021 12:29:03 +0200 Subject: [PATCH 04/19] Update test to include linear polar coefficients in StabilityDerivatives --- tests/uvlm/static/polars/generate_wing.py | 5 +- tests/uvlm/static/polars/test_polars.py | 56 +++++++++++++++++++++++ 2 files changed, 59 insertions(+), 2 deletions(-) diff --git a/tests/uvlm/static/polars/generate_wing.py b/tests/uvlm/static/polars/generate_wing.py index 5cdb1671c..2b77aa224 100644 --- a/tests/uvlm/static/polars/generate_wing.py +++ b/tests/uvlm/static/polars/generate_wing.py @@ -60,6 +60,7 @@ def generate_infinite_wing(case_name, alpha, **kwargs): alpha_rad, 0.]))} + u_inf_direction = algebra.rotation3d_y(alpha_rad * 0).T.dot(np.array([1, 0, 0])) settings['AerogridLoader'] = {'unsteady': 'on', 'aligned_grid': 'on', 'mstar': int(kwargs.get('wake_length', 100) * m), @@ -70,7 +71,7 @@ def generate_infinite_wing(case_name, alpha, **kwargs): 'wake_shape_generator': 'StraightWake', 'wake_shape_generator_input': { 'u_inf': u_inf, - 'u_inf_direction': [1., 0., 0.], + 'u_inf_direction': u_inf_direction, 'dt': wing.dt, }, } @@ -94,7 +95,7 @@ def generate_infinite_wing(case_name, alpha, **kwargs): 'vortex_radius': 1e-6, 'velocity_field_generator': 'SteadyVelocityField', 'velocity_field_input': {'u_inf': u_inf, - 'u_inf_direction': [1., 0, 0]}, + 'u_inf_direction': u_inf_direction}, 'rho': rho} settings['StaticCoupled'] = {'print_info': 'off', diff --git a/tests/uvlm/static/polars/test_polars.py b/tests/uvlm/static/polars/test_polars.py index ac5b883df..452c9e2f2 100644 --- a/tests/uvlm/static/polars/test_polars.py +++ b/tests/uvlm/static/polars/test_polars.py @@ -7,6 +7,9 @@ from sharpy.aero.utils.utils import local_stability_axes import sharpy.utils.algebra as algebra +import h5py +import sharpy.utils.h5utils as h5utils +import sharpy.aero.utils.airfoilpolars as airfoilpolars class InfiniteWing: @@ -25,6 +28,10 @@ class TestAirfoilPolars(unittest.TestCase): polar_data = np.loadtxt(route_test_dir + '/xf-naca0018-il-50000.txt', skiprows=12) + polar = airfoilpolars.polar() + polar.initialise(np.column_stack((polar_data[:, 0] * np.pi / 180, polar_data[:, 1], polar_data[:, 2], + polar_data[:, 4]))) + def test_infinite_wing(self): """ Infinite wing should yield same results as airfoil polar @@ -68,6 +75,55 @@ def test_infinite_wing(self): cm_polar = np.interp(results[:, 0], self.polar_data[:, 0], self.polar_data[:, 4]) np.testing.assert_array_almost_equal(cm_polar, results[:, 3], decimal=3) + def test_linear_infinite_wing(self): + + cases_route = self.route_test_dir + '/cases/' + output_route = self.route_test_dir + '/output/' + + wing = InfiniteWing() + + # put inside function and test diff alpha + for alpha in [0, 2.5, 5]: + with self.subTest(alpha): + self.run_linear_test(alpha, cases_route, output_route) + + def run_linear_test(self, alpha, cases_route, output_route): + case_name = 'linear_alpha{:04g}'.format(alpha * 100).replace('-', 'M') + + flow = ['BeamLoader', + 'AerogridLoader', + 'StaticCoupled', + 'AeroForcesCalculator', + 'Modal', + 'LinearAssembler', + 'StabilityDerivatives', + 'SaveParametricCase'] + + gw.generate_infinite_wing(case_name, + alpha=alpha, + flow=flow, + case_route=cases_route, + polar_file=self.route_test_dir + '/xf-naca0018-il-50000.txt', + aspect_ratio=1e7, + main_ea=0.25, + output_route=output_route) + + derivatives = self.postprocess_linear(case_name) + + with self.subTest(msg='CL_alpha at {:f}'.format(alpha)): + cla_sharpy = derivatives['force_angle_velocity'][2, 1] + cla_polar = self.polar.get_derivatives_at_aoa(alpha * np.pi / 180)[0] + delta = np.abs(cla_sharpy - cla_polar) / cla_polar + np.testing.assert_array_less(delta, 0.15, f'Difference in polar {cla_polar:0.3f} and ' + f'SHARPy CL_alpha {cla_sharpy:0.3f} ' + f'greater than 15% ({delta * 100:0.3f} %)') + + def postprocess_linear(self, case_name): + with h5py.File(self.route_test_dir + '/output/' + case_name + + '/derivatives/aerodynamic_stability_derivatives.h5', 'r') as f: + aerodynamic_ders = h5utils.load_h5_in_dict(f) + return aerodynamic_ders + def postprocess(self, output_folder): cases = glob.glob(output_folder + '/*') From 16758cb56e4bdf2073fb7c0b34cb6ba9ba504006 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Mon, 10 May 2021 13:42:01 +0100 Subject: [PATCH 05/19] Support for efficiency corrections in linear --- sharpy/generators/polaraeroforces.py | 104 +++++++++++++++++++++--- tests/uvlm/static/polars/test_polars.py | 2 +- 2 files changed, 95 insertions(+), 11 deletions(-) diff --git a/sharpy/generators/polaraeroforces.py b/sharpy/generators/polaraeroforces.py index 00ffb261e..83db8b4d1 100644 --- a/sharpy/generators/polaraeroforces.py +++ b/sharpy/generators/polaraeroforces.py @@ -8,7 +8,7 @@ @generator_interface.generator class PolarCorrection(generator_interface.BaseGenerator): - """ + r""" This generator corrects the aerodynamic forces from UVLM based on the airfoil polars provided by the user in the ``aero.h5`` file. Polars are entered for each airfoil, in a table comprising ``AoA (rad), CL, CD, CM``. @@ -93,7 +93,7 @@ def initialise(self, in_dict, **kwargs): self.vortex_radius = kwargs.get('vortex_radius', 1e-6) def generate(self, **params): - """ + r""" Keyword Args: aero_kstep (:class:`sharpy.utils.datastructures.AeroTimeStepInfo`): Current aerodynamic substep structural_kstep (:class:`sharpy.utils.datastructures.StructTimeStepInfo`): Current structural substep @@ -125,9 +125,6 @@ def generate(self, **params): cga = algebra.quat2rotation(structural_kstep.quat) pos_g = np.array([cga.dot(structural_kstep.pos[inode]) + np.array([0, 0, 0]) for inode in range(nnode)]) - aero_kstep.polar_coeff = np.zeros((nnode, 6)) - aero_kstep.stability_transform = np.zeros((nnode, 3, 3)) - structural_kstep.postproc_node['aoa'] = np.zeros(nnode) for inode in range(nnode): @@ -247,9 +244,6 @@ def generate_linear(self, **params): cga = algebra.quat2rotation(structural_kstep.quat) - aero_kstep.polar_coeff = np.zeros((nnode, 6)) - aero_kstep.stability_transform = np.zeros((nnode, 3, 3)) - # Matrix allocation num_dof_str = beam.sys.num_dof_str num_dof_rig = beam.sys.num_dof_rig @@ -331,7 +325,7 @@ def generate_linear(self, **params): @generator_interface.generator class EfficiencyCorrection(generator_interface.BaseGenerator): - """ + r""" The efficiency and constant terms are introduced by means of the array ``airfoil_efficiency`` in the ``aero.h5`` .. math:: @@ -401,4 +395,94 @@ def generate(self, **params): return new_struct_forces def generate_linear(self, **params): - raise NotImplementedError('Efficiency factors not yet implemented') + beam = params['beam'] + tsstruct0 = params['tsstruct0'] + tsaero0 = params['tsaero0'] + + aerogrid = self.aero + structure = self.structure + + structural_kstep = tsstruct0 + aero_kstep = tsaero0 + nnode = tsstruct0.pos.shape[0] + + aero_dict = aerogrid.aero_dict + + cga = algebra.quat2rotation(structural_kstep.quat) + + # Matrix allocation + num_dof_str = beam.sys.num_dof_str + num_dof_rig = beam.sys.num_dof_rig + use_euler = beam.sys.use_euler + + # GEBM degrees of freedom + jj_for_tra = range(num_dof_str - num_dof_rig, + num_dof_str - num_dof_rig + 3) + jj_for_rot = range(num_dof_str - num_dof_rig + 3, + num_dof_str - num_dof_rig + 6) + + if use_euler: + jj_euler = range(num_dof_str - 3, num_dof_str) + euler = algebra.quat2euler(tsstruct0.quat) + tsstruct0.euler = euler + else: + jj_quat = range(num_dof_str - 4, num_dof_str) + + polar_gain = np.zeros((num_dof_str, num_dof_str)) + + jj = 0 # Global DOF index + for inode in range(nnode): + if aero_dict['aero_node'][inode]: + + ### detect bc at node (and no. of dofs) + bc_here = structure.boundary_conditions[inode] + + if bc_here == 1: # clamp (only rigid-body) + dofs_here = 0 + jj_tra, jj_rot = [], [] + + elif bc_here == -1 or bc_here == 0: # (rigid+flex body) + dofs_here = 6 + jj_tra = 6 * structure.vdof[inode] + np.array([0, 1, 2], dtype=int) + jj_rot = 6 * structure.vdof[inode] + np.array([3, 4, 5], dtype=int) + else: + raise NameError('Invalid boundary condition (%d) at node %d!' \ + % (bc_here, inode)) + + jj += dofs_here + + ielem, inode_in_elem = structure.node_master_elem[inode] + iairfoil = aero_dict['airfoil_distribution'][ielem, inode_in_elem] + isurf = aerogrid.struct2aero_mapping[inode][0]['i_surf'] + i_n = aerogrid.struct2aero_mapping[inode][0]['i_n'] + N = aerogrid.aero_dimensions[isurf, 1] + + airfoil_efficiency = aero_dict['airfoil_efficiency'] + + cab = algebra.crv2rotation(structural_kstep.psi[ielem, inode_in_elem, :]) + cgb = np.dot(cga, cab) + + dir_span, span, dir_chord, chord = span_chord(i_n, aero_kstep.zeta[isurf]) + + # Define the relative velocity and its direction + urel, dir_urel = magnitude_and_direction_of_relative_velocity(structural_kstep.pos[inode, :], + structural_kstep.pos_dot[inode, :], + structural_kstep.for_vel[:], + cga, + aero_kstep.u_ext[isurf][:, :, i_n]) + + # Stability axes - projects forces in B onto S + c_bs = local_stability_axes(cgb.T.dot(dir_urel), cgb.T.dot(dir_chord)) + cas = cab.dot(c_bs) + + local_correction = np.diag(airfoil_efficiency[ielem, inode_in_elem, 0, :]) + + weight = 1 / np.sum(aero_dict['aero_node']) + if bc_here != 1: + polar_gain[np.ix_(jj_tra, jj_tra)] += cas.dot(local_correction[:3, :3]).dot(cas.T) + polar_gain[np.ix_(jj_rot, jj_rot)] += cas.dot(local_correction[3:, 3:]).dot(cas.T) + + polar_gain[np.ix_(jj_for_tra, jj_for_tra)] += weight * cas.dot(local_correction[:3, :3]).dot(cas.T) + polar_gain[np.ix_(jj_for_rot, jj_for_rot)] += weight * cas.dot(local_correction[3:, 3:]).dot(cas.T) + + return polar_gain diff --git a/tests/uvlm/static/polars/test_polars.py b/tests/uvlm/static/polars/test_polars.py index 38b6689ea..6b699ff5b 100644 --- a/tests/uvlm/static/polars/test_polars.py +++ b/tests/uvlm/static/polars/test_polars.py @@ -28,7 +28,7 @@ class TestAirfoilPolars(unittest.TestCase): polar_data = np.loadtxt(route_test_dir + '/xf-naca0018-il-50000.txt', skiprows=12) - polar = airfoilpolars.polar() + polar = airfoilpolars.Polar() polar.initialise(np.column_stack((polar_data[:, 0] * np.pi / 180, polar_data[:, 1], polar_data[:, 2], polar_data[:, 4]))) From 3e885e61214fe2a0fae84b0e40dbb14f35533312 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Mon, 10 May 2021 16:10:10 +0100 Subject: [PATCH 06/19] Fix class attribute initialisation --- sharpy/linear/utils/derivatives.py | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/sharpy/linear/utils/derivatives.py b/sharpy/linear/utils/derivatives.py index ee7ace461..e63149bce 100644 --- a/sharpy/linear/utils/derivatives.py +++ b/sharpy/linear/utils/derivatives.py @@ -54,22 +54,20 @@ def initialise_derivatives(self, state_space, steady_forces, quat, v0, phi=None) """ cls = DerivativeSet - if cls.quat is None: - cls.quat = quat - cls.cga = algebra.quat2rotation(cls.quat) - cls.v0 = v0 - cls.coefficients = self.coefficients - - if phi is not None: - cls.modal = True - cls.phi = phi[-9:-3, :6] - cls.inv_phi_forces = np.linalg.inv(phi[-9:-3, :6].T) - cls.inv_phi_vel = np.linalg.inv(phi[-9:-3, :6]) - else: - cls.modal = False + cls.quat = quat + cls.cga = algebra.quat2rotation(cls.quat) + cls.v0 = v0 + cls.coefficients = self.coefficients + + if phi is not None: + cls.modal = True + cls.phi = phi[-9:-3, :6] + cls.inv_phi_forces = np.linalg.inv(phi[-9:-3, :6].T) + cls.inv_phi_vel = np.linalg.inv(phi[-9:-3, :6]) + else: + cls.modal = False cls.steady_forces = steady_forces - H0 = state_space.freqresp(np.array([1e-5]))[:, :, 0] # A, B, C, D = state_space.get_mats() # H0 = C.dot(np.linalg.inv(np.eye(state_space.states) - A).dot(B)) + D From fca64322e2624ea3aacddf1aa65d74b5d040ffa1 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Fri, 28 May 2021 12:26:45 +0100 Subject: [PATCH 07/19] Update tests removing deprecated settings --- tests/io/generate_pazy_udpout.py | 10 +--------- tests/linear/derivatives/test_stabilityderivatives.py | 3 +-- tests/linear/gusts/test_external_gust.py | 1 - tests/linear/horten/test_horten.py | 1 - tests/uvlm/static/polars/generate_wing.py | 5 ++--- 5 files changed, 4 insertions(+), 16 deletions(-) diff --git a/tests/io/generate_pazy_udpout.py b/tests/io/generate_pazy_udpout.py index 1c0537ef0..a63ce263d 100644 --- a/tests/io/generate_pazy_udpout.py +++ b/tests/io/generate_pazy_udpout.py @@ -110,14 +110,10 @@ def generate_pazy_udp(u_inf, case_name, output_folder='/output/', cases_folder=' 'velocity_field_input': { 'u_inf': ws.u_inf, 'u_inf_direction': ws.u_inf_direction}, - 'rollup_dt': ws.dt, 'print_info': 'on', 'horseshoe': 'on', 'num_cores': 4, - 'n_rollup': 0, - 'rollup_aic_refresh': 0, - 'vortex_radius': 1e-9, - 'rollup_tolerance': 1e-4} + 'vortex_radius': 1e-9} settings = dict() settings['NonLinearStatic'] = {'print_info': 'off', @@ -140,10 +136,6 @@ def generate_pazy_udp(u_inf, case_name, output_folder='/output/', cases_folder=' 'print_info': 'off', 'horseshoe': 'on', 'num_cores': 4, - 'n_rollup': 0, - 'rollup_dt': ws.dt, - 'rollup_aic_refresh': 1, - 'rollup_tolerance': 1e-4, 'velocity_field_generator': 'SteadyVelocityField', 'velocity_field_input': { 'u_inf': ws.u_inf, diff --git a/tests/linear/derivatives/test_stabilityderivatives.py b/tests/linear/derivatives/test_stabilityderivatives.py index 1bde06bac..6a0fb5f22 100644 --- a/tests/linear/derivatives/test_stabilityderivatives.py +++ b/tests/linear/derivatives/test_stabilityderivatives.py @@ -143,9 +143,8 @@ def run_sharpy(self, alpha_deg, flow, target_system, not_run=False): 'minus_m_star': 0} ws.config['AeroForcesCalculator'] = {'write_text_file': 'on', - # 'text_file_name': ws.case_name + '_aeroforces.csv', 'screen_output': 'on', - 'unsteady': 'off'} + } ws.config['BeamPlot'] = {'include_rbm': 'off', 'include_applied_forces': 'on'} diff --git a/tests/linear/gusts/test_external_gust.py b/tests/linear/gusts/test_external_gust.py index ec1577990..9ce4992ac 100644 --- a/tests/linear/gusts/test_external_gust.py +++ b/tests/linear/gusts/test_external_gust.py @@ -106,7 +106,6 @@ def generate_sharpy(alpha=0., case_name='hale_static', case_route='./', **kwargs settings['AeroForcesCalculator'] = {'write_text_file': 'on', 'text_file_name': 'aeroforces.txt', 'screen_output': 'on', - 'unsteady': 'off', 'coefficients': True, 'q_ref': 0.5 * rho * u_inf ** 2, 'S_ref': 12.809, diff --git a/tests/linear/horten/test_horten.py b/tests/linear/horten/test_horten.py index 62de0d655..dc4b43314 100644 --- a/tests/linear/horten/test_horten.py +++ b/tests/linear/horten/test_horten.py @@ -158,7 +158,6 @@ def run_rom_convergence(case_name, case_route='./cases/', output_folder='./outpu settings['AeroForcesCalculator'] = {'write_text_file': 'off', 'text_file_name': ws.case_name + '_aeroforces.csv', 'screen_output': 'on', - 'unsteady': 'off', 'coefficients': True, 'q_ref': 0.5 * ws.rho * ws.u_inf ** 2, 'S_ref': 12.809, diff --git a/tests/uvlm/static/polars/generate_wing.py b/tests/uvlm/static/polars/generate_wing.py index 1b7ec7b4f..fdae26fd9 100644 --- a/tests/uvlm/static/polars/generate_wing.py +++ b/tests/uvlm/static/polars/generate_wing.py @@ -121,11 +121,11 @@ def generate_infinite_wing(case_name, alpha, **kwargs): 'rigid_body_modes': True, 'write_modes_vtk': 'on', 'print_matrices': 'on', - 'write_data': 'on', + 'save_data': 'off', 'continuous_eigenvalues': 'off', 'plot_eigenvalues': False, 'rigid_modes_ppal_axes': 'on', - 'folder': output_route} + } # ROM settings rom_settings = dict() @@ -154,7 +154,6 @@ def generate_infinite_wing(case_name, alpha, **kwargs): 'density': rho, 'remove_predictor': 'off', 'use_sparse': False, - 'rigid_body_motion': True, 'vortex_radius': 1e-7, 'remove_inputs': ['u_gust'], 'convert_to_ct': 'on', From 224a53045fd7977cdcbc22b734a2a1d37bf4238f Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Wed, 20 Oct 2021 18:21:41 +0100 Subject: [PATCH 08/19] Increase verbose for StabilityDerivatives (output moment coeffs) --- sharpy/linear/utils/derivatives.py | 10 ++++++++-- sharpy/postproc/stabilityderivatives.py | 7 ++++--- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/sharpy/linear/utils/derivatives.py b/sharpy/linear/utils/derivatives.py index cbe84d94a..10862af85 100644 --- a/sharpy/linear/utils/derivatives.py +++ b/sharpy/linear/utils/derivatives.py @@ -34,7 +34,11 @@ def __init__(self, reference_dimensions, static_state, target_system=None): 'force_angular_vel': self.dynamic_pressure * s_ref * c_ref / u_inf, 'moment_lon_angular_vel': self.dynamic_pressure * s_ref * c_ref * c_ref / u_inf} # missing rates - self.steady_coefficients = np.array(self.static_state) / self.coefficients['force'] + self.steady_coefficients = np.array(self.static_state) + self.steady_coefficients[:3] /= self.coefficients['force'] + self.steady_coefficients[3] /= self.coefficients['moment_lat'] + self.steady_coefficients[4] /= self.coefficients['moment_lon'] + self.steady_coefficients[5] /= self.coefficients['moment_lat'] self.filename = 'stability_derivatives.txt' if target_system is not None: @@ -70,6 +74,8 @@ def initialise_derivatives(self, state_space, steady_forces, quat, v0, phi=None, cls.modal = False cls.steady_forces = steady_forces + import pdb; pdb.set_trace() + H0 = state_space.freqresp(np.array([1e-5]))[:, :, 0] # A, B, C, D = state_space.get_mats() # H0 = C.dot(np.linalg.inv(np.eye(state_space.states) - A).dot(B)) + D @@ -157,7 +163,7 @@ def savetxt(self, folder): outfile.write('\t{:4f}\t\t\t # Reference span\n'.format(b_ref)) outfile.write(separator) - outfile.write('\nCoefficients:\n') + outfile.write('\nSteady-state Force and Moment Coefficients:\n') for ith, coeff in enumerate(self.steady_coefficients): outfile.write('\t{:4e}\t\t\t # {:s}\n'.format(coeff, labels_out[ith])) diff --git a/sharpy/postproc/stabilityderivatives.py b/sharpy/postproc/stabilityderivatives.py index 1be543d53..bf6190731 100644 --- a/sharpy/postproc/stabilityderivatives.py +++ b/sharpy/postproc/stabilityderivatives.py @@ -260,11 +260,12 @@ def uvlm_steady_state_transfer_function(self): def steady_aero_forces(self): # Find ref forces in G fx, fy, fz = self.data.linear.tsaero0.total_steady_inertial_forces[:3] - return fx, fy, fz + mx, my, mz = self.data.linear.tsaero0.total_steady_inertial_forces[3:] + return fx, fy, fz, mx, my, mz def static_state(self): - fx, fy, fz = self.steady_aero_forces() - force_coeff = 0.5 * self.data.linear.tsaero0.rho * self.settings['u_inf'].value ** 2 * self.settings['S_ref'].value + fx, fy, fz, mx, my, mz = self.steady_aero_forces() + force_coeff = 0.5 * self.data.linear.tsaero0.rho * self.settings['u_inf'].value ** 2 * self.settings['S_ref'] Cfx = fx / force_coeff Cfy = fy / force_coeff Cfz = fz / force_coeff From ea8e2cc469fca4e3599ebf0624bd461d273fb5c6 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Wed, 20 Oct 2021 18:22:36 +0100 Subject: [PATCH 09/19] Fix bug in computing total moments in A frame --- sharpy/aero/utils/mapping.py | 31 +++++++++++++++++++++++++ sharpy/postproc/aeroforcescalculator.py | 11 +++++++-- 2 files changed, 40 insertions(+), 2 deletions(-) diff --git a/sharpy/aero/utils/mapping.py b/sharpy/aero/utils/mapping.py index e3b359294..9bde08983 100644 --- a/sharpy/aero/utils/mapping.py +++ b/sharpy/aero/utils/mapping.py @@ -75,3 +75,34 @@ def aero2struct_force_mapping(aero_forces, struct_forces[i_global_node, 3:6] += np.dot(cbg, algebra.cross3(chi_g, aero_forces[i_surf][0:3, i_m, i_n])) return struct_forces + + +def total_forces_moments(forces_nodes_a, + pos_def, + ref_pos=np.array([0., 0., 0.])): + """ + Performs a summation of the forces and moments expressed at the structural nodes in the A frame of reference. + + Note: + If you need to transform forces and moments at the nodes from B to A, use the + :func:`~sharpy.structure.models.beam.Beam.nodal_b_for_2_a_for()` function. + + Args: + forces_nodes_a (np.array): ``n_node x 6`` vector of forces and moments at the nodes in A + pos_def (np.array): ``n_node x 3`` vector of nodal positions in A + ref_pos (np.array (optional)): Location in A about which to compute moments. Defaults to ``[0, 0, 0]`` + + Returns: + np.array: Vector of length 6 containing the total forces and moments expressed in A at the desired location. + """ + + num_node = pos_def.shape[0] + ra_vec = pos_def - ref_pos + + total_forces = np.zeros(3) + total_moments = np.zeros(3) + for i_global_node in range(num_node): + total_forces += forces_nodes_a[i_global_node, :3] + total_moments += forces_nodes_a[i_global_node, 3:] + algebra.cross3(ra_vec[i_global_node], forces_nodes_a[i_global_node, :3]) + + return np.concatenate((total_forces, total_moments)) diff --git a/sharpy/postproc/aeroforcescalculator.py b/sharpy/postproc/aeroforcescalculator.py index 6e7d2981e..0e74b722f 100644 --- a/sharpy/postproc/aeroforcescalculator.py +++ b/sharpy/postproc/aeroforcescalculator.py @@ -144,8 +144,15 @@ def calculate_forces(self, ts): self.data.structure.timestep_info[ts]) # Express total forces in A frame - self.data.aero.timestep_info[ts].total_steady_body_forces = np.sum(steady_forces_a, axis=0) - self.data.aero.timestep_info[ts].total_unsteady_body_forces = np.sum(unsteady_forces_a, axis=0) + moment_reference_location = np.array([0., 0., 0.]) + 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=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=moment_reference_location) # Express total forces in G frame self.data.aero.timestep_info[ts].total_steady_inertial_forces = \ From ff07fdb808a0c6cd1dfc13e78ea0000a906c32e0 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Fri, 22 Apr 2022 17:46:47 +0200 Subject: [PATCH 10/19] Remove import pdb --- sharpy/linear/utils/derivatives.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/sharpy/linear/utils/derivatives.py b/sharpy/linear/utils/derivatives.py index 10862af85..3af533e6d 100644 --- a/sharpy/linear/utils/derivatives.py +++ b/sharpy/linear/utils/derivatives.py @@ -74,8 +74,6 @@ def initialise_derivatives(self, state_space, steady_forces, quat, v0, phi=None, cls.modal = False cls.steady_forces = steady_forces - import pdb; pdb.set_trace() - H0 = state_space.freqresp(np.array([1e-5]))[:, :, 0] # A, B, C, D = state_space.get_mats() # H0 = C.dot(np.linalg.inv(np.eye(state_space.states) - A).dot(B)) + D @@ -397,7 +395,6 @@ def apply_coefficients(self): phi = data.linear.linear_system.linearisation_vectors['mode_shapes'] coefficients = {'force': 1960., 'moment_lon': 1960., 'moment_lat': 62720.} - import pdb; pdb.set_trace() DerivativeSet.initialise_derivatives(state_space, steady_forces, quat, np.array([-10., 0, 0]), phi) From 715e6d84ae4cb1891a27f27c477c8fe2aaf92ce7 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Fri, 22 Apr 2022 17:47:44 +0200 Subject: [PATCH 11/19] Fix total aero forces and moments in the case of corrected forces --- sharpy/generators/polaraeroforces.py | 81 +++++++++++++++++--- sharpy/linear/assembler/linearaeroelastic.py | 19 +++-- 2 files changed, 82 insertions(+), 18 deletions(-) diff --git a/sharpy/generators/polaraeroforces.py b/sharpy/generators/polaraeroforces.py index 689ddc398..eadb9ef10 100644 --- a/sharpy/generators/polaraeroforces.py +++ b/sharpy/generators/polaraeroforces.py @@ -217,8 +217,42 @@ def generate(self, **params): return new_struct_forces def generate_linear(self, **params): - """ - Generate Linear Gain matrix with correction factors (lift only at the moment) + r""" + Generate Linear Gain matrix with correction factors (lift only at the moment). These corrections are implemented + by means of a Gain that goes from forces and moments in A frame to corrected forces and moments in A frame. + + The nodal corrections are implemented as + + .. math:: K_f^i \leftarrow w_f^i + + where + + .. math:: w_f^i = C^{AS} W_f^i C^{SA} + + and :math:`W_f` is a 3x3 correction matrix expressed in stability axes and is given for each node. + + For the moments, the correction is applied in a similar way, noting that the moments are expressed as + :math:`T^\top \delta m_B` and thus need pre- and post- multiplication + + .. math:: K_m^i \leftarrow T^\top C^{BA} w_m^i C^{AB} T^{-\top} + + where + + .. math:: w_m^i = C^{AS} W_m^i C^{SA} + + and :math:`W_m^i` is a 3x3 correction matrix expressed in stability axes for the moment correction at each node. + + The total forces and moments are then computed as a summation of the above, where the total forces + + are computed simply as + + .. math:: K_F \leftarrow \Sum w_f^i + + and the total moments + + .. math:: K_M \leftarrow \Sum \tilde{R}_{A,i} w_f^i + \Sum w_m^i C^{AB} T^{-\top} + + which leaves the end result expressed in the A frame. Keyword Args: beam (sharpy.linear.assembler.linearbeam.LinearBeam): Beam object @@ -292,8 +326,11 @@ def generate_linear(self, **params): i_n = aerogrid.struct2aero_mapping[inode][0]['i_n'] N = aerogrid.aero_dimensions[isurf, 1] polar = aerogrid.polars[iairfoil] - cab = algebra.crv2rotation(structural_kstep.psi[ielem, inode_in_elem, :]) + psi = structural_kstep.psi[ielem, inode_in_elem, :] + cab = algebra.crv2rotation(psi) cgb = np.dot(cga, cab) + Tan = algebra.crv2tan(psi) + r_a = tsstruct0.pos dir_span, span, dir_chord, chord = span_chord(i_n, aero_kstep.zeta[isurf]) @@ -313,13 +350,24 @@ def generate_linear(self, **params): cla, cda, cma = polar.get_derivatives_at_aoa(tsstruct0.postproc_node['aoa'][inode]) local_correction[2, 2] = cla / 2 / np.pi - weight = 1 / np.sum(aero_dict['aero_node']) + mom_b2a = cab.dot(np.linalg.inv(Tan.T)) + mom_a2b = Tan.T.dot(cab) + + wfi = cas.dot(local_correction[:3, :3]).dot(cas.T) + wmi = cas.dot(local_correction[3:, 3:]).dot(cas.T) + if bc_here != 1: - polar_gain[np.ix_(jj_tra, jj_tra)] += cas.dot(local_correction[:3, :3]).dot(cas.T) - polar_gain[np.ix_(jj_rot, jj_rot)] += cas.dot(local_correction[3:, 3:]).dot(cas.T) + polar_gain[np.ix_(jj_tra, jj_tra)] += wfi + polar_gain[np.ix_(jj_rot, jj_rot)] += mom_a2b.dot(wmi.dot(mom_b2a)) + + # Total forces and moments + # total forces + polar_gain[np.ix_(jj_for_tra, jj_tra)] += wfi - polar_gain[np.ix_(jj_for_tra, jj_for_tra)] += weight * cas.dot(local_correction[:3, :3]).dot(cas.T) - polar_gain[np.ix_(jj_for_rot, jj_for_rot)] += weight * cas.dot(local_correction[3:, 3:]).dot(cas.T) + # forces contribution to total moments + polar_gain[np.ix_(jj_for_rot, jj_tra)] += algebra.skew(r_a).dot(wfi) + # moments contribution to total moments + polar_gain[np.ix_(jj_for_rot, jj_rot)] += wmi.dot(mom_b2a) return polar_gain @@ -457,10 +505,13 @@ def generate_linear(self, **params): isurf = aerogrid.struct2aero_mapping[inode][0]['i_surf'] i_n = aerogrid.struct2aero_mapping[inode][0]['i_n'] N = aerogrid.aero_dimensions[isurf, 1] + psi = structural_kstep.psi[ielem, inode_in_elem, :] + cab = algebra.crv2rotation(psi) + Tan = algebra.crv2tan(psi) + r_a = tsstruct0.pos airfoil_efficiency = aero_dict['airfoil_efficiency'] - cab = algebra.crv2rotation(structural_kstep.psi[ielem, inode_in_elem, :]) cgb = np.dot(cga, cab) dir_span, span, dir_chord, chord = span_chord(i_n, aero_kstep.zeta[isurf]) @@ -478,12 +529,18 @@ def generate_linear(self, **params): local_correction = np.diag(airfoil_efficiency[ielem, inode_in_elem, 0, :]) - weight = 1 / np.sum(aero_dict['aero_node']) if bc_here != 1: polar_gain[np.ix_(jj_tra, jj_tra)] += cas.dot(local_correction[:3, :3]).dot(cas.T) polar_gain[np.ix_(jj_rot, jj_rot)] += cas.dot(local_correction[3:, 3:]).dot(cas.T) - polar_gain[np.ix_(jj_for_tra, jj_for_tra)] += weight * cas.dot(local_correction[:3, :3]).dot(cas.T) - polar_gain[np.ix_(jj_for_rot, jj_for_rot)] += weight * cas.dot(local_correction[3:, 3:]).dot(cas.T) + # Total forces and moments + # total forces + polar_gain[np.ix_(jj_for_tra, jj_tra)] += cas.dot(local_correction[:3, :3]).dot(cas.T) + # forces contribution to total moments + polar_gain[np.ix_(jj_for_rot, jj_tra)] += \ + algebra.skew(r_a) * cas.dot(local_correction[:3, :3]).dot(cas.T) + # moments contribution to total moments + polar_gain[np.ix_(jj_for_rot, jj_rot)] += \ + cas.dot(local_correction[3:, 3:]).dot(cas.T) * cab * np.linalg.inv(Tan.T) return polar_gain diff --git a/sharpy/linear/assembler/linearaeroelastic.py b/sharpy/linear/assembler/linearaeroelastic.py index d1372da2b..9843464ae 100644 --- a/sharpy/linear/assembler/linearaeroelastic.py +++ b/sharpy/linear/assembler/linearaeroelastic.py @@ -845,12 +845,19 @@ def get_gebm2uvlm_gains(self, data): np.dot(Tan.T, np.dot(Cbg, algebra.skew(Xg))) # or, equivalently, np.dot( algebra.skew(Xb),Cbg) - # total forces - Kforces[np.ix_(jj_for_tra, ii_vert)] += Cag - - # total moments - Kforces[np.ix_(jj_for_rot, ii_vert)] += \ - np.dot(Cag, algebra.skew(zetag)) + # Total forces + if not self.correct_forces: + # total forces + Kforces[np.ix_(jj_for_tra, ii_vert)] += Cag + + # total moments + Kforces[np.ix_(jj_for_rot, ii_vert)] += \ + np.dot(Cag, algebra.skew(zetag)) + else: + pass + # when forces are corrected via polar these terms are included in the polar_gain, + # as the coefficient for each nodal force cannot be added to the total component + # a posteriori # quaternion equation # null, as not dep. on external forces From ee2b8f15ae644822d5bf65de21707e0fd5bdfd72 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Mon, 25 Apr 2022 11:44:26 +0200 Subject: [PATCH 12/19] Fix getting position of current node --- sharpy/generators/polaraeroforces.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sharpy/generators/polaraeroforces.py b/sharpy/generators/polaraeroforces.py index eadb9ef10..98054b590 100644 --- a/sharpy/generators/polaraeroforces.py +++ b/sharpy/generators/polaraeroforces.py @@ -330,7 +330,7 @@ def generate_linear(self, **params): cab = algebra.crv2rotation(psi) cgb = np.dot(cga, cab) Tan = algebra.crv2tan(psi) - r_a = tsstruct0.pos + r_a = tsstruct0.pos[inode] dir_span, span, dir_chord, chord = span_chord(i_n, aero_kstep.zeta[isurf]) @@ -508,7 +508,7 @@ def generate_linear(self, **params): psi = structural_kstep.psi[ielem, inode_in_elem, :] cab = algebra.crv2rotation(psi) Tan = algebra.crv2tan(psi) - r_a = tsstruct0.pos + r_a = tsstruct0.pos[inode] airfoil_efficiency = aero_dict['airfoil_efficiency'] From c38dd12b1be18df4a4a5eef91d636fe479d16460 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Tue, 26 Apr 2022 13:49:14 +0200 Subject: [PATCH 13/19] Fix process for computing the corrected forces at the A frame node --- sharpy/generators/polaraeroforces.py | 50 +++++++++++++++---- sharpy/linear/assembler/linearaeroelastic.py | 20 +++----- .../derivatives/test_stabilityderivatives.py | 24 ++++++++- tests/uvlm/static/polars/test_polars.py | 8 +++ 4 files changed, 78 insertions(+), 24 deletions(-) diff --git a/sharpy/generators/polaraeroforces.py b/sharpy/generators/polaraeroforces.py index 98054b590..20368eb3f 100644 --- a/sharpy/generators/polaraeroforces.py +++ b/sharpy/generators/polaraeroforces.py @@ -254,6 +254,11 @@ def generate_linear(self, **params): which leaves the end result expressed in the A frame. + An additional gain, referred to as ``forces_at_ref_gain`` is needed to obtain the nodal force at the + :math:`A` frame in order to correct the forces at this point with the local polar. Thus, this gain + subtracts from the total forces and moments the sum of those from the remaining nodes. Then the ``polar_gain`` + corrects the local force at the A frame and sums the corrected forces and moments at the other nodes. + Keyword Args: beam (sharpy.linear.assembler.linearbeam.LinearBeam): Beam object tsstruct0 (sharpy.utils.datastructures.StructTimestepInfo): Ref structural time step @@ -299,6 +304,10 @@ def generate_linear(self, **params): polar_gain = np.zeros((num_dof_str, num_dof_str)) + # this gain computes the aerodynamic forces and moments at the node at the A frame, which else + # would be included in the total force and total moment. This is needed to correct these nodal forces too + forces_at_ref_gain = np.zeros_like(polar_gain) + jj = 0 # Global DOF index for inode in range(nnode): if aero_dict['aero_node'][inode]: @@ -351,7 +360,7 @@ def generate_linear(self, **params): local_correction[2, 2] = cla / 2 / np.pi mom_b2a = cab.dot(np.linalg.inv(Tan.T)) - mom_a2b = Tan.T.dot(cab) + mom_a2b = Tan.T.dot(cab.T) wfi = cas.dot(local_correction[:3, :3]).dot(cas.T) wmi = cas.dot(local_correction[3:, 3:]).dot(cas.T) @@ -360,16 +369,39 @@ def generate_linear(self, **params): polar_gain[np.ix_(jj_tra, jj_tra)] += wfi polar_gain[np.ix_(jj_rot, jj_rot)] += mom_a2b.dot(wmi.dot(mom_b2a)) - # Total forces and moments - # total forces - polar_gain[np.ix_(jj_for_tra, jj_tra)] += wfi + # Total forces and moments + # total forces + polar_gain[np.ix_(jj_for_tra, jj_tra)] += wfi - # forces contribution to total moments - polar_gain[np.ix_(jj_for_rot, jj_tra)] += algebra.skew(r_a).dot(wfi) - # moments contribution to total moments - polar_gain[np.ix_(jj_for_rot, jj_rot)] += wmi.dot(mom_b2a) + # forces contribution to total moments + polar_gain[np.ix_(jj_for_rot, jj_tra)] += algebra.skew(r_a).dot(wfi) + # moments contribution to total moments + polar_gain[np.ix_(jj_for_rot, jj_rot)] += wmi.dot(mom_b2a) - return polar_gain + # keep nodal forces and moments + forces_at_ref_gain[np.ix_(jj_tra, jj_tra)] += np.eye(3) + forces_at_ref_gain[np.ix_(jj_rot, jj_rot)] += np.eye(3) + + # obtain the nodal force at the A frame + # total force minus the sum of all other nodes + # f_node_at_A = F_A - sum(f) + forces_at_ref_gain[np.ix_(jj_for_tra, jj_tra)] -= np.eye(3) + + # obtain the nodal moment in A frame at the A frame + # m_node_at_A = M_A - sum(C^AB T^-\top m_b) - sum(r_a.cross(f_a)) + forces_at_ref_gain[np.ix_(jj_for_rot, jj_tra)] -= algebra.skew(r_a) + forces_at_ref_gain[np.ix_(jj_for_rot, jj_rot)] -= mom_b2a + + else: + # local force and moment at the A frame + forces_at_ref_gain[np.ix_(jj_for_tra, jj_for_tra)] = np.eye(3) + forces_at_ref_gain[np.ix_(jj_for_rot, jj_for_rot)] = np.eye(3) + + # scale nodal force at A frame node + polar_gain[np.ix_(jj_for_tra, jj_for_tra)] = wfi + polar_gain[np.ix_(jj_for_rot, jj_for_rot)] = wmi + + return polar_gain.dot(forces_at_ref_gain) @generator_interface.generator diff --git a/sharpy/linear/assembler/linearaeroelastic.py b/sharpy/linear/assembler/linearaeroelastic.py index 9843464ae..c4aa2fe85 100644 --- a/sharpy/linear/assembler/linearaeroelastic.py +++ b/sharpy/linear/assembler/linearaeroelastic.py @@ -845,19 +845,13 @@ def get_gebm2uvlm_gains(self, data): np.dot(Tan.T, np.dot(Cbg, algebra.skew(Xg))) # or, equivalently, np.dot( algebra.skew(Xb),Cbg) - # Total forces - if not self.correct_forces: - # total forces - Kforces[np.ix_(jj_for_tra, ii_vert)] += Cag - - # total moments - Kforces[np.ix_(jj_for_rot, ii_vert)] += \ - np.dot(Cag, algebra.skew(zetag)) - else: - pass - # when forces are corrected via polar these terms are included in the polar_gain, - # as the coefficient for each nodal force cannot be added to the total component - # a posteriori + # Total forces and Total moments + # total forces + Kforces[np.ix_(jj_for_tra, ii_vert)] += Cag + + # total moments + Kforces[np.ix_(jj_for_rot, ii_vert)] += \ + np.dot(Cag, algebra.skew(zetag)) # quaternion equation # null, as not dep. on external forces diff --git a/tests/linear/derivatives/test_stabilityderivatives.py b/tests/linear/derivatives/test_stabilityderivatives.py index 0874ec40b..548615b56 100644 --- a/tests/linear/derivatives/test_stabilityderivatives.py +++ b/tests/linear/derivatives/test_stabilityderivatives.py @@ -301,7 +301,7 @@ def run_case(self, target_system): NameError('Unrecognised system') case_name_db = [] - not_run = False # for debugging + not_run = False # for debugging # Reference Case at 4 degrees # Run nonlinear simulation and save linerised ROM alpha_deg_ref = 0 @@ -325,7 +325,7 @@ def run_case(self, target_system): qS = 0.5 * ref.rho * u_inf ** 2 * ref.c_ref * ref.wing_span if self.print_info: print(f'Reference span: {ref.wing_span} m') - print(f'Reference chord: {ref.main_chord} m') + print(f'Reference chord: {ref.c_ref} m') print(f'Aspect ratio: {ref.aspect_ratio}') # Run nonlinear cases in the vicinity @@ -337,6 +337,7 @@ def run_case(self, target_system): alpha_vec = np.linspace(alpha_deg_min, alpha_deg_max, n_evals) nlin_forces_g = np.zeros((n_evals, 3)) nlin_forces_a = np.zeros((n_evals, 3)) + nlin_moments_g = np.zeros((n_evals, 3)) for ith, alpha in enumerate(alpha_vec): if alpha == alpha_deg_ref: case_name = ref_case_name @@ -349,19 +350,38 @@ def run_case(self, target_system): '/output/{:s}/forces/forces_aeroforces.txt'.format(case_name), skiprows=1, delimiter=',') + nlin_moments = np.loadtxt(self.route_test_dir + + '/output/{:s}/forces/moments_aeroforces.txt'.format(case_name), + skiprows=1, + delimiter=',') nlin_forces_g[ith, :3] = nlin_forces[1:4] nlin_forces_a[ith, :3] = nlin_forces[7:10] + nlin_moments_g[ith, :3] = nlin_moments[1:4] + + if self.print_info: + np.savetxt(self.route_test_dir + '/output/summary.txt', + np.column_stack((alpha_vec, nlin_forces_g, nlin_moments_g))) nlin_forces_g /= qS nlin_forces_a /= qS + nlin_moments_g /= (qS * ref.c_ref) + + if self.print_info: + np.savetxt(self.route_test_dir + '/output/summary_coeff.txt', + np.column_stack((alpha_vec, nlin_forces_g, nlin_moments_g))) lift_poly = np.polyfit(alpha_vec * np.pi/180, nlin_forces_g[:, 2], deg=1) nonlin_cla = lift_poly[0] drag_poly = np.polyfit(alpha_vec * np.pi/180, nlin_forces_g[:, 0], deg=2) nonlin_cda = 2 * drag_poly[0] * alpha0 + drag_poly[1] + + moment_poly = np.polyfit(alpha_vec * np.pi/180, nlin_moments_g[:, 1], deg=1) + nonlin_cma = moment_poly[0] + if self.print_info: print('Nonlinear coefficients') print('CLa', nonlin_cla) print('CDa', nonlin_cda) + print('CMa', nonlin_cma) lift_poly = np.polyfit(alpha_vec * np.pi/180, nlin_forces_a[:, 2], deg=1) nonlin_cza = lift_poly[0] diff --git a/tests/uvlm/static/polars/test_polars.py b/tests/uvlm/static/polars/test_polars.py index 6b699ff5b..b0933b906 100644 --- a/tests/uvlm/static/polars/test_polars.py +++ b/tests/uvlm/static/polars/test_polars.py @@ -32,6 +32,8 @@ class TestAirfoilPolars(unittest.TestCase): polar.initialise(np.column_stack((polar_data[:, 0] * np.pi / 180, polar_data[:, 1], polar_data[:, 2], polar_data[:, 4]))) + print_info = False + def test_infinite_wing(self): """ Infinite wing should yield same results as airfoil polar @@ -111,9 +113,15 @@ def run_linear_test(self, alpha, cases_route, output_route): derivatives = self.postprocess_linear(case_name) with self.subTest(msg='CL_alpha at {:f}'.format(alpha)): + # Corrections are only implemented in the lift cla_sharpy = derivatives['force_angle_velocity'][2, 1] cla_polar = self.polar.get_derivatives_at_aoa(alpha * np.pi / 180)[0] delta = np.abs(cla_sharpy - cla_polar) / cla_polar + if self.print_info: + print(f'Alpha = {alpha}') + print(f'Polar CL_alpha = {cla_polar:0.3f}') + print(f'SHARPy CL_alpha = {cla_sharpy:0.3f}') + print(f'Relative difference = ({delta * 100:0.3f} %)') np.testing.assert_array_less(delta, 0.15, f'Difference in polar {cla_polar:0.3f} and ' f'SHARPy CL_alpha {cla_sharpy:0.3f} ' f'greater than 15% ({delta * 100:0.3f} %)') From de0e42c0f180475c54b984b48c395c403a3e417d Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Wed, 15 Jun 2022 15:52:54 +0200 Subject: [PATCH 14/19] Mix issues in merge --- sharpy/linear/utils/derivatives.py | 2 +- sharpy/postproc/stabilityderivatives.py | 27 ++++++++++++++++--------- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/sharpy/linear/utils/derivatives.py b/sharpy/linear/utils/derivatives.py index bd1f2f6ad..2faa84a10 100644 --- a/sharpy/linear/utils/derivatives.py +++ b/sharpy/linear/utils/derivatives.py @@ -59,7 +59,7 @@ def initialise_derivatives(self, state_space, steady_forces, quat, v0, phi=None, tpa (np.array (optional)): Transformation matrix onto principal axes """ - cls = DerivativeSet # explain what is the DerivativeSet class + cls = DerivativeSet # type: DerivativeSet if cls.quat is None: cls.quat = quat cls.cga = algebra.quat2rotation(cls.quat) diff --git a/sharpy/postproc/stabilityderivatives.py b/sharpy/postproc/stabilityderivatives.py index 0d8bb5fc6..5ce783bca 100644 --- a/sharpy/postproc/stabilityderivatives.py +++ b/sharpy/postproc/stabilityderivatives.py @@ -194,6 +194,15 @@ def get_freestream_velocity(self): return v0 def get_state_space(self, target_system): + """ + Returns the target state-space ``target_system`` either ``aeroelastic`` or ``aerodynamic`` + + Returns: + libss.StateSpace: relevant state-space + + Raises: + NameError: if the target system is not ``aeroelastic`` or ``aerodynamic`` + """ if target_system == 'aerodynamic': ss = self.data.linear.linear_system.uvlm.ss elif target_system == 'aeroelastic': @@ -204,15 +213,13 @@ def get_state_space(self, target_system): return ss def steady_aero_forces(self): - # Find ref forces in G - fx, fy, fz = self.data.linear.tsaero0.total_steady_inertial_forces[:3] - mx, my, mz = self.data.linear.tsaero0.total_steady_inertial_forces[3:] - return fx, fy, fz, mx, my, mz - - fy = np.sum(self.data.aero.timestep_info[0].inertial_steady_forces[:, 1], 0) + \ - np.sum(self.data.aero.timestep_info[0].inertial_unsteady_forces[:, 1], 0) - - fz = np.sum(self.data.aero.timestep_info[0].inertial_steady_forces[:, 2], 0) + \ - np.sum(self.data.aero.timestep_info[0].inertial_unsteady_forces[:, 2], 0) + """Retrieve steady aerodynamic forces at the linearisation reference + + Returns: + tuple: (fx, fy, fz) in the inertial G frame + """ + fx = np.sum(self.data.linear.tsaero0.inertial_steady_forces[:, 0], 0) + fy = np.sum(self.data.linear.tsaero0.inertial_steady_forces[:, 1], 0) + fz = np.sum(self.data.linear.tsaero0.inertial_steady_forces[:, 2], 0) return fx, fy, fz From eceff9c70d77e100c358b1b9309eb88363cdebee Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Wed, 15 Jun 2022 20:06:27 +0200 Subject: [PATCH 15/19] Recover moments to output force/mom coefficients at ref state --- sharpy/linear/utils/derivatives.py | 10 +++++++++- sharpy/postproc/stabilityderivatives.py | 20 +++++++++++++------- 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/sharpy/linear/utils/derivatives.py b/sharpy/linear/utils/derivatives.py index 2faa84a10..e1ea3d838 100644 --- a/sharpy/linear/utils/derivatives.py +++ b/sharpy/linear/utils/derivatives.py @@ -7,13 +7,19 @@ class Derivatives: """ Class containing the derivatives set for a given state-space system (i.e. aeroelastic or aerodynamic) + + Args: + reference_dimensions (dict): Info on the reference dimensions ``S_ref`` (area) ``b_ref`` (span) ``c_ref`` chord + ``u_inf`` velocity and ``rho`` density + static_state (tuple): Force and moment coefficient at the reference state (dim6) (inertial frame) + target_system (str): Name of target system ``aeroelastic`` or ``aerodynamic`` """ def __init__(self, reference_dimensions, static_state, target_system=None): self.target_system = target_system # type: str # name of target system (aerodynamic/aeroelastic) self.transfer_function = None # type: np.array # matrix of steady-state TF for target system - self.static_state = static_state # type: tuple # [fx, fy, fz] at ref state + self.static_state = static_state # type: tuple # [fx, fy, fz, mx, my, mz] at ref state self.reference_dimensions = reference_dimensions # type: dict # name: ref_dimension_value dictionary self.separator = '\n' + 80 * '#' + '\n' @@ -33,6 +39,8 @@ def __init__(self, reference_dimensions, static_state, target_system=None): 'force_angular_vel': self.dynamic_pressure * s_ref * c_ref / u_inf, 'moment_lon_angular_vel': self.dynamic_pressure * s_ref * c_ref * c_ref / u_inf} # missing rates + # only used to show the forces/moment coefficient at ref state + # for derivatives calculation the vectors at the linearisation ref are used self.steady_coefficients = np.array(self.static_state) self.steady_coefficients[:3] /= self.coefficients['force'] self.steady_coefficients[3] /= self.coefficients['moment_lat'] diff --git a/sharpy/postproc/stabilityderivatives.py b/sharpy/postproc/stabilityderivatives.py index 5ce783bca..c4b12312c 100644 --- a/sharpy/postproc/stabilityderivatives.py +++ b/sharpy/postproc/stabilityderivatives.py @@ -9,7 +9,16 @@ @solver_interface.solver class StabilityDerivatives(solver_interface.BaseSolver): """ - Outputs the stability derivatives of a free-flying aircraft + Outputs the stability derivatives of a free-flying aircraft. + + The simulation set-up to obtain the stability derivatives is not standard, and requires specific settings + in the solvers ran prior to this post-processor. Please see the tutorial at: + https://github.com/ngoiz/hale-ders/blob/main/Delivery/01_StabilityDerivatives/01_StabilityDerivatives.ipynb + + In the future, a routine will be available where these required solvers' settings are pre-populated + + Note: + Requires the AeroForcesCalculator post-processor to have been run before. """ solver_id = 'StabilityDerivatives' @@ -213,13 +222,10 @@ def get_state_space(self, target_system): return ss def steady_aero_forces(self): - """Retrieve steady aerodynamic forces at the linearisation reference + """Retrieve steady aerodynamic forces and moments at the linearisation reference at the Returns: - tuple: (fx, fy, fz) in the inertial G frame + tuple: (fx, fy, fz, mx, my, mz) in the inertial G frame """ - fx = np.sum(self.data.linear.tsaero0.inertial_steady_forces[:, 0], 0) - fy = np.sum(self.data.linear.tsaero0.inertial_steady_forces[:, 1], 0) - fz = np.sum(self.data.linear.tsaero0.inertial_steady_forces[:, 2], 0) - return fx, fy, fz + return self.data.linear.tsaero0.total_steady_inertial_forces From 684e5c30ff22a5068f867a255504c85dce7ba0b2 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Thu, 16 Jun 2022 08:44:05 +0100 Subject: [PATCH 16/19] Add verbose to test failing in CI but pass in local --- tests/uvlm/static/polars/generate_wing.py | 2 +- tests/uvlm/static/polars/test_polars.py | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/uvlm/static/polars/generate_wing.py b/tests/uvlm/static/polars/generate_wing.py index 65bf21fd5..00141e103 100644 --- a/tests/uvlm/static/polars/generate_wing.py +++ b/tests/uvlm/static/polars/generate_wing.py @@ -50,7 +50,7 @@ def generate_infinite_wing(case_name, alpha, **kwargs): settings['SHARPy'] = {'case': case_name, 'route': case_route, 'flow': kwargs.get('flow', []), - 'write_screen': 'off', + 'write_screen': kwargs.get('write_screen', 'off'), 'write_log': 'on', 'log_folder': output_route, 'log_file': case_name + '.log'} diff --git a/tests/uvlm/static/polars/test_polars.py b/tests/uvlm/static/polars/test_polars.py index b0933b906..4824f7df6 100644 --- a/tests/uvlm/static/polars/test_polars.py +++ b/tests/uvlm/static/polars/test_polars.py @@ -32,7 +32,7 @@ class TestAirfoilPolars(unittest.TestCase): polar.initialise(np.column_stack((polar_data[:, 0] * np.pi / 180, polar_data[:, 1], polar_data[:, 2], polar_data[:, 4]))) - print_info = False + print_info = True def test_infinite_wing(self): """ @@ -108,7 +108,8 @@ def run_linear_test(self, alpha, cases_route, output_route): polar_file=self.route_test_dir + '/xf-naca0018-il-50000.txt', aspect_ratio=1e7, main_ea=0.25, - output_route=output_route) + output_route=output_route, + write_screen=self.print_info) derivatives = self.postprocess_linear(case_name) From 9a62d9e1de71d0a34a2b2651db198289980c4a56 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Thu, 30 Jun 2022 14:30:06 +0200 Subject: [PATCH 17/19] Remove legacy code Apply revisions from code review --- sharpy/linear/assembler/linearaeroelastic.py | 128 ------------------- 1 file changed, 128 deletions(-) diff --git a/sharpy/linear/assembler/linearaeroelastic.py b/sharpy/linear/assembler/linearaeroelastic.py index c4aa2fe85..e70ac1de5 100644 --- a/sharpy/linear/assembler/linearaeroelastic.py +++ b/sharpy/linear/assembler/linearaeroelastic.py @@ -1025,131 +1025,3 @@ def load_uvlm(filename): # uvlm_ss_read = read_data.linear.linear_system.uvlm.ss uvlm_ss_read = read_data return libss.StateSpace(uvlm_ss_read.A, uvlm_ss_read.B, uvlm_ss_read.C, uvlm_ss_read.D, dt=uvlm_ss_read.dt) - - - def polar_gains(self, data, ss): - - import sharpy.aero.utils.utils as aeroutils - aero = data.aero - structure = data.structure - tsaero = self.uvlm.tsaero0 - tsstr = self.beam.tsstruct0 - - Kzeta = self.uvlm.sys.Kzeta - num_dof_str = self.beam.sys.num_dof_str - num_dof_rig = self.beam.sys.num_dof_rig - num_dof_flex = self.beam.sys.num_dof_flex - use_euler = self.beam.sys.use_euler - - # allocate output - Kdisp = np.zeros((3 * Kzeta, num_dof_str)) - Kdisp_vel = np.zeros((3 * Kzeta, num_dof_str)) # Orientation is in velocity DOFs - Kvel_disp = np.zeros((3 * Kzeta, num_dof_str)) - Kvel_vel = np.zeros((3 * Kzeta, num_dof_str)) - Kforces = np.zeros((num_dof_str, 3 * Kzeta)) - - Kss = np.zeros((num_dof_flex, num_dof_flex)) - Csr = np.zeros((num_dof_flex, num_dof_rig)) - Crs = np.zeros((num_dof_rig, num_dof_flex)) - Crr = np.zeros((num_dof_rig, num_dof_rig)) - Krs = np.zeros((num_dof_rig, num_dof_flex)) - - # get projection matrix A->G - # (and other quantities indep. from nodal position) - Cga = algebra.quat2rotation(tsstr.quat) # NG 6-8-19 removing .T - Cag = Cga.T - - # for_pos=tsstr.for_pos - for_vel = tsstr.for_vel[:3] - for_rot = tsstr.for_vel[3:] - skew_for_rot = algebra.skew(for_rot) - Der_vel_Ra = np.dot(Cga, skew_for_rot) - - Faero = np.zeros(3) - FaeroA = np.zeros(3) - - # GEBM degrees of freedom - jj_for_tra = range(num_dof_str - num_dof_rig, - num_dof_str - num_dof_rig + 3) - jj_for_rot = range(num_dof_str - num_dof_rig + 3, - num_dof_str - num_dof_rig + 6) - - if use_euler: - jj_euler = range(num_dof_str - 3, num_dof_str) - euler = algebra.quat2euler(tsstr.quat) - tsstr.euler = euler - else: - jj_quat = range(num_dof_str - 4, num_dof_str) - - polar_matrix_vel = np.zeros((ss.outputs, num_dof_str)) - polar_matrix_disp = np.zeros((ss.outputs, num_dof_str)) - - jj = 0 # nodal dof index - for node_glob in range(structure.num_node): - - ### detect bc at node (and no. of dofs) - bc_here = structure.boundary_conditions[node_glob] - - if bc_here == 1: # clamp (only rigid-body) - dofs_here = 0 - jj_tra, jj_rot = [], [] - # continue - - elif bc_here == -1 or bc_here == 0: # (rigid+flex body) - dofs_here = 6 - jj_tra = 6 * structure.vdof[node_glob] + np.array([0, 1, 2], dtype=int) - jj_rot = 6 * structure.vdof[node_glob] + np.array([3, 4, 5], dtype=int) - else: - raise NameError('Invalid boundary condition (%d) at node %d!' \ - % (bc_here, node_glob)) - - jj += dofs_here - - # retrieve element and local index - ee, node_loc = structure.node_master_elem[node_glob, :] - - # get position, crv and rotation matrix - Ra = tsstr.pos[node_glob, :] # in A FoR, w.r.t. origin A-G - Rg = np.dot(Cag.T, Ra) # in G FoR, w.r.t. origin A-G - psi = tsstr.psi[ee, node_loc, :] - psi_dot = tsstr.psi_dot[ee, node_loc, :] - Cab = algebra.crv2rotation(psi) - Cba = Cab.T - Cbg = np.dot(Cab.T, Cag) - Tan = algebra.crv2tan(psi) - - track_body = self.settings['track_body'] - isurf = data.aero.struct2aero_mapping[node_glob][0]['i_surf'] - i_n = data.aero.struct2aero_mapping[node_glob][0]['i_n'] - - cas = Cab.dot(tsaero.stability_transform[node_glob]) - - dir_span, span, dir_chord, chord = aeroutils.span_chord(i_n, tsaero.zeta[isurf]) - - # Define the relative velocity and its direction - urel, dir_urel = aeroutils.magnitude_and_direction_of_relative_velocity(tsstr.pos[node_glob, :], - tsstr.pos_dot[node_glob, :], - tsstr.for_vel[:], - Cga, - tsaero.u_ext[isurf][:, :, i_n]) - coef = 0.5 * 1.225 * np.linalg.norm(urel) - cla = np.zeros((3, 3)) - cla[0, :] = 0.5 - cla[2, :] = 18 - u_rel_norm = np.linalg.norm(urel) - ldir = np.zeros((3, 3)) - ldir[:, 2] = 1 - if bc_here != 1: - polar_matrix_vel[np.ix_(jj_tra, jj_tra)] -= cas.dot(coef * cla * ldir).dot(cas.T) # eta_dot - polar_matrix_vel[np.ix_(jj_tra, jj_for_tra)] -= cas.dot(coef * cla * ldir).dot(cas.T) # va - polar_matrix_vel[np.ix_(jj_tra, jj_for_rot)] += cas.dot(coef * cla * ldir).dot(cas.T).dot(algebra.skew(Ra)) # eta * Ra - - polar_matrix_vel[np.ix_(jj_for_tra, jj_tra)] -= cas.dot(coef * cla * ldir).dot(cas.T) # eta_dot - polar_matrix_vel[np.ix_(jj_for_tra, jj_for_tra)] -= cas.dot(coef * cla * ldir).dot(cas.T) # va - polar_matrix_vel[np.ix_(jj_for_tra, jj_for_rot)] += cas.dot(coef * cla * ldir).dot(cas.T).dot(algebra.skew(Ra)) # eta * Ra - # polar_matrix_disp[np.ix_(jj_for_tra, jj_euler)] += cas.dot(coef * cla * ldir).dot(cas.T) - - - - import pdb; pdb.set_trace() - return polar_matrix_disp, polar_matrix_vel From 852b7e656d4cfa672bd756acc46116979610e5c6 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Thu, 30 Jun 2022 14:30:28 +0200 Subject: [PATCH 18/19] Fix differencing scheme --- sharpy/aero/utils/airfoilpolars.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sharpy/aero/utils/airfoilpolars.py b/sharpy/aero/utils/airfoilpolars.py index df61cff8a..d64aa340d 100644 --- a/sharpy/aero/utils/airfoilpolars.py +++ b/sharpy/aero/utils/airfoilpolars.py @@ -75,7 +75,7 @@ def get_derivatives_at_aoa(self, aoa): p = np.array(self.get_coefs(aoa)) m1 = np.array(self.get_coefs(aoa - delta_aoa)) - der = (3 * p1 - 4 * p + m1) / 2 / delta_aoa + der = (p1 - m1) / 2 / delta_aoa return der From a8db20a26884039f019bec74b0f4dc4ad9028ee9 Mon Sep 17 00:00:00 2001 From: Norberto Goizueta Date: Thu, 30 Jun 2022 16:00:08 +0200 Subject: [PATCH 19/19] Try test with improved precision --- tests/uvlm/static/polars/generate_wing.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/tests/uvlm/static/polars/generate_wing.py b/tests/uvlm/static/polars/generate_wing.py index 00141e103..79d942054 100644 --- a/tests/uvlm/static/polars/generate_wing.py +++ b/tests/uvlm/static/polars/generate_wing.py @@ -3,6 +3,7 @@ import sharpy.utils.algebra as algebra import sharpy.sharpy_main +np.set_printoptions(precision=16) def generate_infinite_wing(case_name, alpha, **kwargs): @@ -84,7 +85,7 @@ def generate_infinite_wing(case_name, alpha, **kwargs): settings['StaticUvlm'] = {'print_info': 'on', 'horseshoe': 'on', 'num_cores': 4, - 'vortex_radius': 1e-6, + 'vortex_radius': 1e-9, 'velocity_field_generator': 'SteadyVelocityField', 'velocity_field_input': {'u_inf': u_inf, 'u_inf_direction': u_inf_direction}, @@ -124,7 +125,7 @@ def generate_infinite_wing(case_name, alpha, **kwargs): 'save_data': 'off', 'continuous_eigenvalues': 'off', 'plot_eigenvalues': False, - 'rigid_modes_ppal_axes': 'on', + 'rigid_modes_ppal_axes': 'off', } # ROM settings @@ -137,10 +138,8 @@ def generate_infinite_wing(case_name, alpha, **kwargs): settings['LinearAssembler'] = {'linear_system': 'LinearAeroelastic', 'linear_system_settings': { 'beam_settings': {'modal_projection': 'off', - 'inout_coords': 'modes', + 'inout_coords': 'nodes', 'discrete_time': 'off', - 'newmark_damp': 0.5e-3, - 'discr_method': 'newmark', 'dt': wing.dt, 'proj_modes': 'undamped', 'use_euler': 'on', @@ -154,7 +153,7 @@ def generate_infinite_wing(case_name, alpha, **kwargs): 'density': rho, 'remove_predictor': 'off', 'use_sparse': False, - 'vortex_radius': 1e-7, + 'vortex_radius': 1e-9, 'remove_inputs': ['u_gust'], 'convert_to_ct': 'on', },