¿Por qué se comba así el puente? El método de rigidez directa para cerchas#

En 1956, el ingeniero de Boeing M. J. Turner se atascó calculando a mano las tensiones de un ala en flecha. Con cientos de barras no se persigue el equilibrio ecuación por ecuación. La salida que publicó con sus colegas fue el método de rigidez directa. Se forma la rigidez de un elemento, se suma donde los grados de libertad coinciden y se resuelve $Ku = f$ una vez.

El ángulo va primero. Una barra tumbada apaga las casillas verticales. Al bajar el control del ángulo a 0, las filas y columnas verticales de la rigidez 4×4 pasan a cero juntas.

Lo que conviene mirar es qué casillas se vacían con $\theta = 0$. La barra solo es rígida a lo largo de su eje. Al subir el ángulo, esos ceros se reparten entre $x$ y $y$.

La rigidez empieza con un solo resorte#

Un elemento barra es un resorte que solo lleva fuerza en su eje. Longitud $L$, área $A$, módulo $E$, y la rigidez axial es $k = EA/L$. Fuerza por unidad de alargamiento. Con desplazamientos de extremo $u_1, u_2$, las fuerzas nodales quedan así:

$f_1, f_2$ son las fuerzas de extremo, y el corchete es la rigidez local $k_e$. Cada fila suma cero. Un desplazamiento rígido de los dos nudos no cuesta fuerza. Sin apoyos, la matriz es singular por eso.

De lo local a lo global: la rotación#

Las barras apuntan a todos lados. Para ensamblar, el desplazamiento axial local tiene que escribirse en $x, y$ globales. Sea $\theta$ el ángulo con el eje $x$, $c = \cos\theta$, $s = \sin\theta$. El desplazamiento axial es una proyección: $u_{\text{axial}} = c\,u_x + s\,u_y$. Esa transformación, aplicada por ambos lados, da la rigidez de elemento global 4×4.

Cada entrada une el $x, y$ del nudo 1 con el $x, y$ del nudo 2. Con $\theta = 0$, $s = 0$ y $c = 1$, la matriz por $EA/L$ tiene una primera fila 1, 0, −1, 0 y una segunda fila de ceros. La tercera fila cambia el signo de la primera. La cuarta también es cero. Esos ceros explican que cada bloque 2×2 tenga rango uno. Si las coordenadas se inclinan, el tensor de rigidez gira con ellas. Ese giro está en rotar el tensor constitutivo con la métrica.

Sumar los GDL solapados: ensamblaje global#

El ensamblaje es sumar rigidez en los grados de libertad compartidos. El nudo $i$ vive en los índices globales $2i, 2i+1$, primero $x$ y luego $y$. Se dispersan los cuatro GDL locales del elemento sobre esos índices y se acumula.

$\mathbf{L}_e$ es la matriz de selección que lleva los GDL del elemento a los globales, y $K$ es la rigidez global. El código no forma $\mathbf{L}_e$. Suma directo por la lista de índices. Varias barras en un nudo se apilan en ese bloque diagonal.

Condiciones de contorno: borrar los apoyos#

Una $K$ recién ensamblada es singular. La estructura todavía flota y el movimiento rígido queda libre. Los desplazamientos en los apoyos tienen que quedar en cero. Los GDL se parten en un conjunto libre $f$ y uno restringido $c$.

Con $u_c = 0$ solo queda el bloque superior: $K_{ff}\,u_f = F_f$. Ese sistema reducido no es singular. Las reacciones vuelven después como $R_c = K_{cf}\,u_f$. $R_c$ es la fuerza que el apoyo devuelve a la estructura.

Al subir a la vez el control de carga y EA, se ve cómo baja la flecha de la punta.

Un EA mayor recorta la flecha bajo la misma carga. La matriz de rigidez es el mapa de carga a desplazamiento. Con la carga fija y el EA al doble, la flecha de la punta debería ir hacia la mitad.

Python: resolver una cercha de 12 barras#

Una ménsula de tres paneles empotrada en el muro. Ocho nudos, doce barras. Entran coordenadas, conectividad y cargas. Salen desplazamientos y fuerzas axiales.

import numpy as np
 
nodes = np.array([[0,0],[0,1],[1,0],[1,1],[2,0],[2,1],[3,0],[3,1]], float)
members = [(0,2),(2,4),(4,6),(1,3),(3,5),(5,7),
           (2,3),(4,5),(6,7),(1,2),(3,4),(5,6)]
EA = 8.0e6           # rigidez axial EA [N]
fixed = [0, 1]       # nodos anclados al muro
 
def bar_stiffness(p1, p2, EA):
    d = p2 - p1
    L = np.hypot(*d)
    c, s = d / L
    k = EA / L * np.array([[ c*c,  c*s, -c*c, -c*s],
                           [ c*s,  s*s, -c*s, -s*s],
                           [-c*c, -c*s,  c*c,  c*s],
                           [-c*s, -s*s,  c*s,  s*s]])
    return k, L, (c, s)
 
ndof = nodes.shape[0] * 2
K = np.zeros((ndof, ndof))
geom = []
for a, b in members:                       # ensamblar la rigidez global
    k, L, cs = bar_stiffness(nodes[a], nodes[b], EA)
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    K[np.ix_(dof, dof)] += k
    geom.append((L, cs))
 
F = np.zeros(ndof)                          # carga de 24 kN en el extremo libre (nodos 6,7)
for n in (6, 7):
    F[2*n+1] -= 12.0e3
 
fixed_dof = [d for n in fixed for d in (2*n, 2*n+1)]
free_dof = [d for d in range(ndof) if d not in fixed_dof]
 
u = np.zeros(ndof)                          # sistema reducido K_ff u_f = F_f
u[free_dof] = np.linalg.solve(K[np.ix_(free_dof, free_dof)],
                              F[free_dof])
 
for m, (a, b) in enumerate(members):        # fuerza axial N = (EA/L)[-c,-s,c,s]·u
    L, (c, s) = geom[m]
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    N = EA / L * np.array([-c, -s, c, s]) @ u[dof]
    print(f"member {a}-{b}: N = {N/1e3:+7.2f} kN "
          f"({'tension' if N > 0 else 'compression'})")
 
print(f"free-end drop = {u[2*6+1]*1e3:.3f} mm")

np.ix_ coloca la rigidez de cada elemento en las casillas globales correctas. Ese bucle tiene la misma forma que el ensamblaje de un FEA comercial.

La separación entre nudos es 1 m. Nueve barras horizontales y verticales tienen $L = 1$. Con el EA del código, $8.0\times 10^6$ N, queda $k = EA/L = 8.0\times 10^6$ N/m. La carga vertical de un nudo de punta es $12\times 10^3$ N, así que esa barra sola se mueve $1.5\times 10^{-3}$ m, o sea 1.5 mm. La punta de la cercha no es esa división. Las barras se reparten la carga, y la matriz que de verdad se resuelve es la $K_{ff}$ reducida.

Hay tres diagonales. La que une (0, 1) con (1, 0) mide $\sqrt{2}$, así que su rigidez es $1/\sqrt{2}$ veces la de una barra horizontal. Mismo EA, diagonal más blanda. Esa barra está a 45°. Entonces $c^2$ y $s^2$ valen $1/2$, y el valor absoluto de $cs$ también es $1/2$. Toda casilla distinta de cero en la 4×4 es la mitad de $EA/L$. Con longitud $\sqrt{2}$, el tamaño de la casilla es $EA/(2\sqrt{2})$. Una barra horizontal deja vacías las casillas verticales. Una diagonal llena las cuatro direcciones con ese tamaño. El mismo EA se parte así solo por el ángulo. Ocho nudos dan 16 GDL. Los nudos de muro 0 y 1 sujetan 4, de modo que el sistema es 12×12. En 2D hacen falta al menos 3 restricciones. Este ejemplo usa 4. Si se suelta la $y$ de un nudo de muro, la cuenta queda en ese mínimo. Si se sueltan las dos, $K_{ff}$ vuelve a ser singular.

Cuando la matriz de rigidez se vuelve singular#

El fallo que aparece en la práctica es una matriz singular. Sale por tres lados.

Si falta un apoyo, queda un modo rígido y $K_{ff}$ sigue singular. En dos dimensiones hay que sujetar al menos 3 GDL. En tres, 6.

Un panel cuadrado que nunca se partió en triángulos es un mecanismo. Una cercha tiene que llenarse de triángulos o se afloja.

Longitud nula o nudos duplicados hacen $L = 0$ y la división estalla. Las coordenadas coincidentes se filtran antes de unir la malla.

Si falta un nudo del muro, el montaje se detiene ahí#

El 1.5 mm pertenece a una sola barra horizontal. La flecha de la ménsula de 12 barras sale de resolver $K_{ff}\,u_f = F_f$, y esa ecuación no está definida hasta quitar los cuatro GDL del muro. Un apoyo de menos deja un modo rígido, y el desplazamiento es un espacio nulo, no un número.

Al cambiar las barras por elementos de superficie, el mismo montaje puede agarrotarse por otra razón. Cómo aparece ese agarrotamiento en una cáscara delgada está en el bloqueo a cortante y cómo MITC devuelve la rigidez. En la cercha se comprueban primero los cuatro apoyos. Esa página solo hace falta cuando el modelo se adelgaza.