Skip to content

Commit

Permalink
Merge branch 'main' into dev_pypi
Browse files Browse the repository at this point in the history
  • Loading branch information
@kccwing
kccwing authored Oct 8, 2024
2 parents 904e247 + e231805 commit 923bbc6
Show file tree
Hide file tree
Showing 3 changed files with 76 additions and 159 deletions.
105 changes: 63 additions & 42 deletions sharpy/postproc/aeroforcescalculator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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 = \
Expand All @@ -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
"""
Expand All @@ -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,
Expand All @@ -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:
Expand All @@ -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]
Expand All @@ -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, '
Expand Down
113 changes: 4 additions & 109 deletions sharpy/utils/datastructures.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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()

Expand Down Expand Up @@ -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

Expand All @@ -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.
Expand Down
Loading

0 comments on commit 923bbc6

Please sign in to comment.