Modelos Dinámicos en Ingeniería Aeroespacial
Modelos Dinámicos en Ingeniería Aeroespacial
IO
C
C
RU
ST
N
O
C
Sistemas Dinámicos
EN
IO
C
C
1. Modelos Dinámicos 1
1.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
RU
1.2. Modelos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1. Modelo Conceptual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.2. Cómputos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.3. Modelo Matemático . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3. Análisis del Modelo Matemático . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3.1. Variables de Entrada e Invarianza Temporal del Modelo . . . . . . . . . . . . . 7
ST
2. Respuesta a Perturbaciones 49
2.1. Funciones de Transferencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.1.1. Entradas y Salida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.1.2. Matriz de Transferencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
EN
N
2.4.2. Transformada de Fourier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
2.5. Señales Aleatorias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
2.5.1. Análisis Estático . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
IO
2.5.2. Análisis Dinámico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
2.5.3. Respuesta a Señales Aleatorias . . . . . . . . . . . . . . . . . . . . . . . . . . 87
2.6. Procesamiento Digital de Señales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
2.6.1. Muestreo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
C
2.6.2. Transformada Discreta de Fourier . . . . . . . . . . . . . . . . . . . . . . . . . 91
2.6.3. Procedimiento de Welch para Señales Aleatorias . . . . . . . . . . . . . . . . . 91
C
3. Modelado de Sistemas Mecánicos 93
3.1. Modelos de Parámetros Concentrados . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
3.1.1. Ecuación de Movimiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
RU
3.1.2. Ecuación Matricial de Movimiento . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.1.3. Sistemas No Amortiguados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
3.2. Modelos de Parámetros Distribuidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
3.2.1. Deformaciones Unidireccionales . . . . . . . . . . . . . . . . . . . . . . . . . . 114
3.2.2. Vibraciones Transversales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
3.3. Métodos Aproximados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
ST
A. Apéndices 145
A.1. Dinámica del Avión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
O
IO
C
Modelos Dinámicos
C
RU
ST
N
O
C
EN
1
Introducción
Toda nuestra existencia se desarrolla en un contexto dinámico. Nuestro entorno
evoluciona en mayor o menor medida con el transcurso del tiempo, aunque no todo es
perceptible de forma cotidiana. Hay procesos que duran años, siglos o milenios; mientras
que otros evolucionan en fracciones de segundo.
Podemos notar que a la larga todo proceso dinámico tiende a equilibrarse en ciertas
condiciones que dependen de su propia naturaleza y de las interacciones del entorno,
N
a las que llamaremos perturbaciones. Es oportuno aclarar que el equilibrio solo puede
darse si las perturbaciones se mantienen constantes, o si se neutralizan entre sı́.
IO
Como ejemplo simple consideremos lo que ocurre al calentar agua en un recipiente.
El proceso dinámico en este caso es de naturaleza térmica, y se manifiesta a través de
un cambio en la temperatura del contenido.
C
La perturbación es el flujo de calor que entregamos al recipiente, y la pérdida del mismo
por conducción con el aire circundante. La condición de equilibrio se alcanza cuando el
calor perdido es igual al calor entregado por el calentador, lo que da como resultado que
C
la temperatura se mantenga constante.
Modelos
O
Por ejemplo, si en nuestro hogar tenemos un problema con el agua caliente analizamos la
relación entre elementos tales como el calefón, el flujo de gas y agua, tuberı́as, válvulas,
etc. En ese contexto en nuestra mente el calefón es simplemente una entidad en la cual
entra un caudal de agua frı́a y un caudal de gas (elementos que podemos reducir a
N
Si, a modo de ejemplo, el problema resulta ser un mal funcionamiento del calefón,
probablemente dejaremos de lado elementos del primer modelo conceptual (tuberı́as,
IO
suministro de agua, válvulas, etc.) y comenzaremos a analizar ese dispositivo en base a
un modelo más detallado, pero que seguirá siendo un modelo conceptual por mas detalle
que se incluya en el mismo. Pensaremos en el intercambiador de calor (serpentina), la
válvula de corte por bajo caudal, la válvula de gas, sus acoples, etc.
C
Nuestra mente nunca será capaz de manejar la realidad si no es a través de modelos.
C
Pero si el modelo incluye todos los elementos relevantes de la problemática bajo
análisis una captura completa de la realidad, aunque fuese posible, resultarı́a totalmente
innecesaria e incluso no recomendable por resultar ineficiente para realizar cualquier tipo
RU
de análisis concreto.
1.2.2. Cómputos
Los cómputos asociados a los ejemplos que se presentarán pueden realizarse mediante
MatlabTM o GNU Octave. Para los ejemplos del capı́tulo 2 se requiere la librerı́a control
ST
toolbox.
permite construir un modelo más definido que un cierto conjunto de ideas dispersas a
nivel mental, y permite aportar conceptos al modelo de forma más sistemática.
O
Continuando con el ejemplo del calefón, en este punto se asume como modelo
conceptual que el calefón es un elemento puntual (ya que su geometrı́a real es
irrelevante) y que cumple con los siguientes principios y relaciones termodinámicas:
C
E = ce m T
1 la temperatura no es un escalar (T ∈ R) sino un campo escalar (T : R3 → R). Por lo tanto en este contexto T se refiere
P = k (Tcomb − T )
N
En esta expresión se plantea que toda la masa de agua se encuentra a una misma
temperatura T , lo cual claramente es falso. Si embargo, el efecto del gradiente de
temperatura dentro del calentador puede absorberse en el coeficiente k .
IO
3 - conservación de energı́a La energı́a transferida se convierte en energı́a térmica
de la masa de agua. Pero esta masa de agua se renueva, dado que en el calefón existe
un caudal másico de agua Qm que ingresa con temperatura T0 , y un caudal igual de
C
agua que sale con temperatura T .
Haciendo un balance de energı́a podemos decir que la diferencia entre la energı́a de
salida y la de entrada es energı́a que se almacena en el agua contenida en el calefón,
C
produciendo un cambio en su energı́a interna:
dT
ce m = P − ce Qm (T − T0 )
RU
dt
= k (Tcomb − T ) + k (T0 − T0 ) − ce Qm (T − T0 )
Ṫ = a (Tcomb − T0 ) + b (Tcomb − T0 )
1 k k
a=− + Qm , b=
m ce ce m
ẋ = a x + b u (1.1)
O
derivadas parciales.
Pero en la mayorı́a de los modelos conceptuales estos campos se sustituyen por
propiedades puntuales (valores medios, posiciones y velocidades del centro de masa
en el caso de cuerpos rı́gidos, coordenadas modales en el caso de estructuras flexibles,
etc.), por lo cual el resultado normalmente puede ser reducido a una o más ecuaciones
2x y u serı́an apartamientos respecto de la temperatura del ambiente en lugar de temperaturas absolutas
N
Un modelo sencillo bastante conocido es el del péndulo plano (un modelo conceptual
que puede aplicarse por ejemplo a una carga suspendida por una grúa):
IO
C
C
RU
ST
distancia entre este y el punto de anclaje del cable en la grúa; a partir del principio de
conservación del momento angular se llega al siguiente modelo matemático:
C
g 1
θ̈ + sin θ = f (1.2)
l ml
Si consideramos el caso en el cual la fuerza f es de origen aerodinámico, debemos
modelar su relación con la velocidad relativa entre la masa y el aire v .
v 1
f (v) = − ρ v(t)2 cx S , v = w + lθ̇ (1.3)
|v| 2
donde ρ es la densidad del aire, w la velocidad del viento, S es el área frontal de la carga
y Cx su coeficiente de resistencia.
Este modelo es lo suficientemente preciso para analizar el movimiento de la carga
N
suspendida, en tanto el cable sea lo “suficientemente” rı́gido (se verá luego como evaluar
esta condición).
IO
El modelo 1.2 es no lineal (ver 1.4.1) e invariante en el tiempo, pero incluye una variable
exógena w, por lo cual el sistema no es autónomo. Tratemos de aclarar estos conceptos.
C
1.3.1. Variables de Entrada e Invarianza Temporal del Modelo
En la ecuación 1.2 la variable es θ. Obtener una solución para dicha ecuación implica
C
encontrar una expresión que describa la evolución de θ a medida que transcurre el
tiempo. Si la solución es analı́tica se tendrá una función θ(t). Si la solución es numérica
se tendrá una sucesión de pares {t, θ}.
RU
Podemos decir que además de θ, en la ecuación 1.2 el resto de los elementos
que la constituyen son operadores algebraicos, operadores diferenciales, funciones
trigonométricas y parámetros constantes respecto del tiempo a excepción de f , que
según la ecuación 1.3 evolucionarı́a si w varı́a en el tiempo.
ST
Que un modelo sea invariante en el tiempo implica que la solución no cambia tomando
diferentes instantes iniciales t0 .
C
Perturbaciones
La variable exógena en general se denomina perturbación, y cuantifica el efecto de
procesos externos independientes del proceso bajo análisis pero que inciden en su
evolución.
C
componente longitudinal de la velocidad al atravesar el vórtice a velocidad constante a diferentes alturas (derecha)
C
particulares se requiere de un modelo matemático para dicha perturbación. Realzamos
el término ”modelo”, dado que como todo fenómeno real nunca será posible describirla
de forma perfecta.
RU
Por otra parte, a la hora de describir la perturbación debemos discriminar entre
perturbaciones determinı́sticas y perturbaciones estocásticas.
para las que se conoce su evolución temporal, y por lo tanto pueden ser descritas de
forma analı́tica o numérica por una función p(t) bien definida; aunque es posible que
dicha función sea solo la descripción ”idealizada” de dicha evolución.
También es posible describir una perturbación mediante un modelo dinámico, cuya salida
al partir de una cierta condición inicial o al recibir otra perturbación sea la perturbación
N
desequilibrio.
Esto último en realidad no es un simple recurso matemático para describir una
perturbación, dado que realmente las perturbaciones son consecuencias de la evolución
C
Para el ejemplo de la carga suspendida por la grúa podrı́amos estudiar el efecto que
EN
Combinando este modelo de distribución espacial de velocidad con una velocidad media
del viendo U0 ; la variable r en la ecuación 1.4 se convierte en la variable tiempo a través
de la relación:
t = r/U0
N
Con esto se obtiene el siguiente modelo para la variación temporal de la perturbación:
U02
Γ
IO
2
w(t) = U0 + 1 − exp − t ,
2πU0 t 4νt + rc (0)2
C
Perturbaciones Armónicas Un tipo de perturbación determinı́stica especial lo consti-
tuyen aquellas que pueden modelarse mediante funciones periódicas. El caso más sim-
ple es la función senoidal, pero sabemos que cualquier función periódica puede descom-
C
ponerse en una superposición de funciones senoidales cuyas frecuencias son múltiplos
enteros de la frecuencia de dicha perturbación.
RU
Esta superposición está definida por la expansión en serie de Fourier para dicha función:
∞
X
p(t) = a0 + an sin nω0 t + bn cos ω0 t (1.5)
n=1
ST
o bien
∞
X
p(t) = a0 + cn sin nω0 t + φn
n=1
donde p an
a2n + b2n φn = tg−1
N
cn = ,
bn
La sucesión de coeficientes {an } y {bn } (o bien {cn } y {φn }) junto con la frecuencia
O
evolución temporal, y por lo tanto debe recurrirse a una descripción estadı́stica para su
caracterización.
Comenzaremos analizando modelos autónomos, que según hemos dicho son aquellos
en donde el tiempo no aparece de forma explı́cita. Esto significa que por el momento no
vamos a considerar perturbaciones (variables exógenas) o parámetros que varı́en con el
N
tiempo.
Uno de los elementos más significativos para caracterizar un proceso dinámico son los
IO
puntos crı́ticos o puntos estacionarios del modelo matemático (valores para los cuales
las derivadas temporales se anulan), que asociaremos a las condiciones de equilibrio del
proceso dinámico. Otro de los aspectos importantes es la estabilidad del comportamiento
en el entorno de estos puntos.
C
Es común confundir los términos equilibrio y estabilidad. Aclaramos desde ahora que
son conceptos diferentes pero relacionados entre si, por cuanto la estabilidad es una de
C
las propiedades que poseen las condiciones de equilibrio.
En otras palabras, una condición de equilibrio es un elemento del modelo, y la estabilidad
es una propiedad de dicho elemento.
RU
Equilibrio
Equilibrio es la condición en la cual nada cambia a medida que transcurre el tiempo.
En notación matemática:
ST
d(·)
=0 (1.6)
dt
sin θ = p
L
Puede notarse que para un dado valor constante de p existen infinitas soluciones para
O
Estabilidad
Las condiciones de equilibrio para el péndulo plano no restringido son infinitas, pero
EN
C
serı́a inestable. Pero, ¿como podrı́amos expresar el concepto de inestabilidad?.
C
equilibrio inestable, este no siempre evoluciona hacia dicha condición. ¿Qué es esta
“condición”?.
RU
Estado
Para el caso del péndulo plano la condición de equilibrio para el sistema no perturbado
(p = 0) implica no solo θ = ±n · π , sino que además θ̇ = 0.
Es decir, la situación en la cual se encuentra el sistema en un momento t1 dado no queda
ST
definida únicamente por la posición angular θ(t1 ), sino que también se debe especificar
la velocidad angular θ̇(t1 ).
Serı́a más realista usar una letra especı́fica para la velocidad angular, dado que se trata
de una variable diferente al ángulo3 .
N
θ̇ = ω
g
ω̇ + b ω + sin θ = 0
l
C
que se puede reescribir dejando a la izquierda solo derivadas y del lado derecho solo
relaciones algebraicas:
EN
θ̇ = ω
g
ω̇ = −b ω − sin θ
l
3 Esto puede resultar confuso dado que entre ángulo y velocidad angular existe una relación cinemática. Sin embargo es tan
cierto decir que la velocidad angular no es mas que la derivada del ángulo, como que el ángulo no es más que la integral de la
velocidad angular; y de hecho ambas variables se pueden medir fı́sicamente de forma independiente
N
x˙1
f1 (x1 , x2 , · · · , xn )
x˙2
f2 (x1 , x2 , · · · , xn )
IO
.. = ..
.
.
x˙n fn (x1 , x2 , · · · , xn )
que en forma compacta puede escribirse como una única ecuación diferencial vectorial
C
de primer orden:
ẋ = f (x) (1.8)
C
La ecuación 1.8 es la expresión general de lo que se denomina ecuación de estados
para el caso invariante en el tiempo; y en ella x ∈ Rn es el vector de estados. Para el
RU
modelo del péndulo plano y rı́gido n = 2.
En caso de aparecer el tiempo de forma explı́cita y/o incluir variables exógenas
(perturbaciones) el modelo de estados mas general serı́a:
ẋ = f (x, u, t) (1.9)
ST
Espacio de Estados
La función vectorial f (·) es la expresión matemática de un campo vectorial f : Rn → Rn
O
Sin necesidad de resolver la ecuación de estados podemos calcular esta velocidad para
cualquier punto del espacio de estados x1 x2 , dado que se trata del lado derecho de
la ecuación de estados, que es puramente algebraico. Hemos mencionado que este
término constituye un campo vectorial, y dado que define es la derivada temporal o
velocidad del estado, es denominado campo de velocidades.
N
IO
C
C
RU
Figura 1.4: campo de velocidades para el péndulo plano
ST
Se puede decir que las funciones x1 (t), x2 (t), · · · xn (t) son las proyecciones de una
trayectoria x(t) en los ejes x̂1 , x̂2 , · · · x̂n que utilizamos para definir el espacio de
estados. Es claro que la elección de los ejes no es única, y que podrı́amos usar cualquier
O
transformación que tenga inversa (un homeomorfismo) para definir un conjunto de ejes
alternativos.
C
Esto equivaldrı́a no solamente cambiar la dirección de los ejes, sino incluso podrı́a
implicar una deformación del espacio de estados; y sin embargo el nuevo modelo no
perderı́a su validez en tanto dicha transformación se aplique a ambos lados de la
ecuación 1.8.
EN
C
RU
ST
Lindelöf, que establece que para la ecuación 1.8 podemos garantizar la existencia de
una solución única para una dada condición inicial x(t0 ) = x0 si el campo f (x) verifica
la denominada condición de Lipschitz; que implica lo siguiente:
O
constante de Lipschitz.
Si para una dada condición inicial la solución es única, por cada punto del espacio de
estados solo puede pasar una única trayectoria, y en consecuencia las trayectorias no se
EN
pueden intersectar. Lo que si podrı́a ocurrir es que la trayectoria sea cerrada, o bien que
muchas trayectorias converjan a o diverjan de un único punto, a los que denominamos
puntos crı́ticos y que se corresponden a condiciones de equilibrio.
Además podemos decir que las trayectorias se recorren siempre en el mismo sentido a
medida que transcurre el tiempo.
N
function x d o t = pendulo ( t , x , g , l , b )
dx1 = x ( 2 ) ;
dx2 = −g / l ∗ sin ( x ( 1 ) ) − b ∗ x ( 2 ) ;
IO
x d o t = [ dx1 ; dx2 ] ;
end
C
clear ; clc ; c l f
g = 9.81;
C
l = 2;
b = 0.1;
RU
[ t , x ] = ode45 (@( t , y ) pendulo ( t , y , g , l , b ) , [ 0 2 0 ] , [ 0 . 2 0 ] )
subplot ( 2 1 1 ) ; p l o t ( t , x ( : , 1 ) ) ; g r i d on ;
subplot ( 2 1 2 ) ; p l o t ( t , x ( : , 2 ) ) ; g r i d on ;
ST
Para los modelos de segundo orden (llamados sistemas planares, por tener un espacio
O
Focos En un foco o punto espiral existe un cı́rculo alrededor del punto crı́tico tal que
toda trayectoria que esté dentro del circulo gira en espiral alrededor del punto crı́tico una
infinidad de veces y se aproxima a dicho punto (ver [O’N94, cap.8, p. 519]).
EN
Nodos En un nodo existe un familia infinita de trayectorias que entran a él (ver [O’N94,
cap.8, p. 519]).
Centros Un punto crı́tico es un centro si está encerrado por una familia infinita de
trayectorias cerradas y hay trayectorias cerradas de la familia arbitrariamente cercanas
al punto crı́tico, pero ninguna de ellas se aproxima al mismo (ver [O’N94, cap.8, p. 518]).
Figura 1.6: puntos crı́ticos aislados - arriba: nodo (izquierda) y foco (derecha); abajo: centro (izquierda) y punto de silla (derecha).
N
Punto de Silla Un punto crı́tico es un punto de silla si hay dos trayectorias semirectas
C
que se aproximan al punto crı́tico cuando aumenta el tiempo, y otras dos que se alejan;
mientras que el resto se aproximan en direcciones similares a las primeras y luego se
alejan siguiendo las direcciones de las segundas.
EN
N
El ejemplo clásico de esta clase de comportamiento es el denominado oscilador de
van der Pol. Se trata de un modelo de segundo orden en el cual el término de
IO
amortiguamiento es no lineal, cuya dinámica responde a la ecuación:
ẍ − µ 1 − x2 ẋ + x = 0
C
en donde µ es un parámetro escalar que gobierna la no linealidad y el amortiguamiento.
Podemos ver en la figura 1.7 que partiendo de diferentes condiciones iniciales la solución
siempre confluye a una misma trayectoria cerrada, o cliclo lı́mite.
C
En 1920 Van der Pol construyó el oscilador usando un trı́odo o tetrodo. Después de que
Reona Esaki inventó el diodo túnel en 1957, la fabricación del oscilador de Van der Pol
RU
con circuitos eléctricos se hizo más simple. La oscilación se da por un intercambio de
energı́a entre elementos de distinta naturaleza, como son de capacitor a inductor o de
inductor a capacitor. El diodo funciona como un elemento activo y ayuda a mantener la
oscilación.
d2 x 2 dx
LC 2 − αL 1 − x dt
+x=0
dt
N
Dinámicas Caóticas
El comportamiento más complejo de reconocer y analizar es el de las dinámicas
ST
∂T ∂T ∂ψ ∂θ ∂ψ ∆T ∂ψ
= − + k∇2 T +
∂t ∂z ∂x ∂x ∂z H ∂x
en donde ψ es una función de corriente, definida de forma tal que las velocidades del
fluido en los ejes x e z son u = ∂ψ/∂z y w = −∂ψ/∂x respectivamente.
N
los valores iniciales.
IO
Lorenz incluyó los términos:
X = ψ11
Y = T11
C
Z = T02
C
entre las corrientes descendentes y ascendentes, y Z a la diferencia en el perfil vertical
de temperatura respecto de una distribución lineal.
Con esto propuso el siguiente sistema simplificado de ecuaciones, conocidas como
RU
ecuaciones de Lorenz.
Ẋ = σ (Y − X)
Ẏ = −XZ + rX − Y
Ż = XY − bZ
ST
en donde además:
ν Ra 4
σ≡ , r≡ , b≡
k Rac 1 + a2
N
El punto
n crı́tico {0, 0, 0} corresponde
o anun movimiento no convectivo, y losopuntos crı́ticos
p p p p
en b (r − 1), b (r − 1), r − 1 y − b (r − 1), − b (r − 1), r − 1 corresponden
C
σ (σ + b + 3)
r=
σ−b−1
EN
En la figura 1.8 se muestra la evolución de los estados del modelo de Lorenz para tres
condiciones iniciales casi idénticas y dos muy diferentes.
C
RU
Se observa que luego de un cierto intervalo de tiempo en todas las soluciones se
mezclan, resultando imposible discriminar cuales se corresponden a las diferentes
condiciones iniciales.
Podemos decir que a partir de la condición inicial el pronóstico sobre la evolución futura a
largo plazo se hace prácticamente imposible. Sin embargo, al observar la 1.9 se observa
que la evolución está claramente estructurada.
ST
Esto se puede generalizar para las dinámicas caóticas, para las cuales estos patrones
de comportamiento se denominan atractores extraños.
Debe notarse que, aunque la evolución precisa del estado sea impredecible, esta es
N
Para cada atractor existirá una región del espacio de estados en la cual todas las
trayectorias que pasan por ella convergen a dicho atractor. Esta región se denomina
dominio de atracción del atractor.
análisis de estabilidad.
IO
En otras palabras la estabilidad en el sentido de Lyapunov implica que las trayectorias
se mantienen arbitrariamente cercanas al punto de equilibrio si la condición inicial está
suficientemente cercana al mismo.
La prueba matemática consiste en demostrar que para alguna una bola BR de radio R
C
hay alguna bola Br de radio posiblemente dependiente del anterior r(R), a partir de la
cual todas las trayectorias se mantendrán dentro de BR para todo t ≥ 0; como se ilustra
en la 1.10.
C
El comportamiento será inestable si no existe ninguna bola Br por más pequeña que
sea, a partir de la cual todas las trayectorias que pasen por ella queden dentro de BR
RU
para todo t ≥ 0; independientemente de lo pequeño que se elija R.
Pruebas de Estabilidad
Lyapunov propuso dos métodos para el análisis de estabilidad. El primero de ellos está
basado en las propiedades de la linealización (ver 1.4.1) de las ecuaciones del sistema
N
evolución temporal.
Como ejemplo analicemos el caso de un modelo conceptual compuesto por una masa
con un solo grado de libertad vinculado a un elemento elástico y otro disipativo. Su
aceleración estará dada por:
ẍ = −fd − fe
fd = ẋ |ẋ| , fe = k0 x + k1 x3
ẍ + ẋ |ẋ| + k0 x + k1 x3 = 0 (1.11)
N
Se puede verificar fácilmente que en equilibrio ẋ = ẍ = 0.
Podemos analizar su estabilidad a partir de la energı́a mecánica del sistema. En este
caso esta se compone de energı́a cinética y potencial elástica:
IO
Z x
1 1 1 1
mẋ2 + k0 x + k1 x3 dx = mẋ2 + k0 x2 + k1 x4
E(x) = (1.12)
2 0 2 2 4
C
Comparando las expresiones de energı́a y equilibrio pueden realizarse las siguientes
observaciones
C
El punto de equilibrio se corresponde a energı́a nula (x = 0, ẋ = 0)
Estabilidad asintótica implica la convergencia de la energı́a al valor nulo (sistema
en reposo).
RU
Inestabilidad se corresponde con un aumento de la energı́a total
Estas conclusiones implican una relación entre la estabilidad del sistema y una magnitud
escalar, la energı́a mecánica en este caso, y permiten inferir además propiedades de la
estabilidad a partir de su evolución temporal.
Esta evolución puede determinarse diferenciando la energı́a respecto del tiempo:
ST
∂E(x) ∂x
= mẍ + k0 x + k1 x3 ẋ = (−bẋ|ẋ|) ẋ
Ė(x) =
∂x ∂t
= −bẋ2 |ẋ| = −b|ẋ|3
N
Puede advertirse que Ė(x) es una magnitud siempre negativa, es decir, la energı́a
disminuirá monótonamente hasta anularse. Por lo tanto, el equilibrio es estable y
O
Función Positiva Definida Una función escalar V (x) es localmente positiva definida si
V (0) = 0, y en una bola BR0 alrededor del origen se verifica que x 6= 0 =⇒ V (x) > 0.
1
V (x) = M R2 x22 + M Rg (1 − cos x1 )
2
que es la energı́a mecánica del péndulo, es localmente positiva definida (en el entorno
de {x1 = ±nπ, x2 = 0}, n ∈ N).
En forma similar, V (x) puede definirse como una función negativa definida si −V (x) es
N
positiva definida.
V (x) es positiva semi-definida si V (0) = 0 y V (x) ≥ 0 cuando x 6= 0, es decir, se
permite que V (x) se anule en otros valores diferentes del origen.
IO
Al ser V (x) una función del estado x, es obviamente una función implı́cita del tiempo.
Por lo tanto su derivada temporal es:
C
∂V (x) ∂x ∂V (x)
V̇ (x) = = f (x) = Lf V (x)
∂x ∂t ∂x
C
Esta derivada representa la derivada de la función escalar V (x) a lo largo de la
trayectoria x, coincide con la derivada de Lie de V (x) en la dirección del campo f .
En particular V̇ (x) = 0 en un punto de equilibrio.
RU
Función de Lyapunov Si en una bola BR0 la función V (x) es positiva definida y tiene
derivadas continuas, y si su derivada temporal a lo largo del cualquier trayectoria del
sistema es negativa semi-definida, es decir, V̇ (x) ≤ 0, entonces V (x) es una función
de Lyapunov para ese sistema.
ST
Estabilidad Local Si en una bola BR0 existe una función escalar V (x) con primeras
derivadas continuas tal que:
mv̇ = −bv − kx , v = ẋ
o bien:
mẍ + bẋ + kx = 0
ẋ = Ax
donde:
x1 x 0 1
x= = , A=
x2 v −k/m −b/m
N
La energı́a mecánica (cinética + potencial elástica) es:
IO
1 1
E(x) = mv 2 + kx2 (1.13)
2 2
Esta función es positiva definida. Su derivada a lo largo de las trayectorias resulta:
C
Ė(x, v) = mv v̇ + kxẋ
= v (−kx − bv) + kxv
C
= −bv 2
la cual se anula en los puntos en donde v = 0 (en todo el eje x), en lugar de hacerlo solo
RU
en el origen. Por lo tanto es negativa semi-definida.
Para considerarla función de Lyapunov para este problema tendrı́amos que demostrar
Ė no se anula en el resto de la trayectoria luego de alcanzar el punto v = 0 (excepto
al llegar al equilibrio). Para este caso es sencillo demostrar que Ë 6= 0 cuando v = 0
(excepto en el origen), lo cual prueba que Ė no se mantiene nula al continuar sobre la
trayectoria.
ST
V (x) = xT Px (1.14)
N
en donde P es una matriz real simétrica arbitraria pero que cumpla la siguiente
desigualdad:
xT Px > 0 si kxk =
O
6 0
En este caso se dice que P es una matriz positiva definida (porque la forma cuadrática
construida con ella resulta positiva definida); y en tal caso se verifica que sus menores
C
T
= p1 x21 + p2 x22 + p̄ x1 x2 > 0
x Px = x1 x2
p̄ p2 x2
Para verificar que esta es efectivamente una función de Lyapunov para nuestro sistema
primero calculamos su derivada:
N
V̇ (x) = −xT Qx
IO
siendo Q una matriz positiva definida, que es lo mismo que decir que la siguiente
ecuación tiene solución:
AT P + PA = −Q (1.16)
C
Esta expresión se denomina ecuación de Lyapunov, y en general se utiliza proponiendo
una matriz Q arbitraria pero positiva definida y calculando P.
El algoritmo para obtener su solución puede obtenerse con el comando lyap().
C
En la figura 1.11 se muestra la función cuadrática para todo el espacio de estados x1 , x2
junto con una trayectoria y su proyección sobre la superficie de energı́a.
RU
En la figura 1.12 se ve a la izquierda la misma trayectoria, y a la derecha la elocución
de la energı́a mecánica y la de la forma cuadrática, en donde se verifica los resultados
obtenidos.
ST
N
O
C
Invariantes de Lasalle
El concepto de estabilidad puede extenderse a regiones que no sean puntos crı́ticos
aislados, por ejemplo a los ciclos lı́mites.
Una región G del espacio de estados es un conjunto invariante de un sistema dinámico si
toda trayectoria que comienza en G permanece en G para todo tiempo futuro. Podemos
citar como ejemplos:
C
un punto de equilibrio.
C
una trayectoria de un sistema dinámico autónomo.
un ciclo lı́mite (caso particular de lo anterior).
una región de atracción.
RU
la totalidad del espacio de estados (una obviedad).
Los teoremas de conjuntos invariantes, debidos principalmente a Joseph P. LaSalle,
permiten obtener resultados sobre estabilidad asintótica cuando la derivada de la función
de Lyapunov es solo negativa semi-definida, y además permiten extender el concepto de
ST
Modelos Lineales
N
Mientras que no se cuenta con una herramienta matemática para obtener soluciones
analı́ticas para un modelo de estados general, si la hay para el caso en el cual el campo
de velocidades es un operador lineal. Aclaremos primero el concepto de linealidad.
O
les, las operaciones algebraicas, etc. El campo f (x) en la ecuación 1.8 establece un
mapeo de Rn → Rn (mapea el espacio de estados x al espacio ẋ ).
Los operadores lineales son aquellos que cumplen las siguientes relaciones:
Ω(x + y) = Ω(x) + Ω(y) (1.17)
1.4.2. Linealización
Podemos aproximar un modelo de estados no lineal por uno lineal desarrollando el
N
campo en una serie de Taylor truncada en el término de primer orden.
Una serie de Taylor es una aproximación de una función mediante un polinomio de orden
IO
infinito. Para una función escalar f (y) asjuando la aproximación para un entorno del valor
ȳ se tiene:
n=∞
d(n)
X
(y − ȳ)n = f0 + fy1 (y − ȳ) + fy2 (y − ȳ)2 + · · ·
f (y) ≈ f (y)
C
dy (n)
n=0 ȳ
C
d2 f
df
f0 = f (ȳ) , fy1 = , fy2 = , ···
dy ȳ dy 2 ȳ
RU
Puede apreciarse que la aproximación se construye primero aportando un valor
constante f0 , luego agregando una variación lineal alrededor del punto ȳ con pendiente
fy1 , en tercer término una variación cuadrática alrededor del punto ȳ con curvatura
fy2 y ası́ sucesivamente. En el caso de funciones de dos variables en lugar de curvas
trabajaremos con planos.
ST
∂f ∂f
ẋ = f (x, u) ≈ f (x̄, ū) + (x − x̄) + (u − ū) + · · · (1.19)
∂x x̄,ū ∂u x̄,ū
O
Debe notarse que en la ecuación 1.19 las derivadas parciales del campo respecto del
estado y respecto de las entradas evaluadas en las condiciones de equilibrio resultan en
matrices constantes de n × n y n × m respectivamente:
C
a11 a12 ··· a1n
a21 a22 ··· a2n
∂f
= .
.. .. ..
∂x x̄,ū ..
EN
. . .
an1 an2 · · · ann
b11 · · · b1m
∂f b21 · · · b2m
= .
.. ..
∂u x̄,ū .. . .
bn1 · · · bnm
N
Con esto el modelo resulta:
ẋ = f (x, u) ≈ A ∆x + B ∆u
IO
Si reciclamos las letras x y u para sustituir a ∆x y ∆u respectivamente llegamos a la
notación habitual de la ecuación de estados lineal:
ẋ = Ax + Bu (1.20)
C
C
Por ejemplo, para el caso del péndulo plano (ecuación 1.7), el modelo es de segundo
orden con un vector de entrada escalar, y tiene la forma:
RU
x˙1 f1 (x2 )
=
x˙2 f2 (x1 , x2 , u)
siendo x1 = θ, x2 = ω, u = p, mientras que las componentes del campo son:
ST
f1 ( ) = x2
g
f2 ( ) = −b x2 − sin x1 + u
L
Las matrices jacobianas A y B serán de la forma:
N
∂f ∂f1 ∂f1
1
∂f ∂x1 ∂x2 ∂f ∂u
A= = B= =
∂x ∂f2 ∂f ∂u ∂f2
O
2
∂x1 ∂x2 ∂u
La linealización alrededor de una cierta condición {ū, θ̄, ω̄} es:
C
ω̄ 0 1
θ − θ̄
0
θ̇
≈ g + g + {u − ū}
ω̇ −b ω̄ − sin θ̄ + ū − cos θ̄ −b ω − ω̄ 1
L L
EN
Tomando como referencia {ū, θ̄, ω̄} la condición de equilibrio habitual ueq = 0, θeq =
0, ωeq = 0 el resultado es el siguiente:
x˙1 0 1 x 0
= g 1 + {u} (1.21)
x˙2 − −b x2 1
L
N
1.4.3. Modos Naturales de Respuesta
IO
Autovectores
Los autovalores λj ∈ C de una matriz A de n × n están asociados a autovectores
v j ∈ Cn que verifican la siguiente relación:
C
Av j = λj v j
C
cambia su dirección; solo cambia su magnitud en un factor λj . En otras palabras, los v j
son las direcciones invariantes del operador A·.
RU
Veremos que ası́ como los autovalores están asociados a la forma de la evolución
temporal de la amplitud del vector de estado x, los autovectores definen su dirección
en el espacio de estados.
Cambio de Coordenadas
La elección de las variables para el vector de estados no es única. Esta lleva implı́cita
ST
ẋ = Ax + Bu
x = Tz → z = T−1 x
Sustituyendo:
EN
Tż = [AT] z + Bu
ż = T−1 AT z + T−1 B u
ż = Āz + B̄u
sI − Ā = s T−1 T − T−1 AT
= T−1 (sI − A) T
N
= |sI − A|
IO
tanto su producto es conmutativo.
Coordenadas Modales
Podemos combinar los autovectores para formar una matriz a la que denominaremos
C
matriz modal.
Λ = [v 1 v 2 · · · v n ]
C
Una propiedad interesante de la matriz modal es que, si los autovalores son todos
distintos, permite obtener una matriz diagonal con los mismos autovalores que la matriz
original. Esto puede entenderse fácilmente si se tiene en cuenta que cualquier autovector
RU
cumple:
Av j = λj v j
en donde ambos lados de la ecuación son vectores columna. Construyendo una matriz
con estas columnas se tiene que:
ST
A [v 1 |v 2 | · · · |v n ] = v 1 |v 2 | . . |v n
..
.
λn
O
es decir:
λ1
λ2
C
AΛ = Λ
..
.
λn
Despejando:
EN
λ1
λ2
Λ−1 AΛ =
..
.
λn
N
En estas coordenadas, que llamaremos coordenadas modales, cada variable de estado
es independiente de las demás; ya que desarrollando los productos matriciales puede
verse claramente que el modelo se descompone en n ecuaciones diferenciales de primer
IO
orden desacopladas:
ż1 = λ1 z1 + B̄1 u
ż2 = λ2 z2 + B̄2 u
C
..
.
żn = λn zn + B̄n u
C
en donde B̄k son los vectores fila de la matriz B̄.
Cada coordenada modal responde a una ecuación diferencial de primer orden. Podemos
RU
verificar por sustitución que la solución del caso no perturbado (u = 0) para una
coordenada modal zj es:
zj (t) = zj (0)eλj t
Y como la relación entre las variables de estado originales y estas variables modales es
ST
x = Λz = [v 1 v 2 · · · v n ] z vemos que:
x1
z1
x2
z2
.. = [v 1 v 2 · · · v n ] .
. ..
N
xn zn
v11
v21
vn1
v12 v22 vn2
O
= .. z1 + .. z2 + · · · + .. zn
.
.
.
v1n v2n vnn
C
N
Veamos el espacio de estados para caso inestable:
IO
0 1 0,54
A= λ1 = 1,58, v 1 =
3,27 −0,5 0,84
−0,43
λ2 = −2,08, v 2 =
0,90
C
En la figura 1.13 podemos ver la proyección del estado en las coordenadas originales x1
y x2 (ángulo y velocidad angular), y en las coordenadas modales z1 y z2 definidas por
C
los autovectores v 1 y v 2 .
% c ómputo de a u t o v e c t o r e s y a u t o v a l o r e s
ST
%V : m a t r i z modal [ v1 v2 ]
%D: v e c t o r de a u t o v a l o r e s [ l 1 l 2 ]
[ V , D] = eig ( A , ’ vector ’ ) ;
N
Esto último significa que algunas coordenadas modales (o todas) podrı́an ser números
complejos; aunque ello no plantea una inconsistencia matemática, porque al aplicar
la transformación Λ sobre el vector de estado modal finalmente arroja solo funciones
C
reales.
Si por ejemplo un autovalor i es complejo: λi = p + jq ; el autovector asociado será
de la forma v i = u + jw donde u, w ∈ Rn ; y la solución para la variable modal es
zi (t) = r(t)+jg(t), con las funciones r, g : R → R. Estos estarán siempre acompañados
EN
C
denominada relación de Euler 4 :
zi (t) = r(t) + jg(t) = |zi (t)|ej∠zi (t) = |zi (t)| (cos ∠zi (t) + j sin ∠zi (t))
RU
de lo cual:
r(t) = |zi (t)| cos ∠zi (t) , g(t) = |zi (t)| sin ∠zi (t)
en donde sin es un vector de funciones senoidales con los desfasajes dados por el
vector φ que resulta de los arcos tangente de los cocientes entre las componentes del
N
vector u y del w.
Por lo tanto, los módulos de cada componente del autovector indica el peso de la
variable modal en cada variable de estado, y sus argumentos dan el desfasaje entre cada
O
variable de estado y la variable modal; mientras que sin ∠zi (t) le da un comportamiento
oscilatorio con la misma frecuencia a todos los estados y |zi (t)| una misma modulación
de amplitud.
C
Volviendo al ejemplo del péndulo (subamortiguado) pero ahora para el caso estable,
tanto los autovalores como los autovectores son complejos conjugados. Por ejemplo:
EN
0 1 −0,07 − j0,48
A= λ1 = −0,25 + j1,79, v 1 =
−3,27 −0,5 0,86
−0,07 + j0,48
λ2 = −0,25 − j1,79, v 2 =
0,86
4 La relación de Euler establece que para un número real θ : ejθ = cos θ + j sin θ
N
IO
C
C
RU
Figura 1.14: Trayectorias para el caso con autovalores complejos conjugados.
. .. .. ..
. ..
. . . .
.
0 0 0 ··· λn
O
en donde todas las variables y coeficientes son números reales; y en los bloques p
e q son las partes reales e imaginaria de los autovalores complejos conjugados que
C
Podemos trabajar de forma aislada con uno de estos bloques, ya que están desacopla-
EN
N
donde M se construye de la siguiente forma:
M = B | AB | A2 B | · · · | An−1 B
(1.26)
IO
donde B es una matriz columna (en este caso arbitraria pero no nula); y la matriz T se
construye de la forma:
an−1 an−2 ··· a1 1
C
an−2
an−3 ··· 1 0
W = ... .. .. .. ..
(1.27)
. . . .
C
a1 1 ··· 0 0
1 0 ··· 0 0
RU
siendo ak los coeficientes del polinomio caracterı́stico de la matriz A:
Usando este método, una posible transformación para llevar el bloque de segundo orden
a la forma estándar serı́a:
ST
−p 1
Tc = (1.29)
−q 0
Esto se da por ejemplo en dinámica estructural, en donde los efectos disipativos suelen
ser despreciables frente a los elásticos.
También encontramos este escenario cuando se busca determinar un conjunto de
C
mẍ + kx = 0
Como esta igualdad se debe cumplir para cualquier valor de t, la única solución posible
N
es que el paréntesis sea nulo. Por lo tanto:
r
k
IO
ωn =
m
C
1.4.5. Solución General del Modelo Lineal
Matrix Exponencial
C
Siguiendo el razonamiento planteado en [Oga93, sección 4-9: Solución de la Ecuación
de Estado Invariante en el Tiempo, p. 326]; vemos que un camino para resolver una
ecuación diferencial es proponer la forma de la solución y ajustar sus parámetros para
RU
que se verifique dicha ecuación. En el caso de ecuaciones diferenciales lineales suele
utilizarse la función exponencial, dado que esta resulta invariante en relación al operador
derivada; es decir, la derivada de una exponencial sigue siendo una exponencial (que es
una autofunción para los operadores diferenciales).
Si tomamos como ejemplo una ecuación diferencial de primer orden:
ST
ẋ = a x (1.30)
b c ect = a b ect → b c = a b → c = a
x(0) = b e0 = x0 → b = x0
C
en una suma de otras funciones con diferentes ”pesos”, lo cual es equivalente a proyectar
un vector ∈ Rn en un sistema de n coordenadas (en este caso el vector a proyectar
serı́a la solución que estamos buscando, y los ejes de la base donde proyectamos el
vector corresponderı́an a los términos de la serie). De hecho, desde el punto de vista
matemático, el conjunto de funciones continuas f (x) : R → R definidas en un intervalo
[a, b] conforman un espacio vectorial.
x(t) = b0 + b1 t + b2 t2 + · · · + bk tk + · · · (1.32)
N
de lo cual se deduce que:
b1 = ab0
IO
1 1 2
b2 = ab1 = a b0
2 2
1 1 3
b3 = ab2 = a b0
3 6
C
..
.
1
C
bk = ak b0
k!
y teniendo en cuenta la condición inicial queda claro que x(0) = x0 = b0 .
RU
De las soluciones dadas por las ecuaciones 1.31 y 1.32 se deduce la siguiente identidad:
∞
1 1 X 1
eat = 1 + at + a2 t2 + · · · + a3 t3 + · · · = 1 + ak tk (1.33)
2 3 k!
k=1
ST
ẋ = A x (1.34)
x(t) = b0 + b1 t + b2 t2 + · · · + bk tk + · · ·
O
b1 = Ab0
C
1 1 2
b2 = Ab1 = A b0
2 2
1 1 3
b3 = Ab2 = A b0
3 6
EN
..
.
1
bk = Ak b0
k!
Teniendo en cuenta la condición inicial:
x(0) = b0
N
Esto provee una solución general para cualquier modelo lineal, pero no da una solución
exacta. Podemos obtener una expresión exacta utilizando la transformada de Laplace.
IO
Solución por Laplace
Toda ecuación diferencial lineal ordinaria puede resolverse usando la transformada de
Laplace. Se trata de un operador matemático (lineal) que convierte una función de
C
variable real y en otra de variable compleja s (ver [Oga93, sección 1-3: La Transformada
de Laplace, p. 14];).
Z ∞
C
L {h(y)} = H(s) = h(y) e−sy dy (1.37)
−∞
Z
1
f (t) = L −1 {F (s)} = lı́m est F (s) ds
2πj T →∞ γ−jT
donde γ , la “abcisa de convergencia”, es una constante real arbitraria pero mayor que
las partes reales de las singularidades de F (s).
La utilidad de esta transformada radica en su linealidad:
N
y en la siguiente propiedad:
C
df (t)
L = s F (s) − f (0) (1.38)
dt
Otra de las propiedades de interés es la transformada de la integral de convolución:
EN
Z t
L −1 {F (s) · G(s)} = f (t − τ ) g(τ ) dτ
0
N
con esto se llega a que:
IO
X(s) = Φ(s) [x(0) + BU (s)] (1.40)
C
dominio del tiempo:
Z t
C
x(t) = Φ(t)x(0) + Φ(t − τ )Bu(τ )dτ (1.41)
0
Análisis de la Solución
ST
−1
Φ(s) = [sI − A] = =
. .. .. ..
|sI − A| |sI − A| .. . . .
adjn1 adjn2 ··· adjnn
O
φ11 (s) φ12 (s) · · · φ1n (s)
φ21 (s) φ22 (s) · · · φ2n (s)
= .
. .. . . ..
C
. . . .
φn1 (s) φn2 (s) · · · φnn (s)
X1 (s)
φ11 (s) φ12 (s) · · · φ1n (s)
x1 (0)
b11 b12 ··· b1m U1 (s)
X2 (s) 21 (s) φ22 (s) · · ·
φ φ2n (s) ···
x2 (0)
b21 b22 b2m U2 (s)
= . + .
.. .. .. .. .. .. .. .. ..
.. ..
.
. . .
.
. . . .
Xn (s) φn1 (s) φn2 (s) · · · φnn (s) xn (0) nn1 nn2 ··· bnm Um (s)
Xj (s) = φj1 (s) [x1 (0) + b11 U1 (s) + b12 U2 (s) + · · · + b1m Um (s)]
+ φj2 (s) [x2 (0) + b21 U1 (s) + b22 U2 (s) + · · · + b2m Um (s)]
+ ···
+ φjn (s) [xn (0) + bn1 U1 (s) + bn2 U2 (s) + · · · + bnm Um (s)]
N
asociados al valor inicial de cada variable de estado y a cada una de las funciones
de perturbación.
IO
Por otra parte, vemos que cada elemento φij (s) de la matriz Φ(s) es una función
racional, cuyo polinomio denominador coincide con el polinomio caracterı́stico de la
matriz A (de orden n, que resulta del determinante |sI − A|), y cuyo polinomio
numerador será a lo sumo 5 de orden (n − 1):
C
adjj,i [sI − A] b1 sn−1 + · · · + bn−1 s + bn
Φ(i,j) (s) = = n
|sI − A| s + a1 sn−1 + · · · + an−1 s + an
C
Las anti-transformadas de funciones racionales propias 6 de primero, segundo y tercer
orden se encuentran tabuladas en la bibliografı́a. Para denominadores de orden superior
RU
es sencillo obtener la anti-transformada desarrollando la función racional en fracciones
parciales 7 (ver [Oga93, sec. 1-4: Transformación Inversa de Laplace, p. 37]); dado que
la anti-transformada también es un operador lineal:
num(s) num(s)
Φ(i,j) (s) = =
sn + a1 sn−1 + · · · + an−1 s + an (s + p1 )(s + p2 ) · · · (s + pn )
ST
r1 r2 rn
= + + ··· + , rj , pj ∈ C (1.43)
(s + p1 ) (s + p2 ) (s + pn )
En la expresión precedente los números complejos rj son los residuos de la fracción
correspondiente; y los pj son los negativos de las raı́ces del denominador , que resultan
N
r r 1
φ(t) = L −1
= k e−t/T , k= , T =
s+p p p
Es evidente que la forma de la evolución temporal solo depende del autovalor −p (o de
EN
de dimensión (n − 1) × (n − 1)
6 una función racional es estrictamente propia si el orden del polinomio numerador es menor que el del denominador, bi-propia
C
Si p > 0 la función converge a cero y corresponde a una respuesta exponencialmente
estable; mientras que si p < 0 el autovalor se encuentra en el semiplano derecho del
RU
plano complejo C, y en ese caso esta función diverge (caso inestable); como se muestra
en la figura 1.15.
t φ(t)
0 1
T 0,368
2T 0,135
3T 0,050
N
4T 0,018
Se observa que cuando t = 4T la amplitud es menor al 2 % del valor inicial, y
decimos que el sistema esta nuevamente en equilibrio; aunque desde el punto de vista
O
En el caso inestable podemos calcular el tiempo Td que toma la respuesta para duplicar
C
la amplitud:
t + Td t
= ln 2 + → Td = ln 2 · T = 0,693 T
T T
En una dinámica de primer orden (n = 1) habrá una única variable de estados (x ∈ R) y
la matriz A será en realidad un escalar. La matriz B será un vector fila con m columnas:
x˙1 = [a] x1 + b1 b2 ··· bm u
k ce ṁ
M Ṫ = 0 = Tcombeq − Teq − ṁT0 → Teq = Tcombeq − T0
ce k
Sustituyendo en el modelo x = T − Teq y u = Tcomb − Tcombeq :
N
d(x + Teq ) k
M = u + Tcombeq − x − Teq − ṁT0
dt ce
IO
Como Teq es una constante:
1 k k
ẋ = (u − x) + Tcombeq − Teq − ṁT0
M ce ce
C
La suma de los últimos dos términos del corchete es nula, por tratarse de valores de
equilibrio, por lo que solo queda:
C
k
ẋ = (u − x)
ce M
RU
Llamando a = −k/ce M y b = ṁ/M llegamos a:
ẋ = ax + bu
a
Φ(s) = → Φ(t) = eat
s−a
Debe notarse que la condición inicial x(0) en este modelo corresponde a la diferencia
entre la temperatura inicial T (0) y la temperatura de equilibrio considerada Teq .
C
Si la matriz A tiene todos sus coeficientes reales, sus autovalores serán números reales
(recordar que R ⊂ C) o formarán pares complejos conjugados; y en tal caso los residuos
asociados a las fracciones de los pares complejos conjugados serán a su vez complejos
conjugados. Por lo tanto las fracciones asociadas a polos complejos conjugados se
pk = p̄h pk = a + jb , ph = a − jb
rk = r̄h rk = c + jd , rh = c − jd
rk rh rk (s + ph ) + rh (s + pk )
+ =
(s + pk ) (s + ph ) (s + pk )(s + ph )
(c + jd)(s + a − jb) + (c − jd)(s + a + jb)
N
=
(s + a + jb)(s + a − jb)
cs + ca − jcb + jds + jda − (j)2 db + cs + ca + jcb − jds − jda − (j)2 db
=
IO
s2 + 2as + jbs − jbs − (jb)2 + a2
cs + ca − db
=2 2
s + 2as + a2 + b2
Por conveniencia usaremos para los términos de segundo orden la siguiente notación:
C
r r̄ c1 s + ωn2
+ = c0 2
(s + p) (s + p̄) s + 2ξωn s + ωn2
C
que en el caso en el cual c1 = 0 y c0 = 1 se reduce a la forma estándar normalizada:
RU
ωn2
Φ(s) = (1.44)
s2 + 2ξωn s + ωn2
ωn p
φ(t) = p · e−ξωn t · sin 1 − ξ 2 ωn t (1.45)
ST
1 − ξ2
p
donde ωd = 2 ω es la frecuencia amortiguada que define el perı́odo de la
1 − ξp n
oscilación, k = ωn / 1 − ξ 2 es un factor de amplificación (ganancia) constante y
O
c1 s + ωn2 ωn2
s
EN
φ(s) = 2 = c1 +
s + 2ξωn s + ωn2 s2 + 2ξωn s + ωn2 s2 + 2ξωn s + ωn2
p
donde β = tan−1 1 − ξ 2 /ξ . Combinando ambos resultados se tiene finalmente que:
RU
ωn h p p i
φ(t) = p · e−ξωn t sin 1 − ξ 2 ωn t − c1 sin 1 − ξ 2 ωn t − β
1 − ξ2
que equivale a la expresión (1.45) con el agregado de un ángulo de fase y un ajuste del
término multiplicativo. Por otra parte, si c0 6= 1 solo se introduce un factor de amplificación
ST
Dado que para una respuesta subamortiguada ωn > 0, solo el signo de ξ definirá el
signo de la constante de tiempo; y con ello si la amplitud crece en el tiempo (en el caso
inestable ξ < 0) o decrece y la oscilación converge a cero (caso exponencialmente
N
el comportamiento es el de un “centro”.
p
λ1,2 = −ξωn ± j ωn 1 − ξ2 (1.46)
Vemos entonces que la parte real de los polos es la que define la constante de tiempo,
EN
Re(λ)
∠λ = cos−1 = cos−1 ξ
|λ|
Vemos entonces que:
todos los polos que se encuentren sobre un cı́rculo centrado en el origen tendrán
el mismo ωn
todos los polos que se encuentren sobre una misma linea que pase por el origen
N
(radial) tendrán el mismo factor de amortiguamiento ξ
todos los polos que se encuentren sobre una misma lı́nea horizontal tendrán el
IO
mismo perı́odo o frecuencia de oscilación.
todos los polos (reales o complejos conjugados) que se encuentren sobre una
misma lı́nea vertical tendrán el mismo tiempo de establecimiento, dado que este
depende solo de la constante de tiempo, y esta queda definida por la parte real
C
todos los polos (reales o complejos conjugados) con parte real negativa
corresponderán a comportamientos exponencialmente estables y viceversa.
C
los polos complejos ubicados en el eje imaginario (imaginarios puros) correspon-
derán a oscilaciones no-amortiguadas
RU
Dinámica de Orden Superior
Dado que cualquier función racional se puede descomponer en fracciones parciales,
y teniendo en cuenta que la anti-transformada de Laplace es un operador lineal;
podemos concluir que independientemente de la complejidad del modelo; sus soluciones
siempre podrán ser descompuestas en una suma de exponenciales y/o senoidales con
ST
N
IO
C
C
RU
ST
N
O
C
EN
IO
C
Respuesta a Perturbaciones
C
RU
ST
N
O
C
EN
49
Funciones de Transferencia
2.1.1. Entradas y Salida
Previamente se han introducido los conceptos de variables de estado y variables de
entradas o perturbaciones. Puede notarse que mientras que las variables de estado son
definiciones absolutamente abstractas, independientemente de que en muchos casos se
elijan con algún sentido fı́sico; las variables de entrada pueden asociarse en general a
efectos observables.
N
Respecto del estado del proceso, en general lo observable son ciertas variables medidas,
a las cuales denominaremos variables de salida.
IO
Las variables de salida a considerar estarán relacionadas de forma algebraica con los
estados (si fuesen relaciones diferenciales, serı́an necesarias mas variable de estado en
el modelo). En general las variables de salida estarán definidas por un mapeo no lineal.
Para el caso general el modelo completo serı́a:
C
ẋ = f (x, u)
y = h (x, u)
C
Al igual que la ecuación de estados, esta expresión se puede linealizar. El modelo
completo linealizado resulta:
RU
ẋ = Ax + Bu (2.1)
y = Cx + Du
Tanto en MatlabTM como en GNU Octave con las librerı́as de control instaladas
ST
g = 9.81; % a c e l e r a c i ó n g r a v i t a t o r i a
l = 2; % longitud
O
b = 0.1; % d i s i p a c i ó n v i s c o s a
A = [ 0 1 ; −g / l b ] ; % m a t r i z din ámica
B = [0 ; 1 ] ; % m a t r i z de e n t r a d a
C
nombres = { ’ \ t h e t a ’ , ’ \omega ’ } ;
EN
Sustituyendo:
N
−1 −1
Y (s) = {C [sI − A] B + D}U (s) + C [sI − A] x(0) (2.2)
IO
el efecto de las entradas sobre las salidas partiendo de condiciones iniciales nulas
(equilibrio):
−1
G(s) = C [sI − A] B+D (2.3)
C
En sı́ntesis podemos decir que:
C
A diferencia de la matriz de transición de estados, la matriz de transferencia tendrá tantas
filas como elementos tenga el vector de salidas y (filas de la matriz C) y tantas columnas
como elementos tenga el vector de entradas u (columnas de la matriz D).
RU
φ (s) · · · φ1n (s) b11 ··· b1m
c1n 11
c11 ··· d ··· b1m
.. .. φ21 (s) · · · φ2n (s) b21 ··· b2m 11
. ..
.. ..
G(s) = . .. + ..
. . . .. .. .. .. . .
.. . . . . .
cl1 ··· cln dl1 ··· dlm
φn1(s) · · · φnn (s) bn1 ··· bnm
ST
G11 (s) · · · G1m (s)
.. .. ..
= .
. .
Gl1(s) · · · Gln (s)
N
Es evidente que por superposición se puede agregar el efecto de las otras entradas
sobre una misma salida:
O
U1 (s)
Y1 (s)
G11 (s) G12 (s) · · · G1m (s)
U2 (s)
.. .. .. .. ..
. = .
. . . .
..
C
Yl (s) Gl1(s) Gl2(s) · · · Gln (s)
Um (s)
Yi (s) = Gi1 (s)U1 (s) + Gi2 (s)U2 (s) + · · · + Gim (s)Um (s)
EN
Para ese ejemplo la matriz de transferencia tiene dos filas y una columna.
N
Uj (s) s + a1 sn−1 + · · · + an−1 s + an
IO
también podemos notar que esta función queda completamente definida por las raı́ces
de estos polinomios y un factor multiplicativo:
N (s) (s + z1 ) (s + z2 ) · · · (s + zm )
G(s) = =k
D(s) (s + p1 ) (s + p2 ) · · · (s + pn )
C
donde k es un valor constante. Llamamos ceros de la función de transferencia a las
raı́ces del numerador −zk , ya que evaluar la función con uno de estos valores arroja un
C
resultado nulo; y polos a las raı́ces del denominador −pk , que son sus valores singulares.
Se verá más adelante que resulta mucho más significativo analizar una función de
transferencia a través de sus los polos y ceros que considerar los coeficientes de los
RU
polinomios correspondientes.
N
Respuesta al Impulso
Dado que la transformada de Laplace de una entrada impulsiva es una constante, puede
notarse que la función de transferencia coincide con la transformada de Laplace de la
IO
respuesta al impulso, lo cual da otra forma de definir la función de transferencia:
∞ si t = 0
si u(t) = δ(t) , δ(t) = , → U (s) = 1
6 0
0 si t =
C
Y (s) = G(s) → y(t) = L −1 {G(s)} = g(t)
C
Con modelos LTI la respuesta al impulse se puede obtener fácilmente con el comando
impulse():
RU
impulse (G ) ;
u = 0.1 ∗ t . ˆ 2 ;
l s i m (G, u , t ) ;
Por otra parte, teniendo en cuenta que la respuesta al impulso de un sistema dinámico
N
es la antitransformada de su función de transferencia; estas respuestas puede usarse
también como modelos de entrada para otros proceso dinámicos. Por ejemplo, una
senoidal puede pensarse como la respuesta al impulso de un modelo de segundo orden
IO
no amortiguado.
De hecho podemos afirmar que fı́sicamente las perturbaciones o entradas al proceso en
estudio son en realidad salidas de otros sistemas dinámicos. Por lo tanto podemos usar
expresiones como la (2.4) para generar familias de modelos de perturbación.
C
Escalón Unitario
Una forma muy generalizada de evaluar las caracterı́sticas de un proceso dinámico es la
C
de computar valores que parametricen la respuesta a una perturbación de tipo escalón.
1 si t = 0 1
RU
si u(t) = 1(t) , 1(t) = , → U (s) =
0 si t 6= 0 s
Matemáticamente el escalón unitario es la integral del impulso unitario, y como tal es una
idealización ya que fı́sicamente es imposible lograr una derivada infinita.
Pero esta función describe de forma razonable perturbaciones fáciles de generar
ST
1
Y (s) = G(s)U (s) = G(s)
s
la cual convergerá al valor de estado estacionario:
1 bn
yss = lı́m s G(s) = G(0) =
s→0 s an
N
donde bn y an son los términos independientes del numerador y denominador de la
función de transferencia. A dicho cociente lo denominamos ganancia a frecuencia cero,
o simplemente ganancia.
IO
Podemos afirmar que una función de transferencia queda completamente definida
mediante sus polos, ceros y ganancia a frecuencia cero (o el valor k ) utilizado
previamente, que es directamente proporcional.
Esta ganancia implica solo un factor de escala constante para la amplitud, y no tiene
C
ninguna incidencia en la forma de la evolución temporal de la salida.
Con modelos LTI la respuesta al escalón se puede obtener fácilmente con el comando
C
step():
s t e p (G ) ;
RU
Para una secuencia de pulsos podrı́amos recurrir al comando lsim(). Por ejemplo:
clear ; clc ; c l f
G = t f (1 , [1 0.2 1 ] ) ;
t = 0:0.01:40;
ST
u = sin ( 0 . 5 ∗ t ) ;
u = ( u > 0) − ( u < 0 ) ;
% con e s t o e l r e s u l t a d o es 1 para l o s s e m i c i c l o s p o s i t i v o s
% y −1 para l o s n e g a t i v o s
l s i m (G, u , t ) ;
N
r1 r2 rn
Y (s) = G(s)U (s) = + + ··· + U (s)
s + p1 s + p2 s + pn
EN
Cada fracción (modo natural) tiene una ganancia a frecuencia cero dada por rk /pk .
Normalmente los residuos rk tienen un mismo orden de magnitud; excepto cuando un
cero zi se encuentra próximo a un polo pj , y en ese caso el residuo es tanto menor
cuando la relación zi /pj se aproxima a 1. En el caso lı́mite los polinomios numerador
y denominador no son coprimos, por lo cual se pueden cancelar estas raı́ces, lo que
equivale a rj = 0 para el polo cancelado pj .
10 10/9 10/9 1
= − ≈ (2.6)
(s + 1)(s + 10) s + 1 s + 10 s+1
ST
En la figura 2.1 se muestra la respuesta al escalón para cada fracción (F1 y F2), la
función de transferencia completa (G) y la aproximación considerando solo la dinámica
dominante (G1).
num = 1 0 ;
O
den = [ 1 11 1 0 ] ;
% r : residuos
C
% p : polos
% k : termino constante
[ r , p , k ] = residue ( num , den )
EN
N
IO
C
C
RU
Figura 2.2: Parámetros de respuesta transitoria al escalón
N
y(t + T2 ) = ea(t+T2 ) = 2eat
IO
ln ea(t+T2 ) = ln 2eat → a(t + T2 ) = ln 2 + at
1
T2 = ln 2 (2.8)
C
a
en donde a = 1/T es la parte real del polo o par complejo conjugado inestable.
C
Respuesta al Escalón de Primer Orden
Para un sistema de primer orden el único parámetro es la constante de tiempo T , que
RU
es la inversa de la distancia a la cual se encuentra el polo respecto del origen:
a 1
G(s) = = (2.9)
s+a Ts + 1
La respuesta al escalón resulta:
ST
En este caso solo tiene sentido hablar del tiempo de establecimiento, que para un
margen del 2 % alrededor del valor final resulta ser aproximadamente 4T , ya que
y(4T ) = 1 − e−4 = 0,9817.
N
ωn2
G(s) = (2.11)
s2 + 2ξωn s + ωn2
EN
π−β
tc = , β = cos−1 ξ
ωd
π
tp =
ST
ωd
ξ
−√ π
Mp = e 1−ξ2
1
ts ≈ 4T , T =
ξωn
N
obvia.
Sin embargo, si y1 e y2 son las amplitudes de dos picos consecutivos y τ es el perı́odo,
se encuentra que:
C
donde los senos se cancelan porque estamos considerando picos (sin(·) = 1).
Llamamos a esta relación δ decremento logarı́tmico. Sustituyendo:
2πξ
δ=p
1 − ξ2
de lo cual se puede determinar el factor de amortiguamiento. Para valores pequeños
δ ∼ 2πξ .
N
Esto puede interpretarse diciendo que en la relación dinámica entre la entrada Uj (s) y
la salida Yi (s) no se observa la dinámica de ciertos modos naturales.
Por ejemplo, en una viga empotrada para pequeñas amplitudes (condición de validez del
IO
modelo lineal) los modos de torsión no se ven en la relación entre tensión y deformación
a flexión y viceversa.
C
de transferencia, en donde el residuo para la fracción del factor común entre numerador
y denominador se anula.
C
Por ejemplo:
que es lo mismo que se obtuvo en la ecuación (2.6) del ejemplo precedente. Esto es
obvio, dado que si eliminamos los términos comunes obtenemos la misma función de
transferencia.
ST
Por ejemplo:
Cuando existen ceros pero no se observa una cancelación definida, podemos plantear
el siguiente enfoque. Asumamos una función de transferencia con denominador A(s) y
numerador B(s) = as + b con un cero en z = −b/a:
EN
Y (s) B(s) as + b b b
= G(s) = = = + Ts = Ḡ(s) + T s Ḡ(s)
U (s) A(s) A(s) A(s) A(s)
C
La respuesta y(t) a una entrada u(t) arbitraria será:
C
b b dȳ(t)
y(t) = L {G(s)U (s)} = L U (s) + T s U (s) = ȳ(t) + T
A(s) A(s) dt
RU
Entonces, la respuesta con el cero se puede pensar como la superposición de la
respuesta sin el cero, con la derivada de esta multiplicada por la constante de tiempo.
Podemos hacer una gráfica de polos y ceros de un modelo LTI con el comando pzmap().
La grilla de radiales correspondientes a factores de amortiguamiento constantes y semi-
circunsferencias correspondientes a frecuencias naturales se construye con el comando
O
sgrid().
Por otra parte es posible hacer un listado de los valores de frecuencia y factor de
C
N
modales zi ):
ż1 λ1 0 ··· 0 z1 ··· b̄1j ··· .
IO
..
ż2 0 λ2 ··· 0 z · · · b̄2j · · ·
2
= . .. + .. u
.. .. .. .. .. .. j
..
.
. . . . .
. . .
..
żn 0 0 ··· λn zn ··· b̄nj ···
C
Podrı́amos decir que la columna bj en B = [b1 b2 · · · bm ] mapea variable de entrada
uj en la derivada del vector de estados ẋ.
En algunos casos ciertas variables de entrada no afectan ciertos modos naturales,
C
lo cual se refleja en coeficientes relativamente chicos en la columna correspondiente.
Por lo tanto en las funciones de transferencia asociadas a dicha entrada los polos
RU
correspondientes a los modos naturales inmunes a dicha perturbación estén cancelados
con ceros. Esta cancelación será parcial cuando exista solo alguna incidencia débil entre
una entrada y un modo natural.
Conceptualmente nos estamos refiriendo a modos naturales no controlables con la
entrada considerada.
ST
yj = cj1 cj2 ··· cjn . = Cj x
..
xn
O
plano.
Cuando el estado se mueve sobre cualquier plano paralelo a C la salida se mantiene
constante, y por lo tanto de su observación podrı́a inferirse erróneamente que el proceso
está en equilibrio. También resulta obvio que no será posible identificar
EN
Cualquier modo natural que esté contenido en un plano paralelo a C no será visible
(“observable”) a través de esta salida en particular.
En ese caso, cualquier función de transferencia para esta salida tendrá ceros que
cancelarán los polos correspondientes a los modos naturales restringidos a estos planos.
N
Las entradas en este modelo son el comando de elevador η y el empuje τ . Al analizar
IO
este modelo tı́picamente encontramos dos modos naturales subamortiguados que se
identifican como perı́odo corto y perı́odo largo o fugoide.
Por ejemplo, para una cierta aeronave de transporte en aproximación final obtenemos:
C
0,032 ± j 0,018
0,036
0,999 0,999
λc = −1,019 ± j 1,503 , vc = , |vc | =
−0,003 ± j 0,014 0,014
C
0,007 ± j 0,003 0,008
0,998
0,998
−0,060 ± j 0,005 0,060
RU
λf = −0,008 ± j 0,112 , vf = , |vf | =
0,001
0,001
−0,001 ± j 0,011 0,012
u
α
y= 0 1 0 0
q
θ
N
ż1
0 1 0 0
z1
0,332 −0,037
ż2 −3,297 −2,038 0 0 z2 + −1,957 0,296 η
=
C
ż3
0 0 0 1
z3 1 −0,071 τ
ż4 0 0 −0,013 −0,016 z4 −0,016 −0,112
z1
z2
y = 33,39 1 −0,49 0,0002 0,0019
EN
z3
z4
Vemos en la ecuación de estado que ambas entradas tiene impacto en ambos modos
naturales, mientras que de la ecuación de salida concluimos que el modo fugoide {z3 , z4 }
tiene muy poco peso en la salida.
C
Las funciones de transferencia resultan:
RU
α(s) 0,0141 (s − 57,21)(s2 + 0,024s + 0,017)
= 2
η(s) (s + 0,016s + 0,013)(s2 + 2,038s + 3,297)
En la figura 2.5 vemos que los polos lentos (modo fugoide) están prácticamente
cancelados por los ceros. Por lo tanto no se observará el modo fugoide observando
el indicador de ángulo de ataque.
N
En sı́ntesis concluimos que las cancelaciones entre ceros y polos en las funciones de
transferencia están relacionadas con la “controlabilidad” y “observabilidad” del sistema.
O
esta tiene una “singularidad”, es decir, su norma es infinito: kG(s)k = ∞. Estos puntos
coinciden con las autovalores de la matriz dinámica A, y por lo tanto con los polos de las
funciones de transferencia que la componen (mı́nimo común denominador entre todas
ellas).
EN
N
está dado por los autovalores asociados.
IO
Para un cuerpo rı́gido la aceleración angular es:
ω̇ = J−1 τ
C
inversa coinciden.
Si el vector de entrada τ está alineado con algunos de los autovectores del tensor de
inercia (ejes principales), la aceleración angular tendrá la misma dirección que el vector
C
momento, y el módulo de la aceleración tendrá una amplificación o atenuación dada por
el autovalor correspondiente (que será la inversa del momento principar de inercia en
ese eje).
RU
Si el vector τ no está alineado con un eje principal, la aceleración angular tendrá una
dirección diferente a la del momento.
Por otra parte, en algunas matrices hay ciertas direcciones de entrada para las cuales
el escalado genera como salida un vector nulo. La unión de todas estas direcciones
ST
N
Por otra parte la forma general de la función de transferencia para un modelo lineal es:
B(s) bn sn + bn−1 sn−1 + · · · + b1 s + b0
G(s) = = n
IO
A(s) s + an−1 sn−1 + · · · + a1 s + a0
B(s)
=
(s + p1 ) (s + p2 ) · · · (s + pn )
C
donde A(s) y B(s) representan los polinomios denominador y numerador respectiva-
mente, y los −pk son las raı́ces del denominador (polos de la función de transferencia).
Si la transformada de Laplace de la entrada es también una función racional:
C
N (s) N (s)
X(s) = =
D(s) (s + q1 ) · · · (s + qm )
RU
la transformada de la salida también será una función racional
B(s) N (s) B(s)N (s)
Y (s) = · =
A(s) D(s) (s + p1 ) (s + p2 ) · · · (s + pn ) (s + q1 ) · · · (s + qm )
Esta puede descomponerse en fracciones parciales:
ST
r1 r2 rn h1 hm
Y (s) = + + ··· + + + ··· +
s + p1 s + p2 s + pn s + q1 s + qm
en donde los coeficientes rk son los residuos correspondientes a los polos de la función
de transferencia y los hk corresponden a los polos de la entrada. Todos estos residuos
dependen el polinomio numerador, resultante del producto B(s) · N (s).
N
y(t) = L −1 {Y (s)}
= r1 e−p1 t + r2 e−p2 t + · · · + rn e−pn t + h1 e−q1 t + · · · + hm e−qm t (2.15)
C
Si el sistema es asintóticamente estable, la parte real de todos los polos −pk será
negativa, y por lo tanto las exponenciales correspondientes en la ecuación 2.15 decaerán
en el tiempo. Por lo tanto, en régimen o estado estacionario:
yss = lı́m y(t) = h1 e−q1 t + · · · + rm e−qm t (2.16)
EN
t→∞
que es una función con una descomposición similar a la función de entrada. Puede
decirse entonces que en régimen estacionario “la salida tendrá la misma forma que la
entrada”.
Ya se vio previamente que cuando la entrada es constante (escalón), la salida es régimen
constante, pero su amplitud depende de la ganancia de la función de transferencia. A
continuación se analiza el caso de la respuesta a una entrada senoidal.
Para una función de entrada senoidal x(t) = A sin ωt con amplitud A y frecuencia
angular ω , la transformada de Laplace es:
N
Aω 2 Aω 2
L {x(t)} = X(s) = =
s2 + ω 2 (s + jω) (s − jω)
IO
La respuesta a esta senoidal para un caso general será:
N (s) Aω 2 Aω 2 N (s)
Y (s) = · 2 =
D(s) s + ω 2 (s + p1 ) (s + p2 ) · · · (s + pn ) (s + jω) (s − jω)
C
y su descomposición en fracciones parciales resulta:
r1 r2 rn rω r̄ω
C
Y (s) = + + ··· + + + (2.17)
s + p1 s + p2 s + pn s + jω s − jω
donde debe notarse que los residuos rω y r̄ω correspondientes a los polos jω y −jω
RU
formarán un par complejo conjugado.
rk = lı́m F (s) · (s − pk )
s → pk
N
Aω 2
O
con lo cual:
C
ω ω
rω = AG(jω) , r̄ω = −AG(−jω)
2j 2j
N
y teniendo en cuenta la relación de Euler:
IO
se tiene que:
C
= 2j sin(θ)
C
Por lo tanto:
A
yss = |G(jω)| · [2j sin {ωt + ∠G(jω)}]
2j
RU
Simplificando:
Esto significa que la respuesta en estado estacionario de cualquier sistema lineal ante
una entrada senoidal será otra senoidal con la misma frecuencia ω pero amplificada o
atenuada en un factor |G(jω)| y desfasada un ángulo ∠G(jω).
N
O
C
EN
Ante todo debe notarse que la respuesta en frecuencia es algo que no puede observarse
en ”tiempo real”. Para relevar la respuesta en frecuencia deberı́amos realizar una
sucesión de experimentos, eligiendo para cada uno una cierta frecuencia.
Cada ensayo consistirı́a en inyectar una entrada senoidal con la frecuencia elegida,
N
esperar un tiempo mayor al de establecimiento para alcanzar, y medir la relación de
amplitudes y desfazajes entre la respuesta y la entrada.
De todos los ensayos para los diversos valores ωk se obtendrá una tabla de la forma:
IO
frecuencia ganancia fase
ω1 ··· ···
ω2 ··· ···
.. .. ..
C
. . .
ωk ··· ···
.. .. ..
. . .
C
Podemos realizar dos gráficos cartesianos tomando como abscisas la columna de
frecuencia y como ordenadas la ganancia y la fase respectivamente.
RU
Escala logarı́tmica Existen diversas razones en diferentes contextos para utilizar
escalas logarı́tmicas; pero invariablemente estas razones están asociadas al hecho de
que un desplazamiento en una escala logarı́tmica no representa un valor absoluto sino
una relación.
ST
Por ejemplo la distancia en la gráfica entre los valores 0.8 y 6.4 es la misma que habrı́a
entre los valores 3.2 y 25.6 (un factor 8).
N
O
C
En música la nota LA es por definición una senoidal con frecuencia 440Hz . Sin embargo,
en un piano de 88 teclas hay ocho afinadas a la nota LA, pero con frecuencias
desde 27,5Hz (LA−1 ) hasta los 3520Hz (LA7 ), pasando por los valores intermedios
EN
N
Tratándose de un sistema mecánico, el modelo dinámico mas simple para esta dinámica
serı́a uno lineal de segundo orden:
IO
a0
G(s) = k
s2 + a1 s + a0
donde k es la ganancia a frecuencia cero.
C
El módulo y la fase en función de la frecuencia resultan de la forma (2.20):
1
|G(ω)| = k q
C
2 2
(1 − ω̄ 2 ) + (2ζ ω̄)
−1 ω̄
∠G(ω) = tg 2ζ
RU
1 − ω̄ 2
√ √
en donde se usó la sustitución ω/ωn = ω̄ , siendo ωn = a0 , y ζ = a1 /(2 a0 ).
De esta expresión puede observarse lo siguiente:
Decibeles Dado que la ganancia establece una relación entre amplitudes de entrada y
C
de salida, resulta común utilizar una escala logarı́tmica también para el eje de ganancia.
En el campo de la electrónica se usan logaritmos con base 10 multiplicados por un factor
20, y con lo cual se define el decibel (db). En otras palabras, para obtener una magnitud
en decibeles aplicamos el operador {db}20 · log10 ({g}), y para extraer la relación original
EN
C
decibeles ganancia
RU
−40db 1/100
−20db 1/10
−6db ≈ 0,5√
−3db ≈ 1/ 2
0db 1 √
3db ≈ 2
ST
6db ≈2
20db 10
40db 100
Puede notarse que los valores positivos en decibeles implican amplificación, mientras
N
Cuando lo que se grafican son relaciones entre energı́as de entrada y salida (en general
asociadas a valores cuadráticos de las amplitudes) el decibel se define como 10·log10 (e).
C
En muchos casos se usan los decibeles para representar cantidades absolutas. En esos
casos se calcula el decibel sobre la relación entre el valor a representar y un valor de
referencia estandarizado.
EN
Por ejemplo, en el caso del sonido el valor de referencia es la presión sonora mı́nima
detectable por el oı́do humano. El nivel máximo tolerable es 120db, que corresponde a
106 veces el mı́nimo audible.
Otro uso similar se da en la electrónica, en donde se define la unidad dbm para niveles
de tensión eléctrica tomando como referencia 1mV .
N
términos simples:
IO
Estos términos simples pueden ser constantes, funciones de primer orden (polos y ceros
reales) o de segundo orden (polos y ceros complejos conjugados):
±1
±1 1 2 2ξ
C
k , (T s + 1) , s + s+1
ωn2 ωn
Como caso particular de lo anterior tenemos los polos y ceros en el origen s±1 .
C
Por ejemplo:
RU
0,1s (s + 0,2) (s + 40)
G(s) =
(s + 2) (s + 0,02) (s2 + 0,5s + 4)
1 1 1
= 5 × s × (5s + 1) × (0,025s + 1) × × × 2
0,5s + 1 50s + 1 s 0,5
+2 s+1
4 2
ST
log |G(jω)| = log |G1 (jω)| + log |G2 (jω)| + · · · log |Gn (jω)|
G(jω) = |G(jω)|e∠G(jω)
C
entonces
1
= G(jω)−1 = |G(jω)|−1 e−∠G(jω)
G(jω)
EN
y por lo tanto
1 1
log = − log |G(jω)| , ∠ = −∠G(jω)
G(jω) G(jω)
Esto significa que el diagrama de Bode para la inversa de una determinada función de
transferencia serán iguale al original “reflejados” respecto del eje de abscisas.
1 1 − jωT 1 1
N
G(jω) = = = (1 − jT ω)
jωT + 1 1 − jωT 1 + jωT 1 + T 2 ω2
De esto surge que:
IO
1
|G(jω)| = √ , ∠G(jω) = − tan−1 T ω (2.19)
1 + T 2 ω2
C
C
RU
ST
N
Ası́ntotas Puede notarse que si ω 1/T entonces log |G(jω)| → 0; es decir, coincide
con el eje de abscisas. Por otra parte, si ω 1/T entonces log |G(jω)| → − log ω ; es
C
Por otra parte la gráfica asintótica para la fase no da una buena aproximación. Solo
podemos decir que:
Polos y Ceros en el Origen Un caso particular de dinámica de primer orden son los
polos y ceros en el origen, ya que:
1 1 T 1
N
G(s) = = lı́m = lı́m , T =
s p→0 s + p T →∞ T s + 1 p
En este caso el quiebre de la curva de ganancia se da a frecuencia cero, con lo cual la
IO
respuesta en frecuencia se convierte en la curva asintótica del caso de primer orden con
ω → ∞ y con ganancia k → ∞.
El resultado es una recta que cruza el eje de frecuencias en ω = 1.
1 1
C
G(jω) = = e−π/2 , log |G(jω)| = − log ω
jω ω
lo que equivale a un diagrama de bode cuya curva de ganancia es una recta con
C
pendiente −1 (−20db) que cruza el eje de abscisas en ω = 1; y una curva de fase
constante de −90o .
RU
ST
N
O
C
1 1
G(jω) = 2 = 2
ω 2ζ ω 2ζ
j2 2 + j +1 1− 2 +j
ωn ωn ωn ωn
ω2 ω
N
1− 2 2ζ
ωn ωn
G(jω) = 2 2 − j 2 2
ω2 ω ω2 ω
1− 2 − 2ζ 1− 2 − 2ζ
IO
ωn ωn ωn ωn
De donde surge que:
C
|G(ω)| = q (2.20)
2 2
(1 − ω̄ 2 ) + (2ζ ω̄)
ω̄
C
−1
∠G(ω) = − tan 2ζ (2.21)
1 − ω̄ 2
RU
en donde se usó la sustitución ω/ωn = ω̄ .
La diferencia entre estas curvas asintóticas y la curva exacta depende del factor de
amortiguamiento. Si ζ = 0,7 la diferencia es máxima en ω = ωn (diferencia de 3db), un
poco menor una octava por debajo ω = ωn /2 o por encima ω = 2ωn (diferencia de 1db).
N
Las diferencias son despreciables con diferencias en frecuencia de más de una década.
Por lo tanto la curva de ganancia se puede aproximar por sus ası́ntotas como en el caso
de primer orden.
O
Y al igual que en el caso de primer orden, la gráfica asintótica para la fase no da una
buena aproximación. Solo podemos decir que:
C
ω 1/T ∠G(jω) → 0o
ω = 0,1/T ∠G(jω) = −5,7o
ω 1/T ∠G(jω) → −180o
EN
Resonancias
Los sistemas con modos naturales oscilatorios poco amortiguados presentan picos en
las curvas de ganancia de su respuesta en frecuencia, como se observa en el caso de
la figura 2.10 alrededor de los 10s−1 .
La amplitud de esos picos tienen una relación inversa con el amortiguamiento de esos
C
modos. Para un modelo de segundo orden:
ωn2
G(s) =
s2 + 2ζωn s + ωn2
RU
la respuesta en frecuencia G(jω) = M ejφ tiene ganancia M y fase φ dadas por la
siguiente expresión:
ω
2ζ
1 −1 ωn
M = s , φ = − tan (2.22)
ST
ω 2 2
ω
2 ω2
1− 2 + 2ζ 1− 2
ωn ωn ωn
√
Si 0 ≤ ζ ≤ 1/ 2 = 0,7071 la respuesta en frecuencia tiene un máximo en:
p
ωr = ωn 1 − 2ζ 2
N
(2.23)
y la magnitud de este máximo o pico resonante resulta ser:
O
1
Mr = p (2.24)
2ζ 1 − 2ζ 2
El factor de amortiguamiento es igual al coseno del ángulo entre los polos y el semieje
C
real negativo del plano s: ζ = cos β . Por lo tanto puede notarse que:
1 p
Mr = , ωr = ωn cos 2β (2.25)
sin 2β
EN
Podemos plantear dos miradas sobre esto. Por un lado los modos poco amortiguados
están asociados a picos de amplificación para las señales armónicas de frecuencias
próximas a la de resonancia. Cada par complejo conjugado con amortiguamiento bajo
agregará un pico resonante a la respuesta en frecuencia.
De forma inversa, una dinámica cuya respuesta en frecuencia presenta picos resonantes
incluye modos naturales poco amortiguados. Cada pico resonante implicará un par
complejo conjugado poco amortiguado en la dinámica del sistema.
2 n=1
4 1 1
y(t) = sin(ω0 t) + sin(3ω0 t) + sin(5ω0 t) + · · · (2.27)
π 3 5
N
La función periódica y las armónicas son elementos de un espacio de funciones definidas
en el intervalo [−T /2, T /2]. El peso de cada armónica en la descomposición de Fourier,
dado por los coeficientes an y bn , se obtiene a través del producto escalar de la función
IO
con la armónica correspondiente, lo cual dará una cierta medida cuantitativa de similitud
entre la forma de onda analizada y dicha senoidal:
Z T /2
2
a0 = y(t)dt (2.29)
T
C
−T /2
Z T /2
2
an = y(t) cos (nω0 t) dt (2.30)
T −T /2
C
Z T /2
2
bn = y(t) sin (nω0 t) dt
T −T /2
RU
El coeficiente a0 es el valor medio de la señal, que puede analizarse por separado y por
lo tanto en lo sucesivo lo consideraremos nulo.
Las funciones seno y coseno pueden combinarse en una misma senoidal desfasada:
∞
ST
X
y(t) = rn sin (nω0 t + φn )
n=1
donde
p bn
rn = a2n + b2n , φn = tan−1 −
an
N
Todo esto puede expresarse de forma mas compacta utilizando números complejos.
Teniendo en cuenta la relación de Euler:
O
1 jφ 1 jφ
e − e−jφ e − e−jφ
cos φ = , sin φ = (2.32)
2 2j
y con esto se pueden combinar los senos y cosenos de una misma frecuencia en una
EN
en donde se incluyen frecuencias positivas y negativas (se suma desde n = −∞) para
cancelar los términos imaginarios (tener en cuenta que sin(−φ) = − sin φ).
1 T /2
Z
cn = y(t) e−jnωt dt (2.34)
T −T /2
N
cn = e = (an − jbn )
2 2
Identidad de Parseval
IO
Para cualquier señal podemos definir “energı́a” (en sentido matemático) en un intervalo
[a, b] como:
Z b
E= y 2 (t)dt (2.35)
C
a
Podemos calcular la “potencia media” de una señal periódica mediante su descomposi-
ción en serie de Fourier:
C
∞ ∞
" #" #
Z T /2 Z T /2
1 2 1 X
jnω0 t
X
jmω0 t
y (t)dt = cn e cm e dt
T −T /2 T −T /2
RU
n=−∞ m=−∞
∞ ∞ Z T /2
1 X X
= |cn ||cm | ej(nω0 t+φn ) ej(mω0 t+φm ) dt
T n=−∞ m=−∞ −T /2
Z T /2 ∞
1 X 2
y 2 (t)dt = |cn | (2.36)
T −T /2 n=1
Espectro de Lı́neas
N
frecuencias especı́ficas para transmitir dı́gitos. Esta modulación se conoce como DTMF
(dual-tone multi-frequency) y se utiliza con la siguiente estandarización:
770 Hz 4 5 6 B
852 Hz 7 8 9 C
941 Hz * 0 # D
Puede notarse que la mayorı́a de las frecuencias de la modulación DTMF son co-primas,
por lo cual la superposición genera una señal periódica de frecuencia 1Hz.
N
IO
C
C
En la próxima figura se observa el espectro de lı́neas de esta señal:
RU
ST
N
O
y sus energı́as.
En la banda definida entre cero y una cierta frecuencia ωmax la cantidad de lı́neas
espectrales de una señal periódica de frecuencia ω0 será N = ωmax /ω0 . A medida
que ω0 → 0 (o bien T → ∞), esta cantidad tiende a infinito y el espectro de lı́neas se
convierte en una función continua.
Podemos considerar que una señal no-periódica equivale a una señal periódica con
N
= ω lı́m = dω lı́m = dω
T 2π T →∞ T 2π T →∞
n=−∞
T 2π −∞
Sustituyendo:
IO
Z ∞ Z ∞
1
y(t) = y(t)e jωt
dt e−jωt dω (2.38)
2π −∞ −∞
C
El término entre corchetes define la Transformada de Fourier de la señal y(t):
Z ∞
Y (ω) = F {y(t)} = y(t)ejωt dt (2.39)
C
−∞
Debe remarcarse que para que una función f : R → R tenga transformada de Fourier,
RU
esta debe ser integrable en el eje real; y que su transformada es una función compleja:
F : R → C.
Esta transformada existe solo si la función y(t) se extingue para |t| → ∞.
Resulta obvio que la expresión 2.38 e la operación que permite recuperar la señal original
a partir de su transformada, y por lo tanto constituye la transformada inversa de Fourier :
ST
Z ∞
1
y(t) = F −1 {Y (ω)} = Y (ω)e−jωt dω (2.40)
2π −∞
Laplace dada por (1.37), en la cual se sustituye el exponente complejo s por uno
imaginario puro jω .
O
Z ∞
Y (f ) = F {y(t)} = y(t)ej2πf t dt (2.41)
−∞
Z ∞
y(t) = F −1 {Y (f )} = Y (f )e−j2πf t df (2.42)
EN
−∞
Puede afirmarse entonces que Y (f ) es una forma alternativa de definir la función y(t).
Mientras que y(t) es la proyección de una señal en el dominio del tiempo, Y (f ) es su
proyección en el dominio de la frecuencia.
Finalmente mencionemos que al igual que la serie de Fourier, esta transformada también
se puede aplicar a funciones vectoriales o matriciales.
N
Z ∞ Z ∞ Z ∞ Z ∞
y 2 (t)dt = Y (f ) y(t)ej2πf t df df = Y (f )Y ∗ (f )df (2.44)
−∞ −∞ −∞ −∞
IO
donde Y ∗ (f ) es el complejo conjugado de Y (f ) (notar el signo de la exponencial). Por
lo tanto:
Z ∞ Z ∞ Z ∞
1
C
2 2
y 2 (t)dt = |Y (f )| df = |Y (ω)| dω (2.45)
−∞ −∞ 2π −∞
Se deduce entonces que el área bajo la curva |Y (f )|2 es la energı́a de la señal; que
C
estará definida solo si la función se extingue cuando |t| → ∞. De no ser ası́ el resultado
serı́a ∞, pero de todas formas por lo dicho anteriormente su transformada no estarı́a
RU
definida.
Señales Aleatorias
2.5.1. Análisis Estático
ST
común en los juegos de azar. Sin embargo, para variables definidas en el dominio de los
números reales, la probabilidad de obtener un determinado valor dentro de los infinitos
O
intervalo seleccionado.
N
donde
1 x
Φ(x) = 1 + erf √ (2.49)
IO
2 2
siendo erf la función de error dada por:
Z x
2 2
erf (x) = √ e−t dt (2.50)
C
π 0
Momentos Estadı́sticos
C
valor medio Es el valor esperado de una muestra:
x̄ = E {x} (2.51)
RU
varianza Es el valor esperado del cuadrado del apartamiento de una muestra respecto
del valor medio: n o
2
var = E (x − x̄) (2.52)
ST
sesgo El sesgo es el valor esperado del apartamiento entre una muestra y su valor
medio elevado al cubo:
n o
3
sesgo = E (x − x̄) (2.53)
N
Un valor positivo indica una distribución inclinada hacia la derecha respecto de su valor
O
n o
4
curtosis = E (x − x̄) (2.54)
N
cov(x) = E x xT
(2.57)
IO
en donde los elementos de la diagonal se corresponden a las varianzas de cada
componente del vector; y los elementos restantes corresponden a covarianzas entre
las componentes del vector.
Si no hay relación estadı́stica entre las diferentes componentes del vector, los elementos
C
no diagonales serán nulos. Es fácil por lo tanto intuir el significado de los autovalores y
autovectores de dicha matriz.
C
Procesos Estacionarios y Procesos Ergódicos
Un proceso estocástico es estacionario si sus propiedades estadı́sticas no dependen del
tiempo.
RU
Un proceso estocástico es ergódico si sus propiedades estadı́sticas pueden obtenerse
analizando un conjunto de muestras de una única realización de dicho proceso.
construir un histograma.
La función de correlación entre dos variables aleatorias es el valor esperado del producto
entre el valor que toma una de ellas en un instante de tiempo t y el valor de la otra en el
instante diferente t + τ para diferentes valores del retardo τ :
O
Z T /2
1
Rxy (τ ) = lı́m x(t)y(t + τ )dt (2.59)
T →∞ T −T /2
EN
Es evidente que para la función de autocorrelación se cumple que R(τ ) = R(−τ ), y que
la autocorrelación para τ = 0 se corresponde con la varianza.
N
2π −∞ T →∞ T
IO
1 1
S(ω) = lı́m Y (ω)Y ∗ (ω) = lı́m |Y (ω)|2 (2.63)
T →∞ T T →∞ T
con lo cual
C
Z ∞
1
var = S(ω)dω (2.64)
2π −∞
C
Densidad Espectral y Autocorrelación
Existe una relación biunı́voca entre las funciones de densidad espectral Y (ω) y
RU
autocorrelación R(τ ) de un proceso estocástico:
Z T /2
1
R(τ ) = lı́m y(t)y(t + τ )dt
T →∞ T −T /2
Z ∞
1 T /2
Z
1
= lı́m y(t) Y (ω)ejω(t+τ ) dω dt
ST
T →∞ T −T /2 2π −∞
Z ∞ Z ∞
1 1
= lı́m y(t)e dt Y (ω)ejωτ dω
jωt
2π −∞ T →∞ T −∞
Z ∞
1 1 ∗
= lı́m Y (ω)Y (ω)ejωτ dω
N
2π −∞ T →∞ T
Z ∞
R(τ ) = S(ω)ejωτ dω = F −1 {Y (ω)} (2.65)
−∞
C
∞
EN
Z
1
S(ω) = R(τ )e−jωτ dω = F {R(τ )} (2.66)
2π −∞
2
Sy (ω) = |G(jω)| Sx (ω) (2.67)
N
Procesamiento Digital de Señales
IO
Desde el punto de vista práctico, las señales se procesan mediante computadores
digitales operando sobre ”muestreos” de las mismas.
Para tomar o registrar datos de sensores analógicos con un sistema digital se
C
utilizan conversores analógico/digitales (ADC: analog digital converter ). Cuando es
necesario generar señales continuas desde la computadora se utilizan conversores
digitales/analógico (DAC: digital analog converter ). Un ejemplo común es el hardware
C
de sonido presente en computadoras de escritorio y teléfonos celulares.
Un ADC es básicamente un circuito de muestreo y retención (sample & hold) seguido
de un codificador, como el que se esquematiza en la figura 2.12. Una llave (transistor
RU
JFET o MOSFET) conecta un capacitor con la señal analógica a muestrear durante un
breve intervalo de tiempo. En general se incluyen “buffers” a la entrada y a la salida del
circuito S/H para compensar las pérdidas. Un vez cargado el capacitor se abre la llave y
comienza el proceso de digitalización.
2.6.1. Muestreo
El muestreo de una señal es una secuencia de valores que se corresponde con los de
la señal en ciertos instantes de tiempo. Si el tiempo entre muestras ts es uniforme:
N
y ∗ = {y0 , y1 , y2 , · · · yN } , y ∗ (n) = y(n · ts) (2.68)
IO
en donde n es un tiempo discreto adimensional y N es el tamaño del muestreo.
EL muestreo se puede concebir como un proceso en el cual una señal continua es
“convolucionada” con un tren de impulsos unitarios para dar la secuencia de muestras
(ver figura 2.13). También se podrı́a decir que el muestreo es un tren de impulsos
C
“modulado” por la señal a muestrear.
Un tren de impulsos unitarios separados entre sı́ a intervalos ts se puede expresar
matemáticamente como:
C
∞
X
s(t) = δ(t − nts ) (2.69)
n=0
RU
La señal muestreada y ∗ en el dominio del tiempo continuo resulta:
∞
X ∞
X
∗
y (t) = y(t)δ(t − nts ) = y(nts )δ(t − nts ) (2.70)
n=0 n=0
ST
Queda claro que todas las senoidales cuyas frecuencias den el mismo resto tendrán el
N
mismo muestreo. La condición para poder reconstruir de forma exacta una senoidal a
partir de su muestreo es que se cumpla la condición:
O
ωs
ωnq = >ω (2.73)
2
Si no se respeta esta relación, no será posible discriminar entre una senoidal de baja
C
frecuencia respecto de sus ”alias” con frecuencias superiores en una cantidad entera
de la frecuencia de Nyquist, como se muestra en la figura 2.14a. Este problema es
denominado comúnmente como aliasing.
EN
IO
con frecuencias mayores a la de muestreo (habitualmente 24Hz). Es notable al ver
hélices y rotores de helicópteros, o ruedas con rayos por ejemplo.
La aparición del aliasing tampoco está asociado únicamente al muestreo temporal. Por
C
ejemplo, en fotografı́a digital se realiza un muestreo “espacial” de una escena. En este
caso el intervalo de muestreo está determinado por el tamaño del pixel. En muchos casos
se observan distorsiones como en la figura 2.14b, donde se muestra una secuencia de
C
anillos con centro en el vértice superior izquierdo que por efecto de muestreo lleva a
nuestro cerebro a reconstruir anillos inexistentes.
RU
Lı́mite Superior para la Frecuencia de Muestreo
Desde el punto de vista teórico no existe lı́mite superior para la frecuencia de muestreo.
Pero en el análisis teórico no se considera la “cuantización” de las muestras implı́cita en
todo proceso de digitalización, en el cual es necesario convertir los valores muestreados,
que matemáticamente se modelan como números reales (y por lo tanto con resolución
ST
Lo mismo ocurre aunque sin digitalizar el muestreo cuando la observación se realiza con
sistemas sensoriales de resolución acotada.
Experimentamos esto a diario cuando observamos procesos lentos de nuestro entorno
O
fı́sico. Por ejemplo, nos resulta imposible detectar los cambios producidos por las mareas
en una zona portuaria mirando de forma continua la superficie del agua.
Tampoco vemos los procesos de desarrollo en una planta y mucho menos el desgaste
C
En sı́ntesis podemos decir que la frecuencia de muestreo tiene que ser lo suficiente-
mente alta como para poder reconstruir correctamente el proceso observado, pero lo
N −1
N
n
X
Yk = yn e−j2π N k (2.74)
n=0
IO
La transformada inversa queda definida por la siguiente expresión:
N −1
1 X k
yn = Yk ej2π N n (2.75)
N
C
k=0
C
cómputo con una carga de orden N log N denominado Transformada Rápida de Fourier,
o FFT por sus siglas en inglés (Fast Fourier Transform). Este algorı́tmo requiere que el
RU
muestreo tenga una dimensión equivalente a N = 2k , k ∈ N (una potencia entera de 2).
N −1 N −1
1 X
ST
X
|xn |2 = |Xk |2 (2.76)
n=0
N
k=0
1 X 1 X
var = |xn |2 = 2 |Xk |2 (2.77)
N n=0 N
k=0
O
Para expresar la densidad espectral en unidades de ingenierı́a hay que aplicar los
factores de conversión para llevar las frecuencias adimensionales a valores reales:
C
1 fs
Φ(fk ) = |Xk |2 , fk = k (2.78)
N · fs N
Este procedimiento se puede utilizar en GNU Octave y Matlab con el comando pwelch:
clear ; c l f ; clc
N
ts = 0.01; % tiempo de muestreo
t = 0: ts :3600; % v e c t o r de tiempo
u = randn ( s i z e ( t ) ) ; % v e c t o r de e n t r a d a
IO
% co n t r u i m o s un f i l t r o de ejemplo
[B, A] = butter ( 3 , 0 . 2 ) ;
F = t f (B,A, ts ) ;
% simulamos l a r es p u e s t a d e l f i l t r o a l r u i d o de e n t r a d a
y = lsim (F, u , t ) ;
C
% hay 360.000 muestras
nw = 2 ˆ 9 ; % tama ño d e l bloque para p a r t i c i o n a r e l muestreo
f s = 1 / t s ; % f r e c u e n c i a de muestreo
C
% calculamos l a s d i s t r i b u c i o n e s e s p e c t r a l e s con e l p r o c e d i m i e n t o de Welch
% dejamos e l v a l o r por omisi ón para e l t i p o ( ’ psd ’ )
Su = pwelch ( u , nw , nw / 2 , nw , f s ) ;
Sy = pwelch ( y , nw , nw / 2 , nw , f s ) ;
RU
Df = f s / nw ; % i n t e r v a l o en e l e j e de f r e c u e n c i a
f = ( 1 : length ( Su ) ) ∗ Df ; % v a l o r e s para e l e j e de f r e c u e n c i a
% dibujamos l a s d i s t r i b u c i o n e s e s p e c t r a l e s
subplot ( 2 1 2 ) ;
semilogx ( f , Su ) ; hold on ;
N
semilogx ( f , Sy ) ; g r i d on ;
y l a b e l ( ’PSD ’ ) ;
x l a b e l ( ’ f r e c u e n c i a [ hz ] ’ ) ;
O
legend ( ’ e n t r a d a ’ , ’ s a l i d a ’ , ’ L o c a t i o n ’ , ’ southwest ’ ) ;
IO
C
Modelado de Sistemas Mecánicos
C
RU
ST
N
O
C
EN
93
Modelos de Parámetros Concentrados
Cuando estudiamos el movimiento de un sólido cuya distribución de masas no cambia
utilizamos modelos de “cuerpo rı́gido”.
Un cuerpo rı́gido se caracterizan completamente con sus propiedades másicas: masa
y mementos de inercia. El resto de sus caracterı́sticas fı́sicas son irrelevantes para
el análisis dinámico, y para este estudio resulta conveniente utilizar las ecuaciones
cardinales de la dinámica.
Si además el cuerpo no rota (o no nos interesa estudiar ese aspecto), usamos un modelo
N
de “masa puntual” considerando únicamente su masa concentrada en su baricentro. En
este caso podemos realizar el estudio de forma sencilla a través de la ecuaciones de
Newton.
IO
La situación se torna un poco más compleja cuando se trata de sistemas de cuerpos
rı́gidos vinculados entre sı́, como veremos a continuación.
C
Frecuentemente encontramos problemas en los cuales elementos “rı́gidos” se conectan
entre sı́ a través de vı́nculos elásticos y/o amortiguadores. Los elementos rı́gidos
C
concentran las propiedades másicas, mientras que los vı́nculos proveen rigidez y
disipación.
También es frecuente encontrar modelos de este tipo como aproximación de estructuras
RU
con parámetros distribuidos, como veremos más adelante (ver [Mei97, capı́tulo 8:
Approximate Methods for Distributed-Parameter Systems, p. 501])
El estado de los elementos rı́gidos queda definido por sus desplazamientos r (incluyendo
giros) y las velocidades ṙ correspondientes.
Sumando a las propiedades másicas las correspondientes a los elementos disipativos
ST
y elásticos de los vı́nculos que los conectan obtenemos una definición de las
caracterı́sticas del sistema mecánico.
Un ejemplo clásico es el modelo de parámetros concentrados para la suspensión de un
automóvil, el cual se ilustra en las figuras 3.1 y 3.2.
N
Systems, p. 553]).
Formulación de Newton-Euler
C
j k
IO
C
C
RU
Figura 3.2: Modelo conceptual para la suspensión de un automóvil
1 = zc − zf − l1 θ
2 = zc − zr + l2 θ
3 = zf − hf
4 = zr − hr
N
Ternas no inerciales Debe tenerse presente que los principios de conservación de
IO
cantidad de movimiento en general se formulan tomando como referencia un sistema de
coordenadas inerciales.
Sin embargo, en el caso de cuerpos rı́gidos es más práctico considerar un sistema de
coordenadas solidario al cuerpo, y por lo tanto no inercial.
C
En ese caso a las derivadas de cantidades vectoriales definidas en esta terna del cuerpo
hay que agregar los cambios debidos al movimiento de dicha terna:
dai dab
C
= + ω b × ab
dt dt
En el caso de cantidades compuestas del tipo a = b + c:
RU
dai dbb dcb
= + + ω b × b b + cb
dt dt dt
C
das generalizadas (ver fig.3.4), pero en el caso general serán funciones de las coorde-
nadas generalizadas rj = Rj (q1 , q2 , · · · , qN ).
RU
En otras palabras, las coordenadas generalizadas no necesitan ser desplazamientos
geométricos.
ri = cij qj (3.4)
j=1
N
X N
X N
X
mij q̈i + kij (qj − qi ) + bij (q̇j − q̇i ) = fi (3.6)
j=1 j=1 j=1
C
Formulación de Lagrange
Es posible obtener las ecuaciones dinámicas para un sistema mecánico utilizando la
EN
formulación de Lagrange. Esta surge del principio de mı́nima acción; según el cual las
trayectorias de un sistema mecánico son aquellas de mı́nima energı́a.
Para obtener las ecuaciones del movimiento se construye la función de Lagrange
L = T − V , donde T es la energı́a cinética total del sistema y V la potencial, equivalente
al trabajo realizado por las fuerzas conservativas (las que surgen de un campo potencial
como el gravitatorio, el electromagnético o la deformación elástica).
N
coordenada generalizada qi :
δ W̄ = Qi · δqi
IO
en donde los δqi son desplazamientos virtuales, es decir, desplazamientos infinitesima-
les compatibles con las condiciones de borde. Las contribuciones serán positivas para
aquellas fuerzas en la misma dirección que el desplazamiento y viceversa.
C
En el caso general el lagrangiano será una función L(qi , q̇i , t), porque el modelo podrı́a
incluir parámetros que dependan explı́citamente del tiempo (masa variable por ejemplo).
C
Como ejemplo vamos a considerar el caso de un péndulo doble.
Podemos describir la configuración geométrica con solo dos variables θ1 y θ2 , por lo que
pueden usarse estas como como coordenadas generalizadas.
RU
ST
N
O
Asumimos que las barras son rı́gidas y consideramos masas “equivalentes” localizadas
en los extremos de las mismas.
C
2 2
v22 = (l1 ω1 cos θ1 + l2 ω2 cos θ2 ) + (l1 ω1 sin θ1 + l2 ω2 sin θ2 )
= l12 ω12 cos2 θ1 + l22 ω22 cos2 θ2 + 2l1 l2 ω1 ω2 cos θ1 cos θ2
N
+ l12 ω12 sin2 θ1 + l22 ω22 sin2 θ2 + 2l1 l2 ω1 ω2 sin θ1 sin θ2
v22 = l12 ω12 + l22 ω22 + 2l1 l2 ω1 ω2 cos (θ1 − θ2 )
IO
Para las variaciones de altura tenemos que:
δh1 = l1 (1 − cos θ1 )
δh2 = δh1 + l2 (1 − cos θ2 ) = l1 (1 − cos θ1 ) + l2 (1 − cos θ2 )
C
Sustituyendo:
C
1
m1 v12 + m2 v22 − g (m1 δh1 + m2 δh2 )
L=
2
1 1
RU
= m1 l12 ω12 + m2 l12 ω12 + l22 ω22 + 2l1 l2 ω1 ω2 cos (θ1 − θ2 )
2 2
− m1 g l1 (1 − cos θ1 ) − m2 g [l1 (1 − cos θ1 ) + l2 (1 − cos θ2 )]
1 1
L= (m1 + m2 ) l12 ω12 + m2 l22 ω22 + m2 l1 l2 ω1 ω2 cos (θ1 − θ2 )
2 2
− (m1 + m2 ) g l1 (1 − cos θ1 ) − m2 g l2 (1 − cos θ2 )
ST
− =0
dt ∂ q̇1 ∂q1
d ∂L ∂L
− =0
O
dt ∂ q̇2 ∂q2
∂L
= (m1 + m2 ) l12 ω1 + m2 l1 l2 ω2 cos (θ1 − θ2 )
∂ q̇1
∂L
= −m2 l1 l2 ω1 ω2 sin (θ1 − θ2 ) − (m1 + m2) g l1 sin θ1
EN
∂q1
∂L
= m2 l22 ω2 + m2 l1 l2 ω1 cos (θ1 − θ2 )
∂ q̇2
∂L
= m2 l1 l2 ω1 ω2 sin (θ1 − θ2 ) − m2 g l2 sin θ2
∂q2
N
m1 + m2 l1
ω̇1 + [ω̇2 cos (θ1 − θ2 ) − ω2 (ω1 − ω2 ) sin (θ1 − θ2 )]
m2 l2
IO
m1 + m2 g
−ω1 ω2 sin (θ1 − θ2 ) − sin θ1 = 0
m2 l2
l2
ω̇2 + [ω̇1 cos (θ1 − θ2 ) − ω1 (ω1 − ω2 ) sin (θ1 − θ2 )]
l1
C
g
−ω1 ω2 sin (θ1 − θ2 ) + sin θ2 = 0
l1
C
m 1 + m 2 l1 l1
Obtenemos expresiones más concisas con las sustituciones = λ, = ¯l y
m2 l2 l2
g
= ḡ :
RU
l1
λ ω̇1 + [ω̇2 cos (θ1 − θ2 ) − ω2 (ω1 − ω2 ) sin (θ1 − θ2 )] − ω1 ω2 sin (θ1 − θ2 ) − λḡ sin θ1 = 0
¯l ω̇2 + [ω̇1 cos (θ1 − θ2 ) − ω1 (ω1 − ω2 ) sin (θ1 − θ2 )] − ω1 ω2 sin (θ1 − θ2 ) + ḡ sin θ2 = 0
ST
2 i
donde kij es rigidez del elemento elástico que vincula el desplazamiento ri con el rj
(que será nulo si tal vı́nculo no existe).
∂L ∂L X
= mi q̇i , = kij (qi − qj )
∂ q̇i ∂qi j
N
Reemplazando en la ecuación 3.7:
X
IO
mi q̈i + kij (qi − qj ) = fi (3.9)
j
C
En el caso más general, cuando los elemento de masa mi tiene una posición ri que no
coinciden con alguna coordenada generalizada qj , dicha posición puede expresarse en
términos de estas mediante la ecuación (3.4). La velocidad de este elemento será:
C
N
X ∂ri
vi = q̇j
∂q
RU
i=1 j
2 2 i=1
∂qi j=1
∂qj
k=1 k=1
N N N
!
1 XX X ∂rk ∂rk
= mk q̇i q̇j
2 i=1 j=1 ∂qi ∂qj
k=1
N
X
mij = mi
i=1
∂qi ∂qj
N N
1 XX
T = mij q̇i q̇j
2 i=1 j=1
EN
1 T
T = q̇ Mq̇ (3.10)
2
En notación matricial:
1 T
V = q Kq (3.11)
2
N
La matriz de rigidez K también es son positivas definida, y se demostrará que también
es simétrica.
IO
Entre las fuerzas no conservativas en los sistemas mecánicos es habitual encontrar
disipación.
Si esta disipación es “viscosa”, estas fuerzas disipativas son proporcionales a la
velocidad pero de dirección opuesta. En esos casos para derivar las ecuaciones del
C
movimiento a partir de la formulación de Lagrange es conveniente introducir un término
R conocido como función de disipación de Rayleigh [Rao11, p. 610], definida como:
C
1 T
R= q̇ Cq̇ (3.12)
2
RU
donde C la “matriz de amortiguamiento”, positiva definita como las matrices de masa y
rigidez.
Con esto la ecuación de Lagrange se reescribe de la siguiente forma:
d ∂T ∂T ∂R ∂V
− + + = Qi (3.13)
dt ∂ q̇i ∂qi ∂ q̇i ∂qi
ST
q2 a21 a12 ··· a2N p2
= .
.. .. .. .. ..
. ..
. . . .
qN aN 1 aN 2 ··· aN N pN
o bien:
q =F·p
N
1/2ajj p2j , pero esto también producirı́a un desplazamiento en el nodo i igual a aij pj ,
produciendo un trabajo adicional aij pj pi . El trabajo total serı́a:
IO
1 1
W = aii p2i + ajj p2j + aij pj pi
2 2
Si invertimos el orden en el que aplicamos las cargas tendremos
C
1 1
W = aii p2i + ajj p2j + aji pi pj
2 2
C
Pero como el orden no altera el trabajo neto realizado, se deduce que aij = aji .
Esto es útil para realizar una primer verificación al calcular la matriz de rigidez.
RU
Formas Prácticas
Una forma alternativa para obtener la matriz de rigidez es la de considerar las fuerzas
ejercidas por cada elemento de vı́nculo en las ecuaciones de conservación, y luego
relacionar estas con las coordenadas generalizadas estableciendo desplazamientos
unitarios en cada una por separado. Ası́ primero planteamos:
ST
Mq̈ = F1 f
f = F2 q
mz̈c = f1 + f2
EN
mf z̈f = f3 − f1
mf z̈r = f4 − f2
J θ̈ = f1 l1 − f2 l2
N
Entonces en este caso
1 1 0 0
IO
−1 0 1 0
F1 =
0 −1 0 1
l1 −l2 0 0
C
Luego encontramos la matriz de rigidez entre estas fuerzas y las coordenadas
generalizadas, considerando un desplazamiento unitario en cada una de ellas por
separado.
C
Si por ejemplo consideramos un desplazamiento en zc , aparecerá una fuerza f1 = −k1 zc
y una f2 = −k2 zc , con lo cual completamos la primer columna.
Realizando lo mismo con los otros desplazamientos y con las velocidades tenemos que:
RU
f1
−k1 k1 0 −k1 l1 zc
f2 −k 0 k k l zf
2 2 2 2
=
0 −k3
f3
0 0
zr
f4 0 0 −k4 0 θ
ST
−b1 b1 0 −b1 l1 żc
−b2 0 b 2 b 2 l 2
żf
+
0 −b3
= F2 q + F3 q̇
0 0 żr
0 0 −b4 0 θ̇
N
Mq̈ + Bq̇ + Kq = 0
Por otra parte, de las ecuaciones 3.10 y 3.11 vemos que la energı́a cinética y potencial
elástica surgen de formas cuadráticas construidas con las matrices de inercia y rigidez
respectivamente. Por lo tanto, expresando estas energı́as en función de las coordenadas
EN
Veamos un ejemplo:
N
2V = k1 (x1 − x2 )2 + k2 (x2 − x3 )2
= k1 x21 + k1 x22 − 2k1 x1 x2 + k2 x22 + k2 x23 − 2k2 x2 x3
IO
k1 −k1 0 x 1
= x1 x2 x3 −k1 k1 + k2 −k2 x2
0 −k2 k2 x3
C
Para el armado de la forma matricial tenemos en cuenta que en una forma cuadrática el
coeficiente de x2k ocupa el lugar k en la diagonal principal, mientras que el coeficiente de
xi xj corresponde a los elementos ij y ji (la matriz es simétrica).
C
Como 2V = xT K x y 2T = ẋT M ẋ se deduce que
m1 0 0
RU
M = 0 m2 0
0 0 m3
k1 −k1 0
K = −k1 k1 + k2 −k2
0 −k2 k2
ST
M q̈ + B q̇ + K q = F u (3.14)
C
o como:
[A M] q̈ + [A B] q̇ + q = [A F] u (3.16)
N
q̇ 0 I q 0
ẋ = = + u (3.18)
q̈ −K̄ −B̄ q̇ F̄
IO
donde
K̄ = M−1 K , B̄ = M−1 B , F̄ = M−1 F (3.19)
C
En estructuras elásticas normalmente los amortiguamientos son muy bajos (1 % o
menos). El caso no amortiguado surge de despreciar estos efectos, es decir en la
C
ecuación 3.14 tomar B = 0.
q̈ + K̄ q = F̄ u (3.20)
RU
La ecuación de estados queda:
q̇ 0 I q 0
= + u (3.21)
q̈ −K̄ 0 q̇ F̄
Frecuencias Naturales
La ecuación caracterı́stica para la matriz dinámica en ese caso resulta:
N
I 0 0 I sI −I
0 I − −K̄
s = =0
0 K̄ sI
A B
M=
C D
C
y si D es invertible:
|M| = |D||A − BD−1 C|
Tendiendo en cuenta que si α ∈ R, A ∈ Rn×n se verifica que |αA| = αn |A|, y teniendo
en cuenta que en este caso:
A = sI , B = −I , C = K̄ , D = sI
= |s2 I + K̄| = 0
donde en este caso N son los grados de libertad.
Si sustituimos s2 = −λ:
N
|λI − K̄| = 0
√
En sı́ntesis se tendrán autovalores imaginarios puros s = ±j λi .
IO
Entonces las respuestas asociadas a los modos naturales para sistemas el
Çsticos no
amortiguados serán movimientos armónicos simples con frecuencias ωi = λi , siendo
λi los autovalores de la matriz M−1 K.
C
Formas Modales
Es evidente que para sistemas no-amortiguados no es necesario trabajar con el modelo
de estados, dado que toda la información dinámica puede deducirse a partir de la matriz
C
K̄ = M−1 K.
es decir
ST
q = v 1 η1 + v 2 η2 + · · · + v n ηn
Esto equivale a proyectar el vector de coordenadas generalizadas q en un sistema
de coordenadas definidos por los autovectores v i , los cuales se corresponden con
la configuración de desplazamientos compatibles con el movimiento armónico de
N
frecuencia ωi .
Decimos que el autovector v i da la forma modal para el movimiento armónico de
frecuencia ωi .
O
Podemos ver que a diferencia de un caso dinámico más general, para sistemas
mecánicos no amortiguados la forma modal tiene una representación geométrica real.
q(t) = v 1 0 + v 2 0 + · · · + v i ηi (t) + · · · + v n 0
N
evolucionará en el tiempo según la siguiente expresión:
IO
mientras que las velocidades serán:
C
Puede notarse que cuando la amplitud del desplazamiento alcanza un valor máximo,
la velocidad es nula y viceversa. Con amplitud de desplazamiento máxima la energı́a
potencial elástica alcanza su valor máximo Vimax al tiempo que la energı́a cinética resulta
C
nula. Esta a su vez tendrá una máximo Timax cuando la elástica sea nula, y como se trata
de sistemas conservativos estos valores máximos serán iguales.
RU
Timax = Vimax (3.23)
Debe notarse que esto se corresponde con las energı́as cinética y potencial asociadas
a un determinado modo natural i.
ST
l1 = 1;
N
l2 = 0.5;
g = 10/ l 1 ;
l = l1 / l2 ;
O
M = [m 1 ; 1 l ] ; % m a t r i z de i n e r c i a
K = [m∗g 0 ; 0 g ] ; % m a t r i z de ” r i g i d e z ”
C
% e l de l a p r i m e r b a r r a i g u a l a 0.15 rad
L = repmat ( [ l 1 ; l 2 ] , 1 , 2 ) ;
X = [ 0 0 ; sin (Q) . ∗ L ] ; % coordenadas X de l o s nodos
X(3 ,:) = X(2 ,:) + X(3 ,:)
Y = [ 0 0 ; cos (Q) . ∗ L ] ; % coordenadas Y de l o s nodos
Y(3 ,:) = Y(2 ,:) + Y(3 ,:)
l = l1+l2 ; % escala v e r t i c a l
clf ;
Podemos describir la configuración geométrica con solo dos variables θ1 y θ2 , por lo que
pueden usarse estas como como coordenadas generalizadas.
N
IO
C
C
RU
Cociente de Rayleigh
ST
−1
Los autovalores λi y autovectores v i de K̄ = M K por definición verifican la
siguiente relación:
M−1 K v i = λi v i
Kv i = λi Mv i (3.24)
En esta expresión ambos lados de la igualdad son vectores de orden n, por lo cual
O
v Ti Kv i = λi v Ti Mv i
C
v Ti Kv i
λi = ωi2 = (3.25)
v Ti Mv i
EN
Por lo tanto:
Vimax
ωi2 = (3.27)
Ti∗
Si ahora elegimos un vector v arbitrario en lugar de una forma modal, obtendremos una
cantidad escalar caracterı́stica del sistema que solo depende de esta elección:
N
v T Kv
λ = ω 2 = R(v) = (3.28)
v T Mv
IO
Esta cantidad R(v) se denomina cociente de Rayleigh, y coincide con el cuadrado de
una frecuencia natural en caso que el vector v elegido corresponda a una forma modal.
C
En el ejemplo precedente se consideró un péndulo doble con:
6 1 60 0
M= , K=
C
1 2 0 10
Continuando con el código del ejemplo tomamos el primer autovector (primera columna
de la matriz modal) y computamos el cociente de Rayleigh:
ST
v1 = V ( : , 1 ) ;
R1 = v1 ’ ∗ K∗ v1 / ( v1 ’ ∗M∗ v1 ) ;
w1 = s q r t ( R1 ) / ( 2 ∗ p i ) ;
Se constata que el resultado coincide con el primer elemento del vector de frecuencias
naturales:
N
>> w1 − W( 1 )
ans =
O
>>
C
siempre podemos expresar un vector arbitrario v como una combinación lineal de los
autovectores del sistema (ya que como hemos dicho, estos forman una “base”):
n
X
v= cr v r = Λc
r=1
1 Un punto estacionario de una función f (x) es un valor del argumento x en donde se anula la derivada, lo que implica que
ΛT MΛ = I , ΛT KΛ = λ
N
c2i
i=1
IO
ci serán pequeños si i 6= r. Entonces puede decirse que ci = i cr para i 6= r, donde
i 1. Luego, dividiendo numerador y denominador por c2r :
n
λi (1 − δir ) 2i
P
λr + n
C
i=1
X
R(v) = n ≈ λr + (λi − λr ) 2i
δir ) 2i
P
1+ (1 − i=1
i=1
C
donde δir es el delta de Kronecker (1 si i 6= r, 0 si i = r). El segundo término del
lado derecho es una cantidad de segundo orden, lo que indica que en proximidad de un
RU
autovector el cociente de Rayleigh tiene un punto estacionario.
Si se considera un vector próximo al primer modo:
n
X
R(v) ≈ λ1 + (λi − λ1 ) 2i
i=2
ST
R(v) ≥ λ1
Es implica que el cociente de Rayleigh nunca es menor que el primer autovalor, y por lo
tanto su mı́nimo se corresponde a este.
N
Esto lleva a considerar que es posible obtener una buena aproximación para la
frecuencia del modo fundamental utilizando como vector de prueba la deformada estática
del sistema bajo la acción de un campo gravitatorio; hecho que puede ser útil cuando la
O
Por otra parte, de la ecuación 3.25 puede observarse que aumentando la rigidez
C
v Tj Kv i = λi v Tj Mv i
v Ti Kv j = λj v Ti Mv j
v Tj Kv i = v Ti Kv j
v Tj Mv i = v Ti Mv j
Por lo tanto, en las dos ecuaciones anteriores el lado izquierdo es el mismo; y en el lado
derecho solo difieren en el valor de λ. Restando ambas ecuaciones:
0 = (λi − λj ) v Ti Mv j
N
Dado que hemos asumido λi 6= λj , necesariamente se cumple que:
v Ti Mv j = 0 (3.29)
IO
y por lo tanto también se cumple que
v Ti Kv j = 0 (3.30)
Podemos afirmar entonces que los autovalores o formas modales son ortogonales
C
respecto de las matrices M y K.
C
En muchos problemas de ingenierı́a resulta crı́tico establecer la frecuencia del primer
modo natural. Existen al menos dos formas de obtener aproximaciones para esta
frecuencia:
RU
Método de Rayleigh El método de Rayleigh consiste simplemente en evaluar el
cociente de Rayleigh usando como vector de prueba el correspondiente a la deformada
estática obtenida aplicando como carga el peso propio de cada elemento.
ST
Método de Dunkerley Este método permite obtener una lı́mite inferior para la
frecuencia del primer modo natural para modelos con mas de tres grados de libertad.
En lo sucesivo seguiremos el planteo de [Rao11, cap.7.2, p. 656].
Sabemos que las frecuencias naturales están dadas por las raı́ces de la siguiente
N
como los autovalores de la inversa de una matriz son la inversa de los autovalores de la
matriz original, lo anterior se puede reescribir como:
C
1
− I + AM = 0
λ
1
− + a11 m1 ···
a12 m2 a1N mN
λ
1
a21 m1 − + a22 m2 ··· a2N mN
λ =0
.. .. ..
.
. .
1
aN 1 m1 aN 2 m2 · · · − + aN N mN
λ
N
Llamando 1/λi a las raı́ces de este determinante, podrı́amos expresar la ecuación
caracterı́stica como:
IO
λ−1 − λ−1 λ−1 − λ−1 · · · λ−1 − λ−1
1 2 N =0
λ−N − λ−1 −1 −1
C
−(N −1)
1 + λ2 + · · · + λN λ − ··· = 0
Igualando el coeficiente de λ−(N −1) con la expansión del determinante resulta que:
C
λ−1 −1 −1
1 + λ2 + · · · + λN = a11 m1 + a22 m2 + · · · + aN N mN
RU
En la mayorı́a de los casos las frecuencias altas son significativamente mayores a la
del primer modo, por lo cual sus valores al cuadrado lo serán mucho más. Entonces se
puede aproximar esta expresión por:
λ−1
1 ≈ a11 m1 + a22 m2 + · · · + aN N mN
ST
1 1 1 1
≈ 2 + 2 + ··· + 2
ω12 ω11 ω12 ω1N
N
Cuerdas
N
IO
C
C
RU
ST
N
∂2w
EN
∂T ∂w ∂w ∂ 2 w
dT = dx , sin θ ≈ , sin (θ + dθ) ≈ + dx
∂x ∂x ∂x ∂x2
∂ 2 w(x, t)
∂ ∂w(x, t)
T + f (x, t) = ρ(x) (3.32)
∂x ∂x ∂t2
Oscilaciones libres Para vibraciones libres f (x, t) = 0 proponemos una solución por
separación entre una variación temporal y una espacial:
N
Sustituyendo:
∂ 2 q(t)
∂ ∂y(x)
q(t) T = ρ(x)y(x) (3.34)
IO
∂x ∂x ∂t2
Si además la tensión es constante en el tiempo, agrupando variaciones temporales a la
derecha y espaciales a la izquierda:
1 d2 q(t)
1 d dy(x)
C
T = (3.35)
ρ(x)y(x) dx dx q(t) dt2
Como ambos las dos de la ecuación dependen de variables diferentes, la igualdad
C
solo puede verificarse para cualquier par {x, t} si ambos miembros son iguales a una
constante que llamaremos λ = ω 2 . Por lo tanto tenemos para las variaciones espaciales
que:
RU
d dh(x)
T + ω 2 ρ(x)y(x) = 0 (3.36)
dx dx
mientras que para la evolución temporal:
d2 q(t)
+ ω 2 q(t) = 0 (3.37)
ST
dt2
Antes de proponer una solución veremos que lo mismo puede decirse para el caso de
deformaciones longitudinales o torsionales.
Vibraciones Longitudinales
N
O
C
EN
∂ 2 u(x, t)
(P (x, t) + dP ) − P (x, t) + f (x, t)dx = ρA(x)dx
∂t2
donde f (·) es una distribuida fuerza externa, y ρ la densidad lineal de masa.
∂P
Con la sustitución dP = dx se llega a que:
∂x
N
∂P (x, t) ∂ 2 u(x, t)
+ f (x, t) = ρA(x)
∂x ∂t2
IO
y teniendo en cuenta la ecuación 3.38:
∂ 2 u(x, t)
∂ ∂u(x, t)
EA(x) + f (x, t) = ρA(x) (3.39)
C
∂x ∂x ∂t2
C
forma:
u(x, t) = v(x)q(t) (3.40)
RU
Entonces:
d2 q(t)
d dv(x)
EA(x) q(t) = ρA(x) v(x)
dx dx dt2
de lo que surge:
1 d2 q(t)
1 d dv(x)
ST
EA(x) =
ρA(x)v(x) dx dx q(t) dt2
Dado que el lado izquierdo es solo función de x y el lado derecho solo depende de t,
ambos deben ser constantes. Por lo tanto:
N
1 d dv(x)
EA(x) = −ω 2
v(x)ρA(x) dx dx
1 d2 q(t)
O
= −ω 2
q(t) dt2
d dv(x)
EA(x) + ω 2 ρA(x)v(x) = 0 (3.41)
dx dx
EN
y
d2 q(t)
+ ω 2 q(t) = 0 (3.42)
dt2
De esta última ecuación puede verse que el movimiento libre será armónico simple:
d2 v(x)
+ α2 v(x) = 0
dx2
en donde:
ρ
α2 = ω 2 (3.43)
E
Proponemos la solución:
N
v(x) = A cos αx + B sin αx (3.44)
IO
en donde las constantes deben ajustarse en función de las condiciones de borde.
Para un una condición inicial arbitraria proponemos como solución una serie infinita de
estas posibles soluciones.
C
∞
X x x
u(x, t) = (ai cos ωi t + bi sin ωi t) · ci cos ωi + di sin ωi (3.45)
i=1min
c c
C
donde c = E/ρ es la velocidad de propagación de las ondas longitudinales en el sólido.
RU
Tomemos como ejemplo el caso de vibraciones longitudinales en el árbol de transmisión
de un motor naval.
El timón de una embarcación interfiere con el volumen de agua empujado por la hélice,
generando fluctuaciones en el empuje cada vez que una de las aspas pasa por el timón.
Por lo tanto resultarı́a de interés conocer cuáles pueden ser las frecuencias naturales del
ST
sistema de transmisión para identificar las velocidades de rotación del eje que lo puedan
inducir a la resonancia.
N
O
C
EN
N
IO
C
El problema de hallar las frecuencias y modos normales requiere identificar las
condiciones de borde adecuadas para el problema.
El vı́nculo en x = 0 restringe el desplazamiento longitudinal, por lo que cualquier forma
C
modal vi (x) debe anularse en ese punto.
Por otro lado, al perderse la condición de extremo libre por la presencia de la masa de
la hélice, el esfuerzo axil N (L) = EA(L) en ese punto debe ser igual a la variación del
RU
impulso lineal de la misma:
ST
u(0, t) = 0
2
−EA ∂u(L, t) = M d u(L, t)
∂x dt2
O
La solución está dada por la ecuación (3.45). Aplicando la primer condición de borde
sobre las funciones modales vi (x), se deduce que ci = 0.
C
EA ωL ωL
=⇒ cos = M ω sen
c c c
La expresión que obtuvimos anteriormente se corresponde con la ecuación de
frecuencias del problema.
Si llamamos G(ω) a la diferencia entre los dos miembros de la igualdad anterior, los
valores ωn que verifican G(ωn ) = 0 representan las frecuencias naturales del sistema.
N
IO
C
Con las frecuencias naturales y la ecuación (3.44) obtenemos las formas modales vi (x):
C
RU
ST
N
Para encontrar los cruces por cero de una serie de puntos {x, y} podemos usar el
siguiente código:
EN
% x : abcisas
% y : ordenadas
% z : coordenadas x de l o s cruces por cero
% k : ı́ n d i c e s de l o s puntos donde se d e t e c t a r o n cambios de s i g n o
function [ z , k ] = z e r o s c r o s s i n g s ( x , y )
c = y > 0; % puntos con y > 0
k = f i n d ( c ( 2 : end ) ˜= c ( 1 : end − 1 ) ) ; % t r a n s i c i o n e s 0/+ y +/0
N
dy = y ( k +1) − y ( k ) ;
z = x ( k ) − y ( k ) . ∗ dx . / dy ;
end
IO
end
Vibraciones Torsionales
C
C
RU
ST
∂ 2 θ(x, t)
(M (x, t) + dM ) − M (x, t) + f (x, t)dx = IO (x)dx
∂t2
O
∂θ
M = GJ (3.46)
∂x
siendo GJ la rigidez torsional por unidad de longitud. Procediendo igual a los casos
EN
∂ 2 θ(x, t)
∂ ∂θ(x, t)
GJ(x) + IO (x) = f (x, t) (3.47)
∂x ∂x ∂t2
N
IO
C
C
Figura 3.8: solicitaciones sobre un elemento diferencial de una barra a flexión
RU
En la figura 3.8 se destaca un elemento diferencial ubicado en la coordenada x de masa
dm = ρA(x)dx. Planteando suma de fuerzas de corte en el elemento diferencial:
∂ 2 w(x, t)
ρA(x)dx = − (Q(x, t) + dQ) + f (x, t)dx + Q(x, t)
∂t2
ST
Si despreciamos la energı́a cinética por rotación del elemento diferencial (válido para
vigas esbeltas) podemos asumir que hay equilibrio de momentos. Tomando momentos
respecto del extremo izquierdo como se ve en la figura 3.8:
dx
(M (x, t) + dM ) − (Q(x, t) + dQ) dx + f (x, t)dx −M =0
N
2
Y escribiendo:
O
∂Q ∂M
dQ = dx , dM = dx
∂x ∂x
C
se obtiene:
∂ 2 w(x, t) ∂Q(x, t)
ρA(x) + = f (x, t)
∂t2 ∂x
∂M (x, t) ∂Q(x, t) 1
EN
− Q(x, t) = + f (x, t) dx ≈ 0
∂x ∂x 2
∂2M ∂Q
=
∂x2 ∂x
∂ 2 w(x, t) ∂ 2 M (x, t)
ρA(x) + = f (x, t)
∂t2 ∂x2
De la teorı́a de Euler-Bernoulli la relación entre momento y deformación (o giro) está
dada por:
∂θ ∂2w
M = EI = EI (3.48)
∂x ∂x2
N
donde E es el módulo de Young e I es el momento de inercia (de área) de la sección
transversal respecto del eje y . Sustituyendo:
IO
∂2 ∂ 2 w(x, t) ∂ 2 w(x, t)
EI(x) + µ(x) = f (x, t) (3.49)
∂x2 ∂x2 ∂t2
C
Oscilaciones libres Para vibraciones libres (f (·) = 0) proponemos una solución de la
forma:
C
w(x, t) = y(x)q(t) (3.50)
Entonces:
RU
d2 d2 y(x) d2 q(t)
EI(x) q(t) + ρA(x) y(x) = 0
dx2 dx2 dt2
de lo que surge:
d2 d2 y(x) 1 d2 q(t)
1
ST
EI(x) = −
y(x)ρA(x) dx2 dx2 q(t) dt2
Dado que el lado izquierdo depende solo de x y el lado derecho solo depende de t,
ambos deben ser constantes. Por lo tanto:
N
d2 d2 y(x)
1
EI(x) = ω2
ρA(x)y(x) dx2 dx2
1 d2 q(t)
O
= −ω 2
q(t) dt2
d2 d2 y(x)
EI(x) − ω 2 ρA(x)y(x) = 0 (3.51)
dx2 dx2
EN
y
d2 q(t)
+ ω 2 q(t) = 0 (3.52)
dt2
De esta última ecuación puede verse que el movimiento libre será armónico simple:
N
y(x) = A cosh βx + B sinh βx + C cos βx + D sin βx (3.54)
en donde las constantes deben ajustarse en función de las condiciones de borde.
IO
Como ejemplo realizaremos un análisis de vibraciones transversales libres de una pala
del rotor principal del Westland Lynx (figura 3.9a).
Estas son de fibra de carbono y presentan las siguientes caracterı́sticas:
C
Parámetro Valor
Largo de pala 4,8 m
C
Cuerda 22 cm
Perfil NACA 0020
E efectivo 151 GP a
RU
Densidad 1550 kg/m3
(para este ejemplo) que la pala es homogénea. El elastómero que vincula la raı́z de
la pala al cubo del rotor se modela como un elemento elástico concentrado de rigidez
K = 63000 N m. El modelo conceptual serı́a el siguiente:
O
C
EN
(a) Westland Lynx (b) Esquema del cubo del rotor semi-rı́gido del Linx
IO
plantear condiciones de borde sobre esta función espacial.
Para nuestro modelo conceptual, en el empotramiento (x = 0) no hay desplazamiento y
el momento flector de la viga coincide con el del elastómero:
C
y(0) = 0 , (3.55)
C
y 000 (L) = 0 , y 00 (L) = 0 (3.56)
Imponiendo las condiciones para el extremo x = 0 sobre nuestra solución y(x), se tiene:
RU
A cosh 0 + B senh 0 + C cos 0 + D sen 0 = 0
EJβ 2 (A cosh 0 + B senh 0 − C cos 0 − D sen 0) =
−Kβ (A senh 0 + B cosh 0 − C sen 0 + D cos 0)
−K
A= (B + D) (3.57)
2EIβ
Por otro lado, las condiciones de borde para x = L aplicadas a la solución (3.54) generan
dos ecuaciones complementarias:
N
( 2
β (A cosh 0 + B senh 0 − C cos 0 − D sen 0) = 0
β 3 (A senh 0 + B cosh 0 + C sen 0 − D cos 0) = 0
O
(
(senh βL − aF1 ) B − (sen βL − aF1 ) D = 0
(3.58)
(cosh βL − aF2 ) B − (cos βL − aF2 ) D = 0
EN
siendo:
K
a=
2EIβ
F1 = cosh βL + cos βL
F2 = senh βL − sen βL
G(β) = (senh βL − aF1 ) (cos βL + aF2 ) − (sen βL + aF1 ) (cosh βL − aF2 ) = 0 (3.59)
N
continuación se grafica dicho término en función de la frecuencia, indicándose en
cı́rculos rojos las frecuencias naturales para la condición de empotrado libre (K → ∞)
[Se ha distorsionado la escala vertical para poder notar con mayor claridad los cruces por cero]:
IO
C
C
RU
Una vez obtenidas las frecuencias, es posible usar la relación (3.57) y las ecuaciones
ST
siendo
O
En la próxima figura se muestra en color azul los primeros cuatro modos naturales,
superpuestos con el modo correspondiente para la condición k → ∞ en rojo.
EN
2 0
Como se trata de un movimiento armónico simple v̄(x) = ω ȳ(x), siendo ȳ(x) la forma
modal. Entonces:
C
Z L
1
Tmax = ω 2 T ∗ , T∗ = µ(x)ȳ 2 (x)dx
2 0
Para el numerador debemos calcular la energı́a potencial elástica, que depende del tipo
de deformación considerada.
N
Z L
1 2
V = GJ(x) [φ0 (x)] dx (3.64)
2 0
IO
El cociente de Rayleigth, considerando que ahora la energı́a cinética depende del
momento de inercia y no de la masa, queda:
Z L 2
GJ(x) φ̄0 (x) dx
C
ω2 = 0
Z L
(3.65)
J(x)φ̄2 (x)dx
0
C
Vibraciones flexionales En este caso consideramos el trabajo realizado por el
momento elástico, para el cual la deformación correspondiente es el giro:
RU
dy
θ=
dx
La energı́a potencial elástica de un elemento diferencial es
∂θ
dV = M · dθ = EI · dθ
∂x
ST
Y como
∂θ ∂2y
= = y 00 → dθ = y 00 dx
∂x ∂x2
resulta que
1 L
Z
2
N
Z L
2
EI(x) [ȳ 00 (x)] dx
ω 2 = 0Z L (3.67)
C
ρA(x)ȳ 2 (x)dx
0
Un primer enfoque obvio para resolver modelos continuos cuando no se conoce una
solución exacta es la de distribuir su masa en elementos puntuales conectados entre
sı́ mediante vı́nculos elásticos sin masa; y analizarlos como sistemas de N grados de
libertad.
Una vez construidas las matrices de masa M y rigidez K se obtienen los modos
naturales (formas
y frecuencias) resolviendo el problema de autovalores y autovectores
para la matriz M−1 K .
Z L
1
Tr = ẏ(x, t)2 µ(x)dx
2 0
C
ẏ(x, t) = f (x)q̇(t)
Sustituyendo:
Z L Z L
1 1 2
Tr = µ(x)f 2 (x)q̇ 2 dx = q̇ µ(x)f 2 (x)dx
2 0 2 0
N
Por ejemplo, para el caso de un sistema masa-resorte fijo en un extremo en donde el
resorte tiene propiedades constantes (µ = m/L) podemos asumir una deformación
lineal en su extensión: f = x/L. Entonces:
IO
L
µL3
Z
1
Me = µ x2 dx = 2
= m
0 3L 3
C
Por lo tanto deberı́amos agregar a la masa en el extremo la tercera parte de la masa del
resorte.
C
Solución del Problema de Autovalores
Existen diferentes métodos numéricos para resolver el problema de autovalores. En
RU
[Rao11] e presentan varias alternativas: Método Matricial Iterativo (sec.7.5), Método de
Jabcobi (sec.7.6) y Problema Estándar (sec. 7.7).
Actualmente se cuenta con programas que resuelven esto rápidamente, como Matlab y
Octave mediante el comando eig().
Pero para el caso de vigas y ejes se dispone de algunas formas ventajosas cuando solo
ST
Método de Holzer
El método de Holzer permite calcular modos naturales para deformaciones unidireccio-
nales usando modelos de parámetros concentrados con acoplamiento adyacente.
N
Método de Myklestad
El método de Myklestad es una extensión del de Holzer para vibraciones flexionales.
C
El tema se puede consultar en [Tho82, sección 10.4: Método de Myklestad para Vigas,
p. 304].
Flexo-Torsión
EN
Existen muchos casos de estructuras tipo viga en donde la deformación libre es una
combinación de flexión y torsión. Esto ocurre por ejemplo en una viga recta cuando el
centro de gravedad de las secciones no coincide con sus centros de torsión (eje elástico),
o en el caso general cuando el eje elástico no es recto. Un ejemplo de ello es el caso del
ala de una aeronave.
Este problema se puede encarar como una extensión del método de Miklestad, como se
plantea en [Tho82, sección 10.6: Vibración Acoplada Flexión-Torsión, p. 311].
N
3.3.3. Discretización con Parámetros Distribuidos
IO
La alternativa al modelo de masas concentradas para modelar sistemas continuos es la
de proponer una aproximación para las funciones de deformación con suficientes grados
de libertad como para obtener la precisión requerida en la determinación de los modos
naturales de respuesta.
C
Los parámetros de ajuste de esta función terminan adoptando el rol de coordenadas
generalizadas para la aproximación discreta del modelo continuo.
C
Método de Rayleigh-Ritz
Sabemos que el cociente de Rayleigh siempre arroja un valor igual o superior a la
frecuencia del primer modo natural para cualquier deformada propuesta que cumpla las
RU
condiciones de borde.
Ritz propuso utilizar como función de deformación w(x) una suma de funciones de
interpolación wi (x) multiplicadas por coeficientes a determinar.
N
X
w(x) = c1 w1 (x) + c2 w2 (x) + · · · + cN wN (x) = ci wi (x) (3.68)
ST
i=1
La idea es determinar los coeficientes ci que minimizan el valor del cociente de Rayleigh.
Teniendo en cuenta que este cociente nunca será menor a la frecuencia del primer
modo, el ajuste resultante será el más próximo posible a la deformada real con la forma
N
propuesta.
Para realizar el proceso de minimización buscamos el conjunto de coeficientes que
verifiquen:
O
∂R(w)
=0 , i = 1, 2, ..., N (3.69)
∂ci
Esto resulta en un sistema de N ecuaciones para determinar los valores de los ci .
C
U (w)
R(w) = = U (w)T (w)−1 = ω 2 (3.70)
T (w)
Al derivar este cociente respecto de uno de los parámetros de w(x) se tendrı́a usando
la regla de la cadena:
N
Rc (ω 2 ) c = 0 (3.72)
IO
que anulan el determinante de Rc (ω 2 ), y para cada uno de ellos obtener el espacio nulo
2
de la matriz resultante.
C
Por ejemplo, para un problema de vibraciones transversales en una viga para el cociente
de Rayleigh (ver 3.3.1) se tiene:
2
C
L
d2 y(x)
Z
U (y) = EJ(x) dx
0 dx2
Z L
RU
2
T (y) = ρ(x) {y(x)} dx
0
0 i=1
Z L N
X N
X
= EJ(x) ci wi00 (x) cj wj00 (x)dx
0 i=1 j=1
N X
N
N
X Z L
= ci cj EJ(x)wi00 (x)wj00 (x)dx
i=1 j=1 0
O
N X
X N
= ci cj ki j
i=1 j=1
C
donde Z L
kij = EJ(x)wi00 (x)wj00 (x)dx
0
EN
2 el subespacio nulo (null space) o núcleo (kernel) de una matriz A ∈ Rn×n es el conjunto de vectores u ∈ Rn tales que
A · u = 0. Para una determinada matriz con determinante nulo, estos vectores pueden obtener en Matlab con el comando
null()
N
N X
X N Z L
= ci cj ρ(x)wi (x)wj (x)dx
i=1 j=1 0
IO
N X
X N
= ci cj mij
i=1 j=1
donde
C
Z L
mij = ρ(x)wi (x)wj (x)dx
0
C
Con estos resultados podemos decir que:
N N
∂U (y) X ∂T (y) X
= cj kij , = cj mij
RU
∂ci j=1
∂ci j=1
j=1 j=1
λI − M−1 K c = 0
K c = λM c →
N
−1
siendo c y λ autovalores y autovectores de la matriz M K que deben ser
determinados.
O
Los autovalores λ dan una aproximación a los valores de ω 2 para los diferentes modos,
mientras que los autovectores c combinados con la ecuación 3.68 definen autofunciones
que aproximan a las formas modales.
EN
N
en todo el dominio espacial, en FEM la estructura se divide en elementos, y se utiliza
una única función elemental o de interpolación en cada uno de ellos, lo “suficientemente
continua” en el elemento y nula en el resto del dominio.
IO
Normalmente esta función de interpolación es la misma para todos los elementos, pero
con diferentes ajustes de sus coeficientes en cada uno. Y cualquiera sea la elección,
finalmente estos coeficientes se determinan en función de las deformaciones en los
nodos del elemento.
C
Se busca que las funciones de interpolación sean del menor orden posible, pero deben
ser compatibles con las deformaciones que experimentará el elemento.
C
Funciones de Interpolación Por ejemplo, para deformaciones longitudinales estáticas
u(x) en una viga se tiene para la elástica que:
RU
d2 u(x)
EA =0
dx2
Integrando
du(x)
= c2
ST
dx
u(x) = c1 + c2 x
Podemos entonces proponer como función de interpolación una recta (figura 3.11).
N
y(x) = c1 + c2 x + c3 x2 + c4 x3 (3.73)
C
EN
Figura 3.11: Funciones de interpolación para una barra sometida a deformaciones unidimiensionales (torsión o tracción)
Resulta conveniente utilizar una coordenada adimensional dentro del elemento ξ = x/l.
Con esto la forma genérica para el desplazamiento en un elemento e serı́a:
en donde los coeficientes cek absorben el factor de escala del cambio de coordenadas.
Esto se puede expresar en forma vectorial:
N
ue (ξ) = φT (ξ) ce = cTe φ(ξ) (3.75)
IO
donde:
φ1 (ξ)
ce1
φ2 (ξ) ce2
φ(ξ) = , ce =
···
· · ·
φn (ξ) cen
C
siendo c un vector de coeficientes que dependen deben ajustarse en función de las
condiciones de borde.
C
Para las funciones de interpolación antes mencionadas las φk (ξ) serı́an monomios de
orden k .
RU
Ajuste de Coeficientes Aunque en todos los elementos se use el mismo conjunto
de funciones φ, en cada tendrá su propio vector de coeficientes c. Para el caso de
deformaciones longitudinales o de torsión, las condiciones de borde están establecidas
por los desplazamientos de los nodos:
lo cual equivale a:
wa 1 0 ce1 ce1 1 0 wa
= → =
N
wb 1 1 ce2 ce2 −1 1 wb
que puede escribirse de forma compacta como:
O
T 1 −1
ce = H w e , H=
0 1
C
T
ue (ξ) = φT (ξ) HT we = [H φ(ξ)] we = wTe [H φ(ξ)] (3.76)
EN
u(0) = φT (0)c
u(1) = φT (1)c
l θ(0) = φ0 T (0)c
l θ(1) = φ0 T (1)c
donde
N
ξ2 ξ3
φ= 1 ξ
φ0 = 0 3ξ 2
1 2ξ
IO
En notación matricial:
1 0 0 0 wa 1 0 0 0 c1
C
0 l 0 0 θa = 0 1 0 0 c2
0 0 1 0 wb 1 1 1 1 c3
0 0 0 l θb 0 1 2 3 c4
C
Finalmente para una viga a flexión usando una función de interpolación cúbica:
1 0 0 0 1 0 −3 2
RU
0 l 0 0 0 1 −2 1
H= = Hl H̄ (3.77)
0 0 1 0 0 0 3 −2
0 0 0 l 0 0 −1 1
2 0
Z l
1
Te = ρA(x)u̇(x)u̇(x)dx
2 0
1 l
Z
T
= ρA(x) ẇTe [H φ(x)] [H φ(x)] ẇe dx
2 0
1
= ẇTe Me ẇe
2
N
0
IO
Z l 2
1 ∂u
Ve = EA(x) (x, t) dx
2 0 ∂x
Sustituyendo:
C
Z L iT h
1 h T T i
Ve = EA(x) Hφ0 (ξ) we Hφ0 (ξ) we dx
2 0
C
1 L
Z
T
EA(x) wTe Hφ0 (ξ) Hφ0 (ξ) we dx
=
2 0
RU
1
= wTe Ke we
2
donde
Z l T
Z l h i
EA(x) Hφ0 (x) Hφ0 (x) dx = H EA(x) φ0 (x)φ0T (x)dx HT
Ke =
ST
0 0
es una matriz de rigidez para el elemento. Como dx = ldξ puede notarse que:
Sustituyendo:
1
O
Z
1
Ke = H EA(x) φ0 (ξ)φ0T (ξ)dξ HT (3.79)
l 0
C
L 2
∂2u
Z
1
Ve = EI(x) (x, t) dx
2 0 ∂x2
N
0 1 ξ
IO
1 1 ξ
φ(ξ)φT (ξ) =
1 ξ =
ξ ξ ξ2
0 0 0
φ(ξ)0 φ0T (ξ) =
0 1 =
1 0 1
C
Si la densidad por unidad de longitud es constante:
Z 1
C
1 −1 1 ξ 1 0 ρAl 2 1
Me = ρAl 2 dξ =
0 1 0 ξ ξ −1 1 6 1 2
RU
Para la rigidez, si EA fuese constante:
Z 1
1 −1 0 0 1 0 EA 1 −1
Ke = EAl dξ =
0 1 0 0 1 −1 1 l −1 1
ST
H̄ = , φ= 2 , ,
0 0 3 −2
ξ 2ξ 2
3 2
0 0 −1 1 ξ 3ξ 6ξ
O
ξ2 ξ3
1 ξ 0 0 0 0
ξ ξ2 ξ3 ξ4 0 1 2ξ 3ξ 2
φφT = φ0 φ0T
, =
C
ξ 2 ξ3 ξ4 ξ5 0 2ξ 4ξ 2 6ξ 3
ξ3 ξ4 ξ5 ξ6 0 3ξ 2 6ξ 3 9ξ 4
0 0 0 0
EN
0 0 0 0
φ00 φ00T =
0
0 4 12ξ
0 0 12ξ 36ξ 2
0 0 0 0 0 0 0 0
N
Z 1
0 0T
0 1 2/2 3/3 1
0 30 30 30
φ (ξ)φ (ξ)dξ =
0
=
0 2/2 4/3 6/4 30 0 30 40 45
IO
0 3/3 6/4 9/5 0 30 45 54
0 0 0 0 0 0 0 0
Z 1 0
00 00T 0 0 0 0 0 0 0
φ (ξ)φ (ξ)dξ = = 360
C
0 0 4 12/2 0 0 4 6
0
0 0 12/2 36/3 0 0 6 12
C
Con esto se obtiene:
156 22 54 −13 12 6 −12 6
RU
ρAl
22 4 13 −3 EI 6 4 −6 2
Me = , Ke = 3
420 54 13 156 −22 l −12 −6 12 −6
−13 −3 −22 4 6 2 −6 4
que son válidas si usamos como coordenadas generalizadas los desplazamientos y los
giros multiplicados por las longitudes de los elementos (θl). Si usamos directamente los
ST
giros θ tenemos que escalar estas matrices usando Hl definida en (3.77); por ejemplo
para la rigidez K0 = Hl KHTl .
Matrices de Masa y Rigidez del Sistema Para obtener las ecuaciones dinámicas a
N
E E
X 1X T
T = Te = ẇ Me ẇe
e=1
2 e=1 e
C
E E
X 1X T
V = Ve = w Ke we
e=1
2 e=1 e
Las energı́as de cada elemento dependen de las deformaciones de los nodos que
EN
definen dicho elemento, y en general esos nodos son compartidos por más de un
elemento.
Podemos coleccionar las deformaciones we de cada nodos en un único vector global w,
y poner las energı́as en función de este.
Entonces para cada elemento el vector de deformaciones será:
w e = Ce w
donde M y K son matrices de rigidez para todo el sistema, mientras que M̄e y K̄e
N
son las que corresponden a cada elemento expandidas para trabajar con el vector de
deformaciones globales:
IO
M̄e = CTe Ke Ce , K̄e = CTe Me Ce
C
..
.
w
C
a
u1 ··· 0 1 0 ··· 0 0 0 ···
..
= .
u2 ··· 0 0 0 ··· 0 1 0 ···
wb
RU
..
.
.. ..
..
. . ··· . · · ·
· · · m̄11 ··· m̄12 · · ·
M̄e = ... .. .. .. ..
. . . .
· · · m̄21 ··· m̄22 · · ·
N
.. .. ..
··· . ··· . .
O
Dado que solo unos pocos elementos de estas matrices expandidas serán no nulos,
es posible resolver esta expansión al momento de construir las matrices M y K
C
mismo sistema de coordenadas es necesario definir uno local para cada elemento y
proyectar los desplazamientos locales en dicho sistema global. Para ello utilizamos la
matriz de rotación correspondiente:
r ei = Re r̄ i (3.81)
N
IO
C
C
RU
Figura 3.12: estructura reticulada discretizada por elementos viga
x cos β sin β x̄
=
y − sin β cos β ȳ
1 1 −1
R3 = R7 = √
2 1 1
O
1 1 −1 0 0
C 3 R3 = √
2 0 0 1 −1
Con esto:
EN
x̄1
1 1
x1 −1 0 0 ȳ1
=√
x2 2 0 0 1 −1
x̄2
ȳ2
1 T 1 T
N
T = ẇ Mẇ , V = w Kw
2 2
El trabajo virtual de las fuerzas no conservativas se puede incluir de forma similar:
IO
F
X
δW = f Tk δwk = FT δw
k=1
C
donde
F
X
F = fk
C
k=1
Con esto:
RU
∂L ∂L ∂W
= Mẇ , = −Kw , =F
∂ ẇi ∂wi ∂wi
Si incluimos fuerzas de roce viscoso en el lado derecho, reemplazando en la ecuación
de Lagrange se obtiene finalmente:
ST
El resultado será equivalente a sumar los valores de masa y rigidez de los elementos
concentrados al elemento de la diagonal principal de las matrices de masa y rigidez
global correspondientes a los nodos asociados a dichos elementos.
N
p. 563].
IO
3.3.4. Diferencias Finitas
Desde un punto de vista puramente matemático podemos resolver numéricamente las
ecuaciones diferenciales discretizando el dominio de integración.
En cada nodo se plantea la ecuación diferencial correspondiente, sustituyendo las
C
derivadas por diferencias entre los valores vecinos para las deformaciones; y asegurando
el cumplimiento de las condiciones de borde correspondientes.
El resultado es una ecuación similar a la 3.82, con coeficientes algo diferentes.
C
El tema se puede consultar en [Rao11, sec.11.6: Finite Difference Method for Continuous
Systems, p. 951]).
RU
3.3.5. Análisis Modal de Respuesta Estructural
Agregando una ecuación de salida a esta ecuación de estados 3.14 es posible encontrar
una matriz de transferencia entre dichas salidas y el vector de entradas u. Sin embargo,
ST
esto puede ser laborioso o inexacto para sistemas con muchos grados de libertad, y en
general innecesario.
resultará de la forma:
donde
N = Λ−1 M−1 F (3.84)
El lado derecho de la ecuación 3.83 expresa las fuerzas generalizadas para las
EN
y = Cx + Df (3.85)
y = Cq (3.86)
N
Como el lado derecho en (3.83) constituye un conjunto de ecuaciones diferenciales
de segundo orden desacopladas, y cualquier variable de salida puede expresarse
como combinación lineal de las coordenadas modales según (3.87); puede obtenerse
IO
un modelo simplificado incluyendo solo aquellos modos naturales que contribuyan
significativamente a la combinación lineal para las variables de salida elegidas. La
cantidad de modos a considerar dependerá de la precisión que se requiera, y del rango
de frecuencias a considerar para la perturbación.
C
En el caso de sistemas amortiguados, el cambio de coordenadas propuesto para el
modelo sin amortiguamiento no siempre resulta en un sistema desacoplado.
C
Pero en el caso donde la matriz B de la ecuación 3.14 puede expresarse como αM+βK,
con α, β ∈ R, el cambio de coordenadas resulta en una matriz de amortiguamiento
diagonal.
RU
MΛη̈ + BΛη̇ + KΛη = Fu → η̈ + Cη̇ + λη = Fu
donde
= αI + βλ
IO
C
Apéndices
C
RU
Dinámica del Avión
ST
Casi todas las aeronaves constituyen sistemas “sub-actuados”, ya que no se cuenta con
actuación directa sobre cada grado de libertad. Su control de vuelo depende fuertemente
O
145
N
IO
C
C
Figura A.1: Superficies de control de vuelo
Respecto de una terna de ejes coordenados solidaria a la aeronave (terna móvil o terna-
b) podemos considerar los movimientos de rolido (roll), cabeceo (pitch) y guiñada (yaw)
que son rotaciones alrededor del los ejes xb , y b y z b respectivamente.
Logramos estos movimientos actuando sobre los alerones, el elevador y el timón de
dirección respectivamente, que constituyen los comandos primarios. Debe advertirse
N
Para evitar la “guiñada adversa” se puede combinar el accionamiento de los alerones con
los spoilers, que son superficies auxiliares concebidas para disminuir la sustentación y
aumentar la resistencia aerodinámica.
C
En algunas aeronaves en lugar del elevador se utiliza toda la superficie del estabilizador
horizontal (o del cannard). Algo similar se observa en configuraciones cannard, mientras
que en aquellas sin planos auxiliares las superficies en el borde de fuga del ala cumplen
funciones múltiples.
EN
C
“Trimado”
C
En general la mayor parte del vuelo de una aeronave se realiza en una condición de
vuelo recto. Por una cuestión de comodidad para el pilotaje es necesario realizar ajustes
para que la aeronave quede equilibrada en dicha condición de vuelo con los mandos en
RU
su posición “neutra”. A esta acción se la denomina trimado del avión.
superficie móvil en el elevador denominada trim tab, que permite establecer su ángulo
de deflexión cuando los mandos se dejan en la posición neutra.
necesario aumentar la resistencia usando flaps, spoilers, frenos aerodinámicos (si están
disponibles) e incluso el tren de aterrizaje si se desea realizar el descenso a velocidad
constante. En aeronaves que no cuentan con estos comandos (aeronaves pequeñas
O
como por ejemplo el Piper J3) es necesario realizar una maniobra de “deslizamiento”
para utilizar el fuselaje como freno.
C
masa de aire.
Actitud
Al hablar de actitud nos referimos a la rotación necesaria para llevar una terna de
referencia que pueda considerarse “inercial” (terna-i) hasta hacerla coincidir con la terna
del cuerpo (terna-b) o viceversa.
C
lineales ({ub , v b , wb }) y angulares ({p, q, r}) proyectadas en dicha terna
C
Esta rotación puede describirse mediante ángulos de Euler, matrices de rotación, ángulo
vector o cuaterniones unitarios.
Utilizando ángulos de Euler, de las 21 combinaciones posibles de rotaciones ordenadas
RU
alrededor de ternas cartesianas vamos a considerar la secuencia (intrı́nseca) de Tait-
Bryan (z y 0 x00 ) que implica una rotación inicial alrededor del eje z inicial, seguida por una
rotación alrededor del eje y 0 (eje y de la terna resultante de la primer rotación) y una
rotación alrededor del eje x00 (eje x de la terna precedente).
Solemos referirnos a estos giros como “rolido”, “cabeceo” y “rumbo”, aunque esto puede
generar una cierta confusión en tanto por “cabeceo” podrı́amos referirnos a la elevación
ST
respecto del horizonte o a una rotación respecto del eje y b como habı́amos mencionado
precedentemente. En otros contextos estos mismos ángulos también se denominan
como “azimut”, “elevación” y “rotación propia”.
mostrado en la figura A.3, con un eje x alineado con el fuselaje y apuntando hacia
adelante, un eje z contenido en plano de simetrı́a apuntando hacia abajo, y un eje y
completando una terna derecha.
O
Como terna “cuasi”-inercial se puede adoptar alguna de las mostradas en la figura A.2a.
Lo normal es usar la terna NED (North East Down), con eje z alineado con la “vertical del
lugar”, de forma tal que cuando las ternas móvil y fija se alinean, el vehı́culo mira hacia
C
Incidencia
EN
La incidencia se describe mediante dos ángulos que llamamos ángulo de ataque (en
el plano vertical de la aeronave) y ángulo de deslizamiento (en su plano horizontal),
identificados en la figura A.2a como α y β respectivamente.
Estos dos ángulos dependen únicamente de las componentes de la velocidad relativa
entre la aeronave y el aire proyectada en la terna móvil, y esta velocidad relativa depende
de del vector velocidad de la aeronave y de la velocidad la masa de aire (que podrı́a estar
wb + ww vb + vw
α = tan−1 , β = tan−1 (A.1)
ub + uw ub + uw
N
proyección de las velocidades en dicha terna, afectando a los ángulos de incidencia.
Pero con la posterior evolución del movimiento estas componentes en general sufrirán
modificaciones adicionales, pudiendo incluso volver a los valores iniciales aun con la
IO
nueva actitud. Por ejemplo, en una condición de vuelo recto en régimen estacionario el
ángulo de ataque depende unicamente de la velocidad, siendo afectada mı́nimamente
por el ángulo que forme esta trayectoria respecto del horizonte.
C
Navegación
Al considerar problemas de navegación de una aeronave se busca describir matemáti-
camente su desplazamiento respecto de una terna de navegación (terna-n).
C
Consideramos a la aeronave como un elemento puntual animado con un vector velocidad
cuasi-constante (es decir, si varı́a lo hace lentamente en relación a sus modos naturales
de respuesta), independientemente de los fenómenos fı́sicos que la definen. Por lo tanto
RU
hablamos de la cinemática del punto en un espacio R3 , y podemos separar el análisis
en un plano horizontal respecto del vertical.
Para el plano vertical definimos un ángulo γ de ascenso o descenso (planeo) como aquel
formado por el vector velocidad respecto del horizonte en dicho plano:
ST
N
Para el plano horizontal definimos el rumbo χ como aquel el formado por la velocidad
respecto de una referencia fija en el terreno (un punto cardinal, la orientación de una
pista, etc.)
C
EN
χ=ψ+β (A.3)
C
A.1.3. Modelo de Cuerpo Rı́gido
Es posible modelar cualquier vehı́culo como un cuerpo rı́gido si sus modos elásticos
RU
tienen frecuencias naturales altas en relación a la escala de tiempo asociada al
movimiento de su centro de gravedad. En tal caso el estado del cuerpo puede describirse
en principio mediante un vector velocidad lineal v = {u, v, w} y un vector velocidad
angular ω = {p, q, r}, que podemos referir a una terna solidaria al vehı́culo (terna móvil,
terna del cuerpo o terna b).
ST
dh
+ω×h=M
dt
O
v̇ = −ω × v + f (A.4)
−1
ω̇ = −J ω × Jω + τ (A.5)
N
Con esto es fácil demostrar que:
qw − rv q h̄z − rh̄y
IO
ω × v = ru − pw , ω × h̄ = rh̄x − ph̄z
pv − qu ph̄y − q h̄x
C
Ix 0 Ixz Ix p + Ixz r
J= 0 Iy 0 → Jω = Iy q
C
Ixz 0 Iz Iz r + Ixz p
Entonces, las ecuaciones para un cuerpo rı́gido simétrico respecto del plano xz con
ST
u̇ + qw − rv X
m v̇ + ru − pw = Y (A.7)
ẇ + pv − qu Z
O
Hasta este punto el vehı́culo puede ser de cualquier naturaleza, en tanto la modelización
como cuerpo rı́gido sea válida según lo explicitado al inicio de esta sección.
C
Hay que tener en cuenta que los ángulos de incidencia se relacionan directamente con
las componentes del vector de velocidad relativa respecto del aire descriptas por las
ecuaciones (A.1), por lo cual estos ángulos de incidencia no agregan nuevos estados.
La densidad ρ y el número de Mach M son cantidades que dependen de la condición de
N
vuelo (velocidad y altura) que cambian lentamente en relación a la respuesta dinámica.
Por lo tanto podemos considerarlos como parámetros lentamente variantes en el tiempo.
IO
Es habitual utilizar un modelo aerodinámico cuasi-estacionario. En este se utilizan las
expresiones de fuerza y momento aerodinámico en estado estacionario con los valores
instantáneos de presión dinámica y ángulos de incidencia.
C
Para las fuerzas la expresión es de la forma:
fa = σ S̄ cf (A.8)
C
donde cf es un coeficiente de fuerza adimensional, S̄ es una superficie de referencia (en
general la superficie alar); mientras que para los momentos:
RU
ma = σ S̄ ¯l cm (A.9)
modelo.
Por ejemplo, para el momento de cabeceo (M en la (A.6)) es común usar la expresión:
c̄ c̄
cm (α, α̇, q, δe ) = cm (α) + cmα̇ (α) α̇ + cmq (α) q + cme δe (A.10)
2V 2V
N
∂cm ∂cm ¯ = c̄ α̇ c̄
O
cmα̇ = ¯ , cmq = , α̇ , q̄ = q
∂ α̇ ∂ q̄ 2V 2V
¯ y q̄ velocidades angulares adimensionalizadas con la velocidad angular de
siendo α̇
C
N
c(α)c(β) s(β) −s(α)c(β)
= −c(α)s(β) c(β) s(α)s(β) (A.11)
s(α) 0 c(α)
IO
Esta matriz permite proyectar cantidades en la terna aerodinámica a la terna móvil.
C
orientación o actitud del vehı́culo.
Usaremos ángulos de Euler con la definición de Tait-Byan partiendo de una terna
terrestre NED, que identificamos como rumbo, cabeceo y rolido: θ = {φ, θ, ψ}.
C
La matriz de rotación para realizar la proyección de terna terrestre a una móvil resulta:
RU
1 0 0 c (θ) 0 −s (θ) c (ψ) s (ψ) 0
Cbi = Cφ Cθ Cψ = 0 c (φ) s (φ) 0 1 0 −s (ψ) c (ψ) 0
0 −s (φ) c (φ) s (θ) 0 c (θ) 0 0 1
c (θ) c (ψ) c (θ) s (ψ) −s (θ)
= s (φ) s (θ) c (ψ) − c (φ) s (ψ) s (φ) s (θ) s (ψ) + c (φ) c (ψ) s (φ) c (θ)
ST
c (φ) s (θ) c (ψ) + s (φ) s (ψ) c (φ) s (θ) s (ψ) − s (φ) c (ψ) c (φ) c (θ)
Dado que la gravedad solo tiene componente vertical, de la matriz de rotación interesa
solo la última columna, que con la secuencia de rotaciones adoptada no depende del
rumbo ψ :
0 − sin θ
N
φ̇ 1 sin (φ) tan (θ) cos (φ) tan (θ) p
θ̇ = 0 cos (φ) − sin (φ) · q = Cωė ω
0 sin (φ) sec (θ) cos (φ) sec (θ) r
ψ̇
EN
es decir:
N
de deflexiones de las superficies de control, f p y mp son fuerzas y momentos asociados
a la propulsión.
IO
A.1.5. Linealización
Equilibrio
Normalmente se analiza la dinámica de una aeronave mediante un modelo linealizado
C
alrededor de una condición de vuelo estacionaria (en el sentido dinámico, no en el
cinemático), lo que implica actitud y velocidad de traslación constante.
Teniendo en cuenta que actitud constante implica ω = 0, las ecuaciones (A.4) y (A.5) se
C
convierten en:
g b (θ) + f a (ρ, v, δ) + f p = 0
RU
τ a (ρ, v, δ) + mp = 0
Como este modelo no depende del rumbo, la condición de vuelo podrı́a ser no
solo la de un vuelo rectilı́neo en cualquier dirección, sino también la de un viraje
sostenido; incluyendo una trayectoria en barrena o un tirabuzón (agregando al modelo la
ST
f (ρ, ue , we , δ) = −g b (0, θe , 0) − f p
N
una condición de vuelo especificada mediante el vector velocidad. Esto, que implica
proponer parte del estado de equilibrio, requiere determinar las variables de entrada para
la condición de equilibrio (comandos aerodinámicos y el empuje del motor), en lugar de
C
N
donde cL (α, M ) es un coeficiente adimensional (coeficiente de sustentación) que
depende solo de la forma de la aeronave (pero no de sus dimensiones), del número
de Mach y del ángulo de ataque. Es posible entonces para un cierto número de
IO
Mach determinar el ángulo de ataque en equilibrio calculando primero el coeficiente de
sustentación, y luego el ángulo de ataque correspondiente:
mg
cL (αe ) = cos γe
0,5ρ|v|2 S
C
en donde hemos despreciado la proyección del empuje sobre el eje z w ; pero la
componente del empuje en la dirección de vuelo puede determinarse calculando la
C
resistencia con una expresión similar a la (A.15) y una primer estimación del ángulo
de ataque, y luego incluirla en el cómputo de cL .
RU
Por otra parte, para mantener la velocidad angular nula se requiere que la suma de
momentos se anule. En general esto requiere ajustar el comando de elevador para
”trimar” la aeronave de forma tal que:
cm (αe , δhe ) = 0 (A.16)
ST
Modelo de Estados
Para linealizar los términos inerciales (lado derecho de ec. 1) usamos los siguientes
desarrollos en serie de Taylor:
1 1
∆α ≈ ∆w , ∆β ≈ ∆v
U0 U0
EN
∂P ∂P
P ≈ P (x̄, ȳ) + (x − x̄) + (y − ȳ) + · · ·
∂x x̄,ȳ ∂y x̄,ȳ
P ≈ P (x̄, ȳ) + Px ∆x + Py ∆y + · · ·
Fa σ, M, α, β, α̇, β̇, δ ≈ F0 + Fσ ∆σ + FM ∆M
+ Fα ∆α + Fβ ∆β + Fα̇ ∆α̇ + Fβ̇ ∆β̇
+ Fδ1 ∆δ1 + Fδ2 ∆δ2 + · · · + Fδm ∆δm
Debe tenerse en cuenta que las superficies sustentadoras separadas del centro de
N
gravedad experimentarán un ángulo de incidencia que dependerá de la composición
entre ángulo de ataque de la aeronave y de los movimientos de rotación. Por lo tanto,
también debe considerarse la dependencia de fuerzas y momentos aerodinámicos con
IO
la velocidad angular.
La derivada del ángulo de ataque (y por lo tanto ẇ) tendrá influencia en las fuerzas
y momentos aerodinámicos, debido a las variaciones en el downwash sobre el plano
C
horizontal de cola.
C
dirección ζ y empuje τ :
η
ξ
δ= (A.17)
RU
ζ
τ
= Mu ∆u + Mv ∆v + Mw ∆w + Mp ∆p + Mq ∆q + Mr ∆r + · · ·
Si ω 0 = 0:
Ix ṗ + Ixz ṙ = Lu ∆u + Lv ∆v + Lw ∆w + Lp ∆p + Lq ∆q + Lr ∆r + · · ·
EN
Iy q̇ = Mu ∆u + Mv ∆v + Mw ∆w + Mp ∆p + Mq ∆q + Mr ∆r + M w∆w + · · ·
Iz ṙ + Ixz ṗ = Nu ∆u + Nv ∆v + Nw ∆w + Np ∆p + Nq ∆q + Nr ∆r + · · ·
mu̇ = Xu ∆u + Xv ∆v + Xw ∆w + Xp ∆p + Xq ∆q + Xr ∆r + ∆Xg + Xw ∆w + · · ·
mv̇ = Yu ∆u + Yv ∆v + Yw ∆w + Yp ∆p + Yq ∆q + Yr ∆r + ∆Yg + · · ·
mẇ = Zu ∆u + Zv ∆v + Zw ∆w + Zp ∆p + Zq ∆q + Zr ∆r + ∆Zg + Zw ∆w + · · ·
Xv = Xp = Xr = Zv = Zp = Zr = Mv = Mp = Mr = 0
Xξ = Xζ = Zξ = Zζ = Lη = Nη = Mξ = Mζ = 0
φ̇ ≈ p , θ̇ ≈ q
N
∆Xg ≈ −mg cos θo , ∆Yg ≈ mg cos θo , ∆Zg ≈ −mg sin θo
Con esto, reutilizando las variables de estado originales u, v, · · · para denotar las
IO
variaciones ∆u, ∆v, · · · :
m −Xẇ 0 0 u̇ Xu Xw Xq Xθ u Xη Xτ
0 m − Zẇ 0 ẇ
0 Zu Zw Zq Zθ w Zη Zτ η
= +
0 −Mẇ Iy 0 q̇ Mu Mw Mq 0 q Mη Mτ τ
C
0 0 0 1 θ̇ 0 0 1 0 θ 0 0
m 0 0 0 v̇ Yv Yp Yr Yφ v Yξ Yζ
C
ṗ
0 I −I 0 L Lp Lr 0 p L Lζ
ξ
x xz v ξ
0 −Ixz
= +
Iz 0 ṙ
Nv Np Nr 0 r
Nξ Nζ ζ
0 0 0 1 φ̇ 0 0 1 0 φ 0 0
RU
Estas ecuaciones tienen la forma:
ẋ = Ax + Bu
en donde:
A = M−1 Ā , B = M−1 B̄
N
u̇ xu xw xq xθ u xη xτ
ẇ
zu zw zq zθ w zη zτ
η
C
= +
q̇
m u mw mq 0 q
m η mτ τ
θ̇ 0 0 1 0 θ 0 0
Si aproximamos (u, w U0 ):
EN
w w
α = tan−1 ≈
u + U0 U0
u̇ xu xα xq xθ u xη xτ
α̇
α αα αq αθ α α ατ
η
u η
= +
q̇
mu mα mq 0 q
mη mτ τ
θ̇ 0 0 1 0 θ 0 0
Debe notarse que las derivativas asociadas con θ surgen fundamentalmente por efecto
gravitatorio, mientras que xq y αq tienen que ver con los términos de acoplamiento que
surgen al usar un sistema de coordenadas no inercial.
Si se quisiera incluir la perturbación aerodinámica {uw , ww , qw } se tendrı́a una matriz de
entrada adicional compuesta por las primeras tres columnas de la matriz A eliminando
xq y zq .
N
Podemos decir que para la dinámica longitudinal el modelo resulta:
u̇ xu xα −we −g cos θe u xη xτ
IO
−g/U
α̇
α α u /U sin θ α zα ατ
η
= u α e 0 0 e η
q̇
mu mα mq 0 q + mη mτ τ
(A.18)
θ̇ 0 0 1 0 θ 0 0
xu xw 0
zα zα /U0 0 uw
C
+ w
mu mw mq w
qw
0 0 0
C
donde los coeficientes de la matriz dinámica para una terna principal de inercia son:
2 σS σS ∂cx 1
RU
xu = cx (αe , 0) xw = (αe , 0)
Ve m m ∂α Ve
2 σS σS ∂cz 1
zu = cz (αe , 0) zw = (αe , 0)
Ve m m ∂α Ve
2 σSc̄ σSc̄ ∂cm 1 σSc̄ ∂cm
mu = cm (αe , 0) mw = (αe , 0) mq = (αe , 0)
ST
Ve Iy Iy ∂α Ve Iy ∂q
q̇ mw mq q mη
El polinomio caracterı́stico es:
C
s2 − (mq + zw ) s + (mq zw − mw zq ) = 0
zq = ≈ Ue
m − Zẇ
s2 − (mq + zw ) s + (mq zw − mw Ue ) = 0
Vemos que la frecuencia queda definida por:
r
p 1 Zw 1 p
ωsp ≈ mq zw − mw Ue = p Mq − Mw Ue = p Mα + zα Mα̇
Iy m Iy
N
cabeceo.
IO
Si se puede asumir ẇ ≈ 0 y que q̇ ≈ 0 (ver [Coo07, p. 152]):
u̇ xu xw xq xθ u xη xτ
0
zu zw zq zθ w zη zτ
η
C
= +
0 mu mw mq 0 q mη mτ τ
θ̇ 0 0 1 0 θ 0 0
C
Operando se llega a que:
mu Ue − mq zu mu zw − mw zu
RU
s2 − xu − xw s+g =0
mw Ue − mq zw mw Ue − mq zw
Sustituyendo:
√ g 1 CD
ωlp ≈ 2 , ζlp ≈ √
Ue 2 CL
Pude notarse que la frecuencia del modo fugoide depende básicamente de la velocidad
N
v v
β = tan−1 ≈
u + U0 U0
EN
β̇
yβ yp yr yφ β yζ yξ
ṗ lβ lp lr 0 p lζ lξ
ζ
= +
ṙ
nβ np nr 0 r
nζ nξ ξ
φ̇ 0 0 1 0 φ 0 0
donde:
N
σS ∂cy 1
yv = (αe , 0)
m ∂β Ve
IO
σS b̄ ∂cl 1 σS b̄ ∂cl σS b̄ ∂cl
lv = (αe , 0) lp = (αe , 0) lr = (αe , 0)
Ix ∂β Ve Ix ∂p Ix ∂r
σS b̄ ∂cn 1 σS b̄ ∂cn σS b̄ ∂cn
nv = (αe , 0) np = (αe , 0) nr = (αe , 0)
Iz ∂β Ve Iz ∂p Iz ∂r
C
C
A.1.8. Trayectoria
Plano Vertical
RU
La relación entre el ángulo de planeo γ y la velocidad vertical es:
ż = V sin γ ≈ V · γ
Asumiendo que se cuenta con alguna clase de control sobre el ángulo de ataque, y que
u̇ es pequeño, se puede plantear que:
q̇ mq 0 q mα
= + α
θ̇ 1 0 θ 0
N
ż q
y =γ= = 0 1 + −1 α
V θ
C
Plano Horizontal
Para la cinemática de la trayectoria horizontal se tiene algo similar a lo planteado para el
desplazamiento vertical:
EN
ẏ = V sin ψ ≈ V · ψ (A.20)
N
IO
C
C
RU
ST
N
O
C
EN
IO
C
C
1.1. Diagrama idealizado de la evolución del azúcar en sangre (rojo) y de la insulina (azul)
en humanos a lo largo del dı́a con tres comidas; incluyendo el efecto de comidas
RU
enriquecidas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3. distribución espacial de velocidades para el modelo de vórtice de Lamb-Oseen
(izquierda) y variación de la componente longitudinal de la velocidad al atravesar el
vórtice a velocidad constante a diferentes alturas (derecha) . . . . . . . . . . . . . . . 8
1.4. campo de velocidades para el péndulo plano . . . . . . . . . . . . . . . . . . . . . . . 13
1.5. campo de velocidades y trayectorias para el péndulo plano . . . . . . . . . . . . . . . 14
ST
1.6. puntos crı́ticos aislados - arriba: nodo (izquierda) y foco (derecha); abajo: centro
(izquierda) y punto de silla (derecha). . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.7. respuesta de un oscilador de van der Pol a diferentes condiciones iniciales . . . . . . . 18
1.8. evolución temporal de los estados del modelo de Lorenz con diferentes condiciones
iniciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
N
[Link] para el caso con autovalores reales, uno positivo y otro negativo. Las lı́neas
marcan la dirección de los autovectores. . . . . . . . . . . . . . . . . . . . . . . . . . . 34
[Link] para el caso con autovalores complejos conjugados. . . . . . . . . . . . . 35
[Link] temporal de una dinámica de primer orden con |p| = 1 para el caso estable
EN
163
2.4. Respuesta al escalón de un modelo de segundo orden con ξ = 0,5, ωn = 1, y con un
cero a diferentes distancias del origen . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
2.5. Polos (cruces) y ceros (cı́rculos) de las funciones de transferencia entre elevador η y
empuje τ y ángulo de ataque α de una aeronave . . . . . . . . . . . . . . . . . . . . . 64
2.6. desplazamiento en frecuencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.7. respuesta en frecuencia de primer orden . . . . . . . . . . . . . . . . . . . . . . . . . 73
2.8. respuesta en frecuencia de un integrador y un derivador . . . . . . . . . . . . . . . . . 74
2.9. respuesta en frecuencia de segundo orden sub-amortiguada . . . . . . . . . . . . . . 76
[Link] de banda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
N
[Link]́n en serie de Fourier de una onda cuadrada truncada en el séptimo
armónico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
[Link] funcional de un conversor A/D . . . . . . . . . . . . . . . . . . . . . . . . . . 87
IO
[Link] y reconstrucción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
C
3.3. Modelo conceptual para la suspensión de un automóvil con GL . . . . . . . . . . . . . 96
3.4. modelo de parámetros concentrados con N grados de libertad . . . . . . . . . . . . . . 97
3.5. solicitaciones sobre un elemento diferencial de una cuerda . . . . . . . . . . . . . . . 114
C
3.6. solicitaciones longitudinales sobre un elemento diferencial de una barra . . . . . . . . 115
3.7. solicitaciones torsionales sobre un elemento diferencial de una barra . . . . . . . . . . 120
3.8. solicitaciones sobre un elemento diferencial de una barra a flexión . . . . . . . . . . . 121
RU
[Link] estructurales de parámetros concentrados . . . . . . . . . . . . . . . . . . . 128
[Link] de interpolación para una barra sometida a deformaciones unidimiensionales
(torsión o tracción) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
[Link] reticulada discretizada por elementos viga . . . . . . . . . . . . . . . . . . . 140
IO
C
C
[Coo07] M.V. Cook. Flight Dynamics Principles. Elvesier, 2007.
[Kel12] S. Graham Kelly. Mechanical Vibrations - Theory and Applications. Cengage Learning, 5
edition, 2012.
RU
[Mei86] Leonard Meirovitch. Elements of Vibration Analysis. Mc Graw Hill, 2 edition, 1986.
[Mei97] Leonard Meirovitch. Principles ans Techniques of Vibrations. Prentice Hall, 2 edition, 1997.
[Oga93] Katsuhiko Ogata. Ingenierı́a de Control Moderna. Prentice Hall, 2 edition, 1993.
[O’N94] Peter V. O’Neil. Matemáticas Avanzadas para Ingenierı́a. Prentice Hall, 3 edition, 1994.
[SP01] Sigurd Skogestad and Ian Postlethwaite. Multivariable Feedback Control - Analysis and
Design. John Wiley & Sons, 2 edition, 2001.
O
165