Source code for sharpy.linear.assembler.lincontrolsurfacedeflector

"""
Control surface deflector for linear systems
"""
import sharpy.linear.utils.ss_interface as ss_interface
import numpy as np
import sharpy.utils.algebra as algebra
import sharpy.linear.src.libsparse as libsp
import scipy.sparse as sp
import sharpy.linear.src.libss as libss


[docs]@ss_interface.linear_system class LinControlSurfaceDeflector(object): """ Subsystem that deflects control surfaces for use with linear state space systems. Note: The control surface sign convention is different to the convention in the nonlinear solver. See https://www.github.com/imperialcollegelondon/sharpy/issues/193 for more details. In the linear implementation, the control surface deflection sign convention follows the :math:`x_B` vector, thus, in order to have symmetric control surface deflections, additional inputs may be required compared to the implementation in the nonlinear solver. The current version supports only deflections. Future work will include standalone state-space systems to model physical actuators. """ sys_id = 'LinControlSurfaceDeflector' def __init__(self): # Has the back bone structure for a future actuator model # As of now, it simply maps a deflection onto the aerodynamic grid by means of Kzeta_delta self.n_control_surfaces = 0 self.Kzeta_delta = None # type: np.ndarray self.Kdzeta_ddelta = None # type: np.ndarray self.linuvlm = None # type: sharpy.linear.src.linuvlm.Dynamic self.aero = None # type: sharpy.aero.models.aerogrid.Aerogrid self.structure = None # type: sharpy.structure.models.beam.Beam self.tsaero0 = None self.tsstruct0 = None self.gain_cs = None # type: np.ndarray # input gain to the UVLM self.print_info = False # used for debugging purposes def initialise(self, data, linuvlm): self.n_control_surfaces = data.aero.n_control_surfaces self.linuvlm = linuvlm self.aero = data.aero self.structure = data.structure self.tsaero0 = data.aero.timestep_info[0] self.tsstruct0 = data.structure.timestep_info[0]
[docs] def generate(self): """ Generates a matrix mapping a linear control surface deflection onto the aerodynamic grid. This generates two matrices: * `Kzeta_delta` maps the deflection angle onto displacements. It has as many columns as independent control surfaces. * `Kdzeta_ddelta` maps the deflection rate onto grid velocities. Again, it has as many columns as independent control surfaces. Returns: tuple: Tuple containing `Kzeta_delta` and `Kdzeta_ddelta`. """ # For future development # In hindsight, building this matrix iterating through structural node was a big mistake that # has led to very messy code. Would rework by element and in B frame aero = self.aero structure = self.structure linuvlm = self.linuvlm tsaero0 = self.tsaero0 tsstruct0 = self.tsstruct0 # Find the vertices corresponding to a control surface from beam coordinates to aerogrid data_dict = aero.data_dict n_surf = tsaero0.n_surf n_control_surfaces = self.n_control_surfaces Kdisp = np.zeros((3 * linuvlm.Kzeta, n_control_surfaces)) Kvel = np.zeros((3 * linuvlm.Kzeta, n_control_surfaces)) zeta0 = np.concatenate([tsaero0.zeta[i_surf].reshape(-1, order='C') for i_surf in range(n_surf)]) Cga = algebra.quat2rotation(tsstruct0.quat).T Cag = Cga.T # Initialise these parameters hinge_axis = None # Will be set once per control surface to the hinge axis with_control_surface = False # Will be set to true if the spanwise node contains a control surface for global_node in range(structure.num_node): # Retrieve elements and local nodes to which a single node is attached for i_elem in range(structure.num_elem): if global_node in structure.connectivities[i_elem, :]: i_local_node = np.where(structure.connectivities[i_elem, :] == global_node)[0][0] for_delta = structure.frame_of_reference_delta[i_elem, i_local_node, :] # CRV to transform from G to B frame psi = tsstruct0.psi[i_elem, i_local_node] Cab = algebra.crv2rotation(psi) Cba = Cab.T Cbg = np.dot(Cab.T, Cag) Cgb = Cbg.T # Map onto aerodynamic coordinates. Some nodes may be part of two aerodynamic surfaces. for structure2aero_node in aero.struct2aero_mapping[global_node]: # Retrieve surface and span-wise coordinate i_surf, i_node_span = structure2aero_node['i_surf'], structure2aero_node['i_n'] # Although a node may be part of 2 aerodynamic surfaces, we need to ensure that the current # element for the given node is indeed part of that surface. elems_in_surf = np.where(data_dict['surface_distribution'] == i_surf)[0] if i_elem not in elems_in_surf: continue # Surface panelling M = aero.dimensions[i_surf][0] N = aero.dimensions[i_surf][1] K_zeta_start = 3 * sum(linuvlm.MS.KKzeta[:i_surf]) shape_zeta = (3, M + 1, N + 1) i_control_surface = data_dict['control_surface'][i_elem, i_local_node] if i_control_surface >= 0: if not with_control_surface: i_start_of_cs = i_node_span.copy() with_control_surface = True control_surface_chord = data_dict['control_surface_chord'][i_control_surface] try: control_surface_hinge_coord = \ data_dict['control_surface_hinge_coord'][i_control_surface] * \ data_dict['chord'][i_elem, i_local_node] except KeyError: control_surface_hinge_coord = None i_node_hinge = M - control_surface_chord i_vertex_hinge = [K_zeta_start + np.ravel_multi_index((i_axis, i_node_hinge, i_node_span), shape_zeta) for i_axis in range(3)] i_vertex_next_hinge = [K_zeta_start + np.ravel_multi_index((i_axis, i_node_hinge, i_start_of_cs + 1), shape_zeta) for i_axis in range(3)] if control_surface_hinge_coord is not None and M == control_surface_chord: # fully articulated control surface zeta_hinge = Cgb.dot(Cba.dot(tsstruct0.pos[global_node]) + for_delta * np.array([0, control_surface_hinge_coord, 0])) zeta_next_hinge = Cgb.dot(Cbg.dot(zeta_hinge) + np.array([1, 0, 0])) # parallel to the x_b vector else: zeta_hinge = zeta0[i_vertex_hinge] zeta_next_hinge = zeta0[i_vertex_next_hinge] if hinge_axis is None: # Hinge axis not yet set for current control surface # Hinge axis is in G frame hinge_axis = zeta_next_hinge - zeta_hinge hinge_axis = hinge_axis / np.linalg.norm(hinge_axis) for i_node_chord in range(M + 1): i_vertex = [K_zeta_start + np.ravel_multi_index((i_axis, i_node_chord, i_node_span), shape_zeta) for i_axis in range(3)] if i_node_chord >= i_node_hinge: # Zeta in G frame zeta_node = zeta0[i_vertex] # Gframe chord_vec = (zeta_node - zeta_hinge) if self.print_info: print(f'i_node = {global_node}') print(f'i_node_chord = {i_node_chord}') print(f'zeta_node = {zeta_node}') print(f'chord_vec = {chord_vec}') print(f'zeta_hinge = {zeta_hinge}') print(f'zeta_next_hinge = {zeta_next_hinge}') print(f'hinge_axis = {hinge_axis}') # Flap displacement Kdisp[i_vertex, i_control_surface] = \ der_R_arbitrary_axis_times_v(hinge_axis, 0, chord_vec) # Flap velocity Kvel[i_vertex, i_control_surface] = -algebra.skew(chord_vec).dot( hinge_axis) if self.print_info: print('Matrix entries:') print('Kdisp:') print(Kdisp[i_vertex, i_control_surface]) print('Kvel:') print(Kvel[i_vertex, i_control_surface]) else: with_control_surface = False hinge_axis = None # Reset for next control surface self.Kzeta_delta = Kdisp self.Kdzeta_ddelta = Kvel return Kdisp, Kvel
[docs] def apply(self, ss): """ Applies the control surface deflection to the UVLM state space Args: ss (libss.StateSpace): UVLM state space Returns: libss.StateSpace: UVLM state-space with control surfaces and control surface deflection rate as inputs """ Kzeta_delta, Kdzeta_ddelta = self.generate() n_zeta, n_ctrl_sfc = Kzeta_delta.shape if type(ss.A) is libsp.csc_matrix: gain_cs = sp.eye(ss.inputs, ss.inputs + 2 * self.n_control_surfaces, format='lil') gain_cs[:n_zeta, ss.inputs: ss.inputs + n_ctrl_sfc] = Kzeta_delta gain_cs[n_zeta: 2*n_zeta, ss.inputs + n_ctrl_sfc: ss.inputs + 2 * n_ctrl_sfc] = Kdzeta_ddelta gain_cs = libsp.csc_matrix(gain_cs) else: gain_cs = np.eye(ss.inputs, ss.inputs + 2 * self.n_control_surfaces) gain_cs[:n_zeta, ss.inputs: ss.inputs + n_ctrl_sfc] = Kzeta_delta gain_cs[n_zeta: 2*n_zeta, ss.inputs + n_ctrl_sfc: ss.inputs + 2 * n_ctrl_sfc] = Kdzeta_ddelta control_surface_gain = libss.Gain(gain_cs) in_vars = ss.input_variables.copy() in_vars.append('control_surface_deflection', size=n_ctrl_sfc) in_vars.append('dot_control_surface_deflection', size=n_ctrl_sfc) control_surface_gain.input_variables = in_vars control_surface_gain.output_variables = ss_interface.LinearVector.transform(ss.input_variables, to_type=ss_interface.OutputVariable) self.gain_cs = control_surface_gain return ss
# def generator(): # future feature idea: instead of defining the inputs for the time domain simulations as the whole input vector # etc, we could add a generate() method to these systems that can be called from the LinDynamicSim to apply # the gust and generate the correct input vector. def der_Cx_by_v(delta, v): sd = np.sin(delta) cd = np.cos(delta) v2 = v[1] v3 = v[2] return np.array([0, -v2 * sd - v3 * cd, v2 * cd - v3 * sd]) def der_Cy_by_v(delta, v): s = np.sin(delta) c = np.cos(delta) v1 = v[0] v3 = v[2] return np.array([-s*v1 + v*v3, 0, -c*v1 - s*v3]) def der_R_arbitrary_axis_times_v(u, theta, v): r""" Linearised rotation vector of the vector ``v`` by angle ``theta`` about an arbitrary axis ``u``. The rotation of a vector :math:`\mathbf{v}` about the axis :math:`\mathbf{u}` by an angle :math:`\boldsymbol{\theta}` can be expressed as .. math:: \mathbf{w} = \mathbf{R}(\mathbf{u}, \theta) \mathbf{v}, where :math:`\mathbf{R}` is a :math:`\mathbb{R}^{3\times 3}` matrix. This expression can be linearised for it to be included in the linear solver as .. math:: \delta\mathbf{w} = \frac{\partial}{\partial\theta}\left(\mathbf{R}(\mathbf{u}, \theta_0)\right)\delta\theta The matrix :math:`\mathbf{R}` is .. math:: \mathbf{R} = \begin{bmatrix}\cos \theta +u_{x}^{2}\left(1-\cos \theta \right) & u_{x}u_{y}\left(1-\cos \theta \right)-u_{z}\sin \theta & u_{x}u_{z}\left(1-\cos \theta \right)+u_{y}\sin \theta \\ u_{y}u_{x}\left(1-\cos \theta \right)+u_{z}\sin \theta & \cos \theta +u_{y}^{2}\left(1-\cos \theta \right)& u_{y}u_{z}\left(1-\cos \theta \right)-u_{x}\sin \theta \\ u_{z}u_{x}\left(1-\cos \theta \right)-u_{y}\sin \theta & u_{z}u_{y}\left(1-\cos \theta \right)+u_{x}\sin \theta & \cos \theta +u_{z}^{2}\left(1-\cos \theta \right)\end{bmatrix}, and its linearised expression becomes .. math:: \frac{\partial}{\partial\theta}\left(\mathbf{R}(\mathbf{u}, \theta_0)\right) = \begin{bmatrix} -\sin \theta +u_{x}^{2}\sin \theta \mathbf{v}_1 + u_{x}u_{y}\sin \theta-u_{z} \cos \theta \mathbf{v}_2 + u_{x}u_{z}\sin \theta +u_{y}\cos \theta \mathbf{v}_3 \\ u_{y}u_{x}\sin \theta+u_{z}\cos \theta\mathbf{v}_1 -\sin \theta +u_{y}^{2}\sin \theta\mathbf{v}_2 + u_{y}u_{z}\sin \theta-u_{x}\cos \theta\mathbf{v}_3 \\ u_{z}u_{x}\sin \theta-u_{y}\cos \theta\mathbf{v}_1 + u_{z}u_{y}\sin \theta+u_{x}\cos \theta\mathbf{v}_2 -\sin \theta +u_{z}^{2}\sin\theta\mathbf{v}_3\end{bmatrix}_{\theta=\theta_0} and is of dimension :math:`\mathbb{R}^{3\times 1}`. Args: u (numpy.ndarray): Arbitrary rotation axis theta (float): Rotation angle (radians) v (numpy.ndarray): Vector to rotate Returns: numpy.ndarray: Linearised rotation vector of dimensions :math:`\mathbb{R}^{3\times 1}`. """ u = u / np.linalg.norm(u) c = np.cos(theta) s = np.sin(theta) ux, uy, uz = u v1, v2, v3 = v dR11 = -s + ux ** 2 * s dR12 = ux * uy * s - uz * c dR13 = ux * uz * s + uy * c dR21 = uy * ux * s + uz * c dR22 = -s + uy ** 2 * s dR23 = uy * uz * s - ux * c dR31 = uz * ux * s - uy * c dR32 = uz * uy * s + ux * c dR33 = -s + uz ** 2 dRv = np.zeros((3, )) dRv[0] = dR11 * v1 + dR12 * v2 + dR13 * v3 dRv[1] = dR21 * v1 + dR22 * v2 + dR23 * v3 dRv[2] = dR31 * v1 + dR32 * v2 + dR33 * v3 return dRv