Approximate linear lift curve slope with polar corrections

sharpy/aero/utils/airfoilpolars.py

@@ -59,6 +59,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

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:
@@ -132,6 +135,22 @@ def initialise(self, data):
 
         self.get_gebm2uvlm_gains(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"""
         Assembly of the linearised aeroelastic system.
@@ -227,6 +246,16 @@ def assemble(self):
             self.couplings['Ksa'] = gain_ksa
             self.couplings['Kas'] = gain_kas
 
+            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':
                 # Project UVLM onto modal space
@@ -816,6 +845,7 @@ 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 and Total moments
                     # total forces
                     Kforces[np.ix_(jj_for_tra, ii_vert)] += Cag
 
@@ -996,3 +1026,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

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,7 +39,13 @@ 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']
+        # 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']
+        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:
@@ -55,20 +67,20 @@ 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)
             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
+        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].real
@@ -154,7 +166,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]))
 

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'
@@ -126,6 +135,7 @@ def run(self, online=False):
             self.data.linear.linear_system.update(self.settings['u_inf'])
 
         for target_system in ['aerodynamic', 'aeroelastic']:
+            cout.cout_wrap('-------- {:s} SYSTEM DERIVATIVES ---------'.format(target_system.upper()))
             state_space = self.get_state_space(target_system)
             target_system_derivatives = derivatives[target_system]
 
@@ -193,6 +203,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':
@@ -203,13 +222,10 @@ def get_state_space(self, target_system):
         return ss
 
     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)
+        """Retrieve steady aerodynamic forces and moments at the linearisation reference at the
 
-        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)
+        Returns:
+            tuple: (fx, fy, fz, mx, my, mz) in the inertial G frame
+        """
 
-        return fx, fy, fz
+        return self.data.linear.tsaero0.total_steady_inertial_forces

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,