Los datos se han obtenido del libro Dynamics of Flight Stability and Control, B. Etkin & L.D. Reid pag 165 y corresponden a un Boeing 747-100. La condición de referencia también se ha tomado del Etkin.
from sympy import *
init_printing()
import numpy as np
#se crean los símbolos para todas las variables:
u, alfa, tetha, de, q = symbols('\\Delta\hat{u}, \\Delta\\alpha, \\Delta\\theta, \\Delta\\delta_e, \\Delta\hat{q}')
mu, Iy, D, landa = symbols('mu, \hat{I}_y, D, \hat{\\lambda}')
C_Xu, C_Xalfa, C_Xde = symbols('C_X_\hat{u}, C_X_alpha, C_X_\\delta_e')
C_Zs, C_Zu, C_Zalfap, C_Zalfa, C_Zq, C_Zde = symbols("C_Z_s, C_Z_\\hat{u}, C_Z_\hat{\dot{\\alpha}}, C_Z_alpha, C_Z_\hat{q}, C_Z_\\delta_e")
C_mu, C_malfap, C_malfa, C_mq, C_mde, C_mdep = symbols("C_m_\hat{u}, C_m_\hat{\dot{\\alpha}}, C_m_alpha, C_m_\hat{q}, C_m_\\delta_e, C_m_\hat{\dot{\\delta_e}}")
C_mu, C_malfap, C_malfa, C_mq, C_mde, C_mdep, C_Zs, C_Zu, C_Zalfap, C_Zalfa, C_Zq, C_Zde, C_Xu, C_Xalfa, C_Xde, u, alfa, tetha, de, q, mu, Iy, D, landa
W = 2.83176e6 #N
Sup_alar = 511. #m^2
cuerda = 8.324 #m
Inercia_y = 0.449e8 #kg*m^2
rho = 0.3045 #Kg/m^3
us = 235.9 #m/s
masa_adim = W/9.81/(rho*Sup_alar*cuerda/2)
Inercia_y_adim = 8.*Inercia_y / (rho*Sup_alar*cuerda**3)
coef_fuerza_zs = -W / (1./2.*rho*us**2*Sup_alar)
num_vals = {C_mu:0.1043, C_malfap:-6.314, C_malfa:-1.023, C_mq:-23.92, C_mde:-1.444, C_mdep:0,\
C_Zs:coef_fuerza_zs, C_Zu: -0.106, C_Zalfap:5.9, C_Zalfa:-4.92, C_Zq:-5.92, C_Zde:-0.3648,\
C_Xu: -0.1080, C_Xalfa: 0.2193, C_Xde:0,\
mu: masa_adim, Iy:Inercia_y_adim }
#se ha tenido en cuenta theta_s=0
#se ha cambiado ¿una errata en el libro? en los términos F_3i
F = Matrix([[(C_Xu) / (2*mu), C_Xalfa/(2*mu), 0, C_Zs/(2*mu)],
[(C_Zu + 2*C_Zs) / (2*mu - C_Zalfap), C_Zalfa / (2*mu - C_Zalfap), (2*mu + C_Zq) / (2*mu - C_Zalfap), 0],
[C_mu/Iy + C_malfap/Iy * (C_Zu + 2*C_Zs)/(2*mu - C_Zalfap), C_malfa/Iy + C_malfap/Iy * C_Zalfa/(2*mu - C_Zalfap), C_mq/Iy + C_malfap/Iy * (2*mu + C_Zq)/(2*mu - C_Zalfap), 0],
[0,0,1,0]])
F
F_num = F.subs(num_vals)
F_num.evalf()
(F - landa * eye(4))
(F - landa * eye(4)).det()
polinomio = (F-landa * eye(4)).det().subs(num_vals).evalf()
polinomio.evalf(5)
roots(polinomio)
F_num = np.array(np.array(F_num), np.float)
autovalores, autovectores = np.linalg.eig(F_num)
autovalores #son los mismos que se han sacado del pol característico para comprobar
# Además salen los mismos (lógico) que cuando había puesto la matriz característica.
array([ -6.55727818e-03+0.01564731j, -6.55727818e-03-0.01564731j, -5.80316373e-05+0.00118576j, -5.80316373e-05-0.00118576j])
Se observa que se han obtenidos dos pares de autovalores complejos conjugados (sistema oscilante). Todos con parte real negativa, lo que indica que el sistema va a ser estable. La parte real de una pareja de autovalores es apreciablemente menor que la de la otra. Esto nos indica que ese modo va a estar débilmente amortiguado y corresponderá al fugoide. La otra pareja corresponderá al corto periodo.
Se reordenan los autovalores y sus autovectores de forma que $\lambda_1$ y $\lambda_2$ correspondan al fugoide y $\lambda_3$ y $\lambda_4$ al corto periodo:
landa = np.zeros([4],dtype=complex)
modos = np.zeros([4,4],dtype=complex)
for i in range(4):
j = argmax(autovalores)
landa[i] = autovalores[j]
modos[i] = autovectores[:,j]
autovalores[j] = -100 #esto de quitar los que ya he cogido debería hacerlo más elegante
Se trabajará ahora sólo con la primera pareja de autovalores y sus atovectores asociados.
Los autovalores se pueden caracterizar por su parte real $\hat{n}$ e imaginaria $\omega$:
n = np.zeros(4)
w = np.zeros(4)
for i in range(4):
n[i] = np.real(landa[i])
w[i] = np.imag(landa[i])
print('n1 = %g \n'
'w1 = %g' %(n[0], w[0]))
n1 = -5.80316e-05 w1 = 0.00118576
En lugar de caracterizar la raiz por su parte real $n$ e imaginaria $\omega$, puede resultar más interesante dar el tiempo para que la perturbación se reduzca a la mitad: $\hat{t}_{1/2}$ y el periodo del modo oscilatorio $T$:
t05 = np.zeros(4)
T = np.zeros(4)
for i in range(4):
t05[i] = -0.693/n[i]
T[i] = abs(2*np.pi/w[i])
print('t0.5 = %g \n'
'T = %g' %(t05[0], T[0]))
t0.5 = 11941.8 T = 5298.88
Estos tiempos son adimensionales, si se quiere recuperar las unidades, no hay más que multiplicar por el adimensionalizador de tiempos:
t05_d = np.zeros(4)
T_d = np.zeros(4)
tc = cuerda/(2*us)
for i in range(4):
t05_d[i] = t05[i] * tc
T_d[i] = T[i] * tc
print('t0.5 = %g s\n'
'T = %g s' %(t05_d[0], T_d[0]))
t0.5 = 210.689 s T = 93.4886 s
Una última forma de expresar los autovalores sería mediante su amortiguamiento $\xi$ y frecuencia natural $\hat{\omega}_n$:
xi = np.zeros(4)
wn = np.zeros(4)
for i in range(4):
xi[i] = -n[i] /(np.sqrt(n[i]**2 + w[i]**2))
wn[i] = np.sqrt(n[i]**2 + w[i]**2)
print('xi = %g \n'
'wn = %g' %(xi[0], wn[0]))
xi = 0.0488821 wn = 0.00118718
la frecuencia natural dimensional, $\omega_n$, es:
wn_d = wn * 1./tc
print('wn = %g s^-1' %(wn_d[0]))
wn = 0.0672885 s^-1
Se pasa a analizar ahora la forma del modo asociado a esta pareja de autovalores complejos:
modos[0]
array([ -2.16226620e-02+0.5243972j , 3.84126916e-03+0.03032094j, -4.93645945e-05+0.00100866j, 8.50649694e-01+0.j ])
modos[1]
array([ -2.16226620e-02-0.5243972j , 3.84126916e-03-0.03032094j, -4.93645945e-05-0.00100866j, 8.50649694e-01+0.j ])
Se puede observar que como los autovalores son conjugados, sus modos también lo son.
Estos autovectores están expresados en la forma:
u, alfa, q, tetha
Si se referencian a $\Delta\theta$ todos los valores y mostramos sólo los módulos:
amplitud_1 = abs(modos[0]) / abs(modos[0][3])
amplitud_1
array([ 0.61699052, 0.03592935, 0.00118718, 1. ])
se observa que:
Referenciando también la fase a $\Delta\theta$ se obtiene:
# np.angle da el argumento del vector en sentido antihorario
desfase_1 = np.angle(modos[0], deg=True) - np.angle(modos[0][3], deg=True)
desfase_1
array([ 92.36116035, 82.77983289, 92.80185514, 0. ])
se observa que:
import matplotlib.pyplot as plt
plt.figure()
#me gustaría ponerle nombre a cada fasor
for i in range(4):
dx = amplitud_1[i]*np.cos(np.deg2rad(desfase_1[i]))
dy = amplitud_1[i]*np.sin(np.deg2rad(desfase_1[i]))
plt.arrow(0,0,dx,dy, length_includes_head = True, width = 0.0015)
plt.xlim(-0.1,1)
plt.ylim(-0.11,1)
plt.xlabel(r'$Re$', fontsize=18)
plt.ylabel(r'$Im$', fontsize=18)
plt.grid(color='g', alpha=0.5, linestyle='dashed', linewidth=0.5)
Se procederá ahroa de manera análoga al fugoide, pero trabajando con la segunda pareja de autovalores y autovectores:
landa[2]
(-0.0065572781755465171+0.01564730889682945j)
landa[3]
(-0.0065572781755465171-0.01564730889682945j)
Los autovalores se pueden caracterizar por su parte real $\hat{n}$ e imaginaria $\omega$:
print('n3 = %g \n'
'w3 = %g' %(n[2], w[2]))
n3 = -0.00655728 w3 = 0.0156473
tiempo para que la perturbación se reduzca a la mitad: $\hat{t}_{1/2}$ y el periodo del modo oscilatorio $T$:
print('t0.5 = %g \n'
'T = %g' %(t05[2], T[2]))
t0.5 = 105.684 T = 401.551
print('t0.5 = %g s\n'
'T = %g s' %(t05_d[2], T_d[2]))
t0.5 = 1.86459 s T = 7.08458 s
amortiguamiento $\xi$ y frecuencia natural $\hat{\omega}_n$:
print('xi = %g \n'
'wn = %g' %(xi[2], wn[2]))
xi = 0.386501 wn = 0.0169657
print('wn = %g s^-1' %(wn_d[0]))
wn = 0.0672885 s^-1
Se pasa a analizar ahora la forma del modo asociado a esta pareja de autovalores complejos:
modos[2]
array([ 1.54797745e-02+0.01216979j, 7.33677679e-01+0.j , -7.10104441e-04+0.01149977j, 6.41326188e-01-0.22337705j])
modos[3]
array([ 1.54797745e-02-0.01216979j, 7.33677679e-01+0.j , -7.10104441e-04-0.01149977j, 6.41326188e-01+0.22337705j])
Si se referencian a $\Delta\theta$ todos los valores y mostramos sólo los módulos:
amplitud_2 = abs(modos[2]) / abs(modos[2][3])
amplitud_2
array([ 0.02899479, 1.0803445 , 0.01696574, 1. ])
se observa que:
Referenciando también la fase a $\Delta\theta$ se obtiene:
desfase_2 = np.angle(modos[2], deg=True) - np.angle(modos[2][3], deg=True)
desfase_2
array([ 57.376961 , 19.20347759, 112.73697243, 0. ])
se observa que:
for i in range(4):
dx = amplitud_2[i]*np.cos(np.deg2rad(desfase_2[i]))
dy = amplitud_2[i]*np.sin(np.deg2rad(desfase_2[i]))
plt.arrow(0,0,dx,dy, length_includes_head = True, width = 0.0015)
plt.xlim(-0.1,1.1)
plt.ylim(-0.11,1)
plt.xlabel(r'$Re$', fontsize=18)
plt.ylabel(r'$Im$', fontsize=18)
plt.grid(color='g', alpha=0.5, linestyle='dashed', linewidth=0.5)
class Modos:
''' n, w son las partes re e im del autovalor
phi_u, phi_alfa son los desfases en grados con respecto a theta del modo
u, alfa son los módulos de las dos primeras componentes del modo referidas a theta'''
def __init__(self, n, w, phi_u, phi_alfa, u, alfa):
self.n = n
self.w = w
self.phi_u = phi_u / 180. * 3.14159
self.phi_alfa = phi_alfa / 180. * 3.14159
self.u = abs(u)
self.alfa = abs(alfa)
def du(self,t):
n = self.n
w = self.w
phi = self.phi_u
u = self.u
return 2 * np.real(u * np.exp(n*t) * np.exp(1j*(phi+w*t)))
def dalfa(self,t):
n = self.n
w = self.w
phi = self.phi_alfa
alfa = self.alfa
return 2 * np.real(alfa * np.exp(n*t) * np.exp(1j*(phi+w*t)))
def dtheta(self,t):
n = self.n
w = self.w
return 2 * np.real(np.exp(n*t) * np.exp(1j*(w*t)))
fugoide = Modos(n[0], w[0], desfase_1[0], desfase_1[1], amplitud_1[0], amplitud_1[1])
t_adim = np.linspace(0,40000,1000)
plt.figure()
plt.plot(t_adim * tc, fugoide.du(t_adim))
plt.xlabel('tiempo(s)', fontsize=18)
plt.ylabel(r'$\Delta u/us$', fontsize=18)
plt.grid()
La variación de la velocidad es muy alta porque se ha tomado $\theta=2$ RADIAN
plt.figure()
plt.plot(t_adim * tc, fugoide.dalfa(t_adim))
plt.xlabel('tiempo(s)', fontsize=18)
plt.ylabel(r'$\alpha$', fontsize=18)
plt.grid()
Casi constante
plt.figure()
plt.plot(t_adim * tc, fugoide.dtheta(t_adim))
plt.xlabel('tiempo(s)', fontsize=18)
plt.ylabel(r'$\theta$', fontsize=18)
plt.grid()
cortoper = Modos(n[2], w[2], desfase_2[0], desfase_2[1], amplitud_2[0], amplitud_2[1])
t_adim = np.linspace(0,1000,200)
plt.figure()
plt.plot(t_adim * tc, cortoper.du(t_adim))
plt.xlabel('tiempo(s)', fontsize=18)
plt.ylabel(r'$\Delta u/us$', fontsize=18)
plt.grid()
plt.figure()
plt.plot(t_adim * tc, cortoper.dalfa(t_adim))
plt.xlabel('tiempo(s)', fontsize=18)
plt.ylabel(r'$\Delta \alpha$', fontsize=18)
plt.grid()
plt.figure()
plt.plot(t_adim * tc, cortoper.dtheta(t_adim))
plt.xlabel('tiempo(s)', fontsize=18)
plt.ylabel(r'$\Delta \theta$', fontsize=18)
plt.grid()
t_adim = np.linspace(0,40000,1000)
plt.figure()
plt.plot(t_adim * tc, cortoper.du(t_adim)+fugoide.du(t_adim))
plt.xlabel('tiempo(s)', fontsize=18)
plt.ylabel(r'$\Delta u/us$', fontsize=18)
plt.grid()
plt.figure()
plt.plot(t_adim * tc, cortoper.dalfa(t_adim)+fugoide.dalfa(t_adim))
plt.xlabel('tiempo(s)', fontsize=18)
plt.ylabel(r'$\Delta \alpha$', fontsize=18)
plt.xlim(0,100)
plt.grid()
plt.figure()
plt.plot(t_adim * tc, cortoper.dtheta(t_adim)+fugoide.dtheta(t_adim))
plt.xlabel('tiempo(s)', fontsize=18)
plt.ylabel(r'$\Delta \theta$', fontsize=18)
plt.xlim(0,200)
plt.grid()