from aerosandbox.geometry.airplane import Airplane
from aerosandbox.performance import OperatingPoint
from aerosandbox.weights import MassProperties
import aerosandbox.numpy as np
[docs]def get_modes(
airplane: Airplane,
op_point: OperatingPoint,
mass_props: MassProperties,
aero,
g=9.81,
):
Q = op_point.dynamic_pressure()
S = airplane.s_ref
c = airplane.c_ref
b = airplane.b_ref
QS = Q * S
m = mass_props.mass
Ixx = mass_props.Ixx
Iyy = mass_props.Iyy
Izz = mass_props.Izz
u0 = op_point.velocity
Cxu = -2 * aero["CD"]
Czu = -2 * aero["CL"]
Cmq = aero["Cmq"]
Cma = aero["Cma"]
X_u = QS / m / u0 * Cxu # Units: 1/sec
# X_w = QS / m / u0 * Cxa
# X_q = QS / m * c / (2 * u0) * Cxq
Z_u = QS / m / u0 * Czu # Units: 1/sec
# Z_w = QS / m / u0 * Cza
# Z_q = QS / m * c / (2 * u0) * Czq
# M_u = QS * c / Iyy / u0 * Cmu
M_w = QS * c / Iyy / u0 * Cma # Units: 1/(meter * sec)
M_q = QS * c / Iyy * c / (2 * u0) * Cmq # Units: 1/sec
modes = {}
def get_mode_info(
sigma,
omega_squared,
):
is_oscillatory = omega_squared > 0
return {
"eigenvalue_real": sigma + np.where(
is_oscillatory,
0,
np.abs(omega_squared + 1e-100) ** 0.5,
),
"eigenvalue_imag": np.where(
is_oscillatory,
np.abs(omega_squared + 1e-100) ** 0.5,
0,
),
}
##### Longitudinal modes
### Phugoid
modes['phugoid'] = { # FVA, Eq. 9.55-9.57
**get_mode_info(
sigma=X_u / 2,
omega_squared=-(X_u ** 2) / 4 - g * Z_u / u0
),
"eigenvalue_imag_approx": 2 ** 0.5 * g / u0,
"damping_ratio_approx" : 2 ** -0.5 * aero["CD"] / aero["CL"],
}
#
# ### Short-period
modes['short_period'] = get_mode_info(
sigma=0.5 * M_q,
omega_squared=-(M_q ** 2) / 4 - u0 * M_w
)
##### Lateral modes
### Roll subsidence
modes['roll_subsidence'] = { # FVA, Eq. 9.63
"eigenvalue_real": (
QS * b ** 2 / (2 * Ixx * u0) * aero["Clp"]
),
"eigenvalue_imag": 0,
"damping_ratio" : 1,
}
### Dutch roll
modes['dutch_roll'] = get_mode_info( # FVA, Eq. 9.68
sigma=(
QS * b ** 2 /
(2 * Izz * u0) *
(aero["Cnr"] + Izz / (m * b ** 2) * aero["CYb"])
),
omega_squared=(
QS * b / Izz *
(
aero["Cnb"] + (
op_point.atmosphere.density() * S * b / (4 * m) *
(aero["CYb"] * aero["Cnr"] - aero["Cnb"] * aero["CYr"])
)
)
)
)
### Spiral
spiral_parameter = (aero["Cnr"] - aero["Cnb"] * aero["Clr"] / aero["Clb"]) # FVA, Eq. 9.66
modes['spiral'] = { # FVA, Eq. 9.66
"eigenvalue_real": (
QS * b ** 2 / (2 * Izz * u0) *
spiral_parameter
),
"eigenvalue_imag": 0,
}
### Compute damping ratios of all modes
for mode_name, mode_data in modes.items():
modes[mode_name]['damping_ratio'] = (
-mode_data['eigenvalue_real'] /
(mode_data['eigenvalue_real'] ** 2 + mode_data['eigenvalue_imag'] ** 2) ** 0.5
)
return modes
if __name__ == '__main__':
import aerosandbox as asb
import aerosandbox.numpy as np
from aerosandbox.tools import units as u
from pprint import pprint
# Numbers below are from:
# Caughey, David A., "Introduction to Aircraft Stability and Control, Course Notes for M&AE 5070", 2011
# https://courses.cit.cornell.edu/mae5070/Caughey_2011_04.pdf
[docs] airplane = asb.Airplane(
name='Boeing 737-800',
s_ref=1260 * u.foot ** 2,
c_ref=11 * u.foot,
b_ref=113 * u.foot,
)
aero = dict(
CL=1.83443,
CD=0.13037,
Cm=0,
CLa=5.542930,
Cma=-2.044696,
CYb=-1.103873,
Clb=-0.374933,
Cnb=0.239877,
CYp=0.800161,
Clp=-0.449404,
Cnp=-0.255028,
CLq=18.973344,
Cmq=-74.997742,
CYr=0.796001,
Clr=0.364638,
Cnr=-0.434410
)
mass_props_TOGW = asb.MassProperties(
mass=77146,
x_cg=65.2686 * u.foot,
y_cg=0.0,
z_cg=1.16559 * u.foot,
Ixx=706684,
Iyy=0.270824e7,
Izz=0.330763e7,
Ixy=0.0,
Iyz=0.0,
Ixz=26994.4,
)
op_point = asb.OperatingPoint(
atmosphere=asb.Atmosphere(
altitude=2438.399975619396,
method='differentiable',
),
velocity=85.64176936131635,
)
assert np.allclose(
aero["CL"],
mass_props_TOGW.mass * 9.81 / op_point.dynamic_pressure() / airplane.s_ref,
rtol=0.001
)
eigenvalues_from_AVL = {
'phugoid' : -0.0171382 + 0.145072j, # Real is wrong (2x)
'short_period' : -0.439841 + 0.842195j, # Pretty close
'roll_subsidence': -1.35132, # get_modes says -1.81
'dutch_roll' : -0.385418 + 1.52695j, # Imag is wrong (1.5x)
'spiral' : -0.0573017, # Too small, get_modes says -0.17
}
pprint(
get_modes(
airplane=airplane,
op_point=op_point,
mass_props=mass_props_TOGW,
aero=aero
)
)