Cessna 172, ISA1976 integrated with Flat Earth (Euler angles).
Evolution of the aircraft after a pitch perturbation (delta doublet applied to the elevator) at t=2, a roll perturbation (delta doublet applied to the ailerons) at t=5 and a yaw perturbation (delta doublet applied to the rudder) at t=10.
Initially trimmed to a stationary, horizontal, symmetric, wings level flight.
In [1]:
# -*- coding: utf-8 -*-
Import python libreries needed.
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
In [3]:
import plotly.offline as offline
import plotly.graph_objs as go
import plotly.tools as tls
# notebook mode - inline
offline.init_notebook_mode()
Import PyFME classes
In [4]:
from pyfme.aircrafts import Cessna172
from pyfme.environment.environment import Environment
from pyfme.environment.atmosphere import ISA1976
from pyfme.environment.gravity import VerticalConstant
from pyfme.environment.wind import NoWind
from pyfme.models.systems import EulerFlatEarth
from pyfme.simulator import BatchSimulation
from pyfme.utils.trimmer import steady_state_flight_trimmer
from pyfme.utils.input_generator import doublet
Initialize variables
In [5]:
aircraft = Cessna172()
atmosphere = ISA1976()
gravity = VerticalConstant()
wind = NoWind()
environment = Environment(atmosphere, gravity, wind)
Initial conditions
In [6]:
TAS = 45 # m/s
h0 = 2000 # m
psi0 = 1 # rad
x0, y0 = 0, 0 # m
turn_rate = 0.0 # rad/s
gamma0 = 0.00 # rad
Define system
In [7]:
system = EulerFlatEarth(lat=0, lon=0, h=h0, psi=psi0, x_earth=x0, y_earth=y0)
not_trimmed_controls = {'delta_elevator': 0.0,
'delta_aileron': 0.01 * np.sign(turn_rate),
'delta_rudder': 0.01 * np.sign(turn_rate),
'delta_t': 0.5}
controls2trim = ['delta_elevator', 'delta_aileron', 'delta_rudder', 'delta_t']
In [8]:
trimmed_ac, trimmed_sys, trimmed_env, results = steady_state_flight_trimmer(
aircraft, system, environment, TAS=TAS, controls_0=not_trimmed_controls,
controls2trim=controls2trim, gamma=gamma0, turn_rate=turn_rate, verbose=1)
Steady state flight trimmer results
In [9]:
print('delta_elevator = ',"%8.4f" % np.rad2deg(results['delta_elevator']), 'deg')
print('delta_aileron = ', "%8.4f" % np.rad2deg(results['delta_aileron']), 'deg')
print('delta_rudder = ', "%8.4f" % np.rad2deg(results['delta_rudder']), 'deg')
print('delta_t = ', "%8.4f" % results['delta_t'], '%')
print()
print('alpha = ', "%8.4f" % np.rad2deg(results['alpha']), 'deg')
print('beta = ', "%8.4f" % np.rad2deg(results['beta']), 'deg')
print()
print('u = ', "%8.4f" % results['u'], 'm/s')
print('v = ', "%8.4f" % results['v'], 'm/s')
print('w = ', "%8.4f" % results['w'], 'm/s')
print()
print('psi = ', "%8.4f" % np.rad2deg(psi0), 'deg')
print('theta = ', "%8.4f" % np.rad2deg(results['theta']), 'deg')
print('phi = ', "%8.4f" % np.rad2deg(results['phi']), 'deg')
print()
print('p =', "%8.4f" % results['p'], 'rad/s')
print('q =', "%8.4f" % results['q'], 'rad/s')
print('r =', "%8.4f" % results['r'], 'rad/s')
Initialise simulation
In [10]:
my_simulation = BatchSimulation(trimmed_ac, trimmed_sys, trimmed_env)
In [11]:
tfin = 15 # seconds
N = tfin * 100 + 1
time = np.linspace(0, tfin, N)
initial_controls = trimmed_ac.controls
In [12]:
controls = {}
for control_name, control_value in initial_controls.items():
controls[control_name] = np.ones_like(time) * control_value
In [13]:
amplitude_elev = np.deg2rad(26)
controls['delta_elevator'] = initial_controls['delta_elevator'] + \
doublet(t_init=2,
T=2,
A=amplitude_elev,
time=time,
offset=np.deg2rad(0.0))
amplitude_aile = np.deg2rad(15)
controls['delta_aileron'] = doublet(t_init=5,
T=2,
A=amplitude_aile,
time=time,
offset=np.deg2rad(0.0))
amplitude_rud = np.deg2rad(32)
controls['delta_rudder'] = doublet(t_init=10,
T=2,
A=amplitude_rud,
time=time,
offset=np.deg2rad(0.0))
In [14]:
my_simulation.set_controls(time, controls)
In [15]:
par_list = ['x_earth', 'y_earth', 'height',
'psi', 'theta', 'phi',
'u', 'v', 'w',
'v_north', 'v_east', 'v_down',
'p', 'q', 'r',
'alpha', 'beta', 'TAS',
'F_xb', 'F_yb', 'F_zb',
'M_xb', 'M_yb', 'M_zb']
Run Simulation
In [16]:
my_simulation.set_par_dict(par_list)
my_simulation.run_simulation()
Plot results
In [17]:
plt.style.use('ggplot')
In [18]:
# for ii in range(len(par_list) // 3):
# three_params = par_list[3 * ii:3 * ii + 3]
# fig, ax = plt.subplots(3, 1, sharex=True)
# for jj, par in enumerate(three_params):
# ax[jj].plot(time, my_simulation.par_dict[par])
# ax[jj].set_ylabel(par)
# ax[jj].set_xlabel('time (s)')
# fig = plt.figure()
# ax = Axes3D(fig)
# ax.plot(my_simulation.par_dict['x_earth'],
# my_simulation.par_dict['y_earth'],
# my_simulation.par_dict['height'])
# ax.plot(my_simulation.par_dict['x_earth'],
# my_simulation.par_dict['y_earth'],
# my_simulation.par_dict['height'] * 0)
# ax.set_xlabel('x_earth')
# ax.set_ylabel('y_earth')
# ax.set_zlabel('z_earth')
# plt.show()
In [19]:
camera = dict(
up=dict(x=2, y=1, z=1),
center=dict(x=0, y=0, z=0),
eye=dict(x=2, y=1, z=0.7)
)
layout = {
'autosize': False,
'width': 900,
'height': 500,
'xaxis': {
'range': [-2, 4.5],
'zeroline': False,
},
'yaxis': {
'range': [-2, 4.5]
},
'width': 800,
'height': 800,
'scene':{
'camera': camera
}
}
trace_with_height = go.Scatter3d(
x=my_simulation.par_dict['x_earth'],
y=my_simulation.par_dict['y_earth'],
z=my_simulation.par_dict['height'],
mode='lines')
trace_without_height = go.Scatter3d(
x=my_simulation.par_dict['x_earth'],
y=my_simulation.par_dict['y_earth'],
z=my_simulation.par_dict['height'] * 0,
mode='lines')
data = [trace_with_height, trace_without_height]
fig = {
'data': data,
'layout': layout,
}
offline.iplot(fig)
In [20]:
for ii in range(len(par_list) // 3):
three_params = par_list[3 * ii:3 * ii + 3]
trace1 = go.Scatter(
x=time,
y=my_simulation.par_dict[three_params[0]]
)
trace2 = go.Scatter(
x=time,
y=my_simulation.par_dict[three_params[1]]
)
trace3 = go.Scatter(
x=time,
y=my_simulation.par_dict[three_params[2]]
)
fig = tls.make_subplots(rows=3, cols=1, subplot_titles=('Plot 1', 'Plot 2', 'Plot 3'))
# add subplots
fig.append_trace(trace1, 1, 1)
fig.append_trace(trace2, 2, 1)
fig.append_trace(trace3, 3, 1)
# Edit the layout
# All of the axes properties here: https://plot.ly/python/reference/#XAxis
fig['layout']['xaxis3'].update(title='time (s)', showgrid=True)
# All of the axes properties here: https://plot.ly/python/reference/#YAxis
fig['layout']['yaxis1'].update(title=three_params[0], showgrid=True)
fig['layout']['yaxis2'].update(title=three_params[1], showgrid=True)
fig['layout']['yaxis3'].update(title=three_params[2], showgrid=True)
fig['layout'].update(title='Output subplots')
offline.iplot(fig)