El Método de la Rigidez (o Método Directo de Rigidez) es la base del análisis matricial moderno de estructuras. Su objetivo principal es convertir cualquier estructura física —ya sea una barra simple, una viga, una armadura o un marco complejo— en un sistema de ecuaciones algebraicas lineales que se puede resolver computacionalmente.
Todo el método gira en torno a esta relación fundamental:
$$\mathbf{K} \cdot \mathbf{u} = \mathbf{f}$$
Desglosemos sus componentes:
\(\mathbf{u}\) (Vector de Desplazamientos): Contiene todos los desplazamientos y rotaciones nodales desconocidos de la estructura. Estas son las incógnitas primarias que necesitamos encontrar. Para un nodo de una armadura 2D, podría ser \((u_x, u_y)\); para un nodo de viga, podría incluir la rotación \((v, \theta)\).
\(\mathbf{f}\) (Vector de Fuerzas): Contiene todas las fuerzas y momentos externos aplicados directamente en los nodos de la estructura. Estas son las cargas conocidas que actúan sobre el sistema.
\(\mathbf{K}\) (Matriz de Rigidez Global): Esta matriz cuadrada representa la rigidez general de toda la estructura. Relaciona linealmente las fuerzas aplicadas con los desplazamientos resultantes. Piénsala como la “huella digital” de la estructura: encapsula su geometría, propiedades del material y conectividad.
Analogía: Esta ecuación matricial es una versión generalizada y multidimensional de la simple Ley de Hooke para un resorte (\(F = kx\)). Así como la rigidez \(k\) de un resorte relaciona la fuerza \(F\) con el desplazamiento \(x\), la matriz \(\mathbf{K}\) relaciona el vector de fuerzas \(\mathbf{f}\) con el vector de desplazamientos \(\mathbf{u}\).
El método asume que la estructura está en equilibrio estático. Esto significa que en cada nodo (unión) de la estructura, la suma de todas las fuerzas externas aplicadas (\(f_i\)) debe estar equilibrada por la suma de todas las fuerzas internas ejercidas por los elementos conectados a ese nodo(\(K_{ij} u_j\)).
Implicación Matemática: Cada fila de la matriz de rigidez global \(\mathbf{K}\) corresponde a un grado de libertad (GDL) específico en un nodo particular. La ecuación representada por esa fila \((\sum K_{ij} u_j = f_i)\) es esencialmente la ecuación de equilibrio para ese GDL.
Visualización Matricial (Superposición de Rigideces):
Considere un sistema con solo dos elementos (1 y 2) conectados en un GDL global (ej., \(u_1\), el desplazamiento horizontal del Nodo 1). La rigidez global \(K_{11}\) en ese GDL es la suma de las contribuciones \(k_{11}^{(1)}\) y \(k_{11}^{(2)}\) de los elementos 1 y 2 a ese GDL:
Ejemplo: Superposición de Rigideces: Imagina tres resortes con rigideces \(k_1, k_2, k_3\) que se unen en un nodo central. Se aplica una fuerza externa \(F\) al nodo, causando un desplazamiento \(u\). El otro extremo de cada resorte está fijo. Las fuerzas internas en los resortes son \(f_1 = k_1 u\), \(f_2 = k_2 u\), \(f_3 = k_3 u\).Por equilibrio:
$$f_1 + f_2 + f_3 = F$$
Sustituyendo las relaciones fuerza-desplazamiento:
$$(k_1 + k_2 + k_3) u = F$$
La rigidez total \(K\) del nodo es simplemente la suma de las rigideces de los elementos conectados: \(K = k_1 + k_2 + k_3\).
Conclusión: La Matriz de Rigidez Global \(\mathbf{K}\) se construye sumando (superponiendo) las contribuciones de rigidez de los elementos individuales que comparten los mismos grados de libertad.
Dado que se asume que la estructura es continua, los elementos conectados en un nodo deben moverse juntos. No puede haber separaciones ni superposiciones en las uniones después de la deformación.
Implicación Matemática: Esta condición se impone asignando un grado de libertad único (por ejemplo, \(u_x\), \(u_y\), \(\theta_z\)) para cada dirección posible de movimiento en un nodo. Todos los elementos conectados a ese nodo comparten estas mismas variables de GDL.
Ejemplo: Si el resorte E1 y el resorte E2 se conectan en el nodo 1, el desplazamiento horizontal \(u_{1}\) debe ser el mismo para el extremo de E1 y el extremo de E2. Comparten la variable de desplazamiento común \(u_{1}\).
Mueve los controles para ver ambos principios en acción. En la parte A, observa cómo la rigidez total del nudo es la suma de las rigideces y cómo el resorte más rígido toma más carga. En la parte B, observa qué términos de la matriz \(\mathbf{K}\) aporta cada elemento y por qué solo se suman en el GDL que comparten.
3. Desarrollo Matemático: La Matriz de Rigidez del Elemento (\(\mathbf{k}'\)) 🧱#
La construcción de la matriz global \(\mathbf{K}\) comienza definiendo la matriz de rigidez para un elemento individual.
Notación: usamos minúscula para la matriz de un elemento y mayúscula solo para la estructura completa. El apóstrofo (\('\)) indica ejes locales, igual que en los desplazamientos \(u'\).
Símbolo
Significado
Tamaño (barra 2D)
En el código
\(\mathbf{k}'\)
Elemento, en ejes locales
\(2 \times 2\)
k_local
\(\bar{\mathbf{k}}\)
Elemento, en ejes globales
\(4 \times 4\)
k_elem
\(\mathbf{K}\)
Estructura completa (ensamblada)
\(n_{GDL} \times n_{GDL}\)
K
Cuando hay varios elementos, se indica cuál con un superíndice: \(\bar{\mathbf{k}}^{(1)}, \bar{\mathbf{k}}^{(2)}, \ldots\)
Consideremos el elemento estructural más simple: una barra uniforme de longitud \(L\), área de sección transversal \(A\), y Módulo de Young \(E\), orientada a lo largo de su propio eje local \(x'\). Solo resiste fuerzas axiales.
Ley Constitutiva (Esfuerzo-Deformación): \(\sigma = E \varepsilon\) (Ley de Hooke para materiales).
Fuerza Interna (Esfuerzo-Área): La fuerza axial interna \(P\) en la barra es \(P = \sigma A = E A \varepsilon\).
Deformación-Desplazamiento: La deformación unitaria \(\varepsilon\) es el cambio de longitud dividido por la longitud original. Si los desplazamientos en los extremos 1 y 2 son \(u'_1\) y \(u'_2\) a lo largo del eje local, el alargamiento total es \(\delta = u'_2 - u'_1\). Por lo tanto, \(\varepsilon = \delta / L = (u'_2 - u'_1) / L\).
Fuerza-Desplazamiento para el Elemento: Sustituyendo la deformación unitaria en la ecuación de la fuerza se obtiene la relación fundamental para la fuerza interna \(P\) (positiva en tracción):
$$P = E A \frac{(u'_2 - u'_1)}{L} = \frac{AE}{L} (u'_2 - u'_1)$$
Ahora, consideremos las fuerzas nodales \(f'_1\) y \(f'_2\) aplicadas a los nodos a lo largo del eje local, necesarias para mantener el equilibrio con la fuerza interna \(P\). Si \(P\) es de tracción, entonces \(f'_1 = -P\) (actúa hacia la izquierda) y \(f'_2 = P\) (actúa hacia la derecha).
\(AE/L\): Es la rigidez axial de la barra. Es la constante de proporcionalidad que relaciona la fuerza axial con el alargamiento axial. Un \(AE/L\) mayor significa que la barra es más rígida.
Elemento Diagonal \(k'_{11} = AE/L\): Representa la fuerza \(f'_1\) requerida en el nodo 1 para causar un desplazamiento unitario \(u'_1 = 1\) manteniendo el nodo 2 fijo (\(u'_2 = 0\)).
Elemento Fuera de Diagonal \(k'_{12} = -AE/L\): Representa la fuerza \(f'_1\) requerida en el nodo 1 para causar un desplazamiento unitario \(u'_2 = 1\) manteniendo el nodo 1 fijo (\(u'_1 = 0\)). El signo negativo indica que es una fuerza de reacción que se opone al alargamiento causado por mover el nodo 2.
4. Desarrollo Matemático: Transformación de Coordenadas 🔄#
Los elementos estructurales rara vez están perfectamente alineados con el sistema de coordenadas global (\(X, Y\)). Para ensamblar la matriz de rigidez global \(\mathbf{K}\), debemos transformar la matriz de rigidez de cada elemento desde su sistema de coordenadas local (\(x', y'\)) al sistema común global (\(X, Y\)).
La Matriz de Transformación \(\mathbf{T}\) utiliza cosenos directores para relacionar los desplazamientos locales del elemento (\(u'_1, u'_2\)) con sus desplazamientos nodales globales. Para un elemento de barra 2D orientado a un ángulo \(\theta\) con el eje X global, cada nodo tiene 2 GDL globales (\(U_x, U_y\)). El elemento tiene 4 GDL globales en total: \((U_{1x}, U_{1y}, U_{2x}, U_{2y})\).
Los desplazamientos axiales locales \(u'_1\) y \(u'_2\) son proyecciones de los desplazamientos globales sobre el eje local \(x'\) del elemento. Sean \(c = \cos\theta\) y \(s = \sin\theta\).
$$u'_1 = U_{1x} c + U_{1y} s$$$$u'_2 = U_{2x} c + U_{2y} s$$
En forma matricial, relacionando los 2 GDL locales con los 4 GDL globales:
$$\begin{pmatrix} u'_1 \\ u'_2 \end{pmatrix}_{\text{local}} = \begin{pmatrix} c & s & 0 & 0 \\ 0 & 0 & c & s \end{pmatrix} \begin{pmatrix} U_{1x} \\ U_{1y} \\ U_{2x} \\ U_{2y} \end{pmatrix}_{\text{global}}$$
Así, la matriz de transformación \(\mathbf{T}\) es:
$$\mathbf{T} = \begin{pmatrix} c & s & 0 & 0 \\ 0 & 0 & c & s \end{pmatrix}$$
B. La Matriz de Rigidez del Elemento en Ejes Globales (\(\bar{\mathbf{k}}\))#
Utilizando principios de equivalencia de trabajo o reglas de transformación de coordenadas, la matriz de rigidez del elemento en coordenadas globales (\(\bar{\mathbf{k}}\), ahora de \(4 \times 4\)) se obtiene mediante la transformación congruente:
Sustituyendo \(\mathbf{T}\) y \(\mathbf{k}'\) y realizando la multiplicación de matrices se obtiene la Matriz de Rigidez Global para un Elemento de Barra Axial 2D:
Esta matriz \(4 \times 4\) es el bloque de construcción fundamental que se “ensamblará” en la matriz de rigidez global \(\mathbf{K}\) de toda la estructura. Relaciona correctamente las fuerzas y desplazamientos aplicados en las direcciones globales \(X, Y\).
5. Ejercicio Demostrativo: \(\bar{\mathbf{k}}\) de una Barra Inclinada 💻#
Calculemos \(\bar{\mathbf{k}}\) para una barra con las siguientes propiedades:
\(A = 0.01 \, \text{m}^2\)
\(E = 200 \times 10^9 \, \text{Pa}\) (Acero)
\(L = 5 \, \text{m}\)
Ángulo \(\theta = 30^\circ\)
importnumpyasnp# Propiedades de la barraA=0.01# m^2E=200e9# PaL=5.0# mangulo_grados=30.0# Rigidez axialAE_L=(A*E)/L# Cosenos directoresangulo_rad=np.radians(angulo_grados)c=np.cos(angulo_rad)s=np.sin(angulo_rad)# Matriz de Rigidez Local (k_local)k_local=AE_L*np.array([[1,-1],[-1,1]])# Matriz de Transformación (T)T=np.array([[c,s,0,0],[0,0,c,s]])# Cálculo de k_elem = T_transpuesta * k_local * Tk_elem=T.T@k_local@Tprint(f"Rigidez Axial (AE/L): {AE_L:.2e} N/m")print("\nMatriz de Transformación T:")print(T)print("\nMatriz de Rigidez del Elemento en Ejes Globales (k_elem):")# Imprimir con formato para mejor lecturanp.set_printoptions(precision=3,suppress=True,linewidth=100)print(k_elem)# Verificación rápida: k_elem debe ser simétricaprint("\n¿Es k_elem simétrica?",np.allclose(k_elem,k_elem.T))
Este código calcula \(\mathbf{k}'\), \(\mathbf{T}\), y finalmente \(\bar{\mathbf{k}}\) usando operaciones matriciales de NumPy, demostrando el proceso matemático.
A lo largo de estas secciones construiremos un módulo de Python propio, rigidez.py, con las funciones del método. Cada sección le agrega funciones nuevas y, al final, tendrás una pequeña herramienta de análisis estructural que podrás reutilizar en tareas y proyectos. El módulo completo también se puede descargar.
Crea un archivo rigidez.py en la misma carpeta de tus scripts con el siguiente contenido. La función k_barra hace lo mismo que el ejercicio demostrativo, pero recibe las coordenadas de los nudos en lugar del largo y el ángulo: así calcula \(L\), \(c\) y \(s\) por sí misma y evita errores al calcular ángulos.
"""rigidez.py: funciones del Método de la Rigidez (MECCOM).
Este módulo se construye a lo largo de las Secciones I a V del curso.
Convenciones:
- Los nudos y los GDL se numeran desde 0 en el código.
- Armaduras 2D: 2 GDL por nudo (ux, uy); el nudo i usa los índices 2*i y 2*i + 1.
- Vigas: 2 GDL por nudo (v, θ). Marcos: 3 GDL por nudo (ux, uy, θ).
- Unidades del SI (m, N, Pa).
"""importnumpyasnpimportmatplotlib.pyplotasplt# ============================ Sección I ============================defgeometria_barra(xy1,xy2):"""Largo y cosenos directores de un elemento que va del nudo 1 al nudo 2.
Retorna (L, c, s), con c = cos(θ) y s = sin(θ).
"""dx,dy=np.subtract(xy2,xy1)L=np.hypot(dx,dy)returnL,dx/L,dy/Ldefk_barra(A,E,xy1,xy2):"""Matriz de rigidez 4x4 de una barra axial en ejes globales: Tᵀ k' T."""L,c,s=geometria_barra(xy1,xy2)k_local=A*E/L*np.array([[1,-1],[-1,1]])T=np.array([[c,s,0,0],[0,0,c,s]])returnT.T@k_local@T
Para usarlo desde otro script de la misma carpeta, se importa como cualquier biblioteca. Verifiquemos que entrega la misma matriz que el ejercicio demostrativo:
importnumpyasnpfromrigidezimportk_barra# Barra de 5 m a 30°: el nudo 2 está en (5 cos 30°, 5 sin 30°)theta=np.radians(30)k_elem=k_barra(0.01,200e9,(0,0),(5*np.cos(theta),5*np.sin(theta)))np.set_printoptions(precision=3,suppress=True,linewidth=100)print(k_elem)
Concepto: Explica con tus propias palabras la diferencia entre la matriz de rigidez de un elemento en ejes locales (\(\mathbf{k}'\)) y en ejes globales (\(\bar{\mathbf{k}}\)). ¿Por qué es necesaria la transformación?
Cálculo: Para una barra de aluminio (\(E = 70 \, \text{GPa}\)) de \(2 \, \text{m}\) de longitud y \(5 \, \text{cm}^2\) de área, calcula \(AE/L\) y escribe su matriz \(\mathbf{k}'\).
Cálculo: Si la barra del ejercicio 2 está a \(45^\circ\), calcula \(c\), \(s\), la matriz \(\mathbf{T}\) y, a mano, la matriz \(\bar{\mathbf{k}}\).
Cálculo: A mano, desarrolla el producto \(\mathbf{T}^T \mathbf{k}' \mathbf{T}\) en forma simbólica y verifica la primera fila de \(\bar{\mathbf{k}}\).
Cálculo: ¿Cómo se simplifica \(\bar{\mathbf{k}}\) para una barra horizontal (\(\theta = 0^\circ\)) y para una vertical (\(\theta = 90^\circ\))? Explica por qué aparecen filas y columnas de ceros.
Programación: Usa k_barra de tu módulo para verificar tu resultado del ejercicio 3 (pista: ubica el nudo 1 en el origen y el nudo 2 en \((L\cos\theta, L\sin\theta)\)). Calcula también \(\bar{\mathbf{k}}\) para una barra a \(120^\circ\) y comprueba que la matriz es simétrica con np.allclose.
Concepto: ¿Qué significan físicamente los términos \(\bar{k}_{11}\) y \(\bar{k}_{13}\) de \(\bar{\mathbf{k}}\)? ¿Por qué \(\bar{k}_{13}\) es negativo o cero?
Concepto: ¿Por qué \(\bar{\mathbf{k}}\) siempre es simétrica? (Pista: teorema de reciprocidad de Maxwell-Betti o energía de deformación).
Unidades y escala: ¿Cuáles son las unidades de los términos de \(\bar{\mathbf{k}}\) si \(A\) está en \(\text{m}^2\), \(E\) en \(\text{Pa}\) y \(L\) en \(\text{m}\)? ¿Cómo cambian sus valores si se duplica \(L\)? ¿Y si se duplica \(A\)?
Aplicación: Si una estructura combina barras normales con una barra muy corta y de área muy grande (un soporte casi rígido), ¿qué esperarías de los valores de su \(\bar{\mathbf{k}}\) comparados con los de las demás barras? ¿Qué problema numérico podría aparecer al resolver \(\mathbf{K}\mathbf{u} = \mathbf{f}\)? (Pista: número de condición).
Metodo de la rigidez -
Este artículo es parte de una serie.
