Merge pull request #301 from ImperialCollegeLondon/fix_aeroforces_for

Fix AeroForcesCalculator giving forces in wrong frame of reference

sharpy/postproc/aeroforcescalculator.py

@@ -63,6 +63,10 @@ class AeroForcesCalculator(BaseSolver):
     settings_default['c_ref'] = 1
     settings_description['c_ref'] = 'Reference chord'
 
+    settings_types['u_inf_dir'] = 'list(float)'
+    settings_default['u_inf_dir'] = [1., 0., 0.]
+    settings_description['u_inf_dir'] = 'Flow direction'
+
     settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
@@ -118,55 +122,73 @@ def run(self, **kwargs):
         if self.settings['write_text_file']:
             self.file_output(self.settings['text_file_name'])
         return self.data
- 
+
     def calculate_forces(self, ts):
-        # Forces per surface in G frame
-        self.rot = algebra.quat2rotation(self.data.structure.timestep_info[ts].quat)
+        self.rot = algebra.quat2rotation(self.data.structure.timestep_info[ts].quat)  # R_GA
+
+        # flow rotation angle from x and z components of flow direction
+        # WARNING: this will give incorrect results when there is sideslip
+        rot_xy = np.arctan2(-self.settings['u_inf_dir'][2], self.settings['u_inf_dir'][0])
+        rmat_xy = algebra.euler2rot((0., rot_xy, 0.))
+
+        # Forces per surface in A frame relative to the flow
         force = self.data.aero.timestep_info[ts].forces
         unsteady_force = self.data.aero.timestep_info[ts].dynamic_forces
-        for i_surf in range(self.data.aero.n_surf):
+
+        n_surf = self.data.aero.n_surf
+
+        rmat = np.block([[self.rot @ rmat_xy.T, np.zeros((3, 3))], [np.zeros((3, 3)), self.rot @ rmat_xy.T]])
+        force_g = [np.moveaxis(np.squeeze(rmat @ np.expand_dims(np.moveaxis(
+            force[i], (0, 1, 2), (2, 0, 1)), -1)), source=(0, 1, 2), destination=(1, 2, 0))
+                   for i in range(n_surf)]
+        unsteady_force_g = [np.moveaxis(np.squeeze(rmat @ np.expand_dims(np.moveaxis(
+            unsteady_force[i], (0, 1, 2), (2, 0, 1)), -1)), source=(0, 1, 2), destination=(1, 2, 0))
+                   for i in range(n_surf)]
+
+        for i_surf in range(n_surf):
             (
-            self.data.aero.timestep_info[ts].inertial_steady_forces[i_surf, 0:3], 
-            self.data.aero.timestep_info[ts].inertial_unsteady_forces[i_surf, 0:3],
-            self.data.aero.timestep_info[ts].body_steady_forces[i_surf, 0:3],
-            self.data.aero.timestep_info[ts].body_unsteady_forces[i_surf, 0:3]
-            ) = self.calculate_forces_for_isurf_in_g_frame(force[i_surf], unsteady_force=unsteady_force[i_surf])
-        
+                self.data.aero.timestep_info[ts].inertial_steady_forces[i_surf, 0:3],
+                self.data.aero.timestep_info[ts].inertial_unsteady_forces[i_surf, 0:3],
+                self.data.aero.timestep_info[ts].body_steady_forces[i_surf, 0:3],
+                self.data.aero.timestep_info[ts].body_unsteady_forces[i_surf, 0:3]
+            ) = self.calculate_forces_for_isurf_in_g_frame(force_g[i_surf], unsteady_force=unsteady_force_g[i_surf])
+
         if self.settings["nonlifting_body"]:
             for i_surf in range(self.data.nonlifting_body.n_surf):
                 (
-                self.data.nonlifting_body.timestep_info[ts].inertial_steady_forces[i_surf, 0:3], 
-                self.data.nonlifting_body.timestep_info[ts].body_steady_forces[i_surf, 0:3],
-                ) = self.calculate_forces_for_isurf_in_g_frame(self.data.nonlifting_body.timestep_info[ts].forces[i_surf], nonlifting=True)
-
+                    self.data.nonlifting_body.timestep_info[ts].inertial_steady_forces[i_surf, 0:3],
+                    self.data.nonlifting_body.timestep_info[ts].body_steady_forces[i_surf, 0:3],
+                ) = self.calculate_forces_for_isurf_in_g_frame(
+                    self.data.nonlifting_body.timestep_info[ts].forces[i_surf], nonlifting=True)
+
         # Convert to forces in B frame
         try:
             steady_forces_b = self.data.structure.timestep_info[ts].postproc_node['aero_steady_forces']
         except KeyError:
             if self.settings["nonlifting_body"]:
-                warnings.warn('Nonlifting forces are not considered in aero forces calculation since forces cannot not be retrieved from postproc node.')
-            steady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, force)
+                warnings.warn(
+                    'Nonlifting forces are not considered in aero forces calculation since forces cannot not be retrieved from postproc node.')
+            steady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, force_g)
 
         try:
             unsteady_forces_b = self.data.structure.timestep_info[ts].postproc_node['aero_unsteady_forces']
         except KeyError:
-            unsteady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, unsteady_force)
+            unsteady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, unsteady_force_g)
         # Convert to forces in A frame
         steady_forces_a = self.data.structure.timestep_info[ts].nodal_b_for_2_a_for(steady_forces_b,
                                                                                     self.data.structure)
 
         unsteady_forces_a = self.data.structure.timestep_info[ts].nodal_b_for_2_a_for(unsteady_forces_b,
                                                                                       self.data.structure)
 
-        # Express total forces in A frame
         self.data.aero.timestep_info[ts].total_steady_body_forces = \
-            mapping.total_forces_moments(steady_forces_a,
-                                         self.data.structure.timestep_info[ts].pos,
-                                         ref_pos=self.moment_reference_location)
+            rmat @ mapping.total_forces_moments(steady_forces_a,
+                                                self.data.structure.timestep_info[ts].pos,
+                                                ref_pos=self.moment_reference_location)
         self.data.aero.timestep_info[ts].total_unsteady_body_forces = \
-            mapping.total_forces_moments(unsteady_forces_a,
-                                         self.data.structure.timestep_info[ts].pos,
-                                         ref_pos=self.moment_reference_location)
+            rmat @ mapping.total_forces_moments(unsteady_forces_a,
+                                                self.data.structure.timestep_info[ts].pos,
+                                                ref_pos=self.moment_reference_location)
 
         # Express total forces in G frame
         self.data.aero.timestep_info[ts].total_steady_inertial_forces = \
@@ -179,7 +201,7 @@ def calculate_forces(self, ts):
                       [np.zeros((3, 3)), self.rot]]).dot(
                 self.data.aero.timestep_info[ts].total_unsteady_body_forces)
 
-    def calculate_forces_for_isurf_in_g_frame(self, force, unsteady_force = None, nonlifting = False):
+    def calculate_forces_for_isurf_in_g_frame(self, force, unsteady_force=None, nonlifting=False):
         """
             Forces for a surface in G frame
         """
@@ -193,11 +215,11 @@ def calculate_forces_for_isurf_in_g_frame(self, force, unsteady_force = None, no
                 if not nonlifting:
                     total_unsteady_force += unsteady_force[0:3, i_m, i_n]
         if not nonlifting:
-            return total_steady_force, total_unsteady_force, np.dot(self.rot.T, total_steady_force), np.dot(self.rot.T, total_unsteady_force)
+            return total_steady_force, total_unsteady_force, np.dot(self.rot.T, total_steady_force), np.dot(self.rot.T,
+                                                                                                            total_unsteady_force)
         else:
             return total_steady_force, np.dot(self.rot.T, total_steady_force)
 
-
     def map_forces_beam_dof(self, aero_data, ts, force):
         struct_tstep = self.data.structure.timestep_info[ts]
         aero_forces_beam_dof = mapping.aero2struct_force_mapping(force,
@@ -212,19 +234,19 @@ def map_forces_beam_dof(self, aero_data, ts, force):
 
     def calculate_coefficients(self, fx, fy, fz, mx, my, mz):
         qS = self.settings['q_ref'] * self.settings['S_ref']
-        return fx/qS, fy/qS, fz/qS, mx/qS/self.settings['b_ref'], my/qS/self.settings['c_ref'], \
-               mz/qS/self.settings['b_ref']
+        return fx / qS, fy / qS, fz / qS, mx / qS / self.settings['b_ref'], my / qS / self.settings['c_ref'], \
+               mz / qS / self.settings['b_ref']
 
     def screen_output(self, ts):
         # print time step total aero forces
-        forces  = np.zeros(3)
-        moments  = np.zeros(3)
+        forces = np.zeros(3)
+        moments = np.zeros(3)
 
         aero_tstep = self.data.aero.timestep_info[ts]
         forces += aero_tstep.total_steady_inertial_forces[:3] + aero_tstep.total_unsteady_inertial_forces[:3]
         moments += aero_tstep.total_steady_inertial_forces[3:] + aero_tstep.total_unsteady_inertial_forces[3:]
 
-        if self.settings['coefficients']: # TODO: Check if coefficients have to be computed differently for fuselages
+        if self.settings['coefficients']:  # TODO: Check if coefficients have to be computed differently for fuselages
             Cfx, Cfy, Cfz, Cmx, Cmy, Cmz = self.calculate_coefficients(*forces, *moments)
             self.table.print_line([ts, Cfx, Cfy, Cfz, Cmx, Cmy, Cmz])
         else:
@@ -233,7 +255,7 @@ def screen_output(self, ts):
     def file_output(self, filename):
         # assemble forces/moments matrix
         # (1 timestep) + (3+3 inertial steady+unsteady) + (3+3 body steady+unsteady)
-        force_matrix = np.zeros((self.ts_max, 1 + 3 + 3 + 3 + 3 ))
+        force_matrix = np.zeros((self.ts_max, 1 + 3 + 3 + 3 + 3))
         moment_matrix = np.zeros((self.ts_max, 1 + 3 + 3 + 3 + 3))
         for ts in range(self.ts_max):
             aero_tstep = self.data.aero.timestep_info[ts]
@@ -243,24 +265,23 @@ def file_output(self, filename):
             i += 1
 
             # Steady forces/moments G
-            force_matrix[ts, i:i+3] = aero_tstep.total_steady_inertial_forces[:3]
-            moment_matrix[ts, i:i+3] = aero_tstep.total_steady_inertial_forces[3:]
+            force_matrix[ts, i:i + 3] = aero_tstep.total_steady_inertial_forces[:3]
+            moment_matrix[ts, i:i + 3] = aero_tstep.total_steady_inertial_forces[3:]
             i += 3
 
             # Unsteady forces/moments G
-            force_matrix[ts, i:i+3] = aero_tstep.total_unsteady_inertial_forces[:3]
-            moment_matrix[ts, i:i+3] = aero_tstep.total_unsteady_inertial_forces[3:]
+            force_matrix[ts, i:i + 3] = aero_tstep.total_unsteady_inertial_forces[:3]
+            moment_matrix[ts, i:i + 3] = aero_tstep.total_unsteady_inertial_forces[3:]
             i += 3
 
             # Steady forces/moments A
-            force_matrix[ts, i:i+3] = aero_tstep.total_steady_body_forces[:3]
-            moment_matrix[ts, i:i+3] = aero_tstep.total_steady_body_forces[3:]
+            force_matrix[ts, i:i + 3] = aero_tstep.total_steady_body_forces[:3]
+            moment_matrix[ts, i:i + 3] = aero_tstep.total_steady_body_forces[3:]
             i += 3
 
             # Unsteady forces/moments A
-            force_matrix[ts, i:i+3] = aero_tstep.total_unsteady_body_forces[:3]
-            moment_matrix[ts, i:i+3] = aero_tstep.total_unsteady_body_forces[3:]
-
+            force_matrix[ts, i:i + 3] = aero_tstep.total_unsteady_body_forces[:3]
+            moment_matrix[ts, i:i + 3] = aero_tstep.total_unsteady_body_forces[3:]
 
         header = ''
         header += 'tstep, '

sharpy/utils/datastructures.py

@@ -431,6 +431,10 @@ class AeroTimeStepInfo(TimeStepInfo):
           ``[num_surf x chordwise panels x spanwise panels]``
         dimensions_star (np.ndarray): Matrix defining the dimensions of the vortex grid on wakes
           ``[num_surf x streamwise panels x spanwise panels]``
+
+    Note, Ben Preston 27/09/24: forces and dynamic forces are stated to be in ``G`` however were actually
+    determined to be in ``A``. This issue may apply to other parameters, however errors due to this were only apparent
+    in the AeroForcesCalculator postprocessor.
     """
     def __init__(self, dimensions, dimensions_star):
         super().__init__(dimensions)
@@ -948,22 +952,6 @@ def change_to_local_AFoR(self, for0_pos, for0_vel, quat0):
         # Modify local rotations
         for ielem in range(self.psi.shape[0]):
             for inode in range(3):
-                # psi_master = self.psi[ielem, inode, :] + np.zeros((3,),)
-                # self.psi[ielem, inode, :] = np.dot(Csm, self.psi[ielem, inode, :])
-                # self.psi_dot[ielem, inode, :] = (np.dot(Csm, self.psi_dot[ielem, inode, :] + algebra.cross3(for0_vel[3:6], psi_master)) -
-                #                                  algebra.multiply_matrices(CAslaveG, algebra.skew(self.for_vel[3:6]), CGAmaster, psi_master))
-
-                # psi_master = self.psi[ielem,inode,:] + np.zeros((3,),)
-                # self.psi[ielem, inode, :] = algebra.rotation2crv(np.dot(Csm,algebra.crv2rotation(self.psi[ielem,inode,:])))
-                # psi_slave = self.psi[ielem, inode, :] + np.zeros((3,),)
-                # cbam = algebra.crv2rotation(psi_master).T
-                # cbas = algebra.crv2rotation(psi_slave).T
-                # tm = algebra.crv2tan(psi_master)
-                # inv_ts = np.linalg.inv(algebra.crv2tan(psi_slave))
-
-                # self.psi_dot[ielem, inode, :] = np.dot(inv_ts, (np.dot(tm, self.psi_dot[ielem, inode, :]) +
-                #                                                 np.dot(cbam, for0_vel[3:6]) -
-                #                                                 np.dot(cbas, self.for_vel[3:6])))
                 self.psi[ielem, inode, :] = self.psi_local[ielem,inode,:].copy()
                 self.psi_dot[ielem, inode, :] = self.psi_dot_local[ielem, inode, :].copy()
 
@@ -997,12 +985,6 @@ def change_to_global_AFoR(self, for0_pos, for0_vel, quat0):
 
         for ielem in range(self.psi.shape[0]):
             for inode in range(3):
-                # psi_slave = self.psi[ielem,inode,:] + np.zeros((3,),)
-                # self.psi[ielem, inode, :] = np.dot(Cms, self.psi[ielem, inode, :])
-                # self.psi_dot[ielem, inode, :] = (np.dot(Cms, self.psi_dot[ielem, inode, :] + algebra.cross3(self.for_vel[3:6], psi_slave)) -
-                #                                  algebra.multiply_matrices(CAmasterG, algebra.skew(self.for0_vel[3:6]), CGAslave, psi_slave))
-
-
                 self.psi_local[ielem, inode, :] = self.psi[ielem, inode, :].copy() # Copy here the result from the structural computation
                 self.psi_dot_local[ielem, inode, :] = self.psi_dot[ielem, inode, :].copy() # Copy here the result from the structural computation
 
@@ -1020,93 +1002,6 @@ def change_to_global_AFoR(self, for0_pos, for0_vel, quat0):
                                                                 np.dot(cbas, self.for_vel[3:6]) -
                                                                 np.dot(cbam, for0_vel[3:6])))
 
-    # def whole_structure_to_local_AFoR(self, beam, compute_psi_local=False):
-    #     """
-    #     Same as change_to_local_AFoR but for a multibody structure
-    #
-    #     Args:
-    #         beam(sharpy.structure.models.beam.Beam): Beam structure of ``PreSharpy``
-    #     """
-    #     if not self.in_global_AFoR:
-    #         raise NotImplementedError("Wrong managing of FoR")
-    #
-    #     self.in_global_AFoR = False
-    #     quat0 = self.quat.astype(dtype=ct.c_double, order='F', copy=True)
-    #     for0_pos = self.for_pos.astype(dtype=ct.c_double, order='F', copy=True)
-    #     for0_vel = self.for_vel.astype(dtype=ct.c_double, order='F', copy=True)
-    #
-    #     MB_beam = [None]*beam.num_bodies
-    #     MB_tstep = [None]*beam.num_bodies
-    #
-    #     for ibody in range(beam.num_bodies):
-    #         MB_beam[ibody] = beam.get_body(ibody = ibody)
-    #         MB_tstep[ibody] = self.get_body(beam, MB_beam[ibody].num_dof, ibody = ibody)
-    #         if compute_psi_local:
-    #             MB_tstep[ibody].compute_psi_local_AFoR(for0_pos, for0_vel, quat0)
-    #         MB_tstep[ibody].change_to_local_AFoR(for0_pos, for0_vel, quat0)
-    #
-    #     first_dof = 0
-    #     for ibody in range(beam.num_bodies):
-    #         # Renaming for clarity
-    #         ibody_elems = MB_beam[ibody].global_elems_num
-    #         ibody_nodes = MB_beam[ibody].global_nodes_num
-    #
-    #         # Merge tstep
-    #         self.pos[ibody_nodes,:] = MB_tstep[ibody].pos.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.psi[ibody_elems,:,:] = MB_tstep[ibody].psi.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.psi_local[ibody_elems,:,:] = MB_tstep[ibody].psi_local.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.gravity_forces[ibody_nodes,:] = MB_tstep[ibody].gravity_forces.astype(dtype=ct.c_double, order='F', copy=True)
-    #
-    #         self.pos_dot[ibody_nodes,:] = MB_tstep[ibody].pos_dot.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.psi_dot[ibody_elems,:,:] = MB_tstep[ibody].psi_dot.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.psi_dot_local[ibody_elems,:,:] = MB_tstep[ibody].psi_dot_local.astype(dtype=ct.c_double, order='F', copy=True)
-    #
-    #         # TODO: Do I need a change in FoR for the following variables? Maybe for the FoR ones.
-    #         # tstep.forces_constraints_nodes[ibody_nodes,:] = MB_tstep[ibody].forces_constraints_nodes.astype(dtype=ct.c_double, order='F', copy=True)
-    #         # tstep.forces_constraints_FoR[ibody, :] = MB_tstep[ibody].forces_constraints_FoR[ibody, :].astype(dtype=ct.c_double, order='F', copy=True)
-
-    # def whole_structure_to_global_AFoR(self, beam):
-    #     """
-    #     Same as change_to_global_AFoR but for a multibody structure
-    #
-    #     Args:
-    #         beam(sharpy.structure.models.beam.Beam): Beam structure of ``PreSharpy``
-    #     """
-    #     if self.in_global_AFoR:
-    #         raise NotImplementedError("Wrong managing of FoR")
-    #
-    #     self.in_global_AFoR = True
-    #
-    #     MB_beam = [None]*beam.num_bodies
-    #     MB_tstep = [None]*beam.num_bodies
-    #     quat0 = self.quat.astype(dtype=ct.c_double, order='F', copy=True)
-    #     for0_pos = self.for_pos.astype(dtype=ct.c_double, order='F', copy=True)
-    #     for0_vel = self.for_vel.astype(dtype=ct.c_double, order='F', copy=True)
-    #
-    #     for ibody in range(beam.num_bodies):
-    #         MB_beam[ibody] = beam.get_body(ibody = ibody)
-    #         MB_tstep[ibody] = self.get_body(beam, MB_beam[ibody].num_dof, ibody = ibody)
-    #         MB_tstep[ibody].change_to_global_AFoR(for0_pos, for0_vel, quat0)
-    #
-    #
-    #     first_dof = 0
-    #     for ibody in range(beam.num_bodies):
-    #         # Renaming for clarity
-    #         ibody_elems = MB_beam[ibody].global_elems_num
-    #         ibody_nodes = MB_beam[ibody].global_nodes_num
-    #
-    #         # Merge tstep
-    #         self.pos[ibody_nodes,:] = MB_tstep[ibody].pos.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.psi[ibody_elems,:,:] = MB_tstep[ibody].psi.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.psi_local[ibody_elems,:,:] = MB_tstep[ibody].psi_local.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.gravity_forces[ibody_nodes,:] = MB_tstep[ibody].gravity_forces.astype(dtype=ct.c_double, order='F',
-    #                                                                                    copy=True)
-    #
-    #         self.pos_dot[ibody_nodes,:] = MB_tstep[ibody].pos_dot.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.psi_dot[ibody_elems,:,:] = MB_tstep[ibody].psi_dot.astype(dtype=ct.c_double, order='F', copy=True)
-    #         self.psi_dot_local[ibody_elems,:,:] = MB_tstep[ibody].psi_dot_local.astype(dtype=ct.c_double, order='F', copy=True)
-
-
     def nodal_b_for_2_a_for(self, nodal, beam, filter=np.array([True]*6), ibody=None):
         """
         Projects a nodal variable from the local, body-attached frame (B) to the reference A frame.