Capítulo 9

Desplazamientos y reacciones

Llegamos al momento que justifica todo el trabajo anterior: resolver el sistema. Con los bloques del capítulo 8, la primera ecuación nos da los desplazamientos de los nodos libres y la segunda, las reacciones en los apoyos. NumPy resuelve el sistema en una línea; nuestro trabajo está en hacerlo bien, en detectar a tiempo una armadura inestable y en comprobar que el resultado respeta el equilibrio.

9.1 La solución del sistema

La ecuación (8.1) tiene una sola incógnita, $\vU_l$ . Pasando al lado izquierdo el término conocido obtenemos un sistema de ecuaciones lineales:

\begin{aligned} \mK_{ll}\,\vU_l = \vF_l - \mK_{lr}\,\vU_r \end{aligned}

Con los desplazamientos de los nodos libres conocidos, la ecuación (8.2) da las fuerzas en los grados de libertad restringidos. Como esas fuerzas incluyen las cargas aplicadas directamente sobre los apoyos, las reacciones son

\vect{R} = \mK_{rl}\,\vU_l + \mK_{rr}\,\vU_r - \vF_r, \tag{9.1}

donde $\vF_r$ contiene solo esas cargas aplicadas, casi siempre nulas. Mientras los apoyos no se desplacen, $\vU_r = \vect{0}$ ; en el capítulo 8 dejamos el término $\mK_{lr}\vU_r$ para poder incorporar más adelante los asentamientos de apoyo.

9.2 Resolver un sistema de ecuaciones con NumPy

La función np.linalg.solve(A, b) resuelve el sistema $\vect{A}\vect{x} = \vect{b}$ . Los desplazamientos resultan en metros, porque usamos kN y m; los mostramos en milímetros:

# Apoyos sin asentamientos: U_r es nulo
U_r = np.zeros(len(restringidos))
lado_derecho = F_l - K_lr @ U_r
U_l = np.linalg.solve(K_ll, lado_derecho)
U_l * 1000  # en mm

array([ 0.4,  0.6, -0.7])

Siempre conviene comprobar que la solución satisface el sistema. El residuo (residual) $\mK_{ll}\vU_l - (\vF_l - \mK_{lr}\vU_r)$ debe ser nulo:

residuo = K_ll @ U_l - lado_derecho
print(f"Mayor residuo: {np.abs(residuo).max():.1e} kN")
print(np.allclose(K_ll @ U_l, lado_derecho))

Mayor residuo: 7.1e-15 kN
True

El residuo no es exactamente cero sino del orden de $10^{-15}$ . Las computadoras guardan los números reales con unas 16 cifras significativas, y cada operación redondea la última. Por eso los resultados numéricos se comparan con una tolerancia (tolerance), como hace np.allclose, y nunca exigiendo igualdad exacta.

Cuidado

En los libros de texto la solución se escribe $\vU_l = \mK_{ll}^{-1}(\vF_l - \mK_{lr}\vU_r)$ , y NumPy tiene la función np.linalg.inv para invertir una matriz. El resultado es el mismo:

U_l_con_inv = np.linalg.inv(K_ll) @ lado_derecho
np.allclose(U_l_con_inv, U_l)

True

Pero calcular la inversa para después multiplicarla es más lento y acumula más errores de redondeo que resolver el sistema directamente, porque solve hace solo las operaciones necesarias para esa $\vect{b}$ . En una armadura de pocas barras la diferencia no se nota; en una de miles de grados de libertad, sí. La regla es simple: para resolver un sistema, solve; la inversa, solo si de verdad se necesita la matriz inversa.

9.3 Armaduras inestables: `raise`, `try` y `except`

En el capítulo 8 vimos que si la armadura es inestable, $\mK_{ll}$ es singular. Resolver un sistema singular no tiene sentido, y no podemos confiar en que NumPy siempre lo detecte: por el redondeo, una matriz singular puede parecer «casi» invertible y producir desplazamientos enormes y sin significado. Por eso nuestra función comprueba el rango antes de resolver y, si es insuficiente, se detiene con un mensaje claro:

def resolver_desplazamientos(K_ll, K_lr, F_l, U_r):
    """Desplazamientos de los gdl libres, U_l."""
    if np.linalg.matrix_rank(K_ll) < len(K_ll):
        raise ValueError("La armadura es inestable: "
                         "K_ll es singular.")
    return np.linalg.solve(K_ll, F_l - K_lr @ U_r)

La instrucción raise produce una excepción (exception): detiene la función y comunica qué salió mal. ValueError es el tipo de excepción que Python usa para «un dato con un valor inadecuado». Quien llama a la función puede atrapar la excepción con try y except y decidir qué hacer, en lugar de que el programa se detenga. Probemos con el ejemplo A sin su barra 1, el mecanismo de los ejercicios 6.3 y 8.3:

barras_sin_1 = {2: barras[2], 3: barras[3]}
K_sin_1 = ensamblar_rigidez(nodos, barras_sin_1,
                            rigideces)
K_ll_sin_1, K_lr_sin_1, _, _ = particionar(
    K_sin_1, libres, restringidos)
try:
    resolver_desplazamientos(K_ll_sin_1, K_lr_sin_1,
                             F_l, U_r)
except ValueError as error:
    print(f"No se puede resolver. {error}")

No se puede resolver. La armadura es inestable: K_ll es singular.

9.4 Desplazamientos y reacciones del ejemplo A

Resolvemos ahora el ejemplo A con la función, y armamos el vector completo de desplazamientos, con los valores calculados en los grados de libertad libres y los impuestos en los restringidos. La asignación con una lista de índices coloca cada valor en su lugar:

U_r = np.zeros(len(restringidos))
U_l = resolver_desplazamientos(K_ll, K_lr, F_l, U_r)

# Vector completo de desplazamientos
U = np.zeros(2 * len(nodos))
U[indices(libres)] = U_l
U[indices(restringidos)] = U_r

print(" gdl    U (mm)")
mostrar_vector(U * 1000, list(range(1, len(U) + 1)),
               decimales=3)

 gdl    U (mm)
   1     0.000
   2     0.000
   3     0.431
   4     0.000
   5     0.563
   6    -0.650

Los resultados tienen sentido físico: el nodo 3 se desplaza 0.563 mm hacia la derecha y 0.650 mm hacia abajo, en la dirección de la carga; el rodillo del nodo 2 se desliza 0.431 mm hacia la derecha, porque la barra 1 se alarga; y el nodo 1 no se mueve. Son desplazamientos pequeños, del orden de los que esperamos en una armadura de acero de estas dimensiones.

Las reacciones se calculan con la ecuación (9.1):

def calcular_reacciones(K_rl, K_rr, U_l, U_r, F_r):
    """Reacciones en los grados de libertad restringidos.

    F_r son las cargas aplicadas directamente sobre los
    apoyos (normalmente ceros)."""
    return K_rl @ U_l + K_rr @ U_r - F_r


F_r = F[indices(restringidos)]
R = calcular_reacciones(K_rl, K_rr, U_l, U_r, F_r)

print(" gdl    R (kN)")
mostrar_vector(R, restringidos)

 gdl    R (kN)
   1    -20.00
   2     11.67
   4     38.33

Como el ejemplo A es isostático, podemos comprobarlas con la estática: la reacción horizontal del nodo 1 equilibra la carga de 20 kN, y tomando momentos respecto del nodo 1, $6 R_{2y} = 20 \cdot 4 + 50 \cdot 3$ , de donde $R_{2y} = 38.33\,\mathrm{kN}$ (ejercicio 9.1).

9.5 Verificación del equilibrio

Una armadura en equilibrio cumple, en conjunto, las tres ecuaciones de la estática: la suma de las fuerzas en $x$ , la suma en $y$ y la suma de momentos de las cargas y las reacciones deben ser nulas. Es la verificación más importante del método de rigidez, porque no depende de cómo calculamos los resultados. El método copy crea una copia del vector de cargas, para sumarle las reacciones sin modificar el original:

def comprobar_equilibrio(nodos, F, R, restringidos):
    """Resultantes de cargas más reacciones: Fx, Fy y el
    momento respecto del origen. Deben ser nulas."""
    total = F.copy()
    total[indices(restringidos)] += R
    momento = 0.0
    for nodo in nodos:
        ix, iy = indices(gdl_del_nodo(nodo))
        x, y = nodos[nodo]["x"], nodos[nodo]["y"]
        momento += x * total[iy] - y * total[ix]
    return total[0::2].sum(), total[1::2].sum(), momento

suma_x, suma_y, momento = comprobar_equilibrio(
    nodos, F, R, restringidos)
print(f"Suma de fuerzas en x: {suma_x:.2e} kN")
print(f"Suma de fuerzas en y: {suma_y:.2e} kN")
print(f"Suma de momentos:     {momento:.2e} kN·m")

Suma de fuerzas en x: 3.55e-15 kN
Suma de fuerzas en y: 1.42e-14 kN
Suma de momentos:     5.68e-14 kN·m

Las tres sumas son del orden de $10^{-14}$ , es decir, nulas dentro del redondeo: la armadura está en equilibrio.

Problema resuelto 9.1 · El ejemplo B

Calcular los desplazamientos y las reacciones del ejemplo B y verificar el equilibrio global.

Solución

Encadenamos todas las funciones de la calculadora:

geometria_b = calcular_geometria(nodos_b, barras_b)
rigideces_b = calcular_rigideces(barras_b, geometria_b)
K_b = ensamblar_rigidez(nodos_b, barras_b, rigideces_b)
F_b = ensamblar_cargas(nodos_b, cargas_b)
libres_b, restringidos_b = clasificar_gdl(nodos_b,
                                          apoyos_b)
K_ll_b, K_lr_b, K_rl_b, K_rr_b = particionar(
    K_b, libres_b, restringidos_b)

U_r_b = np.zeros(len(restringidos_b))
U_l_b = resolver_desplazamientos(
    K_ll_b, K_lr_b, F_b[indices(libres_b)], U_r_b)
R_b = calcular_reacciones(
    K_rl_b, K_rr_b, U_l_b, U_r_b,
    F_b[indices(restringidos_b)])

print(" gdl    U (mm)")
mostrar_vector(U_l_b * 1000, libres_b, decimales=3)
print(" gdl    R (kN)")
mostrar_vector(R_b, restringidos_b)

 gdl    U (mm)
   3     0.179
   5     0.900
   6    -0.081
   7     0.854
   8    -0.481
 gdl    R (kN)
   1    -30.00
   2      0.00
   4     80.00

suma_x, suma_y, momento = comprobar_equilibrio(
    nodos_b, F_b, R_b, restringidos_b)
print(f"ΣFx = {suma_x:.1e}, ΣFy = {suma_y:.1e}, "
      f"ΣM = {momento:.1e}")

ΣFx = -1.8e-14, ΣFy = 0.0e+00, ΣM = 5.7e-14

Verificación. Las sumas son nulas dentro del redondeo. Llama la atención que las reacciones puedan comprobarse con la estática, aunque la armadura sea hiperestática: con un apoyo articulado y un rodillo hay solo tres reacciones, de modo que la armadura es isostática externamente. La hiperestaticidad es interna: sobra una barra, y la estática no basta para repartir las fuerzas entre las barras. Eso lo resolverá el capítulo 10.

Repaso y resumen del capítulo

Los desplazamientos libres se obtienen de $\mK_{ll}\vU_l = \vF_l - \mK_{lr}\vU_r$ , y las reacciones, de $\vect{R} = \mK_{rl}\vU_l + \mK_{rr}\vU_r - \vF_r$ .
np.linalg.solve resuelve un sistema lineal; invertir la matriz es más lento y menos preciso.
Los resultados numéricos se comparan con tolerancias: el redondeo deja residuos del orden de $10^{-15}$ .
raise produce una excepción; try y except la atrapan. La calculadora rechaza las armaduras inestables antes de resolver.
El equilibrio global (fuerzas en $x$ e $y$ y momentos) es la verificación más importante de los resultados.

Ejercicios

Ejercicio 9.1

Calcular a mano, con las ecuaciones de la estática, las reacciones del ejemplo A, y compararlas con las del programa.

Ver la solución

La suma de fuerzas horizontales da $R_{1x} = -20\,\mathrm{kN}$ . Con los momentos respecto del nodo 1 obtenemos $R_{2y}$ , y con la suma de fuerzas verticales, $R_{1y}$ :

# Momentos respecto del nodo 1, positivos antihorarios:
# R2y · 6 - 20 · 4 - 50 · 3 = 0
R2y = (20 * 4 + 50 * 3) / 6
R1y = 50 - R2y
R1x = -20
print(f"R1x = {R1x:.2f} kN, R1y = {R1y:.2f} kN, "
      f"R2y = {R2y:.2f} kN")

R1x = -20.00 kN, R1y = 11.67 kN, R2y = 38.33 kN

Los valores coinciden con los del programa.

Ejercicio 9.2

Resolver el ejemplo A con dos apoyos articulados (como en el ejercicio 6.2). ¿Pueden comprobarse ahora las reacciones con la estática? ¿Por qué?

Ver la solución

apoyos_2 = {
    1: {"x": True, "y": True},
    2: {"x": True, "y": True},
}
libres_2, restringidos_2 = clasificar_gdl(nodos, apoyos_2)
K_ll_2, K_lr_2, K_rl_2, K_rr_2 = particionar(
    K_global, libres_2, restringidos_2)
U_r_2 = np.zeros(len(restringidos_2))
U_l_2 = resolver_desplazamientos(
    K_ll_2, K_lr_2, F[indices(libres_2)], U_r_2)
R_2 = calcular_reacciones(K_rl_2, K_rr_2, U_l_2, U_r_2,
                          F[indices(restringidos_2)])
print(" gdl    U (mm)")
mostrar_vector(U_l_2 * 1000, libres_2, decimales=3)
print(" gdl    R (kN)")
mostrar_vector(R_2, restringidos_2)

 gdl    U (mm)
   5     0.347
   6    -0.488
 gdl    R (kN)
   1      8.75
   2     11.67
   3    -28.75
   4     38.33

Con dos apoyos articulados hay cuatro reacciones y solo tres ecuaciones de equilibrio: la armadura es hiperestática externamente, y las reacciones horizontales no pueden obtenerse con la estática. Su suma sí: $8.75 - 28.75 = -20\,\mathrm{kN}$ , que equilibra la carga horizontal. Las reacciones verticales no cambian, porque los momentos respecto del nodo 1 no dependen de las reacciones horizontales.

Ejercicio 9.3

Duplicar las cargas del ejemplo A y comprobar que los desplazamientos se duplican. ¿Qué hipótesis del método garantiza este resultado?

Ver la solución

U_l_doble = resolver_desplazamientos(K_ll, K_lr,
                                     2 * F_l, U_r)
np.allclose(U_l_doble, 2 * U_l)

True

La hipótesis de comportamiento elástico lineal con desplazamientos pequeños: $\mK$ no depende de las cargas, de modo que los desplazamientos son proporcionales a ellas.

Ejercicio 9.4

Calcular el mayor desplazamiento, en valor absoluto, de los nodos del ejemplo B. La función np.abs da el valor absoluto de cada elemento de un arreglo, y el método max, el mayor de ellos.

Ver la solución

geometria_b = calcular_geometria(nodos_b, barras_b)
rigideces_b = calcular_rigideces(barras_b, geometria_b)
K_b = ensamblar_rigidez(nodos_b, barras_b, rigideces_b)
F_b = ensamblar_cargas(nodos_b, cargas_b)
libres_b, restringidos_b = clasificar_gdl(nodos_b,
                                          apoyos_b)
K_ll_b, K_lr_b, _, _ = particionar(K_b, libres_b,
                                   restringidos_b)
U_l_b = resolver_desplazamientos(
    K_ll_b, K_lr_b, F_b[indices(libres_b)],
    np.zeros(len(restringidos_b)))
maximo = np.abs(U_l_b * 1000).max()
print(f"Mayor desplazamiento: {maximo:.3f} mm")

Mayor desplazamiento: 0.900 mm

Desplazamientos y reacciones

9.1 La solución del sistema

9.2 Resolver un sistema de ecuaciones con NumPy

9.3 Armaduras inestables: raise, try y except

9.4 Desplazamientos y reacciones del ejemplo A

9.5 Verificación del equilibrio

Repaso y resumen del capítulo

Ejercicios

9.3 Armaduras inestables: `raise`, `try` y `except`