0% encontró este documento útil (0 votos)
22 vistas23 páginas

Control de Cabeceo en Aeronaves

Sistema de control de cabeceo de un avión

Cargado por

tellefsen879
Derechos de autor
© © All Rights Reserved
Nos tomamos en serio los derechos de los contenidos. Si sospechas que se trata de tu contenido, reclámalo aquí.
Formatos disponibles
Descarga como ODT, PDF, TXT o lee en línea desde Scribd
0% encontró este documento útil (0 votos)
22 vistas23 páginas

Control de Cabeceo en Aeronaves

Sistema de control de cabeceo de un avión

Cargado por

tellefsen879
Derechos de autor
© © All Rights Reserved
Nos tomamos en serio los derechos de los contenidos. Si sospechas que se trata de tu contenido, reclámalo aquí.
Formatos disponibles
Descarga como ODT, PDF, TXT o lee en línea desde Scribd

INSTITUTO UNIVERSITARIO AERONÁUTICO

SISTEMA DE CONTROL
PARA EL CABECEO DE UN AERONAVE

INFORME TÉCNICO: CCA-001/23


INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 2

1. Introducción

El objetivo general de este trabajo es presentar el diseño de un sistema de control de cabeceo


(rotación sobre eje transversal) de una aeronave, es decir se buscara la modificación del angulo
de ataque mediante el timón de profundidad o elevador.
Tomaremos como referencia los datos del Boeing 747-400 para el desarrollo de las ecuaciones
de movimiento.
El diseño se basara en introducir un angulo de ataque deseado que sera leído por el conversor, el
cual mandara un voltaje proporcional que se usara como señal de referencia y sera comparada
con la señal del sensor de salida como realimentación. El diagrama de bloques dibujado utiliza
esta diferencia de señal y obedeciendo la ley de control que diseñamos actuá sobre la aeronave
mediante el elevador.

Imagen 1: Rotacion respecto al eje transversal

2. Arquitectura del diagrama de bloques

PERTURBACIONES

θref Vref ev θout


Vcontrol
CONVERSOR + CONTROLADOR PLANTA
-

Vsensor
SENSOR
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 3

2.1. Conversor

Traduce el angulo de referencia a la cual se quiere equilibrar el avión en una señal de voltaje de
referencia. Supondremos el uso de un conversor proporcional de relación unitaria. Es decir 1 grado
equivale a 1 volt.

Vref (s) v
Función de transferencia conversor: FTconv= =1 [ ]
θ ref (s ) º

2.2. Sensor

Utilizaremos un giroscopio que mida la inclinación del avión y devuelva señal del voltaje como
realimentación del sistema. Supondremos una relación proporcional unitaria.
El tiempo de retraso en la medición se aproxima a una función de transferencia de primer orden.
Para simplificar los cálculos, supondremos un sensor lo suficientemente bueno que nos genere una
constante de tiempo de retraso Ts casi nulo.

Vsensor (s) 1 v
Función de transferencia sensor: FTsensor= = ≈1 [ ]
θ out(s) Ts∗S+1 º

2.3. Planta

La planta consta de un diagrama de bloques propio:

Vcontrol + Vdamp δ θout


- ACTUADOR DINAMICA

Vsas
SAS

2.3.1.Actuador
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 4

El actuador es el mecanismo que se encarga de que las señales que genera el sistema de control se
transformen, en este caso, en movimientos que produzcan las correspondientes deflexiones.

Se utiliza un actuador electrohidraulico para el elevador, de relación proporcional unitaria, es decir un


grado de deflexión de elevador por voltio, y un retardo en la actuación de primer orden.
Supondremos un valor de retardo Te= 0.1s

δ ( s) 1 v
Función de transferencia de actuador: FTactuador= = [ ]
Vdamp (s) TeS +1 º

2.3.2. SAS (Sistema de estabilidad aumentada)

Se encarga de que el avión sea estable en todos sus modos de oscilación, es decir su respuesta ante la
excitación es convergente.
El SAS se encarga de que exista un amortiguamiento de cabeceo para que la aeronave se comporte de
acuerdo a como el piloto lo esperaría. Para aumentar el amortiguamiento el actuador se realimenta con
una señal proporcional a la velocidad de cabeceo. Un sensor de velocidad angular genera la señal que
luego sera amplificada por un factor de proporcionalidad Ksas.
Supondremos que analizamos un Boeing 747-400 con un factor tipico de Ksas=0.5, para lo que el
bloque resulta

Vdamp( s)
SAS= =0.5 S
θ out ( s)

2.3.3. Dinámica del avión

Relaciona el angulo de cabeceo de la aeronave con la deflexión del elevador.


La función de transferencia de este bloque es un modelo matematico deducido con las ecuaciones de
movimiento longitudinal.

2.3.3.1. Modelo Matemático

Se plantean las ecuaciones de movimiento que nos permitirán hacer nuestra función de transferencia,
θ (s)
que relaciona la deflexión del elevador con el angulo de actitud.
δ (s)
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 5

Bajo ciertas hipótesis simplificamos las ecuaciones de movimiento y obtenemos 3 ecuaciones


longitudinales que gobiernan el angulo de cabeceo.
Asumimos que la aeronave se encuentra en vuelo crucero a velocidad y altitud constantes, que su masa
permanece constante, que es un cuerpo rígido y que los ejes se encuentran fijos al cuerpo del avión con
origen en el centro de gravedad.

Imagen 2: Diagrama de fuerzas y momentos actuantes sobre la aeronave

Ecuaciones de movimiento longitudinal:

1 CL
C D u+T u̇+ (C D α −C L ) α + θ = 0
2 2

C Lα
C L u+ α +T α̇ −T θ̇ = 0
2
. 2
Cmu u+C m α α +C m D α T α + Cm D θ T θ̇ −h T θ̈ = −C m δ δ −CmD δ Ṫ δ̇
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 6

δ: Deflexión del elevador


θ: Angulo de cabeceo
α: Angulo de ataque
CL: Coeficiente de sustentación
CDα: Pendiente del coeficiente de resistencia
Clα: Pendiente del coeficiente de sustentación
Cmα: Momento de cabeceo por angulo de ataque
CmDα: Momento de cabeceo por angulo de ataque
CmDθ: Momento de cabeceo por velocidad de cambio
Cmδ: Potencia de elevador
CmDδ: Momento de cabeceo por cambio en la velocidad de elevador
2
2k y
h= : Constante del avión que expresa una relación de inercias
μ c2
m
T= : Constante de tiempo para la condición de equilibrio del avión
ρ SV

Puesto que esperamos que el sistema de control estabilice rápidamente, solo nos interesa el movimiento
en las primeras oscilaciones. En los primeros segundos de movimiento predomina el “short mode”
(corto periodo) caracterizado por oscilaciones fuertemente amortiguadas y alta frecuencia. En este
periodo es despreciable la variacion de la velocidad “u” frente a las variaciones del angulo de cabeceo
y de ataque.
Despreciaremos el valor de CmDδ por ser generalmente pequeño.
Con estas proposiciones las ecuaciones de movimiento de cabeceo se ven reducidas a:

C Lα
α +T α̇ −T θ̇ = 0
2
. 2
Cm α α +C m D α T α +C m D θ T θ̇−hT θ̈ = −C m δ δ

α (s)
Aplicamos la transformada de Laplace a las ecuaciones y la escribimos en terminos de y
δ ( s)
θ ( s)
. En forma matricial me quedan:
δ ( s)
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 7

A=

C Lα
α +T s +Ts -Ts
2
2 2
Cm α +Cm D α T s Cm D θ T s−h T s

( )
α (s)
B= δ (s) C= ( 0 )
θ (s) −C M δ
δ (s)

Siendo la ecuación matricial : A * B = C

θ ( s)
La función de transferencia para se despeja con la regla de Cramer:
δ ( s)

C Lα
+Ts 0
2
θ ( s) C m α +C mD α Ts −Cm δ
=
δ ( s) C Lα
+Ts −Ts
2
2 2
C m α +C mD α Ts C mD θ S−hT S

Para el desarrollo de este trabajo analizamos un Boeing 747-400.


A continuación presentamos todos los datos de esta aeronave necesarios para realizar el modelo
matemático.
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 8

V= 205,13 m/s ρ= .653kg/m³ Ky=12.47m h=.0431


μ= 104.07 T=4.22s CL= 0.4 CDα=0.2
Clα= 4.4 Cmα =-1 CmDα = 0.017 CmDθ = -0.085
Cmδ =-1.3 CmDδ = -0.005

Reemplazando estos valores obtenemos la función de transferencia del bloque DINAMICA

4.22 s +2.2 0
θ ( s) −0.072 s−1 −1.3 1.69 s+0.84
= =
δ ( s) 4.22 s+ 2.2 −4.22 s 2
s (s + 1.18 s+1.59)
−0.072 s−1 0.768 s2−0.359 s
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 9

2.4. Diagrama de bloque de planta

Con las funciones de transferencia encontradas dibujamos el diagrama de bloque de planta

Vdamp 1 δ 1.69s+0.84 θout


Vcontrol +
0.1s+1 s(s²+1.18s+1.59)
-
actuador dinamica

Vsas
0.5s
SAS

θ out(s) actuador∗dinamica 16.90( s+ 0.50)


= =
Vcontrol ( s) 1+ actuador∗dinamica∗SAS s (s +9)(s ²+2.17 s+2.23)

2.5.Diagrama de bloque final

Reemplazando las respectivas funciones de transferencia y utilizando la reducción del diagrama de


bloque de planta obtenemos

θref Vref Vcontrol θout


ev 16.9(s+0.50)
1 v/º + CONTROLADOR
- s(s+9)(s²+2.17s+2.23)

Vsensor

1 v/º
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 10

3. Análisis del sistema

Calcularemos las distintas funciones de transferencia del sistema y posteriormente su estabilidad.

*R(s): Entrada del sistema


*E(s): Señal de error
*B(s): Señal de retroalimentación
*C(s): Señal de salida del sistema
*G(s):Función de transferencia de la trayectoria directa
*H(s): Función de transferencia de retroalimentación

3.1. Función de transferencia de lazo abierto

Los sistemas de control de lazo abierto no utilizan una señal de retroalimentación.


La función de transferencia es la relación entre la señal del sensor B(s) y la señal de error actuante E(s).
B(s)
FTLAZO ABIERTO= = G(s)H(s)
E (s )

Para este análisis tomamos nuestro diagrama de bloques sin tomar en cuenta el controlador.

G(s)= 16.9 (s+0.50)


s(s+9)(s²+2.17s+2.23) ( Función transferencia de trayectoria directa)

H(s)= 1 (realimentación unitaria)

𝐹𝑇|𝐿azo𝐴𝑏𝑖𝑒𝑟𝑡𝑜 = G(s)H(s)= 16.9 (s+0.50)


s(s+9)(s²+2.17s+2.23)

Los ceros de la función de transferencia de lazo abierto ocurren en aquellos valores de la variable
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 11

compleja s para los cuales el numerador se anula. Por otro lado, los polos de la función de
transferencia de lazo abierto son los valores de s que hacen cero el denominador.

Una señal de entrada de tipo impulso (Vref(s)=1) simula una perturbación del sistema y una señal de
entrada de tipo escalón (Vref(s)=1/s) simula el comportamiento del piloto automático.

Encontramos polos y ceros con matlab:

cero

polos

Observamos que hay un polo en 0, por lo que estoy tratando con un sistema tipo 1.

Lo cual quiere decir que tengo un error de estado estable.

3.2. Función de transferencia lazo cerrado

C (s) G(s) G( s)
FT lazo cerrado= = =
R (s ) 1+G(s)∗H (s ) 1+ FTlazo abierto
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 12

3.2.1 Error de régimen

Debido a que la función de transferencia de lazo abierto es de tipo 1, el error de régimen para una
entrada de escalón es cero, para una polinómica de grado dos o mas es infinito.

Para una entrada de tipo rampa: e reg =lim s ¹∗E(s)


s →0

1
Donde: E(s)= R(s)* con R(s)= 1/s²
1+ FT lazo abierto

En la tabla podemos observar que para un sistema de tipo 1, el error de régimen para una entrada de
tipo escalón es nulo, para una entrada en rampa es del tipo 1/k y para una entrada del tipo R(s)=1/s³ es
infinito.

R(s)=1/s R(s)=1/s² R(s)=1/s³

Para verificar utilizamos Matlab, presentamos la salida del sistema lazo cerrado respuesta a una entrada
escalón.

Imagen 3: Sistema lazo cerrado en Simulink


INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 13

Imagen 4: Respuesta grafica a la señal escalon

Observamos que el sistema tiene un tiempo de asentamiento muy tardío, por lo que realizaremos un
compensador que pueda generar un tiempo de establecimiento menor o igual a 10 segundos.

3.2.2. Función de transferencia con compensador

La función de transferencia en lazo cerrado con un compensador proporcional K.

θout
+- K 16.9(s+0.50)
s(s+9)(s²+2.17s+2.23)
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 14

K (16.9 s+8.45)
FT lazocerrado= 4
s +11.17 s ³+21.76 s ²+(16.9 K + 20.11)s +8.45 K

Siendo K la ganancia variable del sistema.

3.2.3. Estabilidad de Routh (Estabilidad Transistoria)

Determina si existe alguna raíz o polo en el semiplano derecho del plano complejo “s”, donde si al
menos existe una raíz el sistema es inestable.

Caso contrario, si NO hay ninguna raíz en el semiplano derecho el sistema es estable.

Para determinar la estabilidad mediante Routh, se toma la ecuación característica o denominador de la


función de transferencia de lazo cerrado.

Ecuación Característica: D(s)= s⁴+11.17s³+21.76s²+(16.9K+20.11)s+8.45K

El criterio de Routh se resuelve siguiendo el procedimiento:

S⁴ 1 21.76 8.45K

S³ 11.17 16.9K+20.11

S² 20+20.24K 8.45K
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 15

342.1 K ²+839.43 k + 402.2


S
20+20.24 k

S⁰ 8.45K

El número de cambios de signo de la primera columna resultante del criterio de Routh – Hürwitz, nos
da la cantidad de elementos que están en el semiplano derecho.

Si encontramos cambios de signo, el sistema será inestable.


Si todos los elementos tienen el mismo signo, el sistema será asintóticamente estable.
Tendremos tantos polos en el semiplano positivo como variaciones de signo en la primera columna.

Las ecuaciones encontradas nos ayudan a determinar para qué valores de K el sistema pasará de estable
a inestable (ganancia K límite).

Resolvemos: 20+20.24K > 0 y 342.1K²+839.43K+402.2 > 0

K > -0.988 K> -0.653 ; K> -1.8

La ganancia critica esta definida por el limite inferior:

Kcr = -1.8
3.2.4. Polos Objetivo

Con Matlab buscamos los valores que deberia tener mi polo dominante para cumplir un Mp menor al
5% y un tiempo menor a 5 segundos.

Imagen 5: comandos y resultados


matlab para busqueda de polos
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 16

Observamos que:

* Para un sobrepico de 5% se necesita una razón de amortiguamiento de 7%

* Frecuencia natural aproximadamente 2

* Polos deseados en -1.4 ±1.4 i

Si los polos del sistema contienen a los polos objetivo entonces se cumplen los requerimientos.

3.3. Lugar de raices

El estudio de esto conlleva hallar los polos y ceros (anteriormente encontrados), asíntotas, lugar en el
eje real, en el eje imaginario, puntos de ruptura, ángulos de entrada y salida y margen de ganancia.

3.3.1. Asíntotas

El número de asíntotas se determina haciendo la diferencia entre la cantidad de


polos y ceros de la función de transferencia de lazo abierto.

NroAsintotas = NroPolos − NroCeros = 4 – 1 = 3

Ángulos que forman las asíntotas:

180 º
θasintotas = (2k+1) = 60º ; 180º ; 300º
nroAsintotas

3.3.2. Lugar de raíces en eje real

El lugar de raíces en el eje real se ubica siguiendo la regla de número impar de raíces y/o polos al
desplazarnos de derecha a izquierda a lo largo del eje real.

En este caso de 0 a -0.5 y de -9 a -infinito.

Re
-infinito -9 -0.5 0
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 17

3.3.3 Puntos de ruptura

En un punto del lugar de las raíces se pueden juntar varios polos del sistema. Son los puntos de ruptura
dGla(s)
y se buscan entre las raíces de la ecuación =0
ds

Derivamos la función de transferencia con respecto a s y tomamos el numerador.

s⁴+8.11s³+12.85s²+7.26s+3.35=0

Las raíces son: -1.39 ; -6.22 ; -0.25+0.57i ; -0.25-0.57i

Las soluciones reales no pertenecen al lugar de raíces por lo que no serán puntos de ruptura.

3.3.4. Grafico de raices

Obtenemos el siguiente grafico de raices con Matlab


INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 18

4. Diseño del compensador

Existen dos tipos de compensadores que se pueden agregar al sistema para modificar el
comportamiento del mismo, de adelanto y de atraso.

La red de adelanto mejora la respuesta transitoria y cumple con los requerimientos de respuesta
temporal, para esto debemos modificar el lugar geométrico de raíces de modo que el mismo contenga
al polo deseado y luego a través de un ajuste de ganancia ubicamos las raíces de la ecuación
característica en los puntos que nos permitan tener el comportamiento transitorio en términos de
frecuencia natural y coeficiente de amortiguamiento que nos proporcione la respuesta temporal
deseada. La red de atraso mejora la respuesta en estado estacionario sin modificar la respuesta
transitoria para cumplir el requerimiento de error en régimen permanente, esto se consigue aumentando
la ganancia directa en estado estacionario sin modificar la ubicación de los polos dominantes de lazo
cerrado.
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 19

Nuestro sistema es de tipo 1, por lo que el error en régimen permanente ante una entrada escalón seria
nulo y solo necesitaremos un compensador de adelanto.

1
s+
Ts+1 T
La forma general de un compensador de adelanto es Gc(s)=Kc α =Kc
α Ts+1 1
s+
αT
con T > 0 y 0 < α < 1

4.1.1. Pertenencia de polos deseados en lugar de raíces

Garantizo la existencia de estos puntos en el nuevo lugar de raíces del sistema. El nuevo polo y
cero se introduce en el compensador.

Tomo el angulo encerrado por las rectas que unen a un polo objetivo con el nuevo polo y cero.

Los polos anteriormente encontrados para un Mp menor al 5% y un tiempo menor a 5 segundos son
-1.4 ±1.4 i

Para calcular el angulo:

Σ(ángulos de los polos) − Σ(ángulos de los ceros) = 180 ±k 360

Φ - (θ1+θ2+θ3+θ4+θ5)= -180º ±k 360

θ1= 90º+tan-1(1.4/1.4)= 135º


θ2= 180º-tan-1((1.4-1)/(1.4-1.1) )= 126.9º
θ3= 180º-tan-1 ((1.4+1)/(1.4-1))= 97º
θ4= tan-1(1.4/0.6) = 66.8º
θ5= tan-1(1.4/7.6)= 10.4º
Φ= 180º-tan-1(1.4/0.9)= 122.7º

Φ – (θ1+θ2+θ3+θ4+θ5) = -301 no es igual a -180º,


por lo tanto busco un θm que satisfaga esta ecuacion.
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 20

θm = -180+ 301 = 121º


θm= 121º
4.1.2. Compensador proporcional derivativo
Utilizando este compensador, el denominador tiende a infinito y queda de la forma
Gc(s)= Kc (s+ 1/T).
1
El valor de 1/T esta dado por la geometría y se calcula como a=−x+ y con x igual a la parte
tan φ
real del polo objetivo e y igual a la parte imaginaria.
En este caso x=-1.4 ; y= 1.4
Resolviendo esto obtengo 1/T aproximadamente igual a 2.

4.2.2. Calculo de la ganancia

Aplicamos la condición del modulo de la ganancia para encontrar la ganancia Kc

Π( distanciadel punto a los polos)


Kc Kla= con Kla= 16.9 y polo= -1.4 + 1.4i
Π( distancia del punto a los ceros)

Resolviendo este obtengo Kc= 1

El bloque del compensador me queda: C= 1 (s+2)

4.2.3. Respuesta Matlab

Con Matlab realizamos la búsqueda de diferentes respuestas a distintos Kc.


INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 21

Podemos observar que el calculo realizado, con el que obtenemos un Kc= 1 (en la grafica 0.97), no nos
da un resultado optimo porque el polo se encuentra muy cerca al origen, por lo que a continuación
rediseñaremos el compensador utilizando otros polos objetivo.

4.2.4. Rediseño del compensador

Realizando algunos cambios en los datos de matlab encontramos los polos objetivos que nos
dan un resultado optimo.

P= −4±5 i
ε= 0.64
ωn= 6.4
Kc= 2.65
Polos= −4±5 i ; -2,6 ; -0.5
Mp= 4.7%
ts= 4.9s

Vcontrol ( s)
El bloque compensador final es C= =2.65( s+2)
ev (s)
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 22

4.3. Diagrama bloque final

5. Conclusiones

Podemos concluir que haciendo un mejor análisis y estudio del problema, se puede obtener un
compensador mas preciso, mejorando el nivel amortiguamiento y el tiempo para optimizar el sistema
de la aeronave.
También podemos observar que con un poco mas de análisis se podría hacer un mejor análisis de Kc.
El uso de un buen compensador puede improvisar mucho la respuesta del sistema, sin embargo puede
ser complicado encontrar el ajuste preciso para que el sistema responda de la mejor forma posible.
Un buen manejo de las teorías de control es indispensable si se desea hacer un sistema de control
optimo.

6. Bibliografia

• Ingenieria de Control moderna – Katsuhiko Ogata


• Ingenieria de control- Universidad Navarra
• Control tutorials Matlab/Simulink
• Analisis Numerico de la Estabilidad y Control de una Aeronave – Emilio Dufeo Delarze

https://dademuch.com/2018/03/26/error-en-estado-estable-de-un-sistema-de-control/

https://zaguan.unizar.es/record/6996/files/TAZ-PFC-2012-150.pdf

file:///home/thais/Descargas/711-3288-1-PB.pdf
INSTITUTO UNIVERSITARIO AERONÁUTICO
INFORME TÉCNICO CCA-001/23 pagina: 23

https://biblus.us.es/bibing/proyectos/abreproy/60320/fichero/Rub%C3%A9n+Gonz
%C3%A1lez+Manzanares+PFC.pdf

https://bibdigital.epn.edu.ec/bitstream/15000/656/1/CD-1596%282008-07-15-01-29-21%29.pdf

https://oa.upm.es/70731/1/TFM_JUAN_DIEGO_CABRERA_MATEOS.pdf

También podría gustarte