Skip to content

Commit

Permalink
Merge branch 'feature-tiltrotor_transition' of https://github.com/sua…
Browse files Browse the repository at this point in the history
…vecode/SUAVE into feature-tiltrotor_transition
  • Loading branch information
rachealerhard committed Apr 13, 2022
2 parents a36f623 + b7d1515 commit 7eee908
Show file tree
Hide file tree
Showing 4 changed files with 139 additions and 37 deletions.
2 changes: 2 additions & 0 deletions regression/scripts/segments/segment_test.py
Expand Up @@ -485,6 +485,8 @@ def mission_setup(analyses):
segment.throttle = 0.6
segment.distance = 500 * Units.km
segment.state.numerics.number_control_points = 16
segment.state.unknowns.accel_x = -1. * ones_row(1)
segment.state.unknowns.time = 10.

# add to misison
mission.append_segment(segment)
Expand Down
Expand Up @@ -5,6 +5,7 @@
# Modified: Feb 2016, Andrew Wendorff
# Mar 2020, M. Clarke
# Aug 2021, R. Erhard
# Apr 2022, A. Blaufox

# ----------------------------------------------------------------------
# Imports
Expand Down Expand Up @@ -78,7 +79,7 @@ def __defaults__(self):
# initials and unknowns
ones_row = self.state.ones_row
self.state.unknowns.body_angle = ones_row(1) * 0.0
self.state.unknowns.velocity_x = ones_row(1) * 0.0
self.state.unknowns.accel_x = ones_row(1) * 1.
self.state.unknowns.time = 100.
self.state.residuals.final_velocity_error = 0.0
self.state.residuals.forces = ones_row(2) * 0.0
Expand Down Expand Up @@ -122,6 +123,7 @@ def __defaults__(self):
# Update Conditions
iterate.conditions = Process()
iterate.conditions.differentials = Methods.Common.Numerics.update_differentials_time
iterate.conditions.velocity = Methods.Cruise.Constant_Throttle_Constant_Altitude.integrate_velocity
iterate.conditions.altitude = Methods.Common.Aerodynamics.update_altitude
iterate.conditions.atmosphere = Methods.Common.Aerodynamics.update_atmosphere
iterate.conditions.gravity = Methods.Common.Weights.update_gravity
Expand Down
110 changes: 98 additions & 12 deletions trunk/SUAVE/Input_Output/OpenVSP/vsp_wing.py
Expand Up @@ -287,7 +287,7 @@ def read_vsp_wing(wing_id, units_type='SI', write_airfoil_file=True, use_scaling
x_sec_1_taper_parm = vsp.GetXSecParm(x_sec_1,'Taper')
x_sec_1_rc_parm = vsp.GetXSecParm(x_sec_1,'Root_Chord')
x_sec_1_tc_parm = vsp.GetXSecParm(x_sec_1,'Tip_Chord')
x_sec_1_t_parm = vsp.GetXSecParm(x_sec_1,'ThickChord')
x_sec_1_t_parm = vsp.GetXSecParm(x_sec_1,'ThickChord')

# Calcs
sweep = vsp.GetParmVal(x_sec_1_sweep_parm) * Units.deg
Expand Down Expand Up @@ -315,9 +315,25 @@ def read_vsp_wing(wing_id, units_type='SI', write_airfoil_file=True, use_scaling
# check if control surface (sub surfaces) are defined
tags = []
LE_flags = []
span_fraction_starts = []
span_fraction_ends = []
U_starts = []
U_ends = []
chord_fractions = []

# determine the number of segments and where the breaks are
if len(wing.Segments.keys())>0:
N = len(wing.Segments.keys())-1

# Do a for loop to find the location of breaks
breaks = []
for seg in wing.Segments:
breaks.append(seg.percent_span_location)

else:
N = 1
# Location of breaks
breaks = [0.,1.]



num_cs = vsp.GetNumSubSurf(wing_id)

Expand All @@ -330,17 +346,17 @@ def read_vsp_wing(wing_id, units_type='SI', write_airfoil_file=True, use_scaling
if 'LE_Flag' == vsp.GetParmName(param_names[p_idx]):
LE_flags.append(vsp.GetParmVal(param_names[p_idx]))
if 'UStart' == vsp.GetParmName(param_names[p_idx]):
span_fraction_starts.append(vsp.GetParmVal(param_names[p_idx]))
U_starts.append(vsp.GetParmVal(param_names[p_idx]))
if 'UEnd' == vsp.GetParmName(param_names[p_idx]):
span_fraction_ends.append(vsp.GetParmVal(param_names[p_idx]))
U_ends.append(vsp.GetParmVal(param_names[p_idx]))
if 'Length_C_Start' == vsp.GetParmName(param_names[p_idx]):
chord_fractions.append(vsp.GetParmVal(param_names[p_idx]))

# assign control surface parameters to wings. Outer most control surface on main/horizontal wing is assigned a aileron
for cs_idx in range(num_cs):
aileron_present = False
if num_cs > 1:
aileron_loc = np.argmax(np.array(span_fraction_starts))
aileron_loc = np.argmax(np.array(U_starts))
if cs_idx == aileron_loc:
aileron_present = True
if LE_flags[cs_idx] == 1.0:
Expand All @@ -354,8 +370,31 @@ def read_vsp_wing(wing_id, units_type='SI', write_airfoil_file=True, use_scaling
else:
CS = SUAVE.Components.Wings.Control_Surfaces.Flap()
CS.tag = tags[cs_idx]
CS.span_fraction_start = np.maximum((span_fraction_starts[cs_idx] * (segment_num + 1) - 1) / (segment_num - 1), 0)
CS.span_fraction_end = np.minimum((span_fraction_ends[cs_idx] * (segment_num + 1) - 1) / (segment_num - 1), 1)

# Do some math to get U into a span fraction
U_scale = 1/(N+2) # U scaling

# Determine which segment the control surfaces begin and end, use floor
segment_start = int(np.floor(U_starts[cs_idx]/U_scale)) - 1
segment_end = int(np.floor(U_ends[cs_idx]/U_scale)) - 1

segment_normalized_start = U_starts[cs_idx]/U_scale - (segment_start+1)
segment_normalized_end = U_ends[cs_idx]/U_scale - (segment_end+1)

# calculate segment spans
start_span = breaks[segment_start+1] - breaks[segment_start]
end_span = breaks[segment_end+1] - breaks[segment_end]

# Calculate the offsets
start_offset = breaks[segment_start]
end_offset = breaks[segment_end]

# Calculate the end points
span_start = segment_normalized_start * start_span + start_offset
span_end = segment_normalized_end * end_span + end_offset

CS.span_fraction_start = np.max([span_start, 0])
CS.span_fraction_end = np.min([span_end,1])
if CS.span_fraction_start > 1 or CS.span_fraction_end < 0:
raise AssertionError("SUAVE import of VSP files does not allow control surfaces defined for the wing caps.")

Expand Down Expand Up @@ -655,7 +694,7 @@ def write_vsp_wing(vehicle,wing, area_tags, fuel_tank_set_ind, OML_set_ind):

if 'control_surfaces' in wing:
for ctrl_surf in wing.control_surfaces:
write_vsp_control_surface(wing_id,ctrl_surf,n_segments)
write_vsp_control_surface(wing,wing_id,ctrl_surf)


if 'Fuel_Tanks' in wing:
Expand All @@ -668,7 +707,7 @@ def write_vsp_wing(vehicle,wing, area_tags, fuel_tank_set_ind, OML_set_ind):


## @ingroup Input_Output-OpenVSP
def write_vsp_control_surface(wing_id,ctrl_surf,n_segments):
def write_vsp_control_surface(wing,wing_id,ctrl_surf):
"""This writes a control surface in a wing.
Assumptions:
Expand All @@ -688,6 +727,27 @@ def write_vsp_control_surface(wing_id,ctrl_surf,n_segments):
Properties Used:
N/A
"""

# determine the number of segments and where the breaks are
if len(wing.Segments.keys())>0:
N = len(wing.Segments.keys())-1

# Do a for loop to find the location of breaks
breaks = []
for seg in wing.Segments:
breaks.append(seg.percent_span_location)

else:
N = 1
# Location of breaks
breaks = [0.,1.]

breaks = np.array(breaks)

span_start = ctrl_surf.span_fraction_start
span_end = ctrl_surf.span_fraction_end


cs_id = vsp.AddSubSurf( wing_id, vsp.SS_CONTROL)
param_names = vsp.GetSubSurfParmIDs(cs_id)
for p_idx in range(len(param_names)):
Expand All @@ -696,10 +756,36 @@ def write_vsp_control_surface(wing_id,ctrl_surf,n_segments):
vsp.SetParmVal(param_names[p_idx], 1.0)
else:
vsp.SetParmVal( param_names[p_idx], 0.0)

# U values are only linear about the section in which they are defined

# Identify the segments they start and end with, if there are no segments there is one segment
segment_start = np.where(breaks<span_start)[0][-1]
segment_end = np.where(breaks>=span_end)[0][0] -1

# Find where in the segment it starts and ends (normalized)
start_span = breaks[segment_start+1] - breaks[segment_start]
end_span = breaks[segment_end+1] - breaks[segment_end]

start_offset = breaks[segment_start]
end_offset = breaks[segment_end]

segment_normalized_start = (span_start - start_offset)/start_span
segment_normalized_end = (span_end - end_offset )/end_span

# Find the scaling for those segments
U_scale = 1/(N + 2)

# Used the normalized positions and the U_scaling to find the U location
U_start = U_scale*segment_normalized_start + (segment_start + 1) * U_scale
U_end = U_scale*segment_normalized_end + (segment_end + 1) * U_scale

# Voila!

if 'UStart' == vsp.GetParmName(param_names[p_idx]):
vsp.SetParmVal(param_names[p_idx], (ctrl_surf.span_fraction_start*(n_segments-1)+1)/(n_segments+1))
vsp.SetParmVal(param_names[p_idx], U_start)
if 'UEnd' ==vsp.GetParmName(param_names[p_idx]):
vsp.SetParmVal(param_names[p_idx], (ctrl_surf.span_fraction_end*(n_segments-1)+1)/(n_segments+1))
vsp.SetParmVal(param_names[p_idx], U_end)
if 'Length_C_Start' == vsp.GetParmName(param_names[p_idx]):
vsp.SetParmVal(param_names[p_idx], ctrl_surf.chord_fraction)
if 'Length_C_End' == vsp.GetParmName(param_names[p_idx]):
Expand Down
Expand Up @@ -5,14 +5,13 @@
# Modified: Jan 2016, E. Botero
# May 2019, T. MacDonald
# Mar 2020, M. Clarke
# Apr 2022, A. Blaufox

# ----------------------------------------------------------------------
# Imports
# ----------------------------------------------------------------------

import numpy as np
from SUAVE.Methods.Geometry.Three_Dimensional \
import angles_to_dcms, orientation_product, orientation_transpose

# ----------------------------------------------------------------------
# Unpack Unknowns
Expand All @@ -23,27 +22,49 @@ def unpack_unknowns(segment):

# unpack unknowns
unknowns = segment.state.unknowns
velocity_x = unknowns.velocity_x
time = unknowns.time
theta = unknowns.body_angle
accel_x = unknowns.accel_x
time = unknowns.time

# unpack givens
v0 = segment.air_speed_start
vf = segment.air_speed_end
# rescale time
t_initial = segment.state.conditions.frames.inertial.time[0,0]
t_nondim = segment.state.numerics.dimensionless.control_points

# time
t_final = t_initial + time
t_nondim = segment.state.numerics.dimensionless.control_points
time = t_nondim * (t_final-t_initial) + t_initial

#apply unknowns
# build acceleration
N = segment.state.numerics.number_control_points
a = np.zeros((N, 3))
a[:, 0] = accel_x[:,0]

# apply unknowns
conditions = segment.state.conditions
conditions.frames.inertial.velocity_vector[:,0] = velocity_x
conditions.frames.inertial.velocity_vector[0,0] = v0
conditions.frames.body.inertial_rotations[:,1] = theta[:,0]
conditions.frames.inertial.acceleration_vector = a
conditions.frames.inertial.time[:,0] = time[:,0]

return

# ----------------------------------------------------------------------
# Integrate Velocity
# ----------------------------------------------------------------------

def integrate_velocity(segment):

# unpack
conditions = segment.state.conditions
v0 = segment.air_speed_start
I = segment.state.numerics.time.integrate
a = conditions.frames.inertial.acceleration_vector

# compute x-velocity
velocity_x = v0 + np.dot(I, a)[:,0]

# pack velocity
conditions.frames.inertial.velocity_vector[:,0] = velocity_x

return

# ----------------------------------------------------------------------
# Initialize Conditions
# ----------------------------------------------------------------------
Expand Down Expand Up @@ -104,9 +125,6 @@ def initialize_conditions(segment):
segment.air_speed_start = v0
segment.air_speed_end = vf

# Initialize the x velocity unknowns to speed convergence:
segment.state.unknowns.velocity_x = np.linspace(v0,vf,N)

# pack conditions
segment.state.conditions.propulsion.throttle[:,0] = throttle
segment.state.conditions.freestream.altitude[:,0] = alt
Expand Down Expand Up @@ -148,17 +166,11 @@ def solve_residuals(segment):
FT = conditions.frames.inertial.total_force_vector
vf = segment.air_speed_end
v = conditions.frames.inertial.velocity_vector
D = segment.state.numerics.time.differentiate
m = conditions.weights.total_mass

# process and pack
acceleration = np.dot(D , v)
conditions.frames.inertial.acceleration_vector = acceleration

a = segment.state.conditions.frames.inertial.acceleration_vector
a = conditions.frames.inertial.acceleration_vector

segment.state.residuals.forces[:,0] = FT[:,0]/m[:,0] - a[:,0]
segment.state.residuals.forces[:,1] = FT[:,2]/m[:,0] #- a[:,2]
segment.state.residuals.forces[:,1] = FT[:,2]/m[:,0]
segment.state.residuals.final_velocity_error = (v[-1,0] - vf)

return

0 comments on commit 7eee908

Please sign in to comment.