Modelo de 6 grados de libertad de un MAV de alas fijas

Para realizar el modelo dinámico del MAV se deben considerar 3 posiciones y 3 velocidades asociadas con el movimiento de traslación del vehículo, así como 3 posiciones y 3 velocidades asociadas con el movimiento de rotación. Esto constituye 12 variables de estado a considerar para el modelo (Beard & McLain, (2012) [1]), las cuales se listan en la Tabla 1.

Tabla 1: Variables de estado.

El control de una aeronave se debe al movimiento de las superficies de control en conjunto con el empuje. Las superficies de control son los alerones, el elevador y el timón. Esto constituye las variables de control del sistema, que se muestran en la Tabla 2.

Tabla 2: Variables de control.

Ahora bien, para realizar el modelo dinámico se requiere conocer la derivada temporal de cada variable de estado. Siendo que las posiciones PN, PE, y PD se encuentran definidas en el marco de referencia inercial, y las velocidades u, v, y w, en el marco de referencia del cuerpo, la derivada temporal de las posiciones se puede obtener al expresar las velocidades en el marco inercial.

$\left(\begin{array}{c}\dot{P_N}\\ \dot{P_E}\\ \dot{P_D}\end{array}\right)=R_b^v\left(\begin{array}{c}u\\ v\\ w\end{array}\right)={(R_v^b)}^{\text{T}}\left(\begin{array}{c}u\\ v\\ w\end{array}\right)={(R_{v2}^b{R}_{v1}^{v2}{R}_v^{v1})}^{\text{T}}\left(\begin{array}{c}u\\ v\\ w\end{array}\right)$  (1)

En donde $R_b^v$ es la matriz de rotación necesaria para pasar del marco del cuerpo (Fb), al marco del vehículo (Fv) y

$R_v^{v1}=\left(\begin{array}{ccc}\cos \mathrm{\psi }& \sin \mathrm{\psi }& 0\\ -\sin \mathrm{\psi }& \cos \mathrm{\psi }& 0\\ 0& 0& 1\end{array}\right)\text{,}{R}_{v1}^{v2}=\left(\begin{array}{ccc}\cos \mathrm{\theta }& 0& -\sin \mathrm{\theta }\\ 0& 1& 0\\ \sin \mathrm{\theta }& 0& \cos \mathrm{\theta }\end{array}\right)\text{,}{R}_{v2}^b=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\varphi }& \sin \mathrm{\varphi }\\ 0& -\sin \mathrm{\varphi }& \cos \mathrm{\varphi }\end{array}\right)$  (2)

son las matrices de rotación necesarias para pasar del marco del vehículo (Fv) al marco del vehículo 1 (Fv1), del marco del vehículo 1 (Fv1) al marco del vehículo 2 (Fv2), y del marco del vehículo 2 (Fv2) al marco del cuerpo (Fb) respectivamente.

Por su parte, la derivada temporal de las posiciones angulares ϕ, θ, y ψ, se puede obtener al expresar las velocidades angulares p, q y r, en los marcos correspondientes de cada posición angular, es decir:

$\left(\begin{array}{c}p\\ q\\ r\end{array}\right)=\left(\begin{array}{c}\dot{\mathrm{\varphi }}\\ 0\\ 0\end{array}\right)+R_{v2}^b\left(\begin{array}{c}0\\ \dot{\mathrm{\theta }}\\ 0\end{array}\right)+R_{v2}^b{R}_{v1}^{v2}\left(\begin{array}{c}0\\ 0\\ \dot{\mathrm{\psi }}\end{array}\right)=\left(\begin{array}{ccc}1& 0& -\sin \mathrm{\theta }\\ 0& \cos \mathrm{\varphi }& \sin \mathrm{\varphi }\cos \mathrm{\theta }\\ 0& -\sin \mathrm{\varphi }& \cos \mathrm{\varphi }\cos \mathrm{\theta }\end{array}\right)\left(\begin{array}{c}\dot{\mathrm{\varphi }}\\ \dot{\mathrm{\theta }}\\ \dot{\mathrm{\psi }}\end{array}\right)$  (3)

De la ecuación anterior se obtiene al despejar la derivada temporal de las posiciones angulares:

$\left(\begin{array}{c}\dot{\mathrm{\varphi }}\\ \dot{\mathrm{\theta }}\\ \dot{\mathrm{\psi }}\end{array}\right)=\left(\begin{array}{ccc}1& \sin \mathrm{\varphi }\tan \mathrm{\theta }& \cos \mathrm{\varphi }\tan \mathrm{\theta }\\ 0& \cos \mathrm{\varphi }& -\sin \mathrm{\varphi }\\ 0& \sin \mathrm{\varphi }\mathit{sec}\mathrm{\theta }& \cos \mathrm{\varphi }\mathit{sec}\mathrm{\theta }\end{array}\right)\left(\begin{array}{c}p\\ q\\ r\end{array}\right)$  (4)

Ahora bien, por segunda ley de Newton, la fuerza que experimenta el MAV es:

$F^b=\left(\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right)=m\frac{d{V}_g^b}{{\mathit{dt}}_i}=m\left(\frac{d{V}_g^b}{{\mathit{dt}}_b}+{\mathrm{\omega }}_{b/i}^b\times V_g^b\right)=m\left(\left(\begin{array}{c}\dot{u}\\ \dot{v}\\ \dot{w}\end{array}\right)+\left(\begin{array}{c}p\\ q\\ r\end{array}\right)\times \left(\begin{array}{c}u\\ v\\ w\end{array}\right)\right)$  (5)

En donde Fb es el vector resultante de fuerzas en el marco del cuerpo, m es la masa del vehículo, Vgb es la velocidad del MAV en el marco del cuerpo, $d/{\mathit{dt}}_i$  es la derivada temporal en Fi, $d/{\mathit{dt}}_b$  es la derivada temporal en Fb, y ${\mathrm{\omega }}_{b/i}^b$  es la velocidad angular del cuerpo con respecto del marco inercial.

Despejando la ecuación (5), se obtiene la derivada temporal para las velocidades:

$\left(\begin{array}{c}\dot{u}\\ \dot{v}\\ \dot{w}\end{array}\right)=\frac{F^b}{m}-\left(\begin{array}{c}p\\ q\\ r\end{array}\right)\times \left(\begin{array}{c}u\\ v\\ w\end{array}\right)$  (6)

También por segunda ley de Newton para momento angular se tiene que:

${\text{m}}^b=\frac{d{h}^b}{{\mathit{dt}}_i}=\frac{d{h}^b}{{\mathit{dt}}_b}+{\mathrm{\omega }}_{b/i}^b\times h^b=\frac{d({\text{J}\mathrm{\omega }}_{b/i}^b)}{{\mathit{dt}}_b}+{\mathrm{\omega }}_{b/i}^b\times {({\text{J}\mathrm{\omega }}_{b/i}^b)}^b$  (7)

en donde $\text{J}=\left(\begin{array}{ccc}{J}_x& 0& -J_{\mathit{xz}}\\ 0& J_y& 0\\ -J_{\mathit{xz}}& 0& J_z\end{array}\right)$  es la matriz de inercia (típicamente calculada mediante software de CAD, y que por simetría de la aeronave toma esa forma), mb es el vector de momentos de inercia resultante en Fb, y $h^b\equiv {(\text{J}{\mathrm{\omega }}_{b/i}^b)}^b$ es el momento angular en Fb.

Puesto que J no cambia con el tiempo se tiene que:

${\text{m}}^b=\text{J}\left(\begin{array}{c}\dot{p}\\ \dot{q}\\ \dot{r}\end{array}\right)+\left(\begin{array}{c}p\\ q\\ r\end{array}\right)\times \text{J}\left(\begin{array}{c}p\\ q\\ r\end{array}\right)$  (8)

Despejando se obtiene la derivada temporal de las velocidades angulares.

$\left(\begin{array}{c}\dot{p}\\ \dot{q}\\ \dot{r}\end{array}\right)={\text{J}}^{-1}\left({\text{m}}^b-\left(\begin{array}{c}p\\ q\\ r\end{array}\right)\times \text{J}\left(\begin{array}{c}p\\ q\\ r\end{array}\right)\right)$  (9)

Las ecuaciones (1), (4), (6) y (9) en conjunto, constituyen el sistema de ecuaciones necesario para el modelo dinámico.

$\{\begin{array}{c}\left(\begin{array}{c}\dot{P_N}\\ \dot{P_E}\\ \dot{P_D}\end{array}\right)={(R_v^b)}^{\text{T}}\left(\begin{array}{c}u\\ v\\ w\end{array}\right)\\ \left(\begin{array}{c}\dot{\mathrm{\varphi }}\\ \dot{\mathrm{\theta }}\\ \dot{\mathrm{\psi }}\end{array}\right)=\left(\begin{array}{ccc}1& \sin \mathrm{\varphi }\tan \mathrm{\theta }& \cos \mathrm{\varphi }\tan \mathrm{\theta }\\ 0& \cos \mathrm{\varphi }& -\sin \mathrm{\varphi }\\ 0& \sin \mathrm{\varphi }\mathit{sec}\mathrm{\theta }& \cos \mathrm{\varphi }\mathit{sec}\mathrm{\theta }\end{array}\right)\left(\begin{array}{c}p\\ q\\ r\end{array}\right)\\ \left(\begin{array}{c}\dot{u}\\ \dot{v}\\ \dot{w}\end{array}\right)=\frac{F^b}{m}-\left(\begin{array}{c}p\\ q\\ r\end{array}\right)\times \left(\begin{array}{c}u\\ v\\ w\end{array}\right)\\ \left(\begin{array}{c}\dot{p}\\ \dot{q}\\ \dot{r}\end{array}\right)={\text{J}}^{-1}\left({\text{m}}^b-\left(\begin{array}{c}p\\ q\\ r\end{array}\right)\times \text{J}\left(\begin{array}{c}p\\ q\\ r\end{array}\right)\right)\end{array}$  (10)

El vector de fuerzas resultante se puede considerar como $F=F_g+F_a+F_p$ , en donde Fg es la fuerza gravitacional, Fa son las fuerzas aerodinámicas y Fp es la fuerza de propulsión. De manera equivalente, el vector resultante de momentos es $m=m_a+m_p$ , en donde ma y mp son los momentos aerodinámicos y de propulsión respectivamente.

Siendo que la fuerza de gravedad actúa únicamente hacia abajo en el marco del vehículo, se tiene que:

$F_g^v=\left(\begin{array}{c}0\\ 0\\ \mathit{mg}\end{array}\right)$  (11)

O bien, en el marco de referencia del cuerpo,

$F_g^b=R_v^b\left(\begin{array}{c}0\\ 0\\ \mathit{mg}\end{array}\right)$  (12)

Por su parte, las fuerzas aerodinámicas son la fuerza de arrastre FDrag, que siempre actúa en sentido opuesto al vector de velocidad, la fuerza de sustentación FLift, y la fuerza aerodinámica lateral Fy.

Cabe recordar que las fuerzas aerodinámicas se generan cuando la aeronave se mueve a través del aire que la rodea. A la velocidad de la aeronave relativa al aire circundante se le denomina vector de velocidad del aire, y se denota como Va. Su magnitud es la velocidad del aire y se denota como Va. Usualmente, para que se genere sustentación, las alas deben mantener un ángulo positivo con respecto al vector de velocidad del aire, conocido como ángulo de ataque α, el cual se define mediante una rotación de mano izquierda (por la definición de ángulo de ataque positivo) en torno a jb. Al rotar el marco del cuerpo en un ángulo de α en torno a jb, se obtiene el marco de estabilidad Fs. De este modo, en el marco de estabilidad Fs, se tiene lo siguiente:

$F_a^s=\left(\begin{array}{c}-F_{\mathit{Drag}}\\ F_y\\ -F_{\mathit{Lift}}\end{array}\right)$  (13)

Que en el marco del cuerpo se representa como:

$F_a^b=R_s^b\left(\begin{array}{c}-F_{\mathit{Drag}}\\ F_y\\ -F_{\mathit{Lift}}\end{array}\right)={(R_b^s)}^{\text{T}}\left(\begin{array}{c}-F_{\mathit{Drag}}\\ F_y\\ -F_{\mathit{Lift}}\end{array}\right)$  (14)

En donde $R_s^b$  es la matriz de rotación (de mano izquierda) requerida para pasar del marco de estabilidad al marco del cuerpo.

$R_b^s=\left(\begin{array}{ccc}\cos \mathrm{\alpha }& 0& \sin \mathrm{\alpha }\\ 0& 1& 0\\ -\sin \mathrm{\alpha }& 0& \cos \mathrm{\alpha }\end{array}\right)$  (15)

Por su parte, las fuerzas y momentos aerodinámicos se crean por la distribución de la presión al rededor de las superficies sustentadoras y son función de la forma y actitud de las mismas, así como de la velocidad relativa del aire y la densidad del mismo. Estas se expresan comúnmente como:

$\begin{array}{c}{F}_{\mathit{Lift}}=\frac{1}{2}\mathrm{\rho }{V_a}^{2}S{C}_L(\mathrm{\alpha },q,{\mathrm{\delta }}_e)\\ F_{\mathit{Drag}}=\frac{1}{2}\mathrm{\rho }{V_a}^{2}S{C}_D(\mathrm{\alpha },q,{\mathrm{\delta }}_e)\\ M=\frac{1}{2}\mathrm{\rho }{V_a}^{2}Sc{C}_M(\mathrm{\alpha },q,{\mathrm{\delta }}_e)\end{array}$  (16)

En donde ρ es la densidad del aire, Va es la velocidad del aire o la velocidad del vehículo a través del aire circundante, S es el área superficial alar del vehículo, c es la cuerda media de las alas, y CLift, CD y CM son coeficientes aerodinámicos adimensionales que se ven influenciados en gran medida por el ángulo de ataque α, la velocidad de cabeceo q, y la deflexión del elevador δe. También son influenciados por el número de Reynolds y el número de Mach, sin embargo, para aeronaves pequeñas, estos permanecen aproximadamente constantes (Beard & McLain, (2012) [1]).

En general estas fuerzas y momentos son ecuaciones no lineales. Sin embargo, para ángulos de ataque pequeños, el flujo sobre las alas permanece laminar y unido. Bajo estas circunstancias, la fuerza de sustentación, de arrastre, y el momento de cabeceo se pueden modelar con precisión aceptable usando aproximaciones lineales.

Utilizando una aproximación de primer orden con una serie de Taylor, la fuerza de sustentación se puede escribir como:

$F_{\mathit{Lift}}=\frac{1}{2}\mathrm{\rho }{V_a}^{2}S\left(C_{L_0}+\frac{\partial C_L}{\partial \mathrm{\alpha }}\mathrm{\alpha }+\frac{\partial C_L}{\partial q}q+\frac{\partial C_L}{\partial {\mathrm{\delta }}_e}{\mathrm{\delta }}_e\right)=\frac{1}{2}\mathrm{\rho }{V_a}^{2}S\left(C_{L_0}+C_{L_{\mathrm{\alpha }}}\mathrm{\alpha }+C_{L_q}\frac{c}{2V_a}q+C_{L_{\mathrm{\delta }e}}{\mathrm{\delta }}_e\right)$  (17)

En donde el coeficiente $C_{L_0}$  es el valor de $C_L$  cuando $\mathrm{\alpha }=q={\mathrm{\delta }}_e=0$ . Por su parte: $C_{L_{\mathrm{\alpha }}}=\frac{\partial C_L}{\partial \mathrm{\alpha }}$ , $C_{L_q}=\frac{\partial C_L}{\partial \frac{qc}{2V_a}}$  y $C_{L_{\mathrm{\delta }e}}=\frac{\partial C_L}{\partial {\mathrm{\delta }}_e}$ ,

son cantidades adimensionales. $C_{L_{\mathrm{\alpha }}}$  y $C_{L_q}$  se llaman comúnmente derivadas de estabilidad, mientras que $C_{L_{\mathrm{\delta }e}}$  es una derivada de control. El término “derivada” proviene del hecho de que los coeficientes se originaron como derivada parcial en la aproximación por serie de Taylor. También, puesto que la única derivada parcial que no es adimensional es $\partial C_L/\partial q$ , es conveniente usar el factor $c/(2V_a)$  para hacerlo adimensional.

Empleando la misma técnica, se tiene también para el arrastre y el momento lo siguiente:

$\begin{array}{c}{F}_{\mathit{Drag}}=\frac{1}{2}\mathrm{\rho }{V_a}^{2}S\left(C_{D_0}+C_{D_{\mathrm{\alpha }}}\mathrm{\alpha }+C_{D_q}\frac{c}{2V_a}q+C_{D_{\mathrm{\delta }e}}{\mathrm{\delta }}_e\right)\\ M=\frac{1}{2}\mathrm{\rho }{V_a}^{2}Sc\left(C_{M_0}+C_{M_{\mathrm{\alpha }}}\mathrm{\alpha }+C_{M_q}\frac{c}{2V_a}q+C_{M_{\mathrm{\delta }e}}{\mathrm{\delta }}_e\right)\end{array}$  (18)

No obstante, si se quiere incorporar el efecto de pérdida aerodinámica, el cual ocurre cuando el ángulo de ataque se incrementa al punto en que existe separación de flujo en la superficie del ala y resulta en una pérdida de sustentación muy grande, se deben modificar las ecuaciones para sustentación y arrastre para una no linealidad en el ángulo de ataque.

$\begin{array}{c}{F}_{\mathit{Lift}}=\frac{1}{2}\mathrm{\rho }{V_a}^{2}S\left(C_L(\mathrm{\alpha })+C_{L_q}\frac{c}{2V_a}q+C_{L_{\mathrm{\delta }e}}{\mathrm{\delta }}_e\right)\\ F_{\mathit{Drag}}=\frac{1}{2}\mathrm{\rho }{V_a}^{2}S\left(C_D(\mathrm{\alpha })+C_{D_q}\frac{c}{2V_a}q+C_{D_{\mathrm{\delta }e}}{\mathrm{\delta }}_e\right)\end{array}$  (19)

En donde $C_L(\mathrm{\alpha })$  y $C_D(\mathrm{\alpha })$  representan funciones no lineales de α.

El coeficiente de sustentación en función del ángulo de ataque α, se puede aproximar al mezclar la función lineal de α con el coeficiente de sustentación de una placa plana, ya que para ángulos de ataque que van más allá de las condiciones de pérdida aerodinámica, las alas actúan aproximadamente como una placa plana, cuyo coeficiente de sustentación se puede modelar como:

$C_{L,\mathit{Placa}\mathit{plana}}=2\mathit{sign}(\mathrm{\alpha }){\sin }^2(\mathrm{\alpha })\cos (\mathrm{\alpha })$  (20)

De este modo, un modelo que incorpora ambos comportamientos está dado por

$C_L(\mathrm{\alpha })=\left(1-\mathrm{\sigma }(\mathrm{\alpha })\right)\left(C_{L_0}+C_{L_{\mathrm{\alpha }}}\mathrm{\alpha }\right)+\mathrm{\sigma }(\mathrm{\alpha })\left(2\mathit{sign}(\mathrm{\alpha }){\sin }^2(\mathrm{\alpha })\cos (\mathrm{\alpha })\right)$  (21)

En donde $\mathrm{\sigma }(\mathrm{\alpha })$  es la función sigmoidal que se emplea para la mezcla de ambas funciones, y tiene un corte en $\pm {\mathrm{\alpha }}_0$  y una taza de transición de MMAV.

$\mathrm{\sigma }(\mathrm{\alpha })=\frac{1+e^{-M_{\mathit{MAV}}(\mathrm{\alpha }-{\mathrm{\alpha }}_0)}+e^{M_{\mathit{MAV}}(\mathrm{\alpha }+{\mathrm{\alpha }}_0)}}{\left(1+e^{-M_{\mathit{MAV}}(\mathrm{\alpha }-{\mathrm{\alpha }}_0)}\right)\left(1+e^{M_{\mathit{MAV}}(\mathrm{\alpha }+{\mathrm{\alpha }}_0)}\right)}$  (22)

A su vez, MMAV y α0 son constantes positivas que dependen de la aeronave. Para una aeronave pequeña, el coeficiente de sustentación lineal se puede aproximar razonablemente como

$C_{L_{\mathrm{\alpha }}}=\frac{\mathrm{\pi }\mathit{AR}}{1+\sqrt{1+\frac{\mathit{AR}}{2}}}$  (23)

En donde $\mathit{AR}=b^2/S$  es la razón de aspecto de las alas, b es la envergadura y S la superficie alar.

El resultado de esta mezcla de funciones, se puede ver en la Figura 1.

Figura 1: Coeficiente de sustentación en función de α

El coeficiente de arrastre es también una función no lineal del ángulo de ataque, y existen dos contribuciones que lo conforman, el arrastre inducido y el arrastre parásito. El arrastre parásito se denota como $C_{D_P}$ , es prácticamente constante y se genera principalmente por el esfuerzo cortante del aire sobre las alas. Por su parte, el arrastre inducido es proporcional al cuadrado de la fuerza de sustentación para ángulos de ataque pequeños. De esta manera, combinando ambos se tiene:

$C_D(\mathrm{\alpha })=C_{D_P}+\frac{{\left(C_{L_0}+C_{L_{\mathrm{\alpha }}}\mathrm{\alpha }\right)}^2}{\mathrm{\pi }\mathit{AR}{e}_{\mathit{Oswald}}}$  (24)

En donde eOswald es el factor de eficiencia de Oswald, que se encuentra en el rango entre 0.8 y 1.0.

El momento de cabeceo de la aeronave es generalmente una función no lineal del ángulo de ataque y se debe determinar por vuelos experimentales o en un túnel de viento para cada aeronave específica de interés.

Para el propósito de la simulación, se usarán los modelos lineales de las 3 funciones, es decir:

$\begin{array}{c}{C}_L(\mathrm{\alpha })=C_{L_0}+C_{L_{\mathrm{\alpha }}}\mathrm{\alpha }\\ C_D(\mathrm{\alpha })=C_{D_0}+C_{D_{\mathrm{\alpha }}}\mathrm{\alpha }\\ C_M(\mathrm{\alpha })=C_{M_0}+C_{M_{\mathrm{\alpha }}}\mathrm{\alpha }\end{array}$  (25)

En donde CM(α) es el momento de cabeceo de la aeronave.

De manera similar, las fuerzas y momentos aerodinámicos laterales causan movimiento de traslación sobre el eje jb así como un movimiento rotacional en alabeo y guiñada que resulta en un cambio de dirección de la trayectoria de vuelo.

La aerodinámica lateral es influenciada significativamente por el ángulo de deriva β, así como de la velocidad de alabeo p, la velocidad de guiñada r y la deflexión de los alerones δa y del timón δr. Denotando la fuerza lateral Fy, y los momentos de alabeo y guiñada como l y N respectivamente, se tiene

$\begin{array}{c}{F}_y=\frac{1}{2}\mathrm{\rho }{V_a}^{2}S{C}_Y(\mathrm{\beta },p,r,{\mathrm{\delta }}_a,{\mathrm{\delta }}_r)\\ l=\frac{1}{2}\mathrm{\rho }{V_a}^{2}Sb{C}_l(\mathrm{\beta },p,r,{\mathrm{\delta }}_a,{\mathrm{\delta }}_r)\\ N=\frac{1}{2}\mathrm{\rho }{V_a}^{2}Sb{C}_N(\mathrm{\beta },p,r,{\mathrm{\delta }}_a,{\mathrm{\delta }}_r)\end{array}$  (26)

En donde CY, Cl y CN son coeficientes aerodinámicos adimensionales y b es la envergadura. Empleando el mismo enfoque para obtener una aproximación de primer orden con una serie de Taylor, y adimensionalizando, se tiene:

$\begin{array}{c}{F}_y=\frac{1}{2}\mathrm{\rho }{V_a}^{2}S\left(C_{Y_0}+C_{Y_{\mathrm{\beta }}}\mathrm{\beta }+C_{Y_p}\frac{b}{2V_a}p+C_{Y_r}\frac{b}{2V_a}r+C_{Y_{\mathrm{\delta }a}}{\mathrm{\delta }}_a+C_{Y_{\mathrm{\delta }r}}{\mathrm{\delta }}_r\right)\\ l=\frac{1}{2}\mathrm{\rho }{V_a}^{2}Sb\left(C_{l_0}+C_{l_{\mathrm{\beta }}}\mathrm{\beta }+C_{l_p}\frac{b}{2V_a}p+C_{l_r}\frac{b}{2V_a}r+C_{l_{\mathrm{\delta }a}}{\mathrm{\delta }}_a+C_{l_{\mathrm{\delta }r}}{\mathrm{\delta }}_r\right)\\ N=\frac{1}{2}\mathrm{\rho }{V_a}^{2}Sb\left(C_{N_0}+C_{N_{\mathrm{\beta }}}\mathrm{\beta }+C_{N_p}\frac{b}{2V_a}p+C_{N_r}\frac{b}{2V_a}r+C_{N_{\mathrm{\delta }a}}{\mathrm{\delta }}_a+C_{N_{\mathrm{\delta }r}}{\mathrm{\delta }}_r\right)\end{array}$  (27)

Estas fuerzas y momentos están alineados con el marco del cuerpo, por lo que no requieren rotaciones para implementarse en las ecuaciones de movimiento. El coeficiente $C_{Y_0}$  es el valor del coeficiente de fuerza lateral CY cuando $\mathrm{\beta }=p=r={\mathrm{\delta }}_a={\mathrm{\delta }}_r=0$ . Para aeronaves que son simétricas en el plano ib - kb, $C_{Y_0}$  es típicamente 0. Los coeficientes $C_{l_0}$  y $C_{N_0}$  se definen de manera similar y también son típicamente cero para aeronaves simétricas.

En cuanto a la fuerza de empuje, se puede desarrollar un modelo simple del empuje generado por la hélice, aplicando el principio de Bernoulli para calcular la presión adelante y atrás de la hélice, y posteriormente aplicando la diferencia de presiones al área de la hélice. Así pues, empleando la ecuación de Bernoulli, la presión adelante y atrás de la hélice se puede escribir como:

$\begin{array}{c}{P}_{\mathit{adelante}}=P_0+\frac{1}{2}\mathrm{\rho }{V}_a^2\\ P_{\mathit{atras}}=P_0+\frac{1}{2}\mathrm{\rho }{V}_{\mathit{Salida}}^2\end{array}$  (28)

En donde P0 es la presión estática y VSalida es la velocidad del aire al salir de la hélice, la cual se puede escribir como $V_{\mathit{Salida}}=k_{\mathit{motor}}{\mathrm{\delta }}_t$  en donde kmotor es una constante que especifíca la eficiencia del motor y δt es el comando de modulación del ancho de pulsos del motor, la cual está linealmente relacionada con la velocidad angular de la hélice.

Sea SProp el área que barre la hélice, entonces el empuje que produce el motor está dado por

$F_{X_p}=S_{\text{Prop}}{C}_{\text{Prop}}\left(P_{\mathit{atras}}-P_{\mathit{adelante}}\right)=\frac{1}{2}\mathrm{\rho }{S}_{\text{Prop}}{C}_{\text{Prop}}[{(k_{\mathit{motor}}{\mathrm{\delta }}_t)}^2-{V_a}^2]$  (29)

Por lo tanto

$F_p=\frac{1}{2}\mathrm{\rho }{S}_{\text{Prop}}{C}_{\text{Prop}}\left(\begin{array}{c}{(k_{\mathit{motor}}{\mathrm{\delta }}_t)}^2-{V_a}^2\\ 0\\ 0\end{array}\right)$  (30)

Al girar la hélice se aplica fuerza al aire que pasa a través de la hélice, incrementando el momento del aire y generando propulsión. En consecuencia, fuerzas iguales pero en sentido opuesto, se aplican a la hélice por el aire. El efecto neto de estas fuerzas resulta en un torque en torno al eje de rotación de la hélice. El torque aplicado por el motor a la hélice (y posteriormente al aire), conlleva un torque igual y opuesto, de la hélice al motor que está fijo en la estructura de la aeronave. Este torque es contrario a la dirección de rotación de la hélice y proporcional al cuadrado de la velocidad angular $T_p=-k_{T_p}{(k_{\mathrm{\Omega }}{\mathrm{\delta }}_t)}^2$ , en donde $\mathrm{\Omega }=k_{\mathrm{\Omega }}{\mathrm{\delta }}_t$  es la velocidad de la hélice y $k_{T_p}$  es una constante determinada experimentalmente. El momento a causa del sistema de propulsión es por consecuencia:

$m_p=\left(\begin{array}{c}-k_{T_p}{(k_{\mathrm{\Omega }}{\mathrm{\delta }}_t)}^2\\ 0\\ 0\end{array}\right)$  (31)

El efecto del torque de la hélice es usualmente bajo y por lo tanto despreciable en este modelo.

Con las ecuaciones (12), (14), (16) y (30) se puede determinar por completo el vector de fuerzas.

$F=F_g+F_a+F_p=R_v^b\left(\begin{array}{c}0\\ 0\\ \mathit{mg}\end{array}\right)+{(R_b^s)}^{\text{T}}\left(\begin{array}{c}-F_{\mathit{Drag}}\\ F_y\\ -F_{\mathit{Lift}}\end{array}\right)+\frac{1}{2}\mathrm{\rho }{S}_{\text{Prop}}{C}_{\text{Prop}}\left(\begin{array}{c}{(k_{\mathit{motor}}{\mathrm{\delta }}_t)}^2-{V_a}^2\\ 0\\ 0\end{array}\right)$  (32)

Sustituyendo en la ecuación (6) se obtiene la derivada temporal de las velocidades.

En cuanto a los momentos, se tiene que:

$m=m_a+m_p=\left(\begin{array}{c}l\\ M\\ N\end{array}\right)+\left(\begin{array}{c}-k_{T_p}{(k_{\mathrm{\Omega }}{\mathrm{\delta }}_t)}^2\\ 0\\ 0\end{array}\right)$  (33)

Y con ello, sustituyendo en (9) se obtiene la derivada temporal de las velocidades angulares.

Con esto quedan definidas todas las ecuaciones requeridas para implementar el modelo dinámico del MAV de alas fijas.

En general, el modelo dinámico de la aeronave se puede representar como:

$\dot{x}=f(x,u)$  (34)

en donde x son las variables de estado, u las variables de control o entradas, y $f\mathrm{:}{R}^n\times R^m\to R^n$  se especifíca por el conjunto de ecuaciones (10).

El modelo dinámico de 6 grados de libertad de una aeronave con alas fijas se desarrolló en Scilab, tomando los valores de coeficientes y parámetros del MAV descritos en el Anexo I. Al aplicar métodos numéricos a la ecuación (34) usando la función ode() de Scilab, para una simulación de 1 minuto con un paso temporal de 1 segundo, y tomando como valores de variables de estado y control a $x={\left(\mathrm{0,0,0,30,0,0,0,0.349,0,0,0,0}\right)}^{\text{T}}$  y $u={\left(\mathrm{0.001745,0}\mathrm{.001745,0,1}.5\right)}^{\text{T}}$ , se obtuvieron los resultados mostrados en la Figura 2.

Figura 2: Resultado de la simulación del modelo dinámico.

Referencias

1. Beard, R. W., McLain, T. W. (2012). Small unmanned aircraft: Theory and practice. Princeton university press.