M2 · M2.4 · Sesión S2
La transformada de Laplace convierte una ecuación diferencial en una ecuación algebraica. Lo que sale de ahí, la función de transferencia, resume cómo un sistema lineal responde a cualquier entrada, y tiene dos tipos de números que lo caracterizan: los polos, que son propiedades del sistema (frecuencias y amortiguamientos de los modos), y los ceros, que dependen de dónde se excita y dónde se mide. En identificación modal, los polos son lo que se quiere estimar, y los ceros explican por qué dos sensores del mismo edificio ven la misma estructura con formas tan distintas.
El sistema de un grado de libertad (polos, residuos, respuesta impulsional) está desarrollado en M5. Esta página aporta lo que falta para varios grados de libertad: la matriz de transferencia, el efecto de los sensores y la conexión con la identificación.
Tiempo orientativo: unas 3 horas. Nivel: intermedio. Tres capas: intuición (mapa de polos y ceros, simulación), fundamento (definición, propiedades y teoremas) y aplicación (el edificio de tres plantas, caso transversal).
Objetivos. Al terminar esta página sabrás:
Prerrequisitos. M2.2 (modos) y M5 (oscilador de un grado de libertad); J3.6 para las ecuaciones diferenciales.
Definición 1 Transformada de Laplace. Para una señal \(x(t)\) con \(x(t)=0\) si \(t<0\), \(X(s)=\mathcal L[x](s)=\displaystyle\int_{0}^{\infty}x(t)\,e^{-st}\,\mathrm dt\), con \(s=\sigma+\mathrm j\omega\) complejo, para los \(s\) en los que la integral converge (región de convergencia).
La transformada «pesa» la señal con una exponencial \(e^{-st}\): la parte real \(\sigma\) amortigua las señales que crecen, y la imaginaria \(\omega\) las proyecta sobre oscilaciones de frecuencia \(\omega\). La transformada de Fourier (M2.5) es el caso \(\sigma=0\) (cuando existe).
| Señal \(x(t)\), \(t\ge0\) | \(X(s)\) |
|---|---|
| \(\delta(t)\) (impulso) | \(1\) |
| \(1\) (escalón) | \(1/s\) |
| \(e^{-at}\) | \(1/(s+a)\) |
| \(\sin\omega t\) | \(\omega/(s^2+\omega^2)\) |
| \(e^{-\sigma_0 t}\sin\omega_dt\) | \(\omega_d/((s+\sigma_0)^2+\omega_d^2)\) |
Propiedades (con condiciones iniciales nulas): linealidad; derivación, \(\mathcal L[\dot x]=sX(s)-x(0)\) y \(\mathcal L[\ddot x]=s^2X(s)-sx(0)-\dot x(0)\); convolución, \(\mathcal L[x*h]=X(s)H(s)\); y valor final, si todos los polos de \(sX(s)\) están en el semiplano izquierdo, \(\lim_{t\to\infty}x(t)=\lim_{s\to0}sX(s)\).
Derivación. Integrando por partes, \(\int_0^\infty\dot x\,e^{-st}\mathrm dt=\big[x\,e^{-st}\big]_0^\infty+s\int_0^\infty x\,e^{-st}\mathrm dt=-x(0)+sX(s)\) (el término en \(\infty\) se anula en la región de convergencia). Aplicando la fórmula dos veces se obtiene la de \(\ddot x\).
Convolución. \(\mathcal L[x*h]=\int_0^\infty\!\!\int_0^t x(\tau)h(t-\tau)\,\mathrm d\tau\,e^{-st}\mathrm dt\). Cambiando el orden de integración y poniendo \(u=t-\tau\), queda \(\int_0^\infty x(\tau)e^{-s\tau}\mathrm d\tau\int_0^\infty h(u)e^{-su}\mathrm du=X(s)H(s)\). \(\square\)
Para \(\mathbf M\ddot{\mathbf u}+\mathbf C\dot{\mathbf u}+\mathbf K\mathbf u=\mathbf B\,\mathbf f(t)\) con condiciones iniciales nulas, la transformada es \((s^2\mathbf M+s\mathbf C+\mathbf K)\mathbf U(s)=\mathbf B\mathbf F(s)\). Si se mide una combinación lineal de desplazamientos, \(\mathbf y=\mathbf L\mathbf u\), la matriz de transferencia es
\[ \mathbf H(s)=\mathbf L\,(s^2\mathbf M+s\mathbf C+\mathbf K)^{-1}\mathbf B. \tag{1}\]
Sus elementos son cocientes de polinomios en \(s\): \(H_{ij}(s)=N_{ij}(s)/D(s)\). Las raíces del denominador común \(D(s)=\det(s^2\mathbf M+s\mathbf C+\mathbf K)\) son los polos; las raíces del numerador de cada elemento, los ceros.
Descomposición modal. Con amortiguamiento proporcional y los modos \(\boldsymbol\varphi_k\) normalizados a masa unidad (M2.2), cada elemento de la receptancia es una suma de \(n\) osciladores:
\[ H_{ij}(s)=\sum_{k=1}^{n}\frac{\varphi_{ik}\,\varphi_{jk}}{s^2+2\zeta_k\omega_k s+\omega_k^2}, \qquad \zeta_k=\frac{\boldsymbol\varphi_k^{\mathsf T}\mathbf C\boldsymbol\varphi_k}{2\omega_k}. \tag{2}\]
Es la fórmula que une toda la unidad: los polos \(s_k=-\zeta_k\omega_k\pm\mathrm j\omega_k\sqrt{1-\zeta_k^2}\) dependen solo de la estructura; la forma modal solo entra como el residuo \(\varphi_{ik}\varphi_{jk}\) de cada polo, que depende de dónde se excita (\(j\)) y de dónde se mide (\(i\)). De ahí sale que, midiendo con sensores en distintos puntos, las frecuencias son las mismas y las formas se obtienen de los residuos.
ω_prueba = 2π * 2.0
Hm = sum(Phi3[:, k] * Phi3[:, k]' / ((2π * f3[k])^2 - ω_prueba^2 + 2im * zeta3[k] * 2π * f3[k] * ω_prueba) for k in 1:3)
err_modal = maximum(abs.(Hm - Hrec(ω_prueba))) / maximum(abs.(Hrec(ω_prueba)))1.054436843031758e-15La suma modal reproduce la inversión directa de \(\mathbf K-\omega^2\mathbf M+\mathrm j\omega\mathbf C\) en \(f=2\) Hz con un error relativo de 1e-15 (redondeo): la ecuación 2 es una identidad exacta cuando el amortiguamiento es proporcional.
Con la forma de estado \(\dot{\mathbf x}=\mathbf A\mathbf x+\mathbf B\mathbf u\) (espacio de estados), los polos de \(\mathbf H(s)\) son valores propios de \(\mathbf A\). Un sistema lineal e invariante es asintóticamente estable si todos los polos están en el semiplano izquierdo abierto, \(\operatorname{Re}(s_k)<0\): las respuestas libres decaen. Un polo en el eje imaginario (sin amortiguamiento) da oscilaciones que no decaen, y un polo con parte real positiva, una respuesta que crece. Una estructura real tiene siempre amortiguamiento positivo, así que si un procedimiento de identificación devuelve polos en el semiplano derecho, lo que se ha identificado no es una estructura estable: es un polo espurio o un error.
Que todos los valores propios de \(\mathbf A\) sean polos de un \(H_{ij}(s)\) concreto no está garantizado: si la entrada no excita un modo o la salida no lo observa, ese polo se cancela con un cero. Es el problema de controlabilidad y observabilidad, que se desarrolla en la página de espacio de estados.
Los ceros de \(H_{ij}(s)\) son las frecuencias a las que, con entrada en \(j\), la salida en \(i\) se anula. Hay un hecho clásico, comprobable con el edificio: la receptancia en un punto (\(i=j\), sensor y excitación en el mismo grado de libertad) tiene un cero entre cada dos polos consecutivos, y esos ceros (sin amortiguamiento) son las frecuencias naturales de la estructura con ese grado de libertad bloqueado. Los ceros de \(H_{ij}\) con \(i\ne j\) son otros y dependen del par. En una cadena tridiagonal como el edificio de cortante, el elemento \((i,j)\) de la inversa se calcula con los menores principales de los extremos; cuando sensor y fuerza están en los dos extremos (plantas 1 y 3), el numerador es una constante y \(H_{13}\) no tiene ceros finitos: solo tiene valles poco profundos entre resonancias, como cualquier función con picos próximos.
f_anti = modos_propios(K3[1:2, 1:2], M3[1:2, 1:2])[1] # planta 3 bloqueada
fs = range(0.5, 5.5; length=200_001)
mag33 = [abs(Hrec(2π * x)[3, 3]) for x in fs]
i_min = [k for k in 2:length(fs)-1 if mag33[k] < mag33[k-1] && mag33[k] < mag33[k+1]]
(frecuencias_con_planta_3_bloqueada = f_anti, minimos_de_la_FRF_H33 = fs[i_min])(frecuencias_con_planta_3_bloqueada = [2.3806999966057445, 4.289061510542444], minimos_de_la_FRF_H33 = [2.38025, 4.2688])Con la planta 3 bloqueada, el edificio tiene frecuencias de 2,381 y 4,289 Hz; la FRF puntual en la cubierta tiene mínimos en 2,380 y 4,269 Hz. La primera coincide dentro del \(0{,}1\,\%\); la segunda se aparta un \(0{,}5\,\%\) porque el amortiguamiento y la cercanía del tercer polo desplazan el mínimo de módulo (con amortiguamiento, el cero exacto sale del eje imaginario).
modulo(i, j) = [abs(Hrec(2π * x)[i, j]) for x in fs]
function valles(i, j)
v = modulo(i, j)
m = [k for k in 2:length(v)-1 if v[k] < v[k-1] && v[k] < v[k+1]]
p = [k for k in 2:length(v)-1 if v[k] > v[k-1] && v[k] > v[k+1]]
[(fs[k], 20log10(v[k] / min(maximum(v[q] for q in p if fs[q] < fs[k]), maximum(v[q] for q in p if fs[q] > fs[k])))) for k in m]
end
sub = Dict(1 => "₁", 2 => "₂", 3 => "₃")
filas = String[]
for (i, j) in [(1, 1), (2, 2), (3, 3), (1, 2), (2, 3), (1, 3)]
vs = valles(i, j)
push!(filas, "| H" * sub[i] * sub[j] * " | " * join(["$(coma(f; d=2)) Hz: $(coma(d; d=1)) dB" for (f, d) in vs], " · ") * " |")
end
Markdown.parse("| Elemento | Valles (frecuencia y profundidad) |\n|:--|:--|\n" * join(filas, "\n"))| Elemento | Valles (frecuencia y profundidad) |
|---|---|
| H₁₁ | 1,50 Hz: -55,6 dB · 3,77 Hz: -38,8 dB |
| H₂₂ | 2,25 Hz: -49,4 dB · 3,56 Hz: -43,2 dB |
| H₃₃ | 2,38 Hz: -49,1 dB · 4,27 Hz: -13,5 dB |
| H₁₂ | 2,25 Hz: -54,9 dB · 3,67 Hz: -15,9 dB |
| H₂₃ | 2,32 Hz: -17,7 dB · 3,56 Hz: -37,7 dB |
| H₁₃ | 2,11 Hz: -23,4 dB · 3,87 Hz: -12,1 dB |
Los valles de las FRF puntuales (\(H_{11}\), \(H_{22}\)) son entre 39 y 56 dB más bajos que los picos vecinos; los de la no puntual de los extremos (\(H_{13}\)), solo 12–23 dB. \(H_{12}\) tiene un cero cerca del eje a 2,25 Hz (la frecuencia \(\sqrt{k_3/m}/2\pi\)) y un valle profundo allí. La segunda antirresonancia de \(H_{33}\) es poco profunda (13 dB) porque su cero queda muy cerca del tercer polo (4,29 Hz frente a 4,39 Hz). La profundidad no es una propiedad de la estructura: depende del par entrada–salida.
A3, _ = espacio_estados(M3, C3, K3, [0.0, 0.0, 1.0])
s_polos = filter(p -> imag(p) > 0, eigvals(A3))
ωz = 2π .* f_anti
fig = Figure(size=(560, 420))
ax = Axis(fig[1, 1]; xlabel="parte real, σ (rad/s)", ylabel="parte imaginaria, ω (rad/s)", limits=(-1.0, 0.4, -32, 32))
vlines!(ax, [0.0]; color=GRIS)
scatter!(ax, vcat(real.(s_polos), real.(s_polos)), vcat(imag.(s_polos), -imag.(s_polos)); marker=:xcross, color=NARANJA, markersize=16, label="polos")
scatter!(ax, zeros(4), vcat(ωz, -ωz); marker=:circle, color=:white, markersize=11, strokewidth=2, strokecolor=AZUL, label="ceros (sin amortiguamiento)")
axislegend(ax; position=:rb)
fig
Los polos que se identifican de un registro de aceleraciones son los mismos con cualquier sensor. Los ceros no: cambian con la posición del sensor y de la excitación, y por eso las antirresonancias no son propiedades de la estructura sino del experimento. Un sensor colocado en un nodo de un modo (residuo nulo) no ve ese polo; y un par polo–cero muy próximo (que casi se cancelan) lleva un residuo pequeño y es difícil de identificar.
Sobre el eje imaginario, \(\mathbf H(\mathrm j\omega)\) es la respuesta en frecuencia (FRF): la amplitud y la fase con que el sistema transforma una excitación armónica de frecuencia \(\omega\) en régimen estacionario. La respuesta impulsional \(\mathbf h(t)=\mathcal L^{-1}[\mathbf H(s)]\) se obtiene por fracciones simples: cada par de polos aporta un término \(\frac{\varphi_{ik}\varphi_{jk}}{\omega_{dk}}e^{-\zeta_k\omega_kt}\sin\omega_{dk}t\). Ambas contienen la misma información; la FRF es la que se mide con excitación controlada, y la respuesta impulsional (o sus covarianzas) es la que alimenta a ERA y a SSI (M2.1.5).
La conexión con la medición ambiental es que, si no se conoce la fuerza, no se puede calcular \(\mathbf H(\mathrm j\omega)\) como cociente salida/entrada. Lo que sí se tiene es la densidad espectral de las salidas, \(\mathbf S_{yy}(\omega)=\mathbf H(\mathrm j\omega)\,\mathbf S_{ff}\,\mathbf H(\mathrm j\omega)^{\mathsf H}\), y si la entrada es ruido blanco (\(\mathbf S_{ff}\) constante), la forma de \(\mathbf S_{yy}\) la fija \(\mathbf H\) (M2.5 y M7). Es la base de FDD.
Simulación 1. FRF de un edificio de tres plantas. Elige la planta donde se excita y la planta donde se mide, y el amortiguamiento modal (el mismo para los tres modos). Se dibuja el módulo de la receptancia \(H_{ij}(\mathrm j\omega)\) en decibelios (con \(1\) m/N de referencia) y su fase, con la frecuencia natural de cada modo marcada (línea vertical gris).
Para explorar. (1) Con \(i=j=3\) cuenta las antirresonancias y compáralas con las frecuencias de la estructura con la planta 3 bloqueada (cálculo de la sección M2.4.4). (2) Con \(i=1\), \(j=3\) observa que los valles son mucho menos profundos y la fase no sube. (3) Sube \(\zeta\) al 5 %: ¿qué pasa con la profundidad de las antirresonancias?
espacio_estados) y con polos_modales; comprueba que dan las mismas frecuencias y los amortiguamientos del 1 %.using Resonancia, LinearAlgebra
M, C, K = edificio_tres_plantas()
A, B = espacio_estados(M, C, K, [0.0, 0.0, 1.0])
f, zeta = polos_modales(A) # polos de la planta 3
H(ω) = inv(K - ω^2 * M + im * ω * C)
f_anti = modos_propios(K[1:2, 1:2], M[1:2, 1:2])[1] # planta 3 bloqueadaimport numpy as np
from scipy import linalg
import matrices as mx # mse/m02-matematicas/python/matrices.py (probado con pytest)
M, K = mx.edificio_tres_plantas("A")
H = lambda w, C: np.linalg.inv(K - w**2 * M + 1j * w * C) # C: matriz de amortiguamiento
f_anti = np.sqrt(linalg.eigh(K[:2, :2], M[:2, :2], eigvals_only=True)) / (2 * np.pi)El código Python no se ejecuta en la web: su versión probada está en python/matrices.py.
«Con Julia y el paquete Resonancia: para el edificio de tres plantas (edificio_tres_plantas), calcula la receptancia \(H_{ij}(\mathrm j\omega)\) por inversión directa para \(i,j\in\{1,2,3\}\), dibuja el módulo de las nueve en dB entre 0,5 y 5,5 Hz y dime en cuáles hay antirresonancias profundas. Verifica con la suma modal. No generalices más allá de este edificio.»
Verifica. (a) Que construye \(\mathbf K-\omega^2\mathbf M+\mathrm j\omega\mathbf C\) y la invierte, en vez de dividir elemento a elemento. (b) Que la matriz es simétrica (\(H_{ij}=H_{ji}\), reciprocidad). (c) Que la suma modal coincide con la inversión (error \(\sim10^{-15}\)). (d) Que cuantifica la profundidad de los valles en lugar de dar un «sí/no» (en este edificio hay valles profundos fuera de la diagonal, como \(H_{12}\)). (e) Que no atribuye las antirresonancias a «daño» ni a propiedades de la estructura.
Todos los valores propios de la matriz de estado \(\mathbf A\) son polos de cualquier función de transferencia \(H_{ij}(s)\) del sistema.
Falso. Si el modo no es controlable desde la entrada \(j\) o no es observable desde la salida \(i\), su polo se cancela con un cero y no aparece en \(H_{ij}\). Un sensor en un nodo de un modo no lo ve.
Un sistema tiene \(H(s)=1/(s^2+2s+26)\). ¿Cuánto vale su amortiguamiento relativo \(\zeta\)?
Compara con \(s^2+2\zeta\omega_ns+\omega_n^2\).
\(\omega_n^2=26\) (\(\omega_n=5{,}099\) rad/s) y \(2\zeta\omega_n=2\), luego \(\zeta=1/\sqrt{26}=0{,}196\). Los polos son \(-1\pm5\mathrm j\).
Un oscilador con \(m=2\) kg, \(c=3\) N·s/m y \(k=50\) N/m recibe un escalón de fuerza de 1 N. ¿A qué valor tiende el desplazamiento?
Teorema del valor final con \(X(s)=H(s)/s\).
\(x(\infty)=\lim_{s\to0}sH(s)/s=H(0)=1/k=0{,}02\) m. Los polos son estables (\(c>0\)), así que el teorema es aplicable.
En el edificio de tres plantas, ¿en qué frecuencia (la primera) tiene un mínimo profundo la FRF puntual de la cubierta \(H_{33}\)?
Es una frecuencia natural del edificio con la planta 3 bloqueada: elimina la fila y la columna 3 de \(\mathbf K\) y \(\mathbf M\).
Con la planta 3 bloqueada quedan dos plantas, con frecuencias de 2,381 Hz y 4,289 Hz; el primer valle de \(H_{33}\) está en 2,381 Hz.
Dos sensores colocados en plantas distintas del mismo edificio identifican las frecuencias naturales. ¿Qué es cierto?
La c. Los polos son propiedades del sistema; los ceros dependen de la entrada y de la salida (ecuación 1).
Un procedimiento de identificación devuelve un polo con parte real positiva para un puente en servicio. ¿Qué puede significar? ¿Qué harías?
Es un polo espurio (o un artefacto numérico): una estructura estable y amortiguada no puede tenerlo. Se descarta, y conviene revisar el orden del modelo, la calidad de los datos y los parámetros del algoritmo (M7).