From 6ac3eea8b4f3716ffd9a9138037880cf1ce84cc8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 18 Sep 2016 21:00:48 +0200 Subject: [PATCH 1/7] Example 1: stationary, horizontal, symmetric. Working! --- examples/example_001.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/example_001.py b/examples/example_001.py index c5595f3..af5f5f1 100644 --- a/examples/example_001.py +++ b/examples/example_001.py @@ -35,12 +35,12 @@ environment = Environment(atmosphere, gravity, wind) # Initial conditions. -TAS = 312.5 * 0.3048 # m/s -h0 = 8000 * 0.3048 # m +TAS = 45 # m/s +h0 = 3000 # m psi0 = 1.0 # rad x0, y0 = 0, 0 # m -turn_rate = 0.1 # rad/s -gamma0 = 0.1 # rad +turn_rate = 0.0 # rad/s +gamma0 = 0.0 # rad system = EulerFlatEarth(lat=0, lon=0, h=h0, psi=psi0, x_earth=x0, y_earth=y0) From 00ba88456abea97cfc45a60ebb043360520db952 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 18 Sep 2016 21:06:58 +0200 Subject: [PATCH 2/7] stationary descent, symmetric, wings level flight. Working! --- examples/example_002.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/examples/example_002.py b/examples/example_002.py index 2eb2ec5..ef82353 100644 --- a/examples/example_002.py +++ b/examples/example_002.py @@ -20,7 +20,7 @@ import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D -from pyfme.aircrafts import Cessna310 +from pyfme.aircrafts import Cessna172 from pyfme.environment.environment import Environment from pyfme.environment.atmosphere import ISA1976 from pyfme.environment.gravity import VerticalConstant @@ -29,24 +29,23 @@ from pyfme.simulator import BatchSimulation from pyfme.utils.trimmer import steady_state_flight_trimmer -aircraft = Cessna310() +aircraft = Cessna172() atmosphere = ISA1976() gravity = VerticalConstant() wind = NoWind() environment = Environment(atmosphere, gravity, wind) # Initial conditions. -TAS = 312.5 * 0.3048 # m/s -h0 = 8000 * 0.3048 # m -psi0 = 1 # rad +TAS = 45 # m/s +h0 = 3000 # m +psi0 = 1.0 # rad x0, y0 = 0, 0 # m turn_rate = 0.0 # rad/s -gamma0 = -0.05 # rad +gamma0 = -0.1 # rad system = EulerFlatEarth(lat=0, lon=0, h=h0, psi=psi0, x_earth=x0, y_earth=y0) not_trimmed_controls = {'delta_elevator': 0.05, - 'hor_tail_incidence': 0.00, 'delta_aileron': 0.01 * np.sign(turn_rate), 'delta_rudder': 0.01 * np.sign(turn_rate), 'delta_t': 0.5} From 8a4ecd1bdeaa239ca5eda6f2472c7268ce8aff8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 18 Sep 2016 21:10:20 +0200 Subject: [PATCH 3/7] stationary ascent, symmetric, wings level flight. Working! --- examples/example_003.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/examples/example_003.py b/examples/example_003.py index 045d30f..a54ca86 100644 --- a/examples/example_003.py +++ b/examples/example_003.py @@ -20,7 +20,7 @@ import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D -from pyfme.aircrafts import Cessna310 +from pyfme.aircrafts import Cessna172 from pyfme.environment.environment import Environment from pyfme.environment.atmosphere import ISA1976 from pyfme.environment.gravity import VerticalConstant @@ -29,24 +29,23 @@ from pyfme.simulator import BatchSimulation from pyfme.utils.trimmer import steady_state_flight_trimmer -aircraft = Cessna310() +aircraft = Cessna172() atmosphere = ISA1976() gravity = VerticalConstant() wind = NoWind() environment = Environment(atmosphere, gravity, wind) # Initial conditions. -TAS = 312.5 * 0.3048 # m/s -h0 = 8000 * 0.3048 # m -psi0 = 1 # rad +TAS = 45 # m/s +h0 = 3000 # m +psi0 = 1.0 # rad x0, y0 = 0, 0 # m turn_rate = 0.0 # rad/s -gamma0 = +0.05 # rad +gamma0 = 0.1 # rad system = EulerFlatEarth(lat=0, lon=0, h=h0, psi=psi0, x_earth=x0, y_earth=y0) not_trimmed_controls = {'delta_elevator': 0.05, - 'hor_tail_incidence': 0.00, 'delta_aileron': 0.01 * np.sign(turn_rate), 'delta_rudder': 0.01 * np.sign(turn_rate), 'delta_t': 0.5} From 9ec8952e2a65036c6bfc74288ea0ebe6a6af8b98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 18 Sep 2016 22:37:44 +0200 Subject: [PATCH 4/7] stationary, horizontal turn. Not fully convinced... --- examples/example_004.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/examples/example_004.py b/examples/example_004.py index d905149..7814942 100644 --- a/examples/example_004.py +++ b/examples/example_004.py @@ -19,7 +19,7 @@ import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D -from pyfme.aircrafts import Cessna310 +from pyfme.aircrafts import Cessna172 from pyfme.environment.environment import Environment from pyfme.environment.atmosphere import ISA1976 from pyfme.environment.gravity import VerticalConstant @@ -28,24 +28,23 @@ from pyfme.simulator import BatchSimulation from pyfme.utils.trimmer import steady_state_flight_trimmer -aircraft = Cessna310() +aircraft = Cessna172() atmosphere = ISA1976() gravity = VerticalConstant() wind = NoWind() environment = Environment(atmosphere, gravity, wind) # Initial conditions. -TAS = 312.5 * 0.3048 # m/s -h0 = 8000 * 0.3048 # m -psi0 = 1 # rad +TAS = 45 # m/s +h0 = 3000 # m +psi0 = 1.0 # rad x0, y0 = 0, 0 # m -turn_rate = 0.1 # rad/s -gamma0 = 0.00 # rad +turn_rate = 0.05 # rad/s +gamma0 = 0.0 # rad system = EulerFlatEarth(lat=0, lon=0, h=h0, psi=psi0, x_earth=x0, y_earth=y0) not_trimmed_controls = {'delta_elevator': 0.05, - 'hor_tail_incidence': 0.00, 'delta_aileron': 0.01 * np.sign(turn_rate), 'delta_rudder': 0.01 * np.sign(turn_rate), 'delta_t': 0.5} @@ -60,7 +59,7 @@ my_simulation = BatchSimulation(trimmed_ac, trimmed_sys, trimmed_env) -tfin = 150 # seconds +tfin = 30 # seconds N = tfin * 100 + 1 time = np.linspace(0, tfin, N) initial_controls = trimmed_ac.controls From 83279ddb5f8321e40a89e423ddac095ec1c176e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 18 Sep 2016 22:42:43 +0200 Subject: [PATCH 5/7] stationary, turn during ascent. Not convinced at all... --- examples/example_005.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/examples/example_005.py b/examples/example_005.py index 4385329..1a1361d 100644 --- a/examples/example_005.py +++ b/examples/example_005.py @@ -19,7 +19,7 @@ import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D -from pyfme.aircrafts import Cessna310 +from pyfme.aircrafts import Cessna172 from pyfme.environment.environment import Environment from pyfme.environment.atmosphere import ISA1976 from pyfme.environment.gravity import VerticalConstant @@ -28,24 +28,23 @@ from pyfme.simulator import BatchSimulation from pyfme.utils.trimmer import steady_state_flight_trimmer -aircraft = Cessna310() +aircraft = Cessna172() atmosphere = ISA1976() gravity = VerticalConstant() wind = NoWind() environment = Environment(atmosphere, gravity, wind) # Initial conditions. -TAS = 312.5 * 0.3048 # m/s -h0 = 8000 * 0.3048 # m -psi0 = 1 # rad +TAS = 45 # m/s +h0 = 3000 # m +psi0 = 1.0 # rad x0, y0 = 0, 0 # m -turn_rate = 0.1 # rad/s +turn_rate = 0.05 # rad/s gamma0 = 0.05 # rad system = EulerFlatEarth(lat=0, lon=0, h=h0, psi=psi0, x_earth=x0, y_earth=y0) not_trimmed_controls = {'delta_elevator': 0.05, - 'hor_tail_incidence': 0.00, 'delta_aileron': 0.01 * np.sign(turn_rate), 'delta_rudder': 0.01 * np.sign(turn_rate), 'delta_t': 0.5} @@ -60,7 +59,7 @@ my_simulation = BatchSimulation(trimmed_ac, trimmed_sys, trimmed_env) -tfin = 150 # seconds +tfin = 30 # seconds N = tfin * 100 + 1 time = np.linspace(0, tfin, N) initial_controls = trimmed_ac.controls From e499a5a7993da73ac476a1e105a41cd68b479ce9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 18 Sep 2016 22:45:29 +0200 Subject: [PATCH 6/7] longitudinal perturbation (delta doublet). ok! --- examples/example_006.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/examples/example_006.py b/examples/example_006.py index f704474..920cb59 100644 --- a/examples/example_006.py +++ b/examples/example_006.py @@ -17,7 +17,7 @@ import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D -from pyfme.aircrafts import Cessna310 +from pyfme.aircrafts import Cessna172 from pyfme.environment.environment import Environment from pyfme.environment.atmosphere import ISA1976 from pyfme.environment.gravity import VerticalConstant @@ -25,15 +25,16 @@ 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 -aircraft = Cessna310() +aircraft = Cessna172() atmosphere = ISA1976() gravity = VerticalConstant() wind = NoWind() environment = Environment(atmosphere, gravity, wind) # Initial conditions. -TAS = 312.5 * 0.3048 # m/s +TAS = 45 # m/s h0 = 8000 * 0.3048 # m psi0 = 1 # rad x0, y0 = 0, 0 # m @@ -43,7 +44,6 @@ system = EulerFlatEarth(lat=0, lon=0, h=h0, psi=psi0, x_earth=x0, y_earth=y0) not_trimmed_controls = {'delta_elevator': 0.05, - 'hor_tail_incidence': 0.00, 'delta_aileron': 0.01 * np.sign(turn_rate), 'delta_rudder': 0.01 * np.sign(turn_rate), 'delta_t': 0.5} @@ -68,11 +68,12 @@ controls[control_name] = np.ones_like(time) * control_value # Elevator doublet -controls['delta_elevator'][np.where(time<2)] = \ - initial_controls['delta_elevator'] * 1.30 - -controls['delta_elevator'][np.where(time<1)] = \ - initial_controls['delta_elevator'] * 0.70 +amplitude = initial_controls['delta_elevator'] * 1.5 +controls['delta_elevator'] = doublet(t_init=1, + T=1, + A=amplitude, + time=time, + offset=initial_controls['delta_elevator']) my_simulation.set_controls(time, controls) From 72df09a39d2fdbc792c5a7a9fdc2d56f2751cabf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 18 Sep 2016 22:56:07 +0200 Subject: [PATCH 7/7] lat-dir perturbation. not convinced... --- examples/example_007.py | 125 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 examples/example_007.py diff --git a/examples/example_007.py b/examples/example_007.py new file mode 100644 index 0000000..a074269 --- /dev/null +++ b/examples/example_007.py @@ -0,0 +1,125 @@ +# -*- coding: utf-8 -*- +""" +Python Flight Mechanics Engine (PyFME). +Copyright (c) AeroPython Development Team. +Distributed under the terms of the MIT License. + +Example +------- + +Cessna 310, ISA1976 integrated with Flat Earth (euler angles). + +Evolution of the aircraft after a lat-dir perturbation. +Trimmed in stationary, horizontal, symmetric, wings level flight. +""" + +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + +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 + +aircraft = Cessna172() +atmosphere = ISA1976() +gravity = VerticalConstant() +wind = NoWind() +environment = Environment(atmosphere, gravity, wind) + +# Initial conditions. +TAS = 45 # m/s +h0 = 8000 * 0.3048 # m +psi0 = 1 # rad +x0, y0 = 0, 0 # m +turn_rate = 0.0 # rad/s +gamma0 = 0.00 # rad + +system = EulerFlatEarth(lat=0, lon=0, h=h0, psi=psi0, x_earth=x0, y_earth=y0) + +not_trimmed_controls = {'delta_elevator': 0.05, + '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'] + +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=2) + +print(results) + +my_simulation = BatchSimulation(trimmed_ac, trimmed_sys, trimmed_env) + +tfin = 15 # seconds +N = tfin * 100 + 1 +time = np.linspace(0, tfin, N) +initial_controls = trimmed_ac.controls + +controls = {} +for control_name, control_value in initial_controls.items(): + controls[control_name] = np.ones_like(time) * control_value + +# Rudder doublet +amplitude = 0.10 +controls['delta_elevator'] = doublet(t_init=1, + T=1, + A=amplitude, + time=time, + offset=initial_controls['delta_rudder']) + +# Rudder aileron +amplitude = 0.15 +controls['delta_elevator'] = doublet(t_init=1.5, + T=1, + A=amplitude, + time=time, + offset=initial_controls['delta_aileron']) + +my_simulation.set_controls(time, controls) + +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'] + +my_simulation.set_par_dict(par_list) +my_simulation.run_simulation() + +# print(my_simulation.par_dict) + +plt.style.use('ggplot') + +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()