from sympy import *
#from sympy.abc import *
init_printing(use_latex=True)
from sympy import symbols
from sympy.physics.mechanics import ReferenceFrame, Point, dynamicsymbols
Se definen los distintos sistemas de referencia:
#ejes tierra
Earth = ReferenceFrame("Earth", latexs=['\mathbf{i_{e}}', '\mathbf{j_{e}}', '\mathbf{k_{e}}'])
#ejes horizonte local
Hor = ReferenceFrame("Hor", latexs=['\mathbf{i_{h}}', '\mathbf{j_{h}}', '\mathbf{k_{h}}'])
#ejes cuerpo
Body = ReferenceFrame("Body", latexs=['\mathbf{i_{b}}', '\mathbf{j_{b}}', '\mathbf{k_{b}}'])
#ejes viento
Wind = ReferenceFrame("Wind", latexs=['\mathbf{i_{w}}', '\mathbf{j_{w}}', '\mathbf{k_{w}}'])
Los sistemas de referencia se orientan mediante los Ángulos de Euler según el convenio de Tait-Bryan $zyx$ ($321$)
$\psi =$ ángulo de guiñada del avión.
$\theta =$ ángulo de cabeceo del avión.
$\phi =$ ángulo de balance del avión.
psi, theta, phi = dynamicsymbols('psi, theta, phi')
Body.orient(Hor, 'Body', (psi, theta, phi), 'ZYX')
L_bh = Body.dcm(Hor)
L_bh
chi, gamma, mu = dynamicsymbols('chi, \gamma, mu')
Wind.orient(Hor, 'Body', (chi, gamma, mu), 'ZYX')
L_wh = Wind.dcm(Hor)
L_wh
beta, alpha = dynamicsymbols('beta, alpha')
Body.orient(Wind, 'Body', (-beta, alpha, 0), 'ZYX')
L_bw = Body.dcm(Wind)
L_bw
Como los ejes horizonte local son paralelos a los ejes tierra con origen en el punto subavión en todo momento por definición:
Hor.orient(Earth, 'Body', (0, 0, 0), 'ZYX')
Hor.dcm(Earth)
Se utilizarán ahora los teoremas de cantidad de movimiento y momento cinético para obtener las ecuaciones del movimiento del avión:
$$\boldsymbol{\mathbf{F}}=m\cdot\frac{d\mathbf{\boldsymbol{V}}}{dt}$$(se ha despreciado la variación de la masa)
$$\boldsymbol{G}=\frac{d\boldsymbol{h}}{dt}=\boldsymbol{I}\boldsymbol{\omega}$$Estas ecuaciones se proyectan sobre el sistema de ejes cuerpo con la doble ventaja de que, en ellos, no varía el tensor de inercia y además debido a la simetría que aparece en todos los aviones convencionales habrá términos de dicho tensor que sean cero.
from sympy.physics.mechanics import inertia
I_x, I_y, I_z, J_xz = symbols('I_x, I_y, I_z, J_xz')
I = inertia(Body, I_x, I_y, I_z, 0, 0, -J_xz)
I
F_x, F_y, F_z = dynamicsymbols('F_x, F_y, F_z')
L, M, N = dynamicsymbols('L, M, N')
u, v, w = dynamicsymbols('u, v, w')
p, q, r = dynamicsymbols('p, q, r')
F = F_x * Body.x + F_y * Body.y + F_z * Body.z
G = L * Body.x + M * Body.y + N * Body.z
V = u * Body.x + v * Body.y + w * Body.z
W = p * Body.x + q * Body.y + r * Body.z
F = (V.dt(Body)) + (W ^ V)
F
G = ((I & W).dt(Body)) + (W ^ (I & W))
G
Body.ang_vel_in(Hor)
ang_vel = Body.ang_vel_in(Hor)-(p * Body.x + q * Body.y + r * Body.z)
solve([ang_vel.args[0][0][i] for i in range(3)],[diff(phi),diff(theta),diff(psi)])