Story Transcript
Departamento de Ingenier´ıa Mec´ anica Facultad de Ciencias F´ısicas y Matem´ aticas Universidad de Chile
Din´ amica Estructural Apuntes para el curso ME706
Viviana Meruane
´Indice general Contenidos
II
1 Introducci´ on
2
1.1
Sistemas de un grado de libertad . . . . . . . . . . . . . . . . . . .
2
1.1.1
Ecuaci´ on de movimiento, funci´on de transferencia . . . . . .
2
1.1.2
Polos del sistema, frecuencias naturales y factores de amortiguamiento . . . . . . . . . . . . . . . . . . . . . . . .
3
1.1.3
Soluci´ on de la ecuaci´on de movimiento . . . . . . . . . . . .
4
1.1.4
Residuos . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.1.5
Funci´ on de respuesta en frecuencia y funci´ on de respuesta a un impulso . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
Influencia de cambios en la masa, rigidez y amortiguaci´on .
8
Sistemas con multiples grados de libertad . . . . . . . . . . . . . .
10
1.2.1
Ecuaci´ on del sistema, funci´on de transferencia . . . . . . . .
10
1.2.2
Polos, frecuencias naturales y factores de amortiguamiento
12
1.2.3
Vectores modales, residuos
. . . . . . . . . . . . . . . . . .
12
1.2.4
Funci´on de respuesta en frecuencia (FRF) y funci´ on de respuesta a un impulso (IRF) . . . . . . . . . . . . . . . . .
14
Sistemas sin amortiguamiento y con amortiguamiento proporcional . . . . . . . . . . . . . . . . . . . . . . . . . .
16
Ortogonalidad y coordenadas modales . . . . . . . . . . . .
18
1.1.6 1.2
1.2.5 1.2.6
ii
´INDICE GENERAL
iii
2 El M´ etodo de Elementos Finitos 2.1
2.2
2.3
2.4
Elemento de barra . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
2.1.1
Coordenadas locales y globales . . . . . . . . . . . . . . . .
24
Elemento de viga . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
2.2.1
Coordenadas locales y globales . . . . . . . . . . . . . . . .
27
Elemento de viga 3D . . . . . . . . . . . . . . . . . . . . . . . . . .
28
2.3.1
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
Ensamble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
Rotaci´ on
3 Integraci´ on num´ erica de las ecuaciones de movimiento 3.1
3.2
4.2
4.3
35
M´etodos de integraci´on directa . . . . . . . . . . . . . . . . . . . .
35
3.1.1
El m´etodo de las diferencias centrales . . . . . . . . . . . .
36
3.1.2
El m´etodo de Wilson θ . . . . . . . . . . . . . . . . . . . . .
37
3.1.3
El m´etodo de Newmark . . . . . . . . . . . . . . . . . . . .
40
3.1.4
Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
Superposici´ on modal . . . . . . . . . . . . . . . . . . . . . . . . . .
43
3.2.1
43
Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4 Procesamiento de se˜ nales 4.1
21
45
La transformada de Fourier . . . . . . . . . . . . . . . . . . . . . .
46
4.1.1
Algunos par´ametros importantes . . . . . . . . . . . . . . .
48
Errores y ventanas . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
4.2.1
Aliasing . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
4.2.2
Leakage . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
4.2.3
Ventanas . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
Funciones de respuesta en frecuencia y funciones de coherencia . .
54
4.3.1
56
Efectos del ruido experimental . . . . . . . . . . . . . . . .
5 Estimaci´ on de par´ ametros modales
60
´INDICE GENERAL
5.1
5.2
5.3
5.4
iv
M´etodos de un grado de libertad . . . . . . . . . . . . . . . . . . .
61
5.1.1
Peak picking . . . . . . . . . . . . . . . . . . . . . . . . . .
61
5.1.2
Circle fitting . . . . . . . . . . . . . . . . . . . . . . . . . .
62
M´etodos con multiples grados de libertad en el dominio de frecuencias 64 5.2.1
M´etodo de m´ınimos cuadrados no lineales, LSFD . . . . . .
64
5.2.2
Rational Fractional Polynomials . . . . . . . . . . . . . . .
65
M´etodos con multiples grados de libertad en el dominio del tiempo
68
5.3.1
El m´etodo Ibrahim . . . . . . . . . . . . . . . . . . . . . . .
68
5.3.2
El m´etodo least-squares complex exponential (LSCE) . . .
70
Diagramas de estabilidad . . . . . . . . . . . . . . . . . . . . . . .
72
6 Medici´ on Experimental 6.1
6.2
6.3
6.4
74
An´ alisis previo a las mediciones . . . . . . . . . . . . . . . . . . . .
74
6.1.1
Rango de frecuencias . . . . . . . . . . . . . . . . . . . . . .
75
6.1.2
Selecci´ on de la ubicaci´on de las respuestas . . . . . . . . . .
75
6.1.3
Selecci´ on de los puntos de excitaci´on . . . . . . . . . . . . .
76
6.1.4
Selecci´ on de los puntos de suspensi´on . . . . . . . . . . . .
77
El montaje experimental . . . . . . . . . . . . . . . . . . . . . . . .
78
6.2.1
Mecanismo de Excitaci´on . . . . . . . . . . . . . . . . . . .
78
6.2.2
Aceler´ ometro . . . . . . . . . . . . . . . . . . . . . . . . . .
79
6.2.3
Sensor de fuerzas . . . . . . . . . . . . . . . . . . . . . . . .
80
Selecci´ on de la fuerza de excitaci´on . . . . . . . . . . . . . . . . . .
82
6.3.1
Excitaci´ on sinusoidal . . . . . . . . . . . . . . . . . . . . . .
82
6.3.2
Excitaci´ on aleatoria . . . . . . . . . . . . . . . . . . . . . .
83
6.3.3
Excitaci´ on pseudo-aleatoria . . . . . . . . . . . . . . . . . .
84
6.3.4
Excitaci´ on aleatoria en trenes (burst random) . . . . . . . .
84
6.3.5
Excitaci´ on de impacto . . . . . . . . . . . . . . . . . . . . .
85
Evaluaci´ on inicial de las FRFs medidas . . . . . . . . . . . . . . . .
86
´INDICE GENERAL
v
6.4.1
Repetibilidad . . . . . . . . . . . . . . . . . . . . . . . . . .
87
6.4.2
Reciprocidad . . . . . . . . . . . . . . . . . . . . . . . . . .
87
6.4.3
Linealidad . . . . . . . . . . . . . . . . . . . . . . . . . . . .
87
6.4.4
Caracter´ısticas especiales de una FRF . . . . . . . . . . . .
87
7 Correlaci´ on Num´ erico-Experimental y Ajuste de Modelos 7.1
7.2
7.3
Parear modelos num´ericos y experimentales . . . . . . . . . . . . .
92
7.1.1
T´ecnicas de reducci´on . . . . . . . . . . . . . . . . . . . . .
92
7.1.2
T´ecnicas de expansi´on . . . . . . . . . . . . . . . . . . . . .
94
T´ecnicas de correlaci´on . . . . . . . . . . . . . . . . . . . . . . . . .
95
7.2.1
Diferencia en las frecuencias de resonancia . . . . . . . . . .
96
7.2.2
Comparaci´on visual de los modos . . . . . . . . . . . . . . .
96
7.2.3
Modal Assurance Criterion . . . . . . . . . . . . . . . . . .
97
7.2.4
Frequency Response Assurance Criterion
. . . . . . . . . .
98
M´etodos iterativos de ajuste de modelos . . . . . . . . . . . . . . .
98
7.3.1
Modos y frecuencias naturales . . . . . . . . . . . . . . . . . 101
7.3.2
Funciones de respuesta en frecuencia . . . . . . . . . . . . . 102
8 Ejemplos de An´ alisis Modal en Estructuras Reales 8.1
Modelo num´erico . . . . . . . . . . . . . . . . . . . . . . . . 114
Tubo de escape de un auto . . . . . . . . . . . . . . . . . . . . . . 116 8.4.1
8.5
Modelo num´erico . . . . . . . . . . . . . . . . . . . . . . . . 111
Viga de concreto reforzado . . . . . . . . . . . . . . . . . . . . . . . 113 8.3.1
8.4
Modelo num´erico . . . . . . . . . . . . . . . . . . . . . . . . 108
Modelo a escala de un avion . . . . . . . . . . . . . . . . . . . . . . 110 8.2.1
8.3
104
Estructura de barras . . . . . . . . . . . . . . . . . . . . . . . . . . 104 8.1.1
8.2
89
Modelo num´erico . . . . . . . . . . . . . . . . . . . . . . . . 117
El puente I-40 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
´INDICE GENERAL
1
8.5.1
Modelo num´erico . . . . . . . . . . . . . . . . . . . . . . . . 121
Bibliograf´ıa
124
Cap´ıtulo 1
Introducci´ on 1.1.
Sistemas de un grado de libertad
1.1.1.
Ecuaci´ on de movimiento, funci´ on de transferencia
Consideremos el sistema de un grado de libertad mostrado en la Figura 1.1. La ecuaci´ on de movimiento de este sistema viene dada por: m¨ x + cx˙ + kx = f (t)
(1.1)
donde, m: masa k: coeficiente de rigidez c: coeficiente de amortiguaci´on x ¨, x, ˙ x: aceleraci´ on, velocidad y desplazamiento f : fuerza de excitaci´ on externa t: tiempo
2
SISTEMAS DE UN GRADO DE LIBERTAD
3
k f(t) m
x(t)
c
Figura 1.1: Transformando la ecuaci´ on 1.1 al dominio de Laplace y asumiendo que el desplazamiento y velocidad inicial son 0, se obtiene: (mp2 + cp + k)X(p) = F (p)
(1.2)
Lo que se puede escribir como: Z(p)X(p) = F (p)
(1.3)
Z(p) se denomina la rigidez din´amica. Invirtiendo la ecuaci´on 1.2 ´o 1.3 se obtiene la funci´on de transferencia H(p) = Z −1 (p): X(p)
=
H(p)
=
1.1.2.
H(p)F (p) p2
1/m + (c/m)p + (k/m)
(1.4) (1.5)
Polos del sistema, frecuencias naturales y factores de amortiguamiento
El denominador de la ecuaci´on 1.5 representa la ecuaci´on caracter´ıstica del sistema. Las ra´ıces de esta ecuaci´ on son: r c c 2 k λ1,2 = − ± − (1.6) 2m 2m m Esta ecuaci´on permite la introducci´on de varios conceptos importantes. Si no hay amortiguamiento, el sistema en estudio es un sistema conservativo (c=0), y la frecuencias natural (rad/seg) del sistema sin amortiguamiento se define como: r k ωn = (1.7) m
SISTEMAS DE UN GRADO DE LIBERTAD
4
Se define como amortiguamiento cr´ıtico, cc , al valor del amortiguamiento que hace que el termino dentro de la ra´ız en la ecuaci´on 1.6 sea cero: r k cc = 2m (1.8) m A partir de esta definici´ on, se define la raz´on de amortiguamiento, ζ, como: c ζ= (1.9) cc Las ra´ıces presentadas en la ecuaci´on 1.6, se pueden escribir como dos ra´ıces complejas conjugadas: λ
=
σ + jωd
(1.10)
λ∗
=
σ − jωd
(1.11)
donde σ es el factor de amortiguamiento, ωd es la frecuencia natural del sistema amortiguado y ∗ representa el complejo conjugado A partir de la definici´ on anterior se pueden definir las siguientes relaciones: λ
=
p −ζ + j 1 − ζ 2 ωn
(1.12)
ωn
=
q ωd2 + σ 2
(1.13)
ζ
=
σ −p 2 ωd + σ 2
(1.14)
σ
=
ωd
=
1.1.3.
−ζωn p ωn 1 − ζ 2
(1.15) (1.16)
Soluci´ on de la ecuaci´ on de movimiento
La soluci´ on de la ecuaci´on 1.1 en el tiempo se puede determinar asumiendo una soluci´ on de la forma: x = Aeλt
(1.17)
donde A es una constante. Por lo tanto: x˙
= λAeλt
(1.18)
x ¨
= λ2 Aeλt
(1.19)
SISTEMAS DE UN GRADO DE LIBERTAD
5
Sustituyendo en la ecuaci´on 1.1 y asumiendo que las fuerzas externas son nulas, se obtiene: mλ2 Aeλt + cλAeλt + kAeλt = 0
(1.20)
Dividiendo por Aeλt se obtiene la ecuaci´on caracter´ıstica del sistema: mλ2 + cλ + k = 0
(1.21)
cuya soluci´on son las ra´ıces de la ecuaci´ on 1.6. Por lo tanto, x = A1 eλ1 t y x = A2 eλ2 t , son soluciones del problema y su suma es la soluci´on general: x = A1 eλ1 t + A2 eλ2 t
(1.22)
Las soluciones particulares dependen de las constantes, A1 y A2 , que se determinan a partir de las condiciones iniciales. La ecuaci´ on 1.22 se puede escribir como: x = e−ζωn t A1 ejωd t + A2 e−jωd t
(1.23)
Utilizando las identidades: ejωd t = cosωd t + jsenωd t y e−jωd t = cosωd t − jsenωd t, la ecuaci´ on 1.23 se puede escribir como: x = e−ζωn t (Bcosωd t + Dsenωd t)
(1.24)
Amplitud
ζ1
0 Tiempo
Figura 1.2:
SISTEMAS DE UN GRADO DE LIBERTAD
6
Dependiendo del valor de la raz´ on de amortiguamiento el sistema se puede clasificar como sobreamortiguado (ξ > 1), con amortiguamiento cr´ıtico (ξ = 1) o con amortiguamiento d´ebil (ξ < 1). Como se muestra en la Figura 1.2, la respuesta de un sistema sobreamortiguado es decreciente sin oscilaciones. En cambio, la respuesta de un sistema con amortiguamiento d´ebil es una respuesta oscilatoria decreciente. Amortiguamiento cr´ıtico es el caso de borde entre un sistema sobreamortiguado y con amortiguamiento d´ebil. Para el caso de una respuesta forzada, en donde f (t) de la ecuaci´ on 1.1 es distinto de cero. La soluci´on viene dada por la soluci´on general cuando f = 0 m´as una soluci´ on particular X: x = A1 eλ1 t + A2 eλ2 t + X
(1.25)
La forma de X depende de la forma de la fuerza de excitaci´on, la tabla a continuaci´on muestra varios casos particulares: Forma de f (t) f (t) = A, constante f (t) = At f (t) = At2 f (t) = Asen(ωt) o Acos(ωt)
Soluci´on particular B Bt + D Bt2 + D Bcos(ωt) + Dsen(ωt)
Tabla 1.1:
1.1.4.
Residuos
Conociendo las ra´ıces del sistema, la ecuaci´on 1.5 se puede escribir como: H(p) =
1/m (p − λ)(p − λ∗ )
(1.26)
Ocupando fracciones parciales se tiene que: H(p) =
A A∗ 1/m + , con A = (p − λ) (p − λ∗ ) 2jωd
(1.27)
En esta formulaci´ on A y A∗ se denominan residuos.
1.1.5.
Funci´ on de respuesta en frecuencia y funci´ on de respuesta a un impulso
En las secciones anteriores se discute la relaci´on entre la fuerza y la respuesta de un sistema de un grado de libertad en el dominio de Laplace. Esta relaci´on
SISTEMAS DE UN GRADO DE LIBERTAD
7
tambi´en se puede expresar en el dominio de frecuencias o del tiempo. La funci´ on de transferencia evaluada en el eje de frecuencias (jω) se denomina funci´ on de respuesta en frecuencia (FRF): H(p)|p=jω = H(ω) =
A A∗ + (jω − λ) (jω − λ∗ )
(1.28)
En general, la contribuci´on de la parte del complejo conjugado (o la parte de la frecuencia negativa) es despreciable cerca de la resonancia, ω = ωd . En consecuencia, la funci´on de respuesta en frecuencia para un grado de libertad frecuentemente se aproxima por: H(ω) '
A (jω − λ)
(1.29)
Amplitud
La Figura 1.3 muestra un gr´afico de la amplitud y fase de una funci´on de respuesta en frecuencia. Notar que el peak en la FRF se obtiene cuando ω = ωd , a esta misma frecuencia la fase cambia en −180◦ .
Fase
0
-pi
ωd
Figura 1.3: Aplicando una transformada inversa de Laplace a la expresi´on de la funci´on de transferencia se obtiene la funci´ on de respuesta a un impulso. La respuesta a un impulso corresponde a la respuesta a un delta de Dirac en t = 0. Para un grado de libertad la funci´ on de respuesta a un impulso es: h(t) = eσt Aeiωd t + A∗ e−iωd t (1.30)
SISTEMAS DE UN GRADO DE LIBERTAD
8
El residuo A define la amplitud inicial, σ la tasa de decrecimiento y ωd la frecuencia de oscilaci´on. Notar que el gr´afico de la funci´on de respuesta a un impulso es equivalente al de la Figura 1.2.
1.1.6.
Influencia de cambios en la masa, rigidez y amortiguaci´ on
Las Figuras 1.4, 1.5 y 1.6 muestran como, la funci´on de respuesta en frecuencia de un sistema de 1 grado de libertad, se ve afectada por cambios en la rigidez, amortiguamiento y masa.
log(FRF)
Un incremento en la rigidez resulta en una frecuencia de resonancia mayor y un valor menor de la FRF en el rango de frecuencias bajas. Debido al efecto dominante de la rigidez a bajas frecuencias, esta regi´on de la FRF se denomina la linea de rigidez o “compliance line”.
frecuencia Figura 1.4: Incremento de la rigidez Un incremento del amortiguamiento produce una peque˜ na reducci´on de la frecuencia de resonancia. Sin embargo, la mayor contribuci´on es una disminuci´on de la amplitud de la funci´on de respuesta en frecuencia en las cercan´ıas a la resonancia. Tambi´en, la fase cambia m´as suavemente. Si el amortiguamiento es cero, la magnitud en la resonancia se vuelve infinito y la fase disminuye repentinamente en 180 grados. Adicionalmente, los polos del sistema, λ, se vuelven puramente imaginarios y en p amplitud iguales a la frecuencia natural sin amortiguamiento ωn = k/m.
9
log(FRF)
SISTEMAS DE UN GRADO DE LIBERTAD
frecuencia Figura 1.5: Incremento del amortiguamiento
log(FRF)
Un incremento de masa disminuye el valor de la frecuencia de resonancia. La amplitud de la FRF a altas frecuencias tambi´en disminuye. Debido al efecto dominante de la masa a altas frecuencias esta a´rea de la FRF se denomina la “linea de masa”.
frecuencia Figura 1.6: Incremento de la masa
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
1.2.
10
Sistemas con multiples grados de libertad
En esta secci´ on vamos a extender los conceptos vistos en la secci´ on anterior para el caso con multiples grados de libertad.
1.2.1.
Ecuaci´ on del sistema, funci´ on de transferencia
En la mayor´ıa de los casos el sistema no se puede describir como un sistema de un grado de libertad. La mayor´ıa de los sistemas consisten en un ensamble de un n´ umero infinito de masas, resortes y amortiguadores. A continuaci´on se mostrar´a como se relaciona la funci´on de transferencia de sistema con multiples grados de libertad con los par´ ametros modales (frecuencias de resonancia y vectores modales). En el caso de un sistema con multiple grados de libertad la ecuaci´on de movimiento equivale a la ecuaci´on de un grado de libertad, pero matricial. Para ilustrar esto, se desarrollar´a el caso de un sistema de dos grados de libertad mostrado en la Figura 1.7.
f2(t) x2 (t)
f1(t) x1 (t) k2
k1
m2
m1 c1
k3
c2
c3
Figura 1.7: Ejemplo de un sistema de dos grados de libertad Las ecuaciones de movimiento para este sistema son: m1 x ¨1 + (c1 + c2 )x˙ 1 − c2 x˙ 2 + (k1 + k2 )x1 − k2 x2
= f1
(1.31)
m2 x ¨2 + (c2 + c3 )x˙ 2 − c2 x˙ 1 + (k2 + k3 )x2 − k2 x1
= f2
(1.32)
En notaci´ on matricial: m1 0 x ¨1 c + c2 −c2 x˙ 1 k + k2 −k2 x1 f1 + 1 + 1 = (1.33) 0 m2 x ¨2 −c2 c2 + c3 x˙ 2 −k2 k2 + k3 x2 f2
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
11
que se puede escribir como: M {¨ x} + C {x} ˙ + K {x} = {f }
(1.34)
donde, M: matriz de masa K: matriz de rigidez C: matriz de amortiguaci´on {x}: vector de respuesta {f }: vector de fuerzas Esta ecuaci´on tambi´en describe el comportamiento de sistemas con m´as grados de libertad. La dimensi´on de las matrices aumenta con el n´ umero de grados de libertad. Transformando esta ecuaci´ on de movimiento en el dominio de Laplace, asumiendo que las velocidades y desplazamientos iniciales son cero, se obtiene: (M p2 + Cp + K)X(p) = F (p)
(1.35)
Lo que se puede escribir como: Z(p)X(p) = F (p)
(1.36)
Z(p) es la matriz de rigidez din´amica. Invirtiendo la ecuaci´on 1.35 ´ o 1.36 se obtiene la matriz funci´on de transferencia H(p): X(p)
=
H(p)F (p)
H(p)
=
Z(p)−1 =
(1.37) adj(Z(p)) |Z(p)|
(1.38)
donde: adj(Z(p)) es la matriz adjunta de Z(p) y |Z(p)| es el determinante de Z(p). T
La matriz adjunta viene dada por adj(Z(p)) = [εij |Zij |] . Donde: |Zij | es el determinante de Z(p), eliminando la fila i y la columna j εij = 1, si i + j es par; = -1, si i + j es impar
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
1.2.2.
12
Polos, frecuencias naturales y factores de amortiguamiento
La ecuaci´ on caracter´ıstica del sistema viene dada por el denominador de la ecuaci´on 1.38, es decir, el determinante de Z(p). Al igual que en el caso de un grado de libertad, las ra´ıces de esta ecuaci´on o polos del sistema, definen las frecuencias naturales. Estas ra´ıces se pueden determinar resolviendo un problema de valores propios. Para transformar la ecuaci´on 1.35 en una ecuaci´on general de valores propios, se a˜ nade la siguiente identidad: (pM − pM ) {x} = 0
(1.39)
combinando las ecuaciones 1.35 y 1.39 se obtiene: (pA + B) {y} = {f 0 }
(1.40)
donde: A=
0 M
M −M , B= C 0
0 px 0 , {y} = y {f 0 } = K x f
(1.41)
Por lo tanto, los polos del sistema son los valores de p que satisfacen: |pA + B| = 0
(1.42)
Notar que las ra´ıces de esta ecuaci´on, tambi´en son las ra´ıces de la ecuaci´on |Z| = 0. Esta ecuaci´on genera 2n (n= n´ umero de grados de libertad) valores propios complejos, que aparecen en pares de complejos conjugados: λ1 .. . 0 λ n [Λ] = λ∗1 .. . 0
λ∗n
σ1 + jω1 .. . 0 σ + jω n n = σ1 − jω1 .. . 0
σn − jωn (1.43)
Como en el caso de un grado de libertad, la parte real del polo, σi , es el factor de amortiguamiento y la parte imaginaria, ωi , es la frecuencia natural amortiguada.
1.2.3.
Vectores modales, residuos
A cada valor propio del sistema le corresponde un vector propio. Para sistemas con multiples grados de libertad estos vectores propios introducen el concepto de
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
13
modos normales, formas modales o vectores modales, φi . Estos tambi´en aparecen en pares de complejos conjugados. Cada vector propio, se relaciona a un valor propio especifico. λ1 φ1 . . . λn φn λ∗1 φ∗1 . . . λ∗n φ∗n Φ= (1.44) φ1 ... φn φ∗1 ... φ∗n En general, los vectores propios contiene desplazamientos modales con valores complejos. Para el valor propio correspondiente, λi , se cumple que: (M λ2i + Cλi + K)φi = 0
(1.45)
De manera similar al caso de un grado de libertad, se puede introducir el concepto de residuos. Dado que λi , λ∗i son ra´ıces de la ecuaci´on caracter´ıstica del sistema, la ecuaci´ on 1.37 se puede escribir como: adj(Z(p)) E(p − λi )(p − λ∗i ) i=1
H(p) = Qn
(1.46)
donde E es una constante. Utilizando fracciones parciales: H(p) =
n X i=1
[A]∗i [A]i + (p − λi ) (p − λ∗i )
(1.47)
Los t´erminos [A]i y [A]∗i se denominan residuos, y vienen dados por: [A]i = Pi adj(Z(λi ))
(1.48)
Pi es una constante que depende del polo. Estudiando estos residuos se puede esclarecer la relaci´ on entre la matriz de funci´ on de transferencia y los vectores modales. Reescribiendo la ecuaci´on 1.38: Z(p)adj(Z(p)) = |Z(p)| [I]
(1.49)
Evaluando para p = λi , se obtiene: Z(λi )adj(Z(λi )) = 0
(1.50)
Si consideramos una columna arbitraria de adj(Z(λi )) se cumple que: Z(λi ) {adj(Z(λi ))}k = 0
(1.51)
Esta ecuaci´on es equivalente a la ecuaci´ on 1.45. Esto quiere decir la columnas de adj(Z(λi )) son proporcionales al i-´esimo vector modal. Considerando que la matriz
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
14
adj(Z(λi )) es sim´etrica, entonces sus filas tambi´en son proporcionales al i-´esimo vector modal. Es decir, la matrix adj(Z(λi )) es de la forma: φi,1 φi,1 φi,1 φi,2 . . . φi,1 φi,n φi,2 φi,1 φi,2 φi,2 . . . φi,2 φi,n adj(Z(λi )) = Ri φi φTi = Ri . (1.52) .. .. .. .. . . . φi,n φi,1
φi,n φi,2
...
φi,n φi,n
donde Ri es una constante asociada al vector φi . A partir de la ecuaci´on anterior, se tiene que cada residuo es de la forma: [A]i = Qi φi φTi
(1.53)
Y se obtiene finalmente que: n X Q∗ φ∗ φ∗T Qi φi φTi + i i ∗i H(p) = (p − λi ) (p − λi ) i=1
1.2.4.
(1.54)
Funci´ on de respuesta en frecuencia (FRF) y funci´ on de respuesta a un impulso (IRF)
La funci´on de respuesta en frecuencia (FRF) es la funci´on de transferencia evaluada en el eje de frecuencias (jω): n X Q∗i φ∗i φ∗T Qi φi φTi i H(jω) = + (jω − λi ) (jω − λ∗i ) i=1
(1.55)
La funci´on de respuesta a un impulso (IRF) viene dada por la transformada inversa de la funci´ on de respuesta en frecuencia:
h(t) =
N X
∗
λr t (Qr φr φTr eλr t + Q∗r φ∗r φ∗T ) r e
(1.56)
r=1
Se debe notar que, en el caso de multiples grados de libertad, para cada valor de ω, H(jω) es una matriz. Hik (jω) corresponde a funci´on de respuesta en frecuencia, cuando se excita la estructura en k y se mide la respuesta en i, o viceversa (H(jω) es sim´etrica):
Hik (jω) = Hki (jω) =
Xi (jω) Xk (jω) = Fk (jω) Fi (jω)
(1.57)
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
15
log(H11)
Las Figuras 1.8,1.9 y 1.10 muestran un ejemplo de las tres FRF para el sistema de dos grados de libertad de la Figura 1.7. Los “peaks” en las curvas son las frecuencias de resonancia, que tienen el mismo valor para todas las FRF, ya que son propiedades globales de la estructura. Por otro lado, las frecuencias donde las FRF tienen “ca´ıdas” se denominan antiresonancias. Los valores de las antiresonancias var´ıan para cada FRF, ya que son propiedades locales de la estructura.
w1
w2
frecuencia
log(H12)
Figura 1.8: Funci´ on de respuesta en frecuencia, excitaci´on en 1, respuesta en 1
w1
w2
frecuencia
Figura 1.9: Funci´on de respuesta en frecuencia, excitaci´on en 1, respuesta en 2 o excitaci´ ´ on en 2, respuesta en 1
16
log(H22)
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
w1
w2
frecuencia
Figura 1.10: Funci´ on de respuesta en frecuencia, excitaci´on en 2, respuesta en 2
Imag(H11)
Se debe recalcar que la funci´on de respuesta en frecuencia posee valores complejos, y se puede representar de distintas maneras: parte real e imaginaria vs frecuencia, amplitud (usualmente en escala logar´ıtmica) y fase vs frecuencia, y por u ´ltimo se puede graficar la parte real vs imaginaria con la frecuencia como un par´ametro ´ variable dentro de la curva. Este u ´ ltimo se denomina diagrama de Nyquist, en la Figura 1.11 se muestra un ejemplo.
Real(H11) Figura 1.11: Diagrama de Nyquist, H11
1.2.5.
Sistemas sin amortiguamiento y con amortiguamiento proporcional
En las secciones anteriores se vio es caso m´as general, en donde, el amortiguamiento es aproximado como un amortiguamiento viscoso general. Estos sistemas entregan polos de valores complejos, vectores modales de valores complejos con fases distintas para cada vector y funciones de respuesta en frecuencia complejas. A continuaci´on
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
17
se estudiaran dos casos particulares: sin amortiguamiento y con amortiguamiento proporcional. Si la matriz de amortiguamiento, C, es cero, el sistema de ecuaciones en el dominio de Laplace es: p2 M + K {x} = {f } (1.58) 2 En este caso, las ra´ıces del polinomio caracter´ıstico p M + K , se encuentran f´ acilmente resolviendo el siguiente problema de valores y vectores propios: p2 M + K {x} = 0 (1.59) La soluci´on a este problema entrega los valores propios, −ωi2 , y los vectores propios correspondientes, φi . Se debe notar que en este caso, los vectores propios tienen valores reales. Es por esto, que en el caso sin amortiguamiento, los vectores modales se denominan vectores modales normales o reales. La matriz de rigidez din´ amica en el dominio de frecuencias (p = jω), Z(jω) = −ω 2 M + K
(1.60)
tiene valores reales. Por lo tanto, su inversa, la matriz de frecuencia en respuesta H(jω), tambi´en tiene valores reales. La expresi´on general 1.55 se puede simplificar a: n X 2jωi Qi φi φTi H(jω) = (1.61) (ωi2 − ω 2 ) i=1 Notar que dado que H(jω) posee valores reales, entonces Qi es puramente imaginario. El caso de amortiguamiento proporcional es una forma hipot´etica de amortiguamiento. En su forma m´as sencilla, la matriz de amortiguamiento, C, toma la forma: C = αM + βK
(1.62)
donde α y β son constantes reales. En este caso el sistema de ecuaciones es: p2 M + p(αM + βK) + K {x}
=
{f }
(1.63)
(p2 + pα)M + (pβ + 1)K {x}
=
{f }
(1.64)
Considerando la versi´ on homog´enea de esta ecuaci´on y dividendo por (pβ + 1): 2 p + pα M + K {x} = 0 (1.65) pβ + 1
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
18
Esta ecuaci´ on es similar al caso sin amortiguamiento. Los sistemas con amortiguamiento proporcional tienen polos complejos, λi , que cumplen: λ2i + λi α = −ωi2 λi β + 1
(1.66)
y modos normales, iguales a los del caso sin amortiguamiento. La complejidad de los c´alculos para el caso de amortiguamiento proporcional es menor al caso de amortiguamiento general. Esta es la raz´ on principal para la introducci´on ´ del amortiguamiento proporcional. Este entrega un intermedio entre el caso sin amortiguamiento y el caso de amortiguamiento viscoso general. Con amortiguamiento proporcional, la funci´on de respuesta en frecuencia posee valores complejos y toma la forma: H(jω) =
n X i=1
1.2.6.
2jωi Qi φi φTi (σi2 + ωi2 − ω 2 ) − 2jωσi
(1.67)
Ortogonalidad y coordenadas modales
Los modos normales satisfacen un set de condiciones, conocidas como condiciones de ortogonalidad. Consideremos el problema de valores y vectores propios para el caso sin amortiguamiento o con amortiguamiento proporcional: −ωi2 M + K φi = 0 Consideremos ahora dos modos arbitrarios r y s: Kφr
=
ωr2 M φr
Kφs
=
ωs2 M φs
pre-multiplicando la primera ecuaci´on por φTs y la segunda ecuaci´on por φTr : φTs Kφr
= ωr2 φTs M φr
φTr Kφs
=
ωs2 φTr M φs
Transponiendo la segunda ecuaci´on y rest´andola a la primera, se tiene que: (ωr2 − ωs2 )φTs M φr = 0
(1.68)
Esta ecuaci´ on nos da dos posibilidades diferentes: φTs M φr = 0
r 6= s
(1.69)
φTs M φr 6= 0
r=s
(1.70)
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
19
La primera ecuaci´on nos dice que dos modos distintos son ortogonales con la matriz de masa (tambi´en con la matriz de rigidez) y que cuando se calcula el producto del mismo modo con la matriz de masa como en la segunda ecuaci´ on, el producto no es cero. El valor de este producto se denomina masa generalizada mi mi = φTi M φi
(1.71)
De manera similar se define la rigidez generalizada ki ki = φTi Kφi
(1.72)
Dado lo anterior, se tiene que: ΦT M Φ
= m
(1.73)
ΦT KΦ
= k
(1.74)
Φ es la matriz de modos normales, m y k son matrices diagonales denominadas matriz de masa modal y matriz de rigidez modal respectivamente. Debido a que la matriz de modos normales esta sujeta a una normalizaci´on arbitraria, los valores de ki y mi no son u ´ nicos. Lo que si se puede demostrar es que la raz´on ki /mi es u ´ nica y tiene un valor igual a ωi2 . Entre los distintos m´etodos de normalizaci´ on de los modos, ´el m´ as relevante en an´ alisis modal es la normalizaci´ on con respecto a la matriz de masa. Los modos normalizados con respecto a la matriz de masa tiene la propiedad: ΦT M Φ
= I
(1.75)
ΦT KΦ
=
(1.76)
Ω
donde Ω es una matriz diagonal que contiene los valores caracter´ısticos ωi2 en la diagonal principal. Introduciendo la transformaci´on modal x = Φq en la ecuaci´on 1.58 y premultiplicando por ΦT , se obtiene que: p2 ΦT M Φq + ΦT KΦq
=
ΦT {f }
(1.77)
(−ω 2 m + k)q
=
ΦT {f }
(1.78)
Para el caso de amortiguamiento proporcional, se define la matriz de amortiguamiento modal: c = ΦT CΦ
(1.79)
SISTEMAS CON MULTIPLES GRADOS DE LIBERTAD
20
Lo que da el siguiente sistema de ecuaciones desacopladas: (p2 m + pc + k)q = ΦT {f }
(1.80)
Cap´ıtulo 2
El M´ etodo de Elementos Finitos Fundamentalmente, existen dos tipos de m´etodos de elementos finitos; el m´etodo de las fuerzas, donde se asumen las fuerzas y se calculan los desplazamientos y el m´etodo de los desplazamientos, donde se asumen los desplazamientos y se calculan las fuerzas. Este u ´ ltimo es el m´as utilizado para an´alisis de vibraciones y ser´a el que se describir´ a a continuaci´on. El primer paso para construir un modelo en elementos finitos es discretizar la estructura en un n´ umero de elementos. Se pueden definir tres familias de elementos. (1) Elementos unidimensionales (linea), (2) Elementos bidimensionales (planos) y (3) elementos tridimensionales (s´olidos). Elementos bidimensionales
Elementos unidimensionales
Resorte, barra, viga, etc.
Membrana, placa, etc.
Elementos tridimensionales
Sólidos, fluidos, temperatura, etc.
Figura 2.1: Tipos de elementos finitos Para cada elemento se define u ˆ como los desplazamientos en los nodos y u(x), ε(x) y σ(x) como los desplazamientos, deformaciones y esfuerzos en el elemento, en 21
´ EL METODO DE ELEMENTOS FINITOS
22
funci´on de las coordenadas del elemento x. u ˆ y u se relacionan por la matriz de funciones de forma N (x), la que es caracter´ıstica para cada tipo de elemento. u(x) = N (x)ˆ u
(2.1)
la deformaci´ on se puede expresar como: ε(x) = Du(x) = DN (x) u ˆ = B(x)ˆ u | {z }
(2.2)
≡B(x)
D es el operador diferencial espacial, que para el caso general viene dado por:
∂ ∂x1
0 ∂ ∂x2
D= 0 0
0
0 0 ∂ ∂x3
∂ ∂x2 ∂ ∂x1
0
0 ∂ ∂x3 ∂ ∂x2
∂ T ∂x3
0
(2.3)
∂ ∂x1
Se puede demostrar que las matrices de rigidez y masa de un elemento viene dadas por: Z K
=
B T HBdV
(2.4)
ρN T N dV
(2.5)
V
Z M
= V
Donde V es el volumen del elemento, ρ es la densidad del material y H es una matriz del material. La tabla a continuaci´on entrega las expresiones de H para cada tipo de elemento. Problema
Matriz H
Barra
E
Viga
EI
Deformaciones planas
E 1−ν 2
1 ν 0
ν 1 0
Esfuerzos planos
E(1−ν) (1+ν)(1−2ν)
0 0
1−ν 2
1
ν (1−ν) 0
ν (1−ν)
1 0
0 0 1−2ν 2(1−ν)
ELEMENTO DE BARRA
23
Axisim´etrico
E(1−ν) (1+ν)(1−2ν)
ν (1−ν)
1
ν (1−ν) 0
1 0
Placa en flexi´ on
Tridimensional
ν (1−ν) ν (1−ν) E(1−ν) (1+ν)(1−2ν) 0 0 0
ν (1−ν)
ν (1−ν) ν (1−ν)
1 ν (1−ν)
1 0 0 0
0 0 0
Notaci´ on: E=Modulo de Young, ν = Coeficiente de Poisson h=Espesor de la placa, I=Momento de inercia
2.1.
Elemento de barra
Consideremos el elemento de barra de Figura 2.2: ûi
ûj A,E
i
j
x L
Figura 2.2: Elemento de barra L: largo del elemento A: area de la secci´ on E: Modulo de Young uˆi : desplazamiento en nodo i u(x): desplazamientos
0 1
1−ν 2
1
0
ν (1−ν) ν (1−ν)
0 0
ν 1 0
1 Eh3 ν 12(1−ν 2 ) 0
1−2ν 2(1−ν)
ν (1−ν)
ν (1−ν)
0 0
0 0 0 1−2ν 2(1−ν)
0 0
0 0 0 0 1−2ν 2(1−ν)
0 0 0 0 0
0
1−2ν 2(1−ν)
ELEMENTO DE BARRA
24
Para el caso de una barra se definen dos funciones de forma: Ni (x)
=
1−
Nj (x)
=
x L
x L
(2.6) (2.7)
Los desplazamientos se pueden escribir como: uˆi u(x) = Ni Nj = Nu ˆ uˆj
(2.8)
Por lo tanto, N
=
B
=
d N = −1/L 1/L dx
1 − x/L
x/L
(2.9) (2.10)
Por ultimo las matrices de rigidez y masa viene dada por: Z L EA 1 −1 −1/L Kl = E −1/L 1/L Adx = 1/L L −1 1 x=0 Z
L
Ml =
ρ x=0
1 − x/L 1 − x/L x/L
L/3 x/L Adx = ρA L/6
ûi
2.1.1.
L/6 L/3
(2.12)
ûj
Coordenadas locales i
(2.11)
yA,E globales j
x
Consideremos una barra como la mostrada en la Figura 2.3. Para determinar las L matrices de rigidez y masa en funci´on de las coordenadas globales, es necesario definir las matrices de rotaci´on. v
y u
y2
x x2 y1 α x1
Figura 2.3: Sistema de coordenadas locales y globales
ELEMENTO DE VIGA
25
De la Figura 2.3 se tiene que:
cos α u1 = 0 u2 | {z } |
sin α 0
0 cos α
0 sin α
{z
ql
}
R
|
x1 y1 x2 y2 {z qg
(2.13)
}
Dado que la energ´ıa es invariante respecto del sistema de referencia: 1 T q Kl ql 2 l
=
1 T q Kg qg 2 g
(2.14)
1 T q˙ Ml q˙l 2 l
=
1 T q˙ Mg q˙g 2 g
(2.15)
De donde se obtiene que las matrices con respecto al sistema global vienen dadas por:
2.2.
Kg
=
RT Kl R
(2.16)
Mg
=
R T Ml R
(2.17)
Elemento de viga
Consideremos elemento de viga como el mostrado en la Figura 2.4.
y vj, θj
vi, θi E,I i
j L Figura 2.4: Elemento de viga
L: Largo I: Momento de inercia de la secci´on E: Modulo de Young
x
ELEMENTO DE VIGA
26
v(x): Desplazamiento vertical θ(x): Rotaci´ on en el eje z Para el caso de un elemento de viga se definen las cuatro funciones de forma mostradas en la Figura 2.5.
ܰଵ = 1-
ܰଷ =
ଷ௫ మ మ
ଷ௫ మ మ
−
+
ଶ௫ య య
ܰଶ = x-
ଶ௫ య య
ܰସ = -
ଶ௫ మ
௫మ
௫య
+ మ
௫య
+ మ
Figura 2.5: Funciones de forma para un elemento de viga Por lo tanto, se tiene que: N (x)
=
B(x)
=
h
2
3
x x 1 − 3L 2 + 2 L3
2
x − 2 xL +
d N (x) = − L62 + 2 d x
12x L3
x3 L2
− L4 +
2
x 3L 2 − 6x L2
6 L2
x3 L3
−
2
− xL + 12x L3
x3 L2
− L2 +
i
6x L2
(2.18) (2.19)
ELEMENTO DE VIGA
27
De donde se obtiene finalmente que: vi Kl
=
θi
EI 12 6L L3 −12 6L
2.2.1.
=
θj
6L −12 6L 4L2 −6L 2L2 −6L 12 −6L 2L2 −6L 4L2
Ml
vj
156 ρA 22L 420 54 −13L
22L 54 4L2 13L 13L 156 −3L2 −22L
(2.20)
−13L −3L2 −22L 4L2
(2.21)
Coordenadas locales y globales
Consideremos una viga como la mostrada en la Figura 2.6. y
v
u θ
y2
x x2 y1 α x1
Figura 2.6: Sistema de coordenadas locales y globales De la Figura 2.6 se tiene que:
v1 − sin α θ1 0 v2 = 0 θ2 0 |
cos α 0 0 0
0 0 1 0 0 − sin α 0 0 {z R
0 0 cos α 0
0 0 0 1 }
x1 y1 θ1 x2 y2 θ2
(2.22)
ELEMENTO DE VIGA 3D
28
Las matrices con respecto al sistema global vienen dadas por:
2.3.
Kg
= RT Kl R
(2.23)
Mg
= R T Ml R
(2.24)
Elemento de viga 3D
Un elemento de viga 3D combina un elemento de viga, barra e incorpora adem´as torsi´ on. La Figura 2.7 muestra un esquema, cada nodo tiene 6 grados de libertad. y x z v2
v1 θv2
θv1 u1
E, A, Ix, Iy, Iz
θw1
u2 θw2
w1
θu1
L
w2
θu2
Figura 2.7: Elemento de viga 3D En este caso, las matrices de rigidez y masa del elemento de viga 3D en el sistema coordenadas local vienen dadas por:
ELEMENTO DE VIGA 3D
Kl
u1
v1
EA L
0
0 12EIz L3 0 0 0 0 EI 0 0 = 6EIy L3 0 L2 −EA 0 L 0 −12EIz L3 0 0 0 0 0 0 6EIz 0 L2
Ml
29
= ρAL
1 3
0 0 0 0 0 1 6
0 0 0 0 0
0 13 35
0 0 0 11L 210
0 9 70
0 0 0 −13L 210
w1 0 0 12EIy L3
0 −6EIy L2
0 0 0 −12EIy L3
0 −6EIy L2
0 0 0 13 35
0 −11L 210
0 0 0 9 70
0 13L 210
0
θu1
θv1
θw1
0 0 0
0 0
0
GIx L
0 0 0 0 0 −GIx L
0 0 0 0 0 r2 3
0 0 0 0 0 r2 6
0 0
−6EIy L2
0 4EIy L
0 0 0 6EIy L2
0 2EIy L2
0 0 0 −11L 210
0 L2 105
6EIy L2
0 0 0 4EIz L
0 −6EIz L2
0 0 0 2EIz L2
0 11L 210
0 0 0 2
0 0 0
L 105
−13 420
0 0 0
0 −L2 140
0
0 13 420
2
−L 140
u2
v2
−EA L
0
0 0 0 0 0 EA L
−12EIz L3
w2
0 0 0 −6EIz L2
0 0 0
0
0 6EIy L2
0 0 0 0 0
0 0 0
12EIy L3
−6EIz L2
0
1 6
0
0 0 0 0 0 1 3
0 0 0 0 0
9 70
0 0 0 13 420
0 13 35
0 0 0 −11L 210
θv2
θw2
0 0 0
0 0
0
0 0 −12EIy L3
12EIz L3
θu2
0 6EIy L2
0 0 9 70
0 −13 420
0 0 0 13 35
0 11L 210
0
−GIx L
0 0 0 0 0 GIx L
0 0 0 0 0 r2 6
0 0 0 0 0 r2 3
0 0
−6EIy L2
0 2EIy L2
0 0 0 6EIy L2
0 4EIy L
0 0 0 13L 210
0 −L2 140
0 0 0 11L 210
0 L2 105
0
6EIz L2
0 0 0 2EIz L2 0 −6EIz 2 L 0 0 0 4EIz L
0 −13L 210
0 0 0 −L2 140
0 −11L 210
0 0 0 L2 105
ELEMENTO DE VIGA 3D
2.3.1.
30
Rotaci´ on
A continuaci´on se construir´a el operador de rotaci´on, que permite expresar las matrices con respecto a un sistema global de coordenadas. Sean,
p1
p2
=
=
x1 y1 z1 x2 y2 z2
las posiciones de los nodos del elemento en el sistema global de coordenadas. − La direcci´ on del eje neutral del elemento, → e se define por: x
p2 − p1 p2 − p1 → − ex = = kp2 − p1 k L
(2.25)
− − ey y → ez , es necesario definir un tercer punto p3 que junto a p1 y p2 Para definir → definen el plano Oxz del elemento. Definiendo p3 se tiene que: (p3 − p1 ) × (p2 − p1 ) k(p3 − p1 ) × (p2 − p1 )k
→ − ey
=
→ − ez
− − = → ex × → ey
(2.26) (2.27)
× es el producto cruz entre dos vectores. Ahora se puede definir el operador de rotaci´ on que relaciona: xg = Rxl
(2.28)
donde xg es la posici´ on de los ejes globales y xl es la posici´on de los ejes locales: → T − − ex → ey → ey R= − (2.29) Se tiene que: u v = w θu θv = θw
x y R z θx θy R θz
(2.30)
(2.31)
ENSAMBLE
31
De donde, se puede expresar el vector de desplazamientos del elemento en los ejes locales y globales como: x1 u1 y1 v1 z w 1 1 θ θ x1 u1 θ θ R 0 0 0 y1 v1 θz1 0 R 0 0 θ w1 = (2.32) x2 0 0 R 0 u2 y2 0 0 0 R v2 | {z } z w 2 2 T θx2 θu2 θ θ y2 v2 θz2 θw2 Las matrices de rigidez y masa del elemento, expresadas en el sistema global de coordenadas quedan,
2.4.
Kg
=
T T Kl T
(2.33)
Mg
=
T T Ml T
(2.34)
Ensamble
La matrices de masa y rigidez de un sistema es un ensamble de las matrices de cada elemento. Para realizar el ensamble, necesitamos saber que grados de libertad est´an asociados que elemento en la estructura. Por ejemplo, consideremos el caso de una viga simple apoyada en ambos extremos que se muestra en la Figura 2.8. 1
2
3
4
5
Figura 2.8: Esta viga esta compuesta por 5 elementos de viga, cada uno con 4 grados de libertad. La estructura en total tiene 5 elementos, 6 nodos y 12 grados de libertad (2 grados por nodo). La matriz de rigidez de la viga completa es un ensamble de las 5 matrices de rigidez, ubicadas en los grados de libertad correspondiente. Esto es, la matriz K1 del primer elemento abarca los grados 1 a 4, la matriz K2 los
ENSAMBLE
32
elementos 3 a 6, etc. La Figura 2.9 ilustra este ensamble. La matriz de masa se construye de manera equivalente.
K1
0 K2
K=
K3 K4 0
K5
Figura 2.9: Ensamble matriz de rigidez La condici´on de borde, es que los grados de libertad de desplazamiento vertical en los nodos 1 y 6 est´an fijos. Para fijar grados de libertad, se eliminan las filas y columnas asociadas a estos grados en las matrices ensambladas. El siguiente c´odigo resuelve el problema en Matlab: %Propiedades de la viga E=3.2e10; %Pa modulo de elasticidad rho=2500; %densidad Kg/m3 A=0.05; %mˆ2 ´ area Iz=1.66e-4; %mˆ4 inercia secci´ on %Definici´ on de los elementos l=6; %largo total en metros n=5; %n´ umero de elementos le=l/n; %m largo de un elemento ne=(n+1)*2; %grados de libertad %Nodos con restricciones de libertad fix=[1 11]; %Nodos con restricciones %matriz de masa de un elemento Me=rho*A*le/420*[156 22*le 54 22*le 4*leˆ2 13*le 54 13*le 156 -13*le -3*leˆ2 -22*le
-13*le; -3*leˆ2; -22*le; 4*leˆ2];
ENSAMBLE
%matriz de rigidez de un elemento Ke=E*Iz/leˆ3*[12 6*le -12 6*le; 6*le 4*leˆ2 -6*le 2*leˆ2; -12 -6*le 12 -6*le; 6*le 2*leˆ2 -6*le 4*leˆ2]; %matrices iniciales K=zeros(ne,ne); M=zeros(ne,ne); %locel: grados de libertad asociados a cada elemento for i=1:n j=2*(i-1)+1; locel(i,1:4)=(j):(j+3); end %ensamble for i=1:n K(locel(i,:),locel(i,:))=K(locel(i,:),locel(i,:)) + Ke; M(locel(i,:),locel(i,:))=M(locel(i,:),locel(i,:)) + Me; end %fijar grados de libertad free=setdiff(1:ne,fix); K=K(free,free); M=M(free,free); %frecuencias naturales y modos [Phif,W]=eig(K,M); w=sqrt(diag(W))/2/pi; Phi=zeros(ne,ne-length(fix)); Phi(free,:)=Phif; %Gr´ afica los primeros tres modos for i=1:3 figure set(gcf,’Position’,[100 100 500 400]) set(gca,’FontSize’,18) set(gca,’Units’,’centimeters’) set(gca,’Position’,[1.0 1.5 11 8]) set(gcf, ’PaperUnits’,’inches’); plot(0:ne/2-1,Phi(1:2:ne,i),’-b’,’LineWidth’,3)
33
ENSAMBLE
34
xlim([0 ne/2-1]) set(gca,’YTick’,[]); set(gca,’XTick’,[0:5]); set(gca,’XTicklabel’,[1:6]); xlabel(’Nodo’,’FontSize’,16) title([’Modo ’ num2str(i) ’, \omega= ’ num2str(w(i),3) ’Hz’],’FontSize’,16) grid on; end
La Figura 2.10 muestra los primeros tres modos de la viga. Modo 1, ω= 9Hz
1
2
3
4 Nodo
Modo 3, ω= 81.6Hz
5
6
1
2
3
4
Modo 2, ω= 36Hz
5
6
Nodo
Figura 2.10: Modos de la viga
1
2
3
4 Nodo
5
6
Cap´ıtulo 3
Integraci´ on num´ erica de las ecuaciones de movimiento Para un problema general, la ecuaci´on de movimiento viene dada por: Mx ¨ + C x˙ + Kx = F
(3.1)
Donde M, C y K son las matrices de masa, amortiguaci´on y rigidez respectivamente. F es un vector de fuerzas externas, y x ¨, x, ˙ x son vectores de aceleraci´ on, velocidad y desplazamiento respectivamente. Matem´aticamente, la ecuaci´on 3.1 representa un sistema de ecuaciones diferenciales lineales de segundo orden, las que en un principio, se pueden resolver por procedimientos est´ andar para ecuaciones diferenciales ordinarias. Sin embargo, estos procedimientos pueden ser muy costosos (en t´erminos de tiempo y cantidad de c´ alculos necesarios) si las matrices son grandes. En el an´alisis de elementos finitos, se utilizan algunos algoritmos eficientes para resolver las ecuaciones de movimiento, los que veremos en las secciones siguientes. Los m´etodos que veremos se dividen en dos categor´ıas: integraci´ on directa y superposici´ on modal. Aunque ambas t´ecnicas puedan parecer muy distintas a primera vista, est´an fuertemente relacionadas, y la selecci´on de una u otra depende de su eficiencia para calcular la soluci´ on num´erica.
3.1.
M´ etodos de integraci´ on directa
En la integraci´on directa las ecuaciones de movimiento son integradas por un procedimiento num´erico por pasos, el t´ermino “directa” significa que no es necesario 35
´ ´ DIRECTA METODOS DE INTEGRACION
36
hacer transformaciones de las ecuaciones previo a la integraci´on num´erica. En los procedimientos de integraci´on temporal, se asume que los vectores de aceleraci´on, velocidad y desplazamiento en el tiempo t = 0 ,¨ x(0), x(0), ˙ x(0), son conocidos. La soluci´ on se divide en pasos de tiempo ∆t.
3.1.1.
El m´ etodo de las diferencias centrales
El m´etodo de las diferencias centrales utiliza la aproximaci´on: x ¨(t) =
1 (x(t − ∆t) − 2x(t) + x(t + ∆t)) ∆t2
(3.2)
El error en la expansi´ on de 3.2 es del orden de ∆t2 , para tener el mismo orden de error en la expansi´ on de la velocidad se utiliza: x(t) ˙ =
1 (x(t + ∆t) − x(t − ∆t)) 2∆t
(3.3)
El desplazamiento en el tiempo t + ∆t se obtiene al considerar la ecuaci´on de movimiento en t: Mx ¨(t) + C x(t) ˙ + Kx(t) = F (t)
(3.4)
Substituyendo, se obtiene:
1 1 2 1 1 M + C x(t+∆t) = F (t)− K − M x(t)− M − C x(t−∆t) ∆t2 2∆t ∆t2 ∆t2 2∆t (3.5)
Se debe notar que para calcular x(t + ∆t), son necesarios x(t) y x(t − ∆t). Por lo tanto, se debe utilizar un procedimiento especial de inicializaci´on. Dado que se conocen x ¨(0), x(0), ˙ x(0), las ecuaciones 3.2 y 3.3 se pueden utilizar para calcular x(−∆t): x(−∆t) = x(0) − ∆tx(0) ˙ +
∆t2 x ¨(0) 2
(3.6)
En general, para llegar a soluciones estables se requiere un paso de tiempo relativamente peque˜ no. De hecho, en el m´etodo de las diferencias centrales es paso de tiempo debe ser menor que un valor cr´ıtico ∆tcr . El paso de tiempo cr´ıtico se determina como: ∆tcr =
Tn π
(3.7)
´ ´ DIRECTA METODOS DE INTEGRACION
37
donde Tn es el menor periodo en el sistema, entre el inverso de la mayor frecuencia natural y la frecuencia de las fuerzas externas. La Tabla 3.1 ilustra el procedimiento paso a paso para resolver las ecuaciones de movimiento mediante el m´etodo de las diferencias centrales. Tabla 3.1: Procedimiento paso a paso, m´etodo de diferencias centrales C´ alculos iniciales: 1. Crear matrices de rigidez K, masa M y amortiguaci´on C. 2. Inicializar x ¨(0), x(0), ˙ x(0). 3. Seleccionar paso de tiempo ∆t, ∆t ≤ ∆tcr . Calcular constantes de integraci´on. a0 =
1 ; ∆t2
a1 =
1 ; 2∆t
a2 = 2a0 ;
a3 =
1 a2
4. Calcular x(−∆t) = x(0) − ∆tx(0) ˙ + a3 x ¨(0) ˆ = a0 M + a1 C. 5. Formar matriz de masa efectiva M Para cada paso de tiempo: 1. Calcular fuerzas efectivas en t: Fˆ (t) = F (t) − (K − a2 M )x(t) − (a0 M − a1 C)x(t − ∆t) 2. Determinar desplazamientos en t + ∆t: ˆ −1 Fˆ (t) x(t + ∆t) = M 3. Si se requiere, evaluar las aceleraciones y velocidades en t: x ¨(t) = a0 (x(t − ∆t) − 2x(t) + x(t + ∆t)) x(t) ˙ = a1 (x(t + ∆t) − x(t − ∆t))
3.1.2.
El m´ etodo de Wilson θ
En el m´etodo de Wilson θ se asume una variaci´on lineal de la aceleraci´on entre un paso de tiempo t y t + θ∆t, donde θ ≥ 1, como se ilustra en la Figura 3.1.
´ ´ DIRECTA METODOS DE INTEGRACION
38
Si θ = 1, el m´etodo se reduce a un esquema de aceleraci´on lineal. Pero se puede demostrar que para asegurar estabilidad incondicional es necesario utilizar θ ≥ 1,37, usualmente se emplea θ = 1,4
ݔሷ ( ݐ+ θ∆)ݐ
ݔሷ ( ݐ+ ∆)ݐ
ݔሷ ()ݐ
t + ∆t
t
t + θ∆t
Figura 3.1: Supuesto de aceleraci´on lineal, m´etodo de Wilson θ Denotemos τ como el incremento en el tiempo, donde 0 ≤ τ ≤ θ∆t. Entonces en el intervalo t a t + θ∆t se tiene que: x ¨(t + τ ) = x ¨(t) +
τ (¨ x(t + θ∆t) − x ¨(t)) θ∆t
(3.8)
Integrando la ecuaci´ on anterior se obtiene que, x(t ˙ + τ ) = x(t) ˙ + τx ¨(t) +
τ2 (¨ x(t + θ∆t) − x ¨(t)) 2θ∆t
(3.9)
y 1 τ3 x(t + τ ) = x(t) + τ x(t) ˙ + τ 2x ¨(t) + (¨ x(t + θ∆t) − x ¨(t)) 2 6θ∆t
(3.10)
Usando 3.9 y 3.10, se tiene en el tiempo t + θ∆t, θ∆t (¨ x(t + θ∆t) + x ¨(t)) 2
x(t ˙ + θ∆t)
=
x(t) ˙ +
x(t + θ∆t)
=
x(t) + θ∆tx(t) ˙ +
θ2 ∆t2 (¨ x(t + θ∆t) + 2¨ x(t)) 6
(3.11) (3.12)
De donde se puede expresar x ¨(t + θ∆t) y x(t ˙ + θ∆t) en t´erminos de x(t + θ∆t): x ¨(t + θ∆t)
=
6 6 (x(t + θ∆t) − x(t)) − x(t) ˙ − 2¨ x(t) θ2 ∆t2 θ∆t
(3.13)
x(t ˙ + θ∆t)
=
θ∆t 3 (x(t + θ∆t) − x(t)) − 2x(t) ˙ − x ¨(t) θ∆t 2
(3.14)
´ ´ DIRECTA METODOS DE INTEGRACION
39
Para obtener la soluci´ on de los desplazamientos, velocidades y aceleraciones en el tiempo t + ∆t, se considera la ecuaci´ on de movimiento en el tiempo t + θ∆t. Sin embargo, dado que la aceleraci´on varia linealmente, se utiliza un vector de fuerzas extrapolado linealmente, i e., la ecuaci´on utilizada es: Mx ¨(t + θ∆t) + C x(t ˙ + θ∆t) + Kx(t + θ∆t) = F¯ (t + θ∆t)
(3.15)
F¯ (t + θ∆t) = F (t) + θ(F (t + ∆t) − F (t))
(3.16)
Sustituyendo las ecuaciones 3.13 y 3.14 en 3.15, se obtiene una ecuaci´on en donde se puede despejar x(t + θ∆t). Luego, sustituyendo x(t + θ∆t) en 3.13, se obtiene x ¨(t + θ∆t) la que es usada en las ecuaciones 3.8 a 3.10, todas evaluadas en τ = ∆t para calcular x(t + ∆t), x(t ˙ + ∆t) y x ¨(t + ∆t). El algoritmo completo se da en la Tabla 3.2 Se debe notar que con este m´etodo no se requiere inicializaciones especiales, dado que los desplazamientos, velocidades y aceleraciones en el tiempo t + ∆t se expresan en t´erminos de las mismas cantidades en el paso de tiempo anterior t. Tabla 3.2: Procedimiento paso a paso, m´etodo de Wilson θ C´ alculos iniciales: 1. Crear matrices de rigidez K, masa M y amortiguaci´on C. 2. Inicializar x ¨(0), x(0), ˙ x(0). 3. Seleccionar paso de tiempo ∆t y θ. Calcular constantes de integraci´on. a0 =
6 ; (θ∆t)2
a5 =
−a2 ; θ
a1 =
3 ; θ∆t
3 a6 = 1 − ; θ
a2 = 2a1 ;
a7 =
∆t ; 2
a3 =
a8 =
θ∆t ; 2
a4 =
a0 θ
∆t2 6
ˆ = K + a0 M + a1 C. 4. Formar matriz de rigidez efectiva K Para cada paso de tiempo: 1. Calcular fuerzas efectivas en t + θ∆t: Fˆ (t + θ∆t) = F (t) + θ(F (t + ∆t) − F (t)) + M (a0 x(t) + a2 x(t) ˙ + 2¨ x(t)) +C(a1 x(t) + 2x(t) ˙ + a3 x ¨(t))
´ ´ DIRECTA METODOS DE INTEGRACION
40
2. Determinar desplazamientos en t + θ∆t: ˆ −1 Fˆ (t + θ∆t) x(t + θ∆t) = K 3. Evaluar aceleraciones, velocidades y desplazamientos en t + ∆t: x ¨(t + ∆t) = a4 (x(t + θ∆t) − x(t)) + a5 x(t) ˙ + a6 x ¨(t) x(t ˙ + ∆t) = x(t) ˙ + a7 (¨ x(t + ∆t) + x ¨(t)) x(t + ∆t) = x(t) + ∆tx(t) ˙ + a8 (¨ x(t + ∆t) + 2¨ x(t))
3.1.3.
El m´ etodo de Newmark
El m´etodo de Newmark tambi´en se puede relacionar con un m´etodo de aceleraciones lineales. Utiliza los siguientes supuestos: x(t ˙ + ∆t)
=
x(t + ∆t)
=
x(t) ˙ + [(1 − δ)¨ x(t) + δ x ¨(t + ∆t)] ∆t 1 x(t) + x(t)∆t ˙ + −α x ¨(t) + α¨ x (t + ∆t) ∆t2 2
(3.17) (3.18)
donde α y δ son par´ametros que determinan la precisi´on y estabilidad del m´etodo de integraci´on. Cuando δ = 12 y α = 16 , las ecuaciones 3.17 y 3.18 corresponden al m´etodo de aceleraciones lineales (que tambi´en se obtiene con θ = 1 en el m´etodo de Wilson θ). Newmark propuso originalmente como condici´on de estabilidad incondicional el esquema de promedio constante de la aceleraci´on (tambi´en denominado regla trapezoidal), en donde δ = 12 y α = 14 (ver Figura 3.2).
ݔሷ ( ݐ+ ∆)ݐ ݔሷ ()ݐ
t
1 ݔሷ ݐ+ ݔሷ ( ݐ+ ∆)ݐ 2
t + ∆t
Figura 3.2: Esquema de promedio de aceleraci´on constante de Newmark
´ ´ DIRECTA METODOS DE INTEGRACION
41
Adicionalmente a las ecuaciones 3.17 y 3.18, para la soluci´ on de los desplazamientos, velocidades y aceleraciones en t + ∆t, se utilizan las ecuaciones de movimiento en t + ∆t: Mx ¨(t + ∆t) + C x(t ˙ + ∆t) + Kx(t + ∆t) = F (t + ∆t)
(3.19)
De la ecuaci´on 3.18 se puede despejar x ¨(t + ∆t) en funci´on de x(t + ∆t) y despu´es se puede reemplazar en la ecuaci´on 3.17. De esta manera se obtienen x ¨(t + ∆t) y x(t ˙ + ∆t) en funci´ on del desplazamiento x(t + ∆t). Sustituyendo estas dos expresiones en la ecuaci´on 3.19, se puede despejar x(t + ∆t), de donde tambi´en se tienen x ¨(t + ∆t) y x(t ˙ + ∆t). El algoritmo completo utilizando el algoritmo de Newmark se ilustra en la tabla 3.3. Se debe notar la cercana relaci´on entre el m´etodo de Wilson θ y el m´etodo de Newmark. Tabla 3.3: Procedimiento paso a paso, m´etodo de Newmark C´ alculos iniciales: 1. Crear matrices de rigidez K, masa M y amortiguaci´on C. 2. Inicializar x ¨(0), x(0), ˙ x(0). 3. Seleccionar ∆t y los par´ametros α y δ. Calcular constantes de integraci´on. δ ≥ 0,50;
α ≥ 0,25(0,5 + δ)2
1 1 1 δ ; a1 = ; a2 = ; a3 = − 1; α∆t2 α∆t α∆t 2α ∆t δ −2 ; a6 = ∆t(1 − δ); a7 = δ∆t a5 = 2 α
a0 =
ˆ = K + a0 M + a1 C. 4. Formar matriz de rigidez efectiva K Para cada paso de tiempo: 1. Calcular fuerzas efectivas en t + θ∆t: Fˆ (t + ∆t) = F (t + ∆t) + M (a0 x(t) + a2 x(t) ˙ + a3 x ¨(t)) +C(a1 x(t) + a4 x(t) ˙ + a5 x ¨(t)) 2. Determinar desplazamientos en t + ∆t:
a4 =
δ −1 α
´ ´ DIRECTA METODOS DE INTEGRACION
42
ˆ −1 Fˆ (t + ∆t) x(t + ∆t) = K 3. Evaluar aceleraciones y velocidades en t + ∆t: x ¨(t + ∆t) = a0 (x(t + ∆t) − x(t)) − a2 x(t) ˙ − a3 x ¨(t) x(t ˙ + ∆t) = x(t) ˙ + a6 x ¨(t) + a7 x ¨(t + ∆t)
3.1.4.
Ejemplo
Consideremos un sistema simple, cuyas ecuaciones de movimiento son las siguientes:
2 0
0 1
x ¨1 x ¨2
0,3 + −0,1
−0,1 0,2
x˙ 1 x˙ 2
6 + −2
−2 4
x1 x2
=
0 10
(3.20)
Las frecuencias naturales del este sistema son 0.23Hz y 0.36Hz, por lo tanto el menor periodo del sistema es 2.78s. El paso de tiempo cr´ıtico es ∆tcr = 2,78 π = 0,88s. Utilizando un paso de tiempo 0.18s se obtienen, de acuerdo a los algoritmos descritos en las secciones anteriores, los resultados mostrados en la Figura 3.3. Los tres algoritmos de integraci´ on dan resultados similares, aunque el m´etodo de las diferencias centrales es que m´etodo que m´ as se acerca a la soluci´on exacta. Se asumieron desplazamientos y velocidades iniciales nulas. Dif. Centrales
Wilson
Newmark
Exacta
3 2.5 2
x1(t)
1.5 1 0.5 0 -0.5 -1 0
1
2 3 Tiempo (s)
4
Figura 3.3: Desplazamientos obtenidos
5
´ MODAL SUPERPOSICION
3.2.
43
Superposici´ on modal
En los m´etodos de integraci´ on directa se deben realizar muchas operaciones en cada paso de tiempo. Si se debe integrar en un numero alto de pasos de tiempo, puede ser m´ as conveniente transformar las ecuaciones de movimiento en una forma que requiera de menos operaciones para cada paso de tiempo. Aprovechando las propiedades ortogonales de los modos propios, se puede definir la transformaci´ on: x(t) = Φy(t)
(3.21)
De donde se obtiene la ecuaci´on de movimiento: ΦT M Φ¨ y + P hiT CΦy˙ + P hiT KΦy = ΦT F
(3.22)
Las matrices ΦT M Φ, ΦT CΦ y ΦT KΦ son matrices diagonales, por lo que el sistema de ecuaciones 3.22 es un sistema de ecuaciones desacoplado. En consecuencia, se tienen n ecuaciones individuales de la forma:
mi y¨i + ci y˙i + ki yi = ri ri = φTi F (t)
i = 1, 2, 3, . . . , n
(3.23)
La soluci´ on para cada una de las ecuaciones se puede obtener utilizando los algoritmos de integraci´on antes descritos o de acuerdo al procedimiento descrito en la secci´ on 1.1.3. Una vez calculada la soluci´on para cada una de las ecuaciones. La soluci´on general se obtiene al superponer las soluciones modales: x(t) =
n X
φi yi (t)
(3.24)
i=1
3.2.1.
Ejemplo
Consideremos el sistema de ecuaciones 3.20. Para este sistema los modos normales son: −0,5774 −0,4082 Φ= (3.25) −0,5774 0,8165 Reemplazando Φ en ecuaci´on 3.22, se obtiene: 1 0 y¨1 0,1 0,0 y˙ 1 2 0 y1 −5,7735 + + = (3.26) 0 1 y¨2 0,0 0,25 y˙ 2 0 5 y2 8,1650
´ MODAL SUPERPOSICION
44
Este sistema de ecuaciones se puede resolver utilizando alg´ un algoritmo de integraci´ on. Utilizando diferencias centrales con el mismo paso de tiempo del ejemplo anterior, se obtiene la soluci´on para y(t). La soluci´on, x(t), viene dada por: x(t) =
−0,5774 −0,5774
y1 (t) +
−0,4082 0,8165
y2 (t)
(3.27)
La soluci´on obtenida se muestra en la Figura 3.4. Como es de esperarse, la soluci´on es la misma que si se utiliza diferencias modales con o sin transformaci´on modal de las ecuaciones. La ventaja de la transformaci´on modal es que se reduce el tiempo requerido para calcular la soluci´on a cada paso de tiempo. 3 2.5 2
1
x (t)
1.5 1 0.5 0 -0.5
Superposicion Modal Dif. Centrales
-1 0
1
2 3 Tiempo (s)
4
Figura 3.4: Desplazamientos obtenidos
5
Cap´ıtulo 4
Procesamiento de se˜ nales El procesamiento de digital se˜ nales es una herramienta muy importante en el an´alisis de sistemas. En esta secci´on se dar´a un resumen de los conceptos b´asicos asociados al procesamiento de se˜ nales. La se˜ nales, en general, se pueden clasificar de acuerdo a la tabla 4.1. Las se˜ nales estacionarias son aquellas cuyas propiedades promedio no var´ıan con el tiempo, pueden ser o deterministicas o aleatorias. El grupo m´as importante de se˜ nales deterministicas son las se˜ nales peri´odicas. Una funci´ on pseudo-aleatoria es una se˜ nal aleatoria que se repite con un cierto periodo. Las se˜ nales no estacionarias se pueden dividir en continuas o transientes. Las se˜ nales transientes se pueden definir como se˜ nales que comienzan y terminar en cero en el periodo de observaci´ on. Se˜ nales Estacionarias No estacionarias Determin´ıstica Aleatoria Continua Transiente Peri´ odica Cuasi-peri´ odica Pseudo-aleatoria Tabla 4.1: Tipos de se˜ nales Dado que el objetivo del procesamiento de se˜ nales es extraer el m´ aximo de informaci´on de las se˜ nales, es en general beneficioso estudiar las se˜ nales en distintos dominios. Las se˜ nales medidas son, claramente, funciones en el dominio del tiempo. Para estudiar su contenido en frecuencias, es m´as f´acil examinar las se˜ nales en el
45
LA TRANSFORMADA DE FOURIER
46
dominio de frecuencias. La transformada (inversa) de Fourier, permite transformar de manera sensilla una se˜ nal en el tiempo a una se˜ nal en frecuencia y viceversa. En consecuencia, la transformada de Fourier es uno de los temas m´as importante en procesamiento de se˜ nales.
4.1.
La transformada de Fourier
J.B. Fourier prob´o que una funci´on peri´odica en el tiempo se puede representar como una suma de componentes sinusoidales a frecuencias equiespaciadas:
g(t) =
+∞ X
G(k∆f )ej2πk∆f t
(4.1)
−∞
Los coeficientes de fourier viene dados por: Z 1 +T /2 G(k∆f ) = g(t)e−j2πk∆f t dt T −T /2
(4.2)
con: t: tiempo k: entero que cuenta los pasos en frecuencia ∆f : espaciado de frecuencias o resoluci´on de frecuencia: ∆f = 1/T √ j = −1 T: periodo de tiempo: T = 1/∆f El set de valores G(k∆f ) se denomina espectro de la funci´ on g(t). En general el espectro posee valores complejos. Al utilizar computadores digitales, es necesario adquirir la se˜ nal continua en intervalos de tiempo. Esto significa que la se˜ nal continua es representada por una se˜ nal discreta con valores a tiempos equidistantes. Considerando esto la transformada de Fourier queda como: Z 1 +fs /2 G(f )ej2πf n∆t df (4.3) g(n∆t) = fs −fs /2 G(f ) =
+∞ X −∞
con:
g(n∆t)e−j2πf n∆t
(4.4)
LA TRANSFORMADA DE FOURIER
47
n: entero contando el numero de pasos de tiempo ∆t: intervalo de muestreo: ∆t = 1/fs fs : frecuencia de muestreo: fs = 1/∆t En condiciones reales de medici´on experimental es imposible medir la se˜ nal temporal hasta un tiempo infinito. Una parte de la se˜ nal debe ser seleccionada. Se asume que la se˜ nal capturada se repite con un periodo T, entregando una funci´on peri´odica. Combinando la hip´ otesis de periodicidad con un muestreo temporal de la se˜ nal, se obtiene la definici´ on de la transformada discreta de Fourier:
g(n∆t) =
Ns −1 1 X G(k∆f )ej2πnk/Ns fs
(4.5)
k=0
G(k∆f ) =
Ns −1 1 X g(n∆t)e−j2πnk/Ns Ns n=0
(4.6)
con: Ns : n´ umero de datos: T = Ns ∆t y fs = Ns ∆f La evaluaci´on directa de la transformada discreta de Fourier requiere Ns2 operaciones. Con la transformada r´apida de Fourier (FFT) se reduce el n´ umero de operaciones a Ns log2 Ns . La transformada r´apida de Fourier es el n´ ucleo de todos los procesadores de se˜ nal modernos. Consideremos como ejemplo una se˜ nal sinusoidal, con periodo T=0.02s y amplitud A=1: x(t) = sen(2π50t)
(4.7)
La se˜ nal es muestreada a una frecuencia de 10kHz y se considera el periodo de tiempo 0-0.2s. En la Figura 4.1 se ilustra el resultado de la transformada r´apida de Fourier de la se˜ nal.
ERRORES Y VENTANAS
48
1
1.5
0.5 Amplitud
Amplitud
1 0
0.5 -0.5
-1 0
0.05
0.1 Tiempo (s)
0.15
0.2
0 0
20
40 60 Frecuencia (Hz)
80
100
Figura 4.1: Espectro de una se˜ nal peri´odica
4.1.1.
Algunos par´ ametros importantes
Aqu´ı se resumen algunos de los par´ametros mencionados en los p´arrafos anteriores: T : Periodo de tiempo en donde se adquieren Ns muestras equiespaciadas de la se˜ nal a ser analizada. En la mayor´ıa de los analizadores de Fourier Ns est´a restringidos a potencias de 2 (por ejemplo, 1024). Relaciones: T = Ns ∆t = 1/∆f fs : Frecuencia de muestreo, frecuencia a la cual se adquieren y digitalizan los datos. Relaciones: fs = 1/∆t = Ns ∆f ∆t : Intervalo de muestreo: intervalo de tiempo al cual la se˜ nal es muestreada. Relaciones: ∆t = T /Ns = 1/fs ∆f : Espaciado de frecuencias en el espectro. Relaciones: ∆f = 1/T = fs /Ns . En consecuencia, mediciones con un espaciado de frecuencias peque˜ no (i.e. una resoluci´ on en frecuencias alta) conllevan tiempos de medici´on mayores. fmax : Frecuencia m´axima: frecuencia mayor contenida o permitida en la se˜ nal temporal. Seg´ un el teorema de Shannon: fmax ≤ fs /2 Ns : N´ umero de muestras en el periodo de tiempo T: Ns = T /∆t = fs /∆f
4.2.
Errores y ventanas
Durante el proceso de an´ alisis digital de una se˜ nal pueden ocurrir muchos errores. Errores t´ıpicos son sobrecargas, ruido digital, errores de cuantificaci´on, limitaciones del rango din´amico. Sin embargo, los dos errores principales son aliasing y leakage.
ERRORES Y VENTANAS
4.2.1.
49
Aliasing
El aliasing se produce debido al hecho que la se˜ nal temporal debe ser muestreada. Componentes de alta frecuencia en la se˜ nal pueden causar errores de amplitud y frecuencia en el espectro. Si la mayor frecuencia contenida en una se˜ nal no cumple con el teorema de Shannon: fmax ≤ fs /2, entonces las frecuencias por sobre fs /2 van a aparecer como frecuencias menores a fs /2. La Figura 4.2 muestra un ejemplo de Aliasing, se muestran tres senos con frecuencias de 1, 4 y 6 Hz muestreados a 5 Hz, a la frecuencia de muestreo los tres senos son id´enticos. 1 1 Hz 4 Hz 6 Hz
Amplitud
0.5
0
-0.5
-1 0
0.2
0.4 0.6 Tiempo (s)
0.8
1
Figura 4.2: Aliasing: Seno a 1, 4 y 6 Hz, muestreados a 5 Hz (c´ırculos) Aliasing se puede evitar removiendo todos los componentes con frecuencias mayores a fs /2. Esto se puede lograr con una se˜ nal de excitaci´ on apropiada, aunque se logra generalmente utilizando un filtro pasa bajas. Dado que no existen filtros que remuevan todas las frecuencias altas a cero sin influenciar en las bajas, los filtros se fijan normalmente a un 40 % de fs .
4.2.2.
Leakage
El Leakage se origina debido a que los datos deben ser adquiridos en un periodo de observaci´on finito T . La transformada discreta de Fourier asume entonces que la se˜ nal es peri´ odica con periodo T. Si esta condici´on no se cumple, se produce un error de “leakage”. La Figura 4.3 ilustra el espectro obtenido de una se˜ nal tipo coseno, cuando la funci´on es peri´ odica en T y cuando no lo es. En el segundo caso,
ERRORES Y VENTANAS
50
el espectro discreto no coincide con el real. El error en la hip´ otesis de periodicidad produce errores importantes de amplitud y frecuencia. 1.5
Amplitud
Amplitud
1
0.5
0
1
2
0 0
3
2
4 6 Frecuencia (Hz)
8
10
2
4 6 Frecuencia (Hz)
8
10
Tiempo (s)
1.5
Amplitud
Amplitud
1
0.5
0
0.8 1.6 Tiempo (s)
2.4
0 0
Figura 4.3: Hip´otesis de periodicidad, Leakage La u ´nica soluci´on al problema de leakage, es asegurarse que la se˜ nal es peri´odica o se observa completamente en el periodo de adquisici´ on, lo que en general, es muy dif´ıcil de lograr. En sistemas perfectamente lineales, se puede lograr al excitarlos con una se˜ nal peri´odica en el intervalo de tiempo considerado. Aumentar el tiempo de adquisici´ on, i.e. aumentar la resoluci´on en frecuencias, ayuda a mejorar la periodicidad de la se˜ nal. El uso de ventanas de tiempo tambi´en ofrecen una soluci´on parcial al problema de leakage.
4.2.3.
Ventanas
El uso de ventanas de tiempo no se puede evitar en el procesamiento digital de una se˜ nal. Al medir una se˜ nal temporal, solo una parte de la se˜ nal total es considerada.
ERRORES Y VENTANAS
51
Esto equivale a multiplicar la se˜ nal actual con una ventana de tiempo rectangular (Figuras 4.4(a) y 4.5(a)). Sin embargo, una mejor selecci´on de la ventana puede reducir considerablemente el error debido a leakage. En general, se buscan ventanas que reduzcan las discontinuidades en los extremos de la se˜ nal, dado que reducen el error por leakage al forzar la se˜ nal a ser peri´odica. La selecci´ on de una ventana de tiempo, es siempre un compromiso entre una buena estimaci´ on de la amplitud y una buena resoluci´on espectral. La figura 4.4 muestra un conjunto de ventanas de tiempo que son usualmente empleadas.
0
Hamming
Amplitud
Amplitud
Hanning
Amplitud
Rectangular
0 Tiempo (s)
0 Tiempo (s)
(a) Rectangular
Tiempo (s)
(b) Hanning
(c) Hamming Flat-top
Amplitud
Amplitud
Gaussian
0
0 Tiempo (s)
(d) Gaussian
Tiempo (s)
(e) Flat-top
Figura 4.4: Diferentes ventanas de tiempo En la Figura 4.5 se muestra el espectro obtenido de una se˜ nal peri´odica utilizando las distintas ventanas de tiempo. La se˜ nal esta compuesta por un coseno de frecuencia f1 y amplitud 1, y un coseno de frecuencia f2 y amplitud 0.01. El coseno de frecuencia f1 no es peri´odico en el intervalo de tiempo considerado. Se observa que solo algunas ventanas son capaces de separar bien ambas se˜ nales. Es importante notar, que dado que la se˜ nal no es peri´odica en el periodo considerado, la amplitud en el dominio de frecuencias depende de la ventana y de la ubicaci´on de la se˜ nal no-peri´odica con respecto a los puntos de frecuencias muestreadas. En consecuencia, este error no puede ser corregido por un factor general. Tambi´en es importante destacar, que el uso de ventanas no rectangulares, reduce el total de energ´ıa en la se˜ nal. Lo que conlleva a una reducci´on en la amplitud mostrada en el dominio de frecuencias. La cantidad de energ´ıa solo depende de la forma de la ventana y por lo tanto, este error puede ser compensado.
52
Amplitud
log(Amplitud)
ERRORES Y VENTANAS
f1 f2
Tiempo (s)
Frecuencia (Hz)
Amplitud
log(Amplitud)
(a) Rectangular
f1 f2
Tiempo (s)
Frecuencia (Hz)
Amplitud
log(Amplitud)
(b) Hanning
f1 f2
Tiempo (s)
Frecuencia (Hz)
Amplitud
log(Amplitud)
(c) Hamming
f1 f2
Tiempo (s)
Frecuencia (Hz)
Amplitud
log(Amplitud)
(d) Gaussian
f1 f2
Tiempo (s)
Frecuencia (Hz)
(e) Flat-top
Figura 4.5: Diferentes ventanas, f (t) = cos(f1 (2πt)) + 0,01cos(f2 (2πt))
ERRORES Y VENTANAS
53
Existen dos ventanas especiales que se utilizan test de impacto. En el caso de test de impacto la se˜ nal de entrada es una se˜ nal tipo pulso y la respuesta es una combinaci´ on de sinusoides que disminuye en el tiempo (Figura 4.6). Para la respuesta se usa una ventana exponencial w(t) = e−αt . En el caso de estructuras con baja amortiguaci´on o si el tiempo de adquisici´on es muy breve, la respuesta no llegar´ a a cero al final del bloque de tiempo, lo que causa discontinuidades y leakage. Al multiplicar esta respuesta con una ventana exponencial da como resultado una se˜ nal que es casi zero al final del bloque de tiempo. En el caso de estructuras con alta amortiguaci´on o con un tiempo de adquisici´on alto, la se˜ nal llega a cero antes del final del bloque de tiempo y el resto de lo medido es b´ asicamente ruido experimental. Al aplicar una ventana exponencial, en este caso, se reduce la contribuci´on del ruido. Dado que la fuerza es de corta duraci´on, cualquier ruido durante el resto de la se˜ nal no es deseable. Aplicar una ventana que sea igual a la ventana exponencial durante el pulso, suavemente desciende a cero justo despu´es del pulso y permanece igual a cero durante el resto del bloque de tiempo, ofrece una soluci´on al problema del ruido (Figura 4.7).
Figura 4.6: Se˜ nal t´ıpica de impacto y respuesta al impacto
Figura 4.7: Ventanas de fuerza y exponencial
FUNCIONES DE RESPUESTA EN FRECUENCIA Y FUNCIONES DE COHERENCIA
54
La combinaci´on de una venta de fuerza y una ventana exponencial, a˜ nade amortiguamiento al sistema. La funciones de respuesta en frecuencia resultantes van a contener el efecto de este amortiguamiento extra. No es f´acil remover este amortiguamiento de las funciones de respuesta en frecuencia, sin embargo, puede ser considerado en la etapa de estimaci´on de par´ametros.
4.3.
Funciones de respuesta en frecuencia y funciones de coherencia
Sea F (f ) el espectro en frecuencia de una se˜ nal de entrada f (t) y X(f ) es espectro en frecuencia de una se˜ nal de salida x(t), la funci´on de respuesta en frecuencia (FRF), H(f ), entre ambas se˜ nales se define como: H(f ) =
X(f ) F (f )
(4.8)
Al calcular H(f ) con la expresi´on anterior se corre el riesgo que existan t´erminos donde F (f ) sea cero. Por lo tanto, en la practica se utilizan maneras alternativas de calcular H(f ), utilizando las potencias espectrales: H1 (f )
=
X(f ) F ∗ (f ) GXF = ∗ F (f ) F (f ) GF F
(4.9)
H2 (f )
=
GXX X(f ) X ∗ (f ) = F (f ) X ∗ (f ) GF X
(4.10)
El principal motivo para estimar las funciones de respuesta en frecuencia con las ecuaciones anteriores es la reducci´on del ruido no correlacionado en las se˜ nales de entrada o salida al promediar.
FUNCIONES DE RESPUESTA EN FRECUENCIA Y FUNCIONES DE COHERENCIA
55
En la pr´actica, la funci´ on de respuesta en frecuencia es estimada con valores promedio de las potencias espectrales,
ˆF F G
=
Na 1 X (GF F )n Na n=1
(4.11)
ˆ XX G
=
Na 1 X (GXX )n Na n=1
(4.12)
ˆF X G
=
Na 1 X (GF X )n Na n=1
(4.13)
ˆ XF G
=
Na 1 X (GXF )n Na n=1
(4.14)
donde Na es el numero de promedios (el ensayo se repite Na veces), lo que entrega una aproximaci´ on de m´ınimos cuadrados de H(f ). Dado que las funciones de respuesta en frecuencia se obtiene de una aproximaci´on de m´ınimos cuadrados, se puede definir un coeficiente de correlaci´on. En este caso, la correlaci´on se denomina funci´on de coherencia y es una medida del error de m´ınimos cuadrados. La coherencia se define por: ˆ 2 GF X H1 (f ) γ2 = (4.15) = ˆ ˆ H2 (f ) GF F GXX La coherencia varia entre 0 y 1. Un valor de 1, indica una relaci´on perfectamente lineal entre las se˜ nales de entrada y salida por sobre todos los promedios. Una coherencia menor a uno, se puede deber a uno de los siguientes motivos: Ruido no correlacionado en las mediciones de f (t) y/o x(t) No-linealidades del sistema en investigaci´on Leakage en el an´ alisis Desfase en las mediciones no compensado en el an´alisis.
FUNCIONES DE RESPUESTA EN FRECUENCIA Y FUNCIONES DE COHERENCIA
4.3.1.
56
Efectos del ruido experimental
Consideremos que las se˜ nales de excitaci´on y respuesta “verdaderas” son e(t) y r(t), respectivamente, las se˜ nales medidas son: f (t)
=
e(t) + m(t)
(4.16)
x(t)
=
r(t) + n(t)
(4.17)
donde m(t) y n(t) representa representan ruido no correlacionado. La ecuaci´ on anterior implica que las potencias espectrales de la se˜ nal son: GF F
=
GEE + GM M + GEM + GM E
GXX
=
GRR + GN N + GRN + GN R
GF X
=
GER + GEN + GM R + GM N
GXF
=
GRE + GRM + GN E + GN M
(4.18)
La funci´ on de respuesta en frecuencia verdadera viene dada por: H(f )
=
GRE GRR = GEE GER
2
=
GRR GEE
|H(f )|
(4.19)
De la ecuaci´ on anterior se tiene: GER GRE = GRR GEE
(4.20)
Bajo el supuesto razonable que las se˜ nales de ruido no est´ an correlacionadas entre ellas ni con las se˜ nales de excitaci´on y de respuesta, entonces GEM = GM E = GRN = GN R = GEN = GN E = GM R = GRM = GM N = GN M = 0, se pueden estudiar tres casos de inter´es.
FUNCIONES DE RESPUESTA EN FRECUENCIA Y FUNCIONES DE COHERENCIA
57
Caso 1. Ruido en la se˜ nal de excitaci´ on, sin ruido en la respuesta Dado que n(t) = 0, ecuaci´ on 4.18 se convierte: GF F
=
GEE + GM M
GXX
=
GRR
GF X
=
GER
GXF
=
GRE
Evaluando H1 de acuerdo a la ecuaci´on 4.9: GRE H H1 = = GEE + GM M 1 + GM M /GEE
(4.21)
(4.22)
Se ve la verdadera funci´on de respuesta en frecuencia es subestimada ya que el denominador es mayor a uno. Evaluando ahora H2 de acuerdo a la ecuaci´on 4.10: GRR GRE = =H (4.23) GER GEE lo que demuestra que, bajo la hip´ otesis de ruido no correlacionado, H2 es insensible a ruido en la se˜ nal de excitaci´on. H2 =
De lo anterior se deduce que si la se˜ nal de la fuerza esta contaminada por ruido no correlacionado, un estimador preferente para H es H2 . La funci´ on de coherencia en este caso viene dada por, 1 γ2 = 1 + GM M /GEE
(4.24)
como es esperado, la coherencia depende de la raz´on entre la se˜ nal y el ruido: mientras menor sea la raz´on entre el ruido y la se˜ nal, m´as cercana a uno es la coherencia. Caso 2. Sin ruido en la se˜ nal de excitaci´ on, ruido en la respuesta Ecuaci´ on 4.18 es en este caso: GF F
=
GEE
GXX
=
GRR + GN N
GF X
=
GER
GXF
=
GRE
(4.25)
FUNCIONES DE RESPUESTA EN FRECUENCIA Y FUNCIONES DE COHERENCIA
58
Evaluando H1 de acuerdo a la ecuaci´on 4.9: H1 =
GRE =H GEE
(4.26)
Por lo tanto, H1 es insensible a ruido no correlacionado en la respuesta. Evaluando ahora H2 de acuerdo a la ecuaci´on 4.10: GRE (1 + GN N /GRR ) GN N GRR + GN N = =H 1+ (4.27) H2 = GER GEE GRR La ecuaci´ on anterior muestra que H2 sobreestima a la verdadera FRF, dado que el valor en el par´entesis es mayor a 1. De lo anterior se deduce que si la respuesta esta contaminada por ruido no correlacionado, un estimador preferente para H es H1 . La coherencia viene dada por, γ2 =
1 1 + GN N GRR
(4.28)
Caso 3. Ruido en ambas se˜ nales En este caso, GF F
=
GEE + GM M
GXX
=
GRR + GN N
GF X
=
GER
GXF
=
GRE
(4.29)
De donde se tiene: H1 =
H GRE = GEE + GM M 1 + GM M /GEE
GRR + GN N GN N H2 = =H 1+ GER GRR
(4.30)
(4.31)
y la funci´ on de coherencia viene dada por, γ2 =
1 (1 + GN N GRR )(1 + GM M /GEE )
(4.32)
FUNCIONES DE RESPUESTA EN FRECUENCIA Y FUNCIONES DE COHERENCIA
59
Es evidente que: H1 < H < H2
(4.33)
por lo tanto una buena estimaci´on para H es la media geom´etrica de estas dos cantidades. Esta estimaci´ on se denomina Hv s p 1 + GM M /GRR Hv = H1 H2 = H (4.34) 1 + GM M /GEE A partir de los resultados anteriores se pueden sacar algunas conclusiones. En las resonancias la se˜ nal de excitaci´on es particularmente susceptible al ruido, ya que se requiere una peque˜ na fuerza para generar desplazamientos significativos, por otro lado, la respuesta tiene una raz´on se˜ nal-ruido baja. Por lo tanto, en las zonas cercanas a la resonancia se puede esperar una mejor estimaci´on si se utiliza el estimador H2 . En contraste, en las zonas lejanas a la resonancia (y espec´ıficamente cerca de las antiresonancias) la respuesta es m´as sensible a contaminaci´on por ruido. Entonces en estas zonas se obtiene una mejor estimaci´on con H1 .
Cap´ıtulo 5
Estimaci´ on de par´ ametros modales Los m´etodos de identificaci´on de par´ametros buscan extraer la informaci´on modal de una estructura a partir de mediciones experimentales. Estos m´etodos se clasifican principalmente en m´etodos en el dominio del tiempo o m´etodos en el dominio de frecuencias. Los m´etodos en el dominio del tiempo siempre se pueden utilizar, ya sea para respuesta libre o forzada (con o sin conocimiento de las fuerzas). Por otro lado, los m´etodos en el dominio de frecuencias se pueden utilizar s´olo en casos de vibraciones forzadas y cuando las fuerzas son conocidas. Para cada dominio hay m´etodos que utilizan informaci´on de un s´ olo punto de medici´on y otros que utilizan la informaci´on de varios puntos simult´aneamente. En cada uno de estos casos pueden haber una o varias fuerzas de excitaci´ on externas. Lo que lleva a la siguiente clasificaci´on: 1. Una respuesta debido a una fuerza (SISO) single-input single-output 2. Varias respuestas debido a una fuerza (SIMO) single-input multiple-output 3. Varias respuestas debido a varias fuerzas (MIMO) multiple-input multipleoutput 4. Una respuesta debido a varias fuerzas (MISO) multiple-input multiple-output En el dominio de tiempo, las respuestas contienen naturalmente informaci´on acerca del contenido en frecuencias, aunque est´a “escondida”, por lo que no es posible definir a priori cuantas resonancias hay presentes en un cierto periodo de tiempo. En
60
´ METODOS DE UN GRADO DE LIBERTAD
61
consecuencias, los m´etodos en el dominio del tiempo deben estimar simult´aneamente varias resonancias de la estructura y para los m´etodos SIMO y MIMO varios modos de vibraci´on. Estos m´etodos son conocidos como m´etodos de multiples grados de libertad (MDOF). En el dominio de frecuencia, dado que los “peaks” de las resonancias son visibles, es posible realizar una identificaci´ on modo por modo. Estos m´etodos se denominan m´etodos de un grado de libertad (SDOF). A continuaci´ on se detallar´an algunos de los m´etodos de identificaci´on de par´ametros m´ as utilizados.
5.1.
M´ etodos de un grado de libertad
En general, la respuesta din´ amica de un sistema es una superposici´on de muchos de sus modos. Sin embargo, si en una cierta banda de frecuencia se asume que un solo modo es importante, lo par´ametros de este modo se pueden determinar separadamente. Los m´etodos basados en este supuesto se denominan “m´etodos de un grado de libertad”, estos m´etodos son una manera r´apida de obtener los par´ ametros modales ya que requieren poco esfuerzo y memoria computacional.
5.1.1.
Peak picking
El m´etodo “peak picking” es talvez el m´etodo m´as sencillo dentro de los de un grado de libertad. Este m´etodo utiliza la curva de la FRF en la cercan´ıa a una resonancia como si fuese la curva de un sistema de un grado de libertad. El procedimiento del m´etodo peak picking es el siguiente: Estimar la frecuencia natural La frecuencia natural se selecciona como la frecuencia del m´aximo en la curva de FRF, ωr = ωm´ax . Estimar el amortiguamiento Para estimar el amortiguamiento, se ubican las frecuencia ubicadas a cada lado del m´ aximo identificado y que corresponden a una amplitud de la FRF αm´ax igual a √ (ver Figura 5.1). La raz´on de amortiguamiento se puede estimar 2 como, ζr =
ωb2 − ωa2 ωb − ωa ≈ 2 4ωr 2ωr
(5.1)
Estimar la constante modal La funci´on de respuesta en frecuencia a un sistema de un grado de libertad
ηr =
ζr =
or (c)
2ω r2
≈
(8.5)
ωr
ω b2 – ω a2 ω b – ω a ≈ 2ω r 4ω r2
(8.6)
Estimating the modal constant
Ar . The η r ω r2 62 ´ METODOS DE UN GRADO DE LIBERTAD 2 modal constant Ar can be estimated from Ar = α max η r ω r . For viscous damping model, this becomes Ar = 2α maxζ r ω r2 .
From the SDoF model, the FRF at the peak is known to be α max =
cerca de la resonancia viene dada por (ecuaci´on 1.29),
Due to its remarkable simplicity, the peak-picking method (Figure 8.4) can derive quick analysis results. However, it is not capable of producing accurate modal data. A A A A ' the peak FRF = value, which is very=difficult − =to measure accurately, (5.2) r ) on This methodH(ω relies (jωr − λr ) (jωr − σr − jωr ) σr ζr ωr to estimate the natural frequency and modal constant. Damping is estimated from Conocido el valor aximo de data la FRF, la are constante A sehalf puede estimar half power points only. del No m´ other FRF points used. The power points como: αm´ax ζr ωr as it is unlikely that they are two of the measured data points. have to be interpolated,
FRF (dB)
| α (ωr) | | α (ω r ) | 2
Frecuencia
ωb
ωr
ωa
Figure 8.4 Peak-picking method
Figura 5.1: Peak picking
5.1.2.
Circle fitting
El m´etodo de “circle fitting” es un de m´etodo de un grado de libertad, puede estimar los polos del sistema y los modos normales (reales o complejos). Este m´etodo se basa en el hecho que la funci´on de respuesta en frecuencia de un sistema de un grado de libertad describe un circulo en el diagrama de Nyquist (ver Figura 1.11). Si la influencia de otros modos se aproxima por una constante compleja R + jI, la funci´ on de respuesta en frecuencia cerca de la resonancia, ωr , se puede escribir como, H(ω) '
U + jV + R + jI −σr + j(ω − ωr )
(5.3)
Para determinar los par´ametros de H, primero se selecciona un conjunto de puntos en las cercan´ıas de la resonancia seleccionada. Luego se ajusta un circulo a estos puntos, como se ilustra en la Figura 5.2.
´ METODOS DE UN GRADO DE LIBERTAD
63
Figura 5.2: Circle fitting La frecuencia natural amortiguada, ωr , corresponde al punto de m´axima raz´on de cambio de a´ngulo entre puntos (m´axima distancia angular entre puntos) o con un ´angulo de fase cercano al ´angulo de fase del centro del circulo, para modos bien separados la diferencia entre estos ´angulos ser´a peque˜ na. La raz´ on de amortiguamiento, ζr , se puede estimar como, ζr =
ω2 − ω1 ωr (tan(θ1 /2) + tan(θ2 /2))
(5.4)
donde, ω1 y ω2 son dos frecuencias a ambos lados de ωr θ1 y θ2 son los ´ angulos entre el radio de ω1 y ω2 y el radio de ωr . El di´ ametro del circulo y la posici´ on angular de la frecuencia natural contiene la informaci´ on del residuo complejo U + jV : √ φ
=
tan(α)
=
U V
U2 + V 2 σr
(5.5) (5.6)
donde φ es el di´ ametro del circulo y α es el ´angulo entre la linea que conecta el centro del circulo con la frecuencia natural y el eje imaginario.
´ METODOS CON MULTIPLES GRADOS DE LIBERTAD EN EL DOMINIO DE FRECUENCIAS
64
5.2.
M´ etodos con multiples grados de libertad en el dominio de frecuencias
5.2.1.
M´ etodo de m´ınimos cuadrados no lineales, LSFD
El m´etodo de m´ınimos cuadrados no lineales en el dominio de frecuencias, es un m´etodo para estimar los polos y modos normales (si se utilizan multiples respuestas). Se basa en el modelo modal en el dominio de frecuencias. La funci´on de respuesta en frecuencia entre una respuesta en el punto i y una excitaci´ on en el punto k se puede aproximar por: Nm X LRik φir Lrk φ∗ir L∗rk Hik (ω) = + U Rik − + ∗ jω − λ jω − λ ω2 r r r=1 | {z }
(5.7)
Gi k(ω)
donde U Rik y LRik son los residuos superiores e inferiores respectivamente. Los residuos aproximan el efecto de modos bajo y sobre el rango de frecuencias en estudio. Hik (ω) es la funci´on de respuesta en frecuencia medida experimentalmente, mientras que el lado derecho de la ecuaci´on es el modelo modal considerando Nu par´ ametros desconocidos λr , φir , Lrk , U Rik , LRik , indicado en la funci´on: Gik (ω) = Gik (ω, λr , φir , Lrk , U Rik , LRik )|r=1,...,Nm
(5.8)
La diferencia entre la funci´on de respuesta en frecuencia medida y estimada viene dada por: eik (ω) = Hik (ω) − Gik (ω)
(5.9)
El error total en el rango de frecuencias de inter´es es: Eik =
Nf X
eik (ωf )e∗ik (ωf )
(5.10)
f =0
Considerando todas las funciones de respuesta en frecuencia entre Ni entradas y No respuestas, el error total es: E=
N0 X Ni X i=1 k=1
Eik
(5.11)
´ METODOS CON MULTIPLES GRADOS DE LIBERTAD EN EL DOMINIO DE FRECUENCIAS
65
Los par´ ametros desconocidos se obtienen al imponer que ´estos minimicen el error total: ∂E ∂λr
=
0
(5.12)
.. . ∂E ∂LRik
=
(5.13) 0
(5.14)
El sistema de ecuaciones anterior es altamente no-lineal, pero puede resolverse iterativamente como un problema linealizado (con una expansi´on de primer orden). Debido a las desventajas conocidas de estos algoritmos, como la necesidad de tener buenos valores de partida, la velocidad de convergencia limitada, el riesgo de divergencia, etc., estos m´etodos nunca se hicieron populares. Sin embargo, pueden ser una herramienta u ´til para mejorar la precisi´on de un modelo modal existente. Si los polos del sistema y los factores de participaci´on modal ya fueron estimados, la ecuaci´on 5.7 se convierte en un set de ecuaciones lineales en funci´on de los par´ametros desconocidos: φir , U Rik y LRik . El set de ecuaciones resultantes es relativamente f´ acil de resolver.
5.2.2.
Rational Fractional Polynomials
La idea detr´as de este m´etodo es expresar la funci´on de respuesta en frecuencia como el cuociente entre dos polinomios. Estableciendo la relaci´on entre los coeficientes de los polinomios y los par´ametros modales, se puede llegar a la identificaci´on de los par´ ametros. La definici´on de la funci´on de respuesta en frecuencia como el cuociente de dos polinomios viene de la definici´ on de la funci´ on de transferencia en el dominio de Laplace, H(s) =
N (s) a0 + a1 s + a2 s2 + . . . + am sm = D(s) b0 + b1 s + b2 s2 + . . . + sn
(5.15)
En este caso el orden del denominador, n, es mayor que el del numerador, m, por 2. Para simplificar se define: p0 (s) = 1,
p1 (s) = s,
p2 (s) = s2 ,
...,
pm (s) = sm
(5.16)
q0 (s) = 1,
q1 (s) = s,
q2 (s) = s2 ,
...,
qm (s) = sm
(5.17)
´ METODOS CON MULTIPLES GRADOS DE LIBERTAD EN EL DOMINIO DE FRECUENCIAS
Entonces la funci´ on de respuesta en frecuencia se puede expresar como: Pm ak pk (jω) H(ω) = Pk=0 n k=0 bk qk (jω)
66
(5.18)
Para el an´alisis subsiguiente se asume que hay mediciones en p frecuencias ω1 , ω2 , . . ., ωp . Tambi´en se asume que hay p frecuencias negativas ωp , ωp−1 , . . ., ω−1 . Se puede demostrar para las frecuencias negativas se cumple que, H(ω−i ) = H(−jωi ) = H ∗ (jωi )
(5.19)
El prop´ osito de las frecuencias negativas quedar´a claro m´as adelante. Para comenzar la identificaci´ on de los coeficientes se define la funci´on de error, ˜ e(ω) = H(ω) − H(ω)
(5.20)
˜ donde H(ω) es la FRF estimada como el cuociente de dos polinomios y H(ω) es la FRF experimental. Lo que lleva a: Pm ak pk (jω) ˜ e(ω) = Pk=0 − H(ω) (5.21) n k=0 bk qk (jω) Este error es una funci´on no-lineal de los coeficientes ak y bk . Para facilitar el an´ alisis el error se re-define como: ! n m n−1 X X X ˜ eˆ(ω) = e(ω) bk qk (jω) = ak pk (jω)− H(ω) bk qk (jω) + qn (jω) (5.22) k=0
k=0
k=0
El error total para todas las frecuencias se puede escribir en forma matricial: E
=
{e(ω−p ), . . . , e(ω−1 ), e(ω1 ), . . . , e(ωp )}
E
=
[P ]{A} − [Q]{B} − {W }
T
(5.23) (5.24)
´ METODOS CON MULTIPLES GRADOS DE LIBERTAD EN EL DOMINIO DE FRECUENCIAS
67
donde,
p0 (jω−p ) p1 (jω−p ) . . . pm (jω−p ) .. . ... ... ... p0 (jω−1 ) p1 (jω−1 ) . . . pm (jω−1 ) p0 (jω1 ) p1 (jω1 ) . . . pm (jω1 ) .. . ... ... ... p0 (jωp ) p1 (jωp ) . . . pm (jωp )
=
[Q]
=
˜ H(jω−p )q0 (jω−p ) ... ˜ H(jω −1 )q0 (jω−1 ) H(jω ˜ )q0 (jω1 ) 1 ... ˜ H(jω p )q0 (jωp )
A
=
{a0 , a1 , . . . , am }T
B
=
{b0 , b1 , . . . , bn−1 }T
W
=
˜ ˜ ˜ {H(jω −p )qn (jω−p ), . . . , H(jω−1 )qn (jω−1 ), H(jω1 )qn (jω1 ), . . .
[P ]
˜ H(jω −p )q1 (jω−p ) . . . .. . ... ˜ H(jω−1 )q1 (jω−1 ) . . . ˜ H(jω ... 1 )q1 (jω1 ) .. . ... ˜ H(jω )q (jω ) . . . p 1 p
˜ H(jω −p )qm (jω−p ) ... ˜ H(jω−1 )qm (jω−1 ) ˜ H(jω 1 )qm (jω1 ) ... ˜ H(jω )q (jω ) p m p
T ˜ , H(jω p )qn (jωp )}
La magnitud total del error se puede definir como: J = {E}H {E}
(5.25)
donde H denomina la transpuesta conjugada. Para determinar los coeficientes {A} y {B} que minimizan la magnitud del error, las siguientes derivadas parciales deben ser iguales a cero: ∂J ∂J = =0 ∂{A} ∂{B}
(5.26)
Lo que lleva al siguiente sistema de ecuaciones para estimar los coeficientes: 1
+ [P ]T [P ]∗ ) (−Re([P ]H [Q]))T
H 2 ([P ] [P ]
−Re([P ]H [Q]) 1 H T ∗ 2 ([Q] [Q] + [Q] [Q] )
{A} Re([P ]H {W }) = {B} Re([Q]H {W })
(5.27)
´ METODOS CON MULTIPLES GRADOS DE LIBERTAD EN EL DOMINIO DEL TIEMPO
68
Los par´ametros modales se puedes determinar a partir de los coeficientes {A} y {B}, al representar el cuociente entre polinomios como fracciones parciales: H(ω) =
n/1 X k=1
rk rk∗ + jω − pk jω − p∗k
(5.28)
donde pk es el k-esimo polo y rk el k-esimo residuo.
5.3.
M´ etodos con multiples grados de libertad en el dominio del tiempo
Los m´etodos en el dominio del tiempo, a diferencias de los algoritmos en el dominio de frecuencias, utilizan informaci´ on de la respuesta en el tiempo. Estos algoritmos ajustan los datos a la funci´on de respuesta a un impulso. La funci´on de respuesta al impulso para un sistema con multiples grados de libertad viene dada por:
h(t)
=
N X λ∗ r t) (Qr φr φTr eλr t + Q∗r φ∗r φ∗T r e
(5.29)
r=1
=
N X ∗ (Ar eλr t + A∗r eλr t )
(5.30)
r=1
=
2N X
Ar eλr t , para r > N : Ar = A∗r ; λr = λ∗r
(5.31)
r=1
5.3.1.
El m´ etodo Ibrahim
El m´etodo de Ibrahim, tambi´en conocido como ITD, construye un problema de valores y vectores propios a partir de la respuesta en el tiempo del sistema. La soluci´on de este problema permite derivar las frecuencias naturales, factores de amortiguamiento y los modos normales. La respuesta en el tiempo en un punto i, para un instante de tiempo tj se puede expresar como la suma de respuestas individuales de cada modo: xi (tj ) =
2N X r=1
φir eλr tj
(5.32)
´ METODOS CON MULTIPLES GRADOS DE LIBERTAD EN EL DOMINIO DEL TIEMPO
Considerando q puntos de medici´on x1 (t1 ) . . . x1 (tL ) φ11 .. . . . .. .. = .. . xq (t1 ) |
... {z
xq (tL )
φq1 }
q×L
|
69
y L instantes de tiempo, se tiene que, λ t . . . φ1L e 1 1 . . . eλ1 tL .. .. .. .. .. (5.33) . . . . . λ2N t1 λ2N tL . . . φq2N e ... e {z }| {z } 2N ×L
q×2N
o bien, X = ΦΛ
(5.34)
donde Λ es una matriz constituida por los elementos eλi tj . Esta ecuaci´on por si sola es insuficiente para determinar los par´ ametros modales. Consideremos entonces un segundo set de L puntos, desfasados en ∆t con respecto al set anterior: xi (tj + ∆t) =
2N X
φir eλr (tj +∆t) =
r=1
2N X
φir eλr ∆t eλr tj
(5.35)
r=1
Definiendo: xi (tj )∗
=
xi (tj + ∆t)
(5.36)
φ∗ir
=
φir eλr ∆t
(5.37)
Siguiendo el mismo procedimiento anterior, se obtiene un segundo set de ecuaciones: X ∗ = Φ∗ Λ
(5.38)
Definiendo una matriz As de q × q con la siguiente propiedad: As Φ = Φ∗
(5.39)
Premultiplicando la ecuaci´on 5.34 por As : As X = As ΦΛ = Φ∗ Λ
(5.40)
Sustituyendo 5.38: As X = X ∗
(5.41)
Si X y X ∗ son conocidos, As se puede despejar como: As = X ∗ X +
(5.42)
donde X + es la pseudoinversa de X y viene dada por: X + = X T (X T X)−1 asumiendo que L > q.
´ METODOS CON MULTIPLES GRADOS DE LIBERTAD EN EL DOMINIO DEL TIEMPO
70
Recordando que As Φr = Φ∗r , se tiene que: A s Φr As − eλr ∆t I Φr
=
Φr eλr ∆t
(5.43)
=
{0}
(5.44)
La ecuaci´on anterior es una ecuaci´on de valores propios. Los modos normales corresponden a los vectores propios y las frecuencias naturales y factores de amortiguamiento se deducen de los valores propios: eλr ∆t
= e(σr +jωdr )∆t = αr + jγr
(5.45)
αr
= eσr ∆t cos(ωr ∆t)
(5.46)
γr
= eσr ∆t sin(ωr ∆t)
(5.47)
σr
=
ωdr
=
ln(αr2 + γr2 ) 2∆t arctan αγrr ∆t
(5.48)
(5.49)
Finalmente, los par´ ametros din´amicos se obtienen como: q
2 + σ2 ωdr r
ωr
=
ζr
σr = −p 2 ωdr + σr2
5.3.2.
(5.50) (5.51)
El m´ etodo least-squares complex exponential (LSCE)
Tenemos que la respuesta a un impulso IRF, se puede expresar como:
h(t) =
2N X r=1
Ar eλr t
(5.52)
´ METODOS CON MULTIPLES GRADOS DE LIBERTAD EN EL DOMINIO DEL TIEMPO
71
Si consideramos que los datos son adquiridos en intervalos de tiempo k∆ (k = 0, 1, . . . , 2N ), se tiene la siguiente serie
h(k∆)
=
2N X
Ar eλr k∆
(k = 0, 1, . . . , 2N )
(5.53)
r=1
hk
=
2N X
Ar zrk
(k = 0, 1, . . . , 2N ),
zrk = eλr k∆
(5.54)
r=1
Cada termino de la suma posee valores reales, aunque los residuos Ar y los polos λr tienen valores complejos. Se puede demostrar que existe un polinomio de coeficientes reales de manera que: β0 + β1 zr + β2 zr2 + . . . + β2N zr2N = 0
(5.55)
Esta ecuaci´on se conoce como la ecuaci´on de Prony. Dado que hay 2N +1 ecuaciones en la expresi´ on 5.58, se pueden multiplicar las ecuaciones por un coeficiente β correspondiente y sumarlas para formar siguiente ecuaci´on: 2N X
βk hk
=
r=1
2N X
βk
r=1
=
2N X
2N X
Ar zrk
(5.56)
βk zrk
(5.57)
r=1
Ar
r=1
2N X r=1
Esta u ´ ltima expresi´ on sabemos que es igual a 0, si zr es una ra´ız de la ecuaci´ on 5.55. Lo que nos lleva a una relaci´ on simple entre los coeficientes β y los datos de la respuesta a un impulso: 2N X
βk hk = 0
(5.58)
r=1
A partir de la ecuaci´on anterior se pueden determinar los valores de βk . Conocidos los coeficientes βk se puede resolver la ecuaci´ on 5.55 y extraer las ra´ıces zr . A partir de zr se calculan las frecuencias naturales y razones de amortiguamiento: ωr
=
1p ln zr ln zr∗ ∆
(5.59)
ζr
=
− ln(zr zr∗ ) 2ωr ∆
(5.60)
DIAGRAMAS DE ESTABILIDAD
72
Para identificar los modos propios del sistema, de la siguiente forma: 1 1 ... 1 A1 z1 A2 z . . . z 2 2N .. .. .. .. .. . . . . . z12N −1
z22N −1
...
2N −1 z2N
podemos reescribir la ecuaci´on 5.54
A2N
=
h0 h1 .. .
(5.61)
h2N −1
La soluci´on de este sistema de ecuaciones entrega los residuos y por lo tanto los modos normales.
5.4.
Diagramas de estabilidad
En todos los algoritmos anteriores se debe definir a priori cual es el numero de polos que se desea estimar. Lo que muchas veces no es una decisi´on simple. Si se utilizan menos polos que lo actuales, el ajuste no ser´ a adecuado. En cambio, si se definen m´as polos que los reales, los algoritmos van a entregar polos “computacionales” que no corresponden a polos reales del sistema, sino que son polos que tratan de modelar el ruido en los datos. Una metodolog´ıa muy u ´til para determinar el numero de polos reales en un sistema son los diagramas de estabilidad. Los diagramas de estabilidad son una herramienta muchas veces imprescindible en un an´ alisis modal experimental, ellos ayudan a separar los polos reales de los polos computacionales. Estos diagramas se obtienen, al repetir el an´alisis modal incrementando el orden del sistema (n´ umero de polos asumidos). Para cada orden se estiman los polos, los resultados son presentados gr´aficamente en el diagrama de estabilidad (Figura 5.3). En el eje vertical se encuentra el orden y el eje horizontal representa la frecuencia natural del polo estimado. En general, los polos reales aparecen a la misma frecuencia en el diagrama, independiente del orden del sistema. En cambio, la frecuencia de los polos computacionales, var´ıa al aumentar el orden del sistema. El diagrama de la Figura 5.3 se representan los polos estables (en frecuencia y amortiguamiento) con un circulo y con una cruz los polos inestables. Normalmente, se define que un polo es estable si la frecuencia no var´ıa m´as del 1 % de su magnitud y la raz´on de amortiguamiento no var´ıa m´as del 5 %.
DIAGRAMAS DE ESTABILIDAD
73
40
Orden
30
20
10
0
20
40 60 80 Frecuencia (Hz)
100
Figura 5.3: Diagrama de estabilidad
120
Cap´ıtulo 6
Medici´ on Experimental El primer paso en un an´alisis modal experimental es adquirir las funciones de respuesta frecuencia experimentales de una estructura. El m´etodo m´ as usual es excitar la estructura con una fuerza conocida y medir de forma simultanea la fuerza y las respuestas en la estructura. Como resultado, se obtiene un grupo de FRFs que pueden ser utilizadas posteriormente para derivar los par´ametros modales de la estructura utilizando algunos de los algoritmos descritos en el capitulo anterior.
6.1.
An´ alisis previo a las mediciones
Al preparar un montaje experimental para un an´alisis modal se debe considerar el prop´osito del experimento, los datos requeridos (FRFs o par´ametros modales), la precisi´on requerida, etc. Para ello se necesita la mayor cantidad de informaci´on posible de la estructura, la que se puede obtener de experimentos anteriores en estructuras similares o de modelos num´ericos de la estructura. A continuaci´on se describir´ an algunas herramientas u ´ tiles para definir el montaje experimental de acuerdo a los requerimientos. Desde un punto de vista pr´actico, los criterios siguientes se deben cumplir en un buen dise˜ no de un montaje experimental: Correspondencia: Los modos medidos experimentalmente deben corresponder a los modos reales, los que desafortunadamente son desconocidos. Sin embargo, experimentos previos en estructuras similares o un modelo en elementos finitos pueden ayudar a estimar los modos. Adicionalmente, el test debe producir modos claramente distinguibles. La independencia del los
74
´ ANALISIS PREVIO A LAS MEDICIONES
75
modos esta directamente relacionada con el rango de la matriz de vectores propios Φ. Excitaci´ on: El montaje debe incluir un sistema de excitaci´on que garantice que todos los modos de inter´es son excitados. Identificaci´ on: Los datos medidos deben contener la informaci´on necesaria para identificar los par´ ametros de inter´es. Por lo tanto, el dise˜ no del montaje depende del prop´ osito del experimento. Visualizaci´ on: En la pr´actica, se requiere visualizar los modos obtenidos, de manera de evaluar de precisi´on de ´estos y para compararlos con modos estimados. La visualizaci´on tambi´en es importante para detectar posibles discrepancias. Robustez: Dado que el montaje est´a basado en experimentos previos o en modelos num´ericos, donde ambos contiene errores, ´este debe ser robusto: No debe ser muy sensible a estos errores. Por lo tanto, es preferible alg´ un grado de redundancia. Accesibilidad: Las ubicaciones seleccionadas para medir la respuesta y para excitar la estructura deben ser accesibles.
6.1.1.
Rango de frecuencias
Es importante definir bien el rango de frecuencias para identificar todos los modos de inter´es. Si se sabe cuantos modos se quieren identificar, el rango de frecuencias se puede definir con un an´alisis previo de un modelo en elementos finitos. Desde luego se debe considerar un margen de error. En caso que no se cuente con un modelo en elementos finitos, el rango de frecuencias se puede definir in situ. Al intentar en una primera instancia con una rango de frecuencias amplio, luego contar los peaks de una de las FRFs para definir la frecuencia m´ axima que se debe medir para identificar n modos.
6.1.2.
Selecci´ on de la ubicaci´ on de las respuestas
En un an´alisis modal se deben definir suficientes puntos de medici´on de manera que los modos identificados sean independientes. A continuaci´on se describir´an algunos de los algoritmos disponibles para determinar la ubicaci´on ´ optima de los puntos de medici´on. En los algoritmos a continuaci´on se asume que se cuenta con los modos normales obtenidos a trav´es de alg´ un m´etodo num´erico. La independencia de los modos se puede determinar a partir del rango de la matriz Φ. Este rango equivale al rango de la matriz ΦT Φ (la matriz de Fisher).
´ ANALISIS PREVIO A LAS MEDICIONES
76
La contribuci´on de cada grado de libertad i al rango de la matriz de Fisher, viene dado por el siguiente indicador (independencia efectiva): −1 T EfIi = diag Φ ΦT Φ Φ (6.1) Al ir eliminando iterativamente los grados de libertad con menor EfIi y re-calculando los nuevos EfI, se llega a un set ´optimo de puntos de medici´ on con respecto a la independencia de los modos de inter´es. Otro indicador de la independencia entre modos es el MAC (Modal Assurance Criterion), el que se define como:
2 φTi φj M ACij = T T φi φi φj φj
(6.2)
Donde φi es el iesimo modo normal. MAC es un factor que mide la correlaci´on entre dos modos. Un valor de 0 indica que no hay correlaci´on, mientras que un valor de 1 indica dos modos perfectamente correlacionados. Si se calcula la correlaci´on para todos los pares de modos, se puede construir la matriz MAC. Si los modos son 100 % independientes, los valores fuera de la diagonal de la matriz son todos cero y en la diagonal son todos 1. Los puntos de medici´on se determinan al encontrar al menor set de puntos de medici´ on que minimize los valores fuera de la diagonal de la matriz MAC. En general, se aceptan valores de hasta un 20 % de correlaci´on cruzada. Este es el m´etodo m´as utilizado en la selecci´ on de los puntos de medici´on.
6.1.3.
Selecci´ on de los puntos de excitaci´ on
Los puntos de excitaci´on de la estructura se deben seleccionar de manera garantizar que todos los modos de inter´es sean excitados adecuadamente. Un modo en particular va a estar bien excitado si la fuerza se aplica en un punto de alto desplazamiento. El m´etodo m´as usual para seleccionar los puntos de excitaci´on es a trav´es de los “driving point residues” (DPR). El residuo Aikr esta definido por la expresi´on de la funci´ on de respuesta en frecuencia en t´erminos de los par´ametros modales:
´ ANALISIS PREVIO A LAS MEDICIONES
Hik (ω)
=
N X Qr φir φkr r=1
= 146
Modal Analysis
jω − λr
77
Q∗ φ∗ φ∗ + r ir kr jω − λr
N X A∗ikr Aikr + jω − λr jω − λr r=1
(6.3)
(6.4)
obtained from del a laboratory measurement to analyse and simulate the De modal donde data se define el DPR punto i como: dynamic behaviour in situ. Other structures cannot fit into a laboratory environment for FRF measurement. They have to be tested in situ. In both cases, it is vital to φ2ir A = (6.5) iir ensure that2jω the test conditions are stable and repeatable so that the measured FRF r data are reliable and representative. If an FRF measurement in the laboratory is desired and feasible, a structure is Al estudiar los DPR para todas candidatos de excitaci´ on subsequent y para todos often prepared to simulate free orpuntos grounded boundary conditions. The use los modos de inter´ e s, se obtiene informaci´ o n u ´ til para la selecci´ o n de los puntos of the modal model to be derived from the measurement determines which boundary de excitaci´ on. Entogeneral, grados detolibertad DPR para el mayor numero conditions use. It islosimpossible imitate con perfect freealtos or grounded conditions. posible modos van a serbebuenos puntosinde on.with reasonable accuracy. Thesede conditions can only approximated theexcitaci´ laboratory The free boundary condition is simulated by supporting the structure with soft materials such as springs or elastic bands. Such an arrangement creates one or more rigid body modesofrom stiffness of these materials and the total mass 6.1.4. Selecci´ n dethelos puntos desupporting suspensi´ on of the structure. If the natural frequencies of these rigid body modes are far from the firstque natural frequency the structure, the measured FRFreales data should be affected en Dado es muy dif´ıcilofimitar condiciones de borde en unnot laboratorio, by this boundary condition. Figure 7.7 shows a simple plate supported by four soft la mayor´ıa de las mediciones experimentales la estructura se monta de manera de springs. For the same principle, a heavy structure can be put on top of an inflated car simular unaorcondici´ on libre (sin condiciones de borde). on se logra tyre tube a few layers of thick porous packaging material.Esta condici´ al suspender la estructura con materiales blandos como resortes o el´ a sticos. Con The grounded boundary condition is more difficult to simulate in the laboratory. esteTheoretically, arreglo se producen de means cuerpothat r´ıgido. lasdegrees frecuencias naturales a groundedmodos condition all theSisix of freedom at the de los boundary modos r´ıgidos son lejanas a lacan primera naturalinde la estructura, are rigidly fixed. This almost frecuencia never be achieved reality. In modal la medici´ on de FRF no se ver´ a afectada condici´ de borde. la Figura testing, thelaarrangement is often to ‘fix’por theesta structure tooan much moreEn rigid and objectun such as a concrete 6.1 heavier se muestra montaje de unafloor. placa suspendida por cuatro resortes.
Figure 7.7 A simple plate supported by four soft springs
Figura 6.1: Una placa suspendida por cuatro resortes. It is only for some special cases that a true grounded condition can be simulated. For example, accurate measurement of the cantilever in Figure 7.8 needs a rigid boundary condition at the built-in end which may be difficult to realize. The alternative is to measure the beam on the right which is freely supported with a doubled length. Any odd number of modes from this beam will be the equivalent modes for the cantilever. This is due to the fact that for these modes, the middle cross-section of the free–free beam is truly ‘fixed’. For the simulation of the free boundary condition, we often know the limitation
EL MONTAJE EXPERIMENTAL
78
Al contrario de los puntos de excitaci´ on, los puntos de suspensi´on se deben seleccionar en puntos con muy bajas amplitud de movimiento. La t´ecnica, por lo tanto es similar a la de la secci´ on anterior, pero ahora se deben seleccionar los grados de libertad con DPR bajos.
6.2.
El montaje experimental
6.2.1.
Mecanismo de Excitaci´ on
La primera parte de un montaje experimental es el mecanismo de excitaci´on que aplica una fuerza de suficiente amplitud y contenido en frecuencia a la estructura. Existen diferentes mecanismo capaces de excitar la estructura. Los dos m´as comunes son el shaker (Figura 6.2) y el martillo (Figura 6.3). Un shaker electromagn´etico, tambi´en conocido como shaker electrodin´ amico, es el tipo de shaker m´as utilizado en an´alisis modal. Consiste en un im´an, un bloque libre y una bobina. Cuando una se˜ nal el´ectrica pasa por la bobina, se genera una fuerza proporcional a la corriente y a la densidad de flujo magn´etico generado, la que mueve al bloque. Un shaker electromagn´etico tiene un rango amplio de frecuencias y de amplitudes. Para frecuencias bajas y altas amplitudes de excitaci´ on, se puede utilizar un shaker electro-hidr´aulico. Generador de la señal
Analisador
Amplificador de poder
Shaker Sensor de fuerzas
Signal flow
Acondicionador de la señal Acelerómetro
Figura 6.2: Set-up experimental con excitaci´on por shaker Un martillo es un aparato que produce una fuerza de excitaci´ on de la forma de un pulso. Consiste en la punta del martillo, un sensor de fuerzas, una masa y el mango. La punta del martillo se puede cambiar para alterar su dureza. Materiales t´ıpicos para puntas de martillo son; goma, pl´ astico y acero. La dureza de la punta, en conjunto con la dureza de la superficie de la estructura, est´a directamente
EL MONTAJE EXPERIMENTAL
79
relacionada con el rango de frecuencias contenida en la se˜ nal de la fuerza. Mientras mayor es la dureza, mayor es el rango de frecuencias contenido. 142 Modal Analysis
Figure 7.2 A measurement set-up with hammer excitation
Figura 6.3: Set-up experimental con excitaci´on por martillo A hammer is a device that produces an excitation force pulse to the test structure. It consists of hammer tip, force transducer, balancing mass and handle. The hammer tip can beAceler´ changedometro to alter the hardness. Typical materials for the tip are rubber, 6.2.2. plastic and steel. The hardness of the tip, together with that of the structure surface to be tested, is directly related to the frequency range of the input pulse force. For a El aceler´ometro es el sensor m´as com´ un utilizado en an´alisis modal. Mide la hard tip striking on a hard surface, we can expect the force pulse to distribute energy aceleraci´on de un punto en la estructura, la se˜ nal de salida viene en la forma de to a wide range of spectrum. This is the only mechanism to control the frequency voltaje. Esta se˜ n al es transformada por un acondicionador de la se˜ nal antes de ser components of excitation in a hammer test. procesada por otro hardware o software. An electromagnetic shaker, also known as an electrodynamic shaker, is the most common type of shaker used inometros modal testing. consists of a magnet, a moving El tipo m´ as com´ un de aceler´ son losItpiezoel´ ectricos, ilustrados en lablock Figura and a coil in the magnet. When an electric current from a signal generator passes 6.4. Estos sensores contienen un cristal piezoel´ectrico en su interior, ´este cristal through the inside shaker, a force proportional to the current the omagnetic produce unacoil carga el´ethe ctrica al ser deformado. Al seleccionar unand aceler´ metro hay flux density is generated which drives the moving block. An electromagnetic shaker distintos factores importantes que se deben considerar, estos son; el rango de has a wide frequency, amplitude and dynamic range. For low frequency and large frecuencias, sensibilidad, masa y estabilidad bajo cambios de temperatura. La amplitude excitation, an electrohydraulic exciter can be used. sensibilidad de un aceler´ometro determina la raz´on entre la se˜ nal medida y el ruido. Mientras m´as alta es la sensibilidad del aceler´ ometro m´as precisas son las mediciones. La masa del aceler´ometro tambi´en es muy importante, ya que puede 7.2.2 Accelerometer modificar las caracter´ısticas de la estructura. Mientras menor es la masa mejor, aunque esto significa veces sensor una menor sensibilidad. An accelerometer is themuchas most common for modal testing. It measures acceleration of a test structure and outputs the signal in the form of voltage. This signal will be transformed by a signal conditioner before it is processed by an analyser, other hardware or software. The accelerometer does not assume the properties of the measured structure such as linearity. An accurate accelerometer only records faithfully the acceleration at the measurement location. There are two aspects in the acceleration measurement that a sensor needs to be capable of dealing with. One is the frequency and the other is the amplitude. Both are reflected in the input–output relationship of an accelerometer. An ideal accelerometer should have a linear input and output relationship in order to ensure that the amplitude content of the acceleration signal at different frequencies is truthfully recorded. The FRF of the accelerometer should be uniformly flat so that the amplitude of no frequency is distorted. The accelerometer should also impose zero phase shift to the signal measured.
shape to ofshow a crystal thus leadingFigure to the7.3change electric charge. At the lowoffrequency their characteristics. shows aof typical frequency response curve end, piezoelectric accelerometers do not respond to DC signal. At the high frequency an accelerometer. end, the accuracy of measurement is degraded by the natural frequency of the sensor. When selecting an accelerometer for modal testing, number of factors need to be Sensitivity Phasea[degrees] [dB] thought through. The main parameters affecting the performance of a piezoelectric 30 accelerometer are: the frequency response property; the sensitivity and its stability under temperature change; cross-axial sensitivity and base strain. 20 EL MONTAJE EXPERIMENTAL 80 Ar
10 10
Amplitude 0
Resorte 0 pre-comprimido –10
Phase
–20
–10 1
2
5
10
20
50
Masa 100 kHz
Cristal + an accelerometer Figure 7.3 A typical chart for q
The characteristics of an accelerometer are its–potential. They can be fully realized Cuerpo if the sensor is connected rigidly on the structure. In reality, this is not to be the case. An accelerometer has to be mounted non-rigidly on a structure for measurement. If of a aceler´ piezoelectric considered Figure as aFigura rigid7.5 mass block, thede accelerometer and piezoel´ its accelerometer mount can be modelled 6.4:Diagram Esquema un ometro ectrico as an SDoF system as shown in Figure 7.4.
The frequency response property determines the linearity of the sensor. The sensitivity Sensor the signal to noise ratio. of an accelerometer dictates Large and stable sensitivity Acelerómetro means accurate Estructura measurement. Cross-axial sensitivity causes inaccuracy in measurement. Base strain is caused by the flexure of the accelerometer base interacting non-rigid Montaje structure surface. Usually, a more sensitive accelerometer is more bulky. The accelerometer mass has the potential to change the characteristics ofEstructura the test structure. This is particularly so if the accelerometer is located at or close to an anti-nodal point 7.4 An SDoF model of an accelerometer and its mounting Figura 6.5: mass Montaje de un aceler´ ometrosignificant natural frequency of a vibration Figure mode where a minute change can cause shift. Thus, when selecting a better sensitivity, we need to be mass conscious. La precisi´on de un aceler´ometro tiene una alta dependencia con la manera en que The accuracy of the acceleration measurement depends largely on the mounting es montado sobre la estructura. En la Figura 6.5 se muestra el montaje t´ıpico de un which is modelled by a spring and a damper. The accelerometer is of course more aceler´ometro. La flexibilidad del montaje afecta las mediciones del sensor, mientras just a mass block. As Figure 7.3 shows, it has its own natural frequency. This 7.2.3than Force m´as r´ıgidotransducer es el montaje mayor es la precisi´on. Existen diferentes mecanismos frequency is usually much higher magn´ than the frequency of the de SDoF Figure de montaje; tornillo, adhesivo, etico, capa delgada cera,system entre in otros. Al 7.4. The best accuracy would arise if the mounting were rigid. The flexibility ofesto the atornillar el sensor a la estructura se obtienen los mejores resultados, aunque A forcemounting transducer is another type of sensor used in modal testing. Like an accelerometer, that the characteristics of the accelerometer ares´ocompromised significa means un montaje experimental m´as complejo. En la pr´actica, lo se utilizan a piezoelectric force transducer generates anthe output charge or voltage that isthat proportional somewhat. Because of it, acceleration from structure may be different sensores atornillados en estructuras que deben ser monitoreadas de forma from continua.
by the However, if the natural frequency of this SDoF though, to the experienced force applied to accelerometer. the transducer (Figure 7.6). Unlike an accelerometer system is five times or more of the frequency of the acceleration signal from the element. a force transducer does not have an inertial mass attached to the transducing Sensor de fuerzas It has to6.2.3. be physically compressed or stretched so that the transducing part can generate output. For a shaker test, a force transducer has to be connected between the El sensor de fuerzas es otro tipo de sensor utilizado en an´alisis modal. Al igual que structureunsurface and the shaker. For a hammer test,produce the transducer located at the aceler´ometro, un sensor de fuerzas piezoel´ectrico una carga ois voltaje hammerproporcional tip and is acompressed whenA impact isdeapplied to. la fuerza aplicada. diferencia un aceler´ ometro, un sensor de fuerzasthe no tiene una masa inercial su interior, si no que debe measurement ser comprimido accuracy The ways characteristics of a enforce transducer affect o estirado f´ısicamente para generar la carga (Figura 6.6). Para un montaje con
EL MONTAJE EXPERIMENTAL
81
un shaker, el sensor de fuerzas se conecta entre la superficie de la estructura en el shaker. En el caso de un martillo, el sensor se ubica en la punta del martillo y es comprimido cuando el martillo impacta la superficie. Las consideraciones para seleccionar un sensor de fuerzas son similares a las de un aceler´ometro: rango de frecuencias, sensibilidad, masa y estabilidad bajo cambios de temperatura.
Frequency response function measurement F
Cristales
+
q –
F
FigureFigura 7.6 6.6: A diagram ofuna sensor piezoelectric forceetransducer Esquema de de fuerzas piezoel´ ctrico
are similar to that for an accelerometer. They include frequency response characteristics, sensitivity and cross-axial sensitivity. The mass of the force transducer can also potentially affect the measurement outcome. For both force and acceleration sensors, this problem is most acute when the structure resonates. The main consideration in selecting a force transducer is to understand how it interacts with the excitation device to which it connects. For example, when a force transducer is used on an impact hammer, variation of the hammer tip and the mass of the hammer handle can cause a different force transducer calibration. When the force transducer is used with an excitation shaker, the presence of its mass may cause significant distortion to the force signal measured at structural resonance. The extent of distortion is dependent on the mass difference between the transducer and the structure. The mass of the transducer may also be responsible for the sensitivity the transducer has to bending moments.
7.3 Preparation of the test structure A real structure is often connected to its surroundings. Therefore, its dynamic characteristics in situ are determined by the boundary conditions as well as by itself. When the FRF measurement is to be carried out, the question to be answered is: ‘Under what conditions do we want to test the structure, stand alone or in situ?’ The answer to this question hinges on two considerations when preparing the structure for test: (1) do we need the modal model of the structure in its working condition or in the laboratory environment; and (2) is it realistic to test the structure in situ or in
145
´ DE LA FUERZA DE EXCITACION ´ SELECCION
6.3.
82
Selecci´ on de la fuerza de excitaci´ on
La precisi´on de las mediciones experimentales depende en gran parte de la forma de la fuerza de excitaci´on utilizada. Aunque te´oricamente una FRF no debiese depender de la excitaci´on utilizada, en la pr´actica la precisi´on y calidad de una FRF depende, entre otros factores, en la elecci´on de la fuerza de excitaci´ on. Es por lo tanto, muy importante en un an´ alisis modal experimental evaluar todos los m´etodos posibles de excitaci´on de manera de elegir el m´as adecuado.
6.3.1.
Excitaci´ on sinusoidal
La excitaci´on sinusoidal es la forma m´as tradicional en an´alisis modal. La fuerza contiene una sola frecuencia en un tiempo y la excitaci´on cambia de una frecuencia a otra con un paso dado, permitiendo a la estructura excitar un modo a la vez. Esta excitaci´ on es efectiva para excitar una estructura con niveles de vibraci´on altos, para caracterizar no-linealidades de una estructura y para excitar modos normales de una estructura amortiguada. La Figura 6.7 presenta una se˜ nal sinusoidal t´ıpica y su contenido en frecuencias.
Figura 6.7: Se˜ nal sinusoidal
´ DE LA FUERZA DE EXCITACION ´ SELECCION
83
En este caso la din´ amica de una estructura se descompone f´ısicamente por frecuencia y posici´on. Al sintonizar cerca de una frecuencia natural, la respuesta de la estructura es dominada por ese modo de vibraci´on. Esta segregaci´on natural proporciona una via para la identificaci´ on directa de par´ametros, la que usualmente tiene una buena raz´on se˜ nal/ruido. Esta es una posibilidad que otros m´etodos de excitaci´ on no ofrecen.
6.3.2.
Excitaci´ on aleatoria
Una se˜ nal de excitaci´on aleatoria es una se˜ nal aleatoria estacionaria que sigue una distribuci´on Gaussiana. Contiene todas las frecuencias en el rango seleccionado. Para una estructura que tiene un comportamiento no-lineal, una excitaci´on aleatoria tiende a linealizar este comportamiento. Por lo tanto la funci´on de respuesta en frecuencia derivada de una excitaci´ on puramente aleatoria va a ser una versi´on linealizada de la FRF. Esta FRF, aunque no contiene informaci´on sobre las nolinealidades, es en realidad una funci´on muy u ´ til y es percibida como la “mejor” estimaci´on para las FRFs. Sin embargo, el hecho que para una excitaci´on aleatoria ni la fuerza ni las respuestas son peri´odicas, puede llevar a errores de leakage. La Figura 6.8 presenta una se˜ nal aleatoria t´ıpica y su contenido en frecuencias.
Figura 6.8: Se˜ nal aleatoria
´ DE LA FUERZA DE EXCITACION ´ SELECCION
6.3.3.
84
Excitaci´ on pseudo-aleatoria
Una se˜ nal de excitaci´on pseudo-aleatoria es una se˜ nal aleatoria estacionaria que consiste en frecuencias discretas formadas m´ ultiplos de la resoluci´ on en frecuencia utilizada en la transformada de Fourier. Es una se˜ nal peri´ odica con amplitud fija y fase aleatoria. Esta excitaci´on elimina el problema de leakage de la se˜ nal aleatoria y usualmente entrega FRFs precisas. La Figura 6.9 presenta una se˜ nal pseudo-aleatoria t´ıpica y su contenido en frecuencias.
Figura 6.9: Se˜ nal pseudo-aleatoria
6.3.4.
Excitaci´ on aleatoria en trenes (burst random)
Una mejor se˜ nal que la pseudo-aleatoria es la denominada se˜ nal aleatoria en trenes. Una se˜ nal aleatoria en trenes se crea al encender y apagar una se˜ nal aleatoria, de esta manera la medici´on comienza con una excitaci´ on aleatoria y continua luego que la excitaci´ on se ha apagado. El espectro de la se˜ nal aleatoria en trenes tiene amplitud y fase aleatorias y contiene suficiente energ´ıa en todo el rango de frecuencias. Al seleccionar adecuadamente el tiempo de apagado de la se˜ nal, esta asemeja una se˜ nal pseudo-aleatoria pero sin la necesidad de esperar el decaimiento de la parte transiente de la respuesta. La proporci´on de tiempo que esta apagada la se˜ nal de excitaci´on se debe seleccionar de manera que la respuesta de la estructura
´ DE LA FUERZA DE EXCITACION ´ SELECCION
85
sea cero al final del periodo de medici´ on, este tiempo depende principalmente del nivel de amortiguamiento en la estructura. La Figura 6.10 presenta una se˜ nal de excitaci´ on aleatoria en trenes, su contenido en frecuencias y la respuesta t´ıpica.
Figura 6.10: Se˜ nal aleatoria en trenes
6.3.5.
Excitaci´ on de impacto
La se˜ nal en el tiempo de una se˜ nal de fuerza debido a un impacto, es un pulso con un contenido en frecuencias no controlable. En t´erminos de hardware involucrado, la excitaci´ on por impacto es relativamente simple comparado con la excitaci´on con un shaker. Es conveniente de usar y muy portable para mediciones en terreno y laboratorio. Debido a que no existe conexi´on f´ısica entre la fuerza de excitaci´on y la estructura, el test de impacto evita el problema de su interacci´on. La principal desventaja del test de impacto es la dificultad de controlar el nivel de la fuerza y su contenido en frecuencias. Esto puede afectar la raz´on se˜ nal/ruido en las mediciones, resultando en datos de baja calidad. En la Figura 6.11 se muestra un t´ıpica se˜ nal de impacto.
´ INICIAL DE LAS FRFS MEDIDAS EVALUACION
86
Figura 6.11: Se˜ nal de impacto
6.4.
Evaluaci´ on inicial de las FRFs medidas
La calidad de un an´ alisis modal recae principalmente en la calidad de las FRFs medidas. Aunque algunos m´etodos de an´alisis modal pueden minimizar ciertos errores en la informaci´on obtenida experimentalmente, ning´ un m´etodo puede rectificar errores fundamentales o mediciones err´oneas. Las propiedades modales obtenidas de FRFs err´ oneas son susceptibles a errores inaceptables. Debido a esto se vuelve esencial verificar la calidad de las funciones de respuesta en frecuencia obtenidas experimentalmente. Las FRFs medidas deben cumplir b´asicamente con dos requerimientos: (1) la estructura satisface los supuestos requeridos en an´alisis modal y (2) los errores de sistema y humanos son minimizados o eliminados. B´asicamente en an´alisis modal la estructura debe cumplir con reciprocidad, invariabilidad en el tiempo y linealidad. Si estos supuestos no son verificados, las propiedades modales obtenidas no ser´an confiables. Aunque existen algunos m´etodos para asistir en la b´ usqueda de posibles errores en las FRFs medidas, algunos de estos errores no pueden ser identificados. Por ejemplo, no es posible detectar a partir de las FRFs si el sistema de medici´ on esta calibrado correctamente.
´ INICIAL DE LAS FRFS MEDIDAS EVALUACION
6.4.1.
87
Repetibilidad
La verificaci´on m´ as simple, pero no menos u ´til, es la repetibilidad de las mediciones. Esto es, asegurar que el comportamiento din´amico de la estructura y de todo el montaje experimental no depende del tiempo. Para una fuerza de excitaci´on y un punto de medici´ on seleccionados, una estructura lineal deber´ıa dar resultados id´enticos en cada medici´on. Para un par de puntos de entrada y salida (fuerza y respuesta), se pueden realizar mediciones a intervalos de tiempo. T´ıpicamente, una FRF seleccionada se mide antes y despu´es que se termine de medir toda la estructura. Este proceso que parece trivial, puede ser muy u ´ til para verificar no s´olo que el comportamiento de la estructura es constante, sino tambi´en que las condiciones de medici´ on no han cambiado durante todo el experimento.
6.4.2.
Reciprocidad
Una estructura lineal e invariante en el tiempo debe tener la propiedad de reciprocidad. Esto significa, que una FRF deber´ıa ser id´entica si se intercambian la ubicaci´on de la fuerza y la respuesta. Te´oricamente, esta propiedad se puede derivar a partir de la simetr´ıa de las matrices de masa, rigidez y amortiguaci´on. Debido a esta simetr´ıa, la matriz FRF, que es la inversa de la matriz din´ amica, tambi´en va a ser sim´etrica. La propiedad de reciprocidad de la FRF se puede utilizar para evaluar la fiabilidad y exactitud de los datos medidos.
6.4.3.
Linealidad
El supuesto m´ as com´ un en an´alisis modal es que la estructura se comporta de manera lineal. Sin este supuesto, el an´alisis modal no tiene sentido. Una manera de verificar la linealidad de una estructura es asegurar que la medici´on de las FRFs es independiente de la amplitud de las fuerzas de excitaci´on. Para verificar esto, se pueden repetir mediciones en una misma ubicaci´on, pero variando la amplitud de la fuerza de excitaci´on.
6.4.4.
Caracter´ısticas especiales de una FRF
A partir de la teor´ıa de an´ alisis modal, se pueden derivar algunas caracter´ısticas de las FRFs, las que se pueden utilizar para evaluar la calidad de los datos medidos. Esta evaluaci´ on puede detectar errores en las mediciones. La primera caracter´ıstica es la de una FRF directa, en donde el punto de excitaci´on es el mismo punto de medici´on. En este caso se espera siempre ver una antiresonancia
Frequency response function measurement
155
response function measurement 7.7.4 Special characteristicsFrequency of an FRF
155
7.7.4 of derive an FRF From theSpecial theory of characteristics modal analysis, we can some characteristics of an FRF and use them to assess the FRF data measured in modal testing. This assessment has the 88 derive some characteristics of an FRF and potential to detect measurement errorscan or mistakes made during the testing. useThe them to assess the FRF data measured in modal testing. This assessment has the first characteristic is that for a point FRF measurement we expect to see an potential to detect measurement errorsresonances. or mistakes made duringofthe anti-resonance between two adjacent The rationale thistesting. characteristic characteristic that for a point measurement wenot expect to see oan entre dos resonancias. loistanto, si esta caracter´ no se observa en la medici´ n hasThe beenfirst explained inPor Chapter 5. Therefore, ifFRF thisıstica characteristic is fully observed between two adjacent resonances. The oftransducers deanti-resonance una FRF measurement, directa, es muy probable quethe losforce sensores de fuerza ythis de characteristic respuesta no on a point it is likely that andrationale response are not has Chapter 5.as Therefore, if this is not observed est´ en been en realidad en in la misma coordinada. actually atexplained the same coordinate, they should be. characteristic Any minor offset offully the two could on a point measurement, it is likely that the force and response transducers are not degenerate some of the anti-resonances. Para unaa estructura empotrada, aasbajas frecuencias, la the caracter´ ıstica actually at the same coordinate, theyfrequency should be. Any minor offset ofpredominante the two could For grounded structure, at very low range, predominant characteristic some of the anti-resonances. dedegenerate la estructura es su rigidez est´ a tica. Por lo tanto, al comienzo de una FRF, of the structure is its static stiffness. Therefore, at the beginning of the FRF, we For aıa grounded structure, at very low frequency range, the predominant characteristic se should deber´ observar una linea de rigidez antes de la primera resonancia, como se see a stiffness line before the first resonance appears. On the other hand, for of the structure is its static stiffness. Therefore, at the beginning of the FRF, we muestra en la Figura 6.13. Por otro lado, para una estructura libre, la caracter´ ıstica a freely supported structure, the prevailing characteristic at very low frequency is the should seeinertia. a astiffness line before thelafirst resonance appears. the other for predominante bajas masa la inercia. significa que se FRF. debe mass and Thisfrecuencias means we es should see ay mass line atEsto the On beginning of hand, the aFigures freelyuna supported structure, prevailing characteristic very lowen frequency is the observar de masa althe inicio la FRF, como seatmuestra la Figura ??. 7.13linea and 7.14 illustrate thesedetwo cases. mass and inertia. This means we should see a mass line at the beginning of the FRF. Figures 7.13 and 7.14 illustrate these two cases.
Receptance FRF(dB) (dB)
´ INICIAL EVALUACI ON FRFS MEDIDASwe From the theory DE of LAS modal analysis,
0
Frecuencia (log)
Figure07.13 A point FRF of a grounded structure with the dotted stiffness line Frequency (log)
Figura 6.12: FRF directa de una estructura empotrada
Figure 7.13 A point FRF of a grounded structure with the dotted stiffness line –150
Receptance (dB) FRF (dB)
–200 –150 –250 –200 –300 –250 –350 –300 –400 –350 –450 –400 –450
0
1
2
3
4
5
Frequency (log)
0 1 2 4 mass line Figure 7.14 A point FRF of a free structure 3with the dotted Frecuencia (log)
Figure 7.14 A point FRF of a free structure with the dotted mass line
Figura 6.13: FRF directa de una estructura libre
5
Cap´ıtulo 7
Correlaci´ on Num´ erico-Experimental y Ajuste de Modelos En el dise˜ no de estructuras mec´anicas, el comportamiento din´ amico es un tema de vital importancia. La vida u ´ til bajo cargas c´ıclicas, niveles de vibraci´on o ruido, interacci´on entre sistemas de control y la vibraci´on de la estructura, . . . son restricciones relevantes en el dise˜ no. Sin embargo, el an´alisis din´amico de una estructura no es directo. Los par´ametros modales se pueden determinar de manera experimental o por m´etodos num´ericos. En general, los resultados experimentales entregan informaci´on de la estructura s´olo para la configuraci´on experimental. Un modelo en elementos finitos, en cambio, permite predecir el comportamiento din´amico de la estructura bajo distintas condiciones de borde y carga, pero la confiabilidad de un modelo en elementos finitos muchas veces no est´a garantizada. Los m´etodos de ajuste de modelos permiten verificar y corregir estos modelos en elementos finitos por medio de los datos experimentales. El resultado de un ajuste de modelo es un modelo en elementos finitos que es m´as confiable para predicciones futuras. El ajuste de un modelo num´erico busca desarrollar un modelo en elementos finitos que entregue predicciones confiables de la din´amica de una estructura mec´ anica. En la Figura 7.1 se ilustra el esquema general de un m´etodo de ajuste de modelos.
89
•
En los objetivos, agrega ajuste del modelo numérico (o matemático)…para tener un modelo validado de la suspensión del auto y a partir de este modelo validado realizar simulaciones y modificaciones.
´ NUMERICO-EXPERIMENTAL ´ CORRELACION Y AJUSTE DE MODELOS
•
90
En la parte de correlación experimental, yo agregaría un esquema como el siguiente:
Funciones de respuesta en frecuencia
Modelo numérico
Estructura real
Análisis modal
Test dinámico
Frecuencias naturales y modos propios
FRFs experimentales
Frecuencias naturales y modos propios
Correlación
buena
Fin
Ajuste del modelo numérico
Figura 7.1: Esquema general de un ajuste de modelo En la parte de avance del trabajo, del modelo matemático también obtuviste las frecuencias y modos de vibración. Pon esas frecuencias, y haz un esquema de los tres El procedimiento comienza con la construcci´ on de un modelo en elementos finitos. modos de vibración. •
La estructura se divide en elementos conectados por nodos. Cada nodo tiene uno o m´as• grados de libertad. Los grados libertad representan Que diferencia existe entre el modelode "computacional" y el modelodesplazamientos matemático? Dado y deformaciones de la estructura en forma discretizada. modelo est´ a representado que el modelo computacional considera al auto comoEste un cuerpo rígido, deberían ser el por tres mismo matrices: de rigidez K, masa M y amortiguaci´on C. A partir de estas modelo... matrices se pueden determinar las propiedades modales de la estructura (ver secci´on 1.2). Por otro lado, se realiza un an´ alisis modal experimental de la estructura. En donde se miden las funciones de respuesta en frecuencia en distintos puntos de la estructura. Luego los par´ametros modales se pueden identificar a partir de un m´etodo de identificaci´ on de par´ametros.
´ NUMERICO-EXPERIMENTAL ´ CORRELACION Y AJUSTE DE MODELOS
91
El ajuste de modelos comienza primero al parear los puntos de medici´on experimentales con nodos del modelo num´erico. Generalmente, los puntos experimentales no coinciden completamente con los nodos num´ericos. Primero, los nodos num´ericos y los puntos experimental pueden diferir en su ubicaci´ on f´ısica. Esto se puede solucionar con una buena comunicaci´ on entre la persona que realiza el modelo num´erico y el que realiza las mediciones experimentales. Segundo, y m´as importante, un modelo en elementos finitos tiene muchos mas grados de libertad que los que se miden experimentalmente. Para algunas t´ecnicas de correlaci´on es necesario que los grados medidos experimentalmente y los grados del modelo num´erico coincidan. Para resolver esta incompatibilidad con las mallas, se puede o reducir el modelo num´erico a los grados de libertad medidos experimentalmente, o expandir los datos experimentales al n´ umero de grados de libertad num´ericos. La reducci´on de matrices o expansi´on datos experimentales es un obst´aculo importante en ajuste de modelos. Luego de coincidir ambos modelos, el procedimiento continua con la correlaci´on entre los datos num´ericos y experimentales. Existen distintas t´ecnicas de correlaci´on que se pueden utilizar. Si la correlaci´on es buena el algoritmo termina, y el modelo en elementos finitos se considera suficientemente bueno. En caso contrario (y m´as probable), si la correlaci´on es mala, el modelo en elementos finitos debe ser corregido por medio de un procedimiento de ajuste. Los datos num´ericos con los experimentales muchas veces no son compatibles por alguna de las razones siguientes: Los grados libertad experimentales no coinciden con los grados de libertad del modelo en elementos finitos: Esto se puede solucionar parcialmente con algoritmos de reduci´on/expansi´on. El conjunto de datos experimentales es incompleto: No s´olo el n´ umero de grados experimentales esta limitado, sino tambi´en el n´ umero de modos que se pueden determinar experimentalmente. La medici´on de las funciones de respuesta en frecuencia se realiza en rango de frecuencias limitado. Por lo tanto, los modos fuera de este rango no pueden ser identificados. Dependiendo de la cantidad de informaci´on experimental, los resultados del ajuste de modelos pueden no ser u ´nicos. Los datos experimentales est´ an contaminados con ruido: En general, se esperan errores del orden del 3 % en las frecuencias naturales y del 10 % en las formas modales. El uso de informaci´ on experimental contaminada con ruido puede llevar a una ajuste de modelos inadecuado. El amortiguamiento no se puede incluir de manera precisa en el modelo num´ erico: La informaci´on de amortiguamiento est´a inherentemente presente en los datos experimentales, pero usualmente no se considera en el
´ PAREAR MODELOS NUMERICOS Y EXPERIMENTALES
92
modelo en elementos finitos. Esta discrepancia puede causar errores en el modelo ajustado.
7.1.
Parear modelos num´ ericos y experimentales
En la mayor´ıa de los casos de n´ umero de grados de libertad del modelo num´erico excede por mucho el n´ umero de grados de libertad medidos. Los modelos en elementos finitos requieren mallas finas para poder entregar resultados precisos. No es practico , y muchas veces no es posible, medir todos los grados de libertad en la estructura real: Muchos de los grados de libertad son internos y no se puede acceder a ellos para realizar mediciones. Los grados de libertad de rotaci´on son dif´ıciles de medir. Para el prop´osito de an´alisis modal, no es necesaria una malla muy fina de puntos de medici´ on. Sin embargo, muchos algoritmos de ajuste de modelos requieren una correspondencia uno-a-uno entre los grados de libertad anal´ıticos y experimentales. Para lograr esto se utilizan algoritmos de reducci´ on de modelos num´ericos y de expansi´on de los datos experimentales.
7.1.1.
T´ ecnicas de reducci´ on
Las t´ecnicas de reducci´on expresan las ecuaciones de movimiento en t´erminos de los grados de libertad medidos experimentalmente. En una primera fase, se identifican los grados de libertad anal´ıticos para cada grado de libertad experimental. Esto define a los grados de libertad “activos”. El resto de los grados de libertad anal´ıticos se denominan grados de libertad “eliminados”:
Xf =
Xa Xd
donde: Xa = grados de libertad activos Xd = grados de libertad eliminados Xf = set completo de grados de libertad
(7.1)
´ PAREAR MODELOS NUMERICOS Y EXPERIMENTALES
93
Las t´ecnicas de reducci´on establecen una relaci´on entre los grados de libertad activos y los grados de libertad eliminados por medio de una matriz de transformaci´on Td o Tf : Xd
=
Td Xa
(7.2)
Xf
=
(7.3)
Tf
=
Tf Xa I Td
(7.4)
Tf se utiliza para construir matrices de masa y rigidez reducidas Mr y Kr respectivamente. Para un sistema sin amortiguamiento se cumple la siguiente relaci´ on: XfT M Xf
=
XaT Mr Xa
(7.5)
XfT KXf
=
XaT Kr Xa
(7.6)
Sustituyendo Xf por Tf Xa se obtiene: Mr
= TfT M Tf
(7.7)
Kr
= TfT KTf
(7.8)
Las matrices del sistema reducido entregan s´ olo una descripci´ on aproximada del comportamiento din´ amico “exacto” del sistema con las matrices completas. Se sabe que todos los grados de libertad deben cumplir con la ecuaci´on de movimiento: K − ω 2 M Xf = Z(ω)Xf = Ff (7.9) Para el sistema reducido no deben haber fuerzas externas en los grados de libertad eliminados. Dividiendo el sistema de ecuaciones en los grados de libertad activos y los eliminados: Zaa Zad Xa Fa = (7.10) Zda Zdd Xd 0 despejando Xd se obtiene que: −1 Xd = −Zdd Zda Xa
(7.11)
´ PAREAR MODELOS NUMERICOS Y EXPERIMENTALES
94
De donde se deduce que: Td = −Zdd (ωr )−1 Zda (ωr )
(7.12)
Se debe notar que Z est´a evaluada a una frecuencia ωr , en consecuencia las matrices Mr y Kr dar´an resultados exactos para esta frecuencia. Para las otras frecuencias Mr y Kr dan s´olo resultados aproximados, siendo peores al alejarse de ωr . Un caso particular es la bien conocida reducci´on de Guyan, en donde se utiliza ωr = 0.
7.1.2.
T´ ecnicas de expansi´ on
las t´ecnicas de expansi´ on expanden los modos experimentales desde el set de grados de libertad activos al set de grados de libertad completo del modelo anal´ıtico. La evaluaci´on de los grados de libertad no medidos en los modos experimentales se lleva a cabo la mayor´ıa de las veces utilizando las relaciones entre los grados activos y eliminados del modelo num´erico. El m´etodo de reducci´ on presentado en la secci´on anterior se puede utilizar tambi´en como un m´etodo de expansi´on ya que define la relaci´on entre los grados de libertad activos y eliminados en el modelo anal´ıtico: Xd = Td Xa
(7.13)
La misma matriz de transformaci´ on se puede utilizar para los modos experimentales: φde = Td φae
(7.14)
Donde φde representa los grados de libertad eliminados de un modo (o matriz de modos) experimental y φae los grados de libertad activos. A continuaci´ on se presentan m´etodos de expansi´on alternativos. Mezcla de vectores modales Este m´etodo simplemente completa los grados de libertad faltantes en los modos experimental con los valores correspondientes de los modos anal´ıticos: a φe f φe = (7.15) φda φfe representa un modo experimental (o matriz de modos) con el set completo de grados de libertad, φda es el modo anal´ıtico en los grados de libertad no medidos y φae es el modo experimental en los grados de libertad medidos.
´ ´ TECNICAS DE CORRELACION
95
Antes de completar el modos experimental con los valores del modo anal´ıtico, ´estos deben tener la misma escala de modo que: kφae k = kφaa k
(7.16)
M´ etodo de coordenadas modales Este m´etodo define los modos experimentales como una combinaci´ on lineal de los modos anal´ıticos. Los coeficientes de la combinaci´on lineal son calculados con los grados de libertad activos: φae
=
φaa q
(7.17)
q
=
a φa+ a φe
(7.18)
La matriz q se utiliza luego para evaluar los grados de libertad experimental faltantes: φde = φda q
(7.19)
Interpolaci´ on Los grados de libertad faltantes tambi´en se pueden aproximar por medio de m´etodos de interpolaci´ on. Una ventaja de este m´etodo es que no se necesitan datos num´ericos para expandir los datos experimentales. Este m´etodo depende fuertemente de la conexi´on entre los distintos grados de libertad. En consecuencia, cada interpolaci´on depende mucho del caso en estudio, lo que dificulta su uso pr´actico. En la pr´actica, estos m´etodos se aplican en estructuras que consisten en componentes unidimensionales. Tambi´en se pueden utilizar para estimar grados de libertad de rotaci´on a partir de los grados de libertad de traslaci´on, pero este procedimiento es num´ericamente poco estable.
7.2.
T´ ecnicas de correlaci´ on
El prop´osito de los m´etodos de ajuste de modelo es lograr que el modelo num´erico se acerque lo m´as posible a los datos experimentales. Para evaluar el nivel de correlaci´on entre los datos num´ericos y experimentales es necesario definir que par´ametros se van a comparar y en que forma. A continuaci´on, se describir´an algunas t´ecnicas usuales de correlaci´on num´erica-experimental.
´ ´ TECNICAS DE CORRELACION
7.2.1.
96
Diferencia en las frecuencias de resonancia
La manera m´as sencilla de verificar la correlaci´on es comparando las frecuencias naturales anal´ıticas y experimentales. La diferencia m´axima permitida depende en la precisi´ on de las frecuencias naturales.
7.2.2.
Comparaci´ on visual de los modos
Modo experimental
Una manera de visualizar la correlaci´ on entre dos modos son los denominados gr´aficos de 45◦ . En estos gr´aficos se gr´ afica un modo experimental versus el modo anal´ıtico correspondiente en un gr´afico x-y. Si ambos modos son id´enticos y tienen la misma escala, entonces todos los puntos se encuentran sobre una linea de 45◦ que parte desde el origen. La distancia entre los puntos y la l´ınea de 45◦ da una indicaci´on de la correlaci´on entre los modos. En la Figura 7.2 se ilustra un gr´afico de 45◦ t´ıpico.
Modo analítico Figura 7.2: Gr´afico de 45◦ Una manera de asegurar que los modos experimentales tengan la misma escala que los modos anal´ıticos es utilizando el “Factor de Escala Modal (MSF)”. El MSF mide el factor de escala entre dos modos, los modos experimentales se pueden escalar a los modos anal´ıticos multiplicados por el MSF correspondiente: φ∗e,i M SFi
= φe,i · M SFi
(7.20)
φTa,i φe,i φTe,i φe,i
(7.21)
=
´ ´ TECNICAS DE CORRELACION
97
Donde φ∗e,i y φe,i es el i´esimo modo experimental escalado y original, φa,i es el i´esimo modo anal´ıtico, y M SFi es el factor de escala modal para el modo i. Al multiplicar el modo experimental por el MSF correspondiente, tambi´en se soluciona el problema que los modos anal´ıtico y experimental pueden estar con un desfase de 180◦ .
7.2.3.
Modal Assurance Criterion
El “Modal assurance criterion (MAC)” expresa la correlaci´on que es visualizada en un gr´ afico de 45◦ en un s´ olo n´ umero. Se define como: 2 φTa,i φe,j M ACij = T φa,i φa,i φTe,j φe,j
(7.22)
Donde φa,i es el iesimo modo anal´ıtico y φe,j es el jesimo modo experimental. Un valor de 0 indica que no hay correlaci´on, mientras que un valor de 1 indica dos modos perfectamente correlacionados. Al ordenar todos los valores M ACij en una matriz, la diagonal deber´ıa tener valores altos (en general mayores a 0.8) para una buena correlaci´ on. Una ventaja del MAC es que la correlaci´on no depende de la escala de los modos, sino que s´olo de forma de ´estos. En la Figura 7.3 se representa gr´ aficamente una matriz de valores MAC t´ıpica. 1
116 0.9
Modos Experimentales (Hz)
58.9 56.6
0.8
44.7
0.7
39
0.6
31.2 0.5
28.4 0.4
26.3 24.2
0.3
19
0.2
15.4
0.1
13 12.7 15.2 19.4 24.2 26.3 28.3 31.4 38.4 43.8
57 58.7 116
Modos numéricos (Hz)
Figura 7.3: Matriz de valores MAC
´ METODOS ITERATIVOS DE AJUSTE DE MODELOS
7.2.4.
98
Frequency Response Assurance Criterion
Tal como el MAC mide la correlaci´on de dos vectores modales, el “Frequency Response Assurance Criterion (FRAC)” mide la correlaci´ on entre dos funciones de respuesta en frecuencia. Se define de manera equivalente al MAC como:
2 a e ω=ω1 Hpq (ω)Hpq (ω) Pω2 Pω2 e e a a ω=ω1 Hpq (ω)Hpq (ω) ω=ω1 Hpq (ω)Hpq (ω)
Pω2
F RACpq =
(7.23)
Donde Hpq es la funci´on de respuesta en frecuencia cuando la excitaci´on es en el punto p y se mide en el punto q. Los indices a y e se refieren a anal´ıtico y experimental respectivamente. Por u ´ ltimo, ω1 y ω2 es el rango de frecuencias en que se quiere correlacionar las FRFs. Como en el caso del MAC, un valor del FRAC igual a 1 equivale a dos FRFs perfectamente correlacionadas, mientras que un valor cercano a 0 indica una mala correlaci´ on.
7.3.
M´ etodos iterativos de ajuste de modelos
Los m´etodos m´as utilizados en ajuste de modelo son m´etodos iterativos basado en gradientes. Estos m´etodos minimizan una funci´ on error ε que relaciona los par´ametros experimentales con los entregados por el modelo num´erico. Esta funci´on de error puede estar compuesta de distintos residuos: La diferencia entre las frecuencias naturales anal´ıticas y experimentales. La diferencia entre los modos anal´ıticos y experimentales. El balance de fuerzas a una frecuencia especifica. La diferencia entre las funciones de respuesta en frecuencia anal´ıticas y experimentales. El problema se resuelve con una aproximaci´on de Taylor de primer orden
ε(β) = ε(β0 ) +
X ∂ε ∆βi ∂βi i
(7.24)
´ METODOS ITERATIVOS DE AJUSTE DE MODELOS
99
Donde β = {β1 , β2 , . . . βm } es un vector de par´ametros a ser actualizados. ∆βi = βij+1 − βij , βij es el valor actual del par´ametro y βij+1 es el valor estimado. Las derivadas parciales de primer orden son evaluadas num´ericamente. Definiendo una matriz S que contenga las derivadas parciales del residuo con respecto a los par´ ametros β: ∂ε1 ∂ε1 ∂ε1 ... ∂β1 ∂β2 ∂βn ∂ε ∂ε2 ∂ε2 2 ... ∂β ∂β ∂β 1 2 n (7.25) S= .. .. ... . ... . ∂εn ∂εn ∂εn ... ∂β1 ∂β2 ∂βn la ecuaci´ on 7.24 se puede escribir como: ∆ε = S∆β,
∆ε = zE − zA (β)
(7.26)
Donde zE and zA (β) son los datos experimentales y anal´ıticos respectivamente, ∆ε es el vector de residuos. El vector de par´ametros β es estimado al minimizar la siguiente funci´on objetivo, T
J(β) = (S∆β − ∆ε) (S∆β − ∆ε)
(7.27)
Este problema de minimizaci´on se resuelve por m´ınimos cuadrados. La soluci´on viene dada por, ∆β = S + ∆ε Donde S + es la pseudoinversa S −1 + T (S S)−1 S T Sn×m = T T −1 S (S S)
(7.28) de la matriz S: n=m n>m n