Capítulo 7

Ensamblaje de la matriz global

Tenemos la matriz de rigidez de cada barra y sabemos a qué grados de libertad de la armadura corresponde cada una de sus filas y columnas. Falta el paso que da nombre al método: ensamblar, es decir, sumar las rigideces de todas las barras en una sola matriz, la matriz de rigidez global (global stiffness matrix) de la armadura. Es una operación sencilla, que a mano resulta tediosa y propensa a errores, y que en Python se escribe en pocas líneas.

7.1 La idea del ensamblaje

Pensemos en el equilibrio de un nodo en la dirección de uno de sus grados de libertad, digamos el $i$ . La fuerza externa $F_i$ que actúa en esa dirección debe ser igual a la suma de las fuerzas que ejercen sobre el nodo todas las barras que llegan a él. Y la fuerza de cada barra depende, a través de su matriz $\mk$ , de los desplazamientos de sus extremos. Al escribir el equilibrio de todos los grados de libertad obtenemos el sistema

\vF = \mK\,\vU,

en el que cada término de $\mK$ es la suma de los aportes de las barras que conectan esos grados de libertad. Si $g_1, g_2, g_3, g_4$ es el vector de correspondencia de una barra, su término $k_{ab}$ se suma en la fila $g_a$ y la columna $g_b$ de la matriz global:

\begin{aligned} K_{g_a g_b} \leftarrow K_{g_a g_b} + k_{ab}, \qquad a, b = 1, \dots, 4. \end{aligned}

La figura 7.1 muestra la idea para la barra 2 del ejemplo A, cuyo vector de correspondencia es $[3, 4, 5, 6]$ : su matriz de $4\times4$ se suma en el bloque de filas y columnas 3 a 6 de la matriz de $6\times6$ . Las demás barras se suman de la misma forma, cada una en sus propias filas y columnas, y los términos que comparten se acumulan.

Figura 7.1:

Ensamblaje de la barra 2 del ejemplo A: cada término de su matriz se suma en la fila y la columna de la matriz global que indica su vector de correspondencia.

7.2 Matrices de ceros e índices

La matriz global empieza llena de ceros, y las barras van sumando sus aportes. np.zeros crea un arreglo de ceros con las dimensiones que se indiquen entre paréntesis:

n_gdl = 2 * len(nodos)
K = np.zeros((n_gdl, n_gdl))
K.shape

(6, 6)

Para leer o modificar un término de una matriz se indican su fila y su columna, separadas por una coma. Los obtenemos siempre con la función indices del capítulo 6, a partir de los grados de libertad. El operador += acumula: cada suma se agrega a lo que el término ya tenía.

# Término de la fila del gdl 5 y la columna del gdl 6
fila, columna = indices([5, 6])
K[fila, columna] += 10.0
K[fila, columna] += 2.5
print(fila, columna, K[fila, columna])

4 5 12.5

Para recorrer los cuatro grados de libertad de una barra usamos range(4), que produce los números 0, 1, 2 y 3. En este caso los necesitamos desde cero, porque son posiciones dentro de la matriz $\mk$ de NumPy:

for a in range(4):
    print(a, end=" ")

0 1 2 3

7.3 El ensamblaje, término a término

La ecuación clave se traduce de forma literal con dos bucles, uno dentro del otro: para cada barra, recorremos sus filas $a$ y, para cada fila, sus columnas $b$ . En total, 16 sumas por barra.

K = np.zeros((n_gdl, n_gdl))
for numero, barra in barras.items():
    g = gdl_de_barra(barra)
    k = rigideces[numero]
    for a in range(4):
        for b in range(4):
            fila, columna = indices([g[a], g[b]])
            K[fila, columna] += k[a, b]

7.4 Mostrar una matriz con sus grados de libertad

Al imprimir un arreglo de NumPy, sus filas y columnas no tienen rótulos, y no es fácil saber qué término corresponde a qué grado de libertad. Escribimos una función que imprime la matriz con los grados de libertad, numerados desde 1, como encabezados de filas y columnas:

def mostrar_matriz(matriz, filas, columnas=None,
                   decimales=1, ancho=7):
    """Imprime una matriz con los grados de libertad
    (desde 1) como encabezados de filas y columnas."""
    if columnas is None:
        columnas = filas
    encabezado = " gdl"
    for g in columnas:
        encabezado += f"{g:>{ancho}}"
    print(encabezado)
    for i, g in enumerate(filas):
        linea = f"{g:>4}"
        for valor in matriz[i]:
            # Sumar 0.0 evita mostrar -0.0
            valor = round(valor, decimales) + 0.0
            linea += f"{valor:>{ancho}.{decimales}f}"
        print(linea)

Si no se indican las columnas, se usan los mismos grados de libertad que para las filas, que es el caso de las matrices cuadradas como $\mK$ . En el capítulo 8 mostraremos también bloques no cuadrados.

Las rigideces de una armadura de acero son del orden de decenas de miles de kN/m. Para leerlas con comodidad las mostramos en kN/mm, dividiendo por 1000:

todos = list(range(1, n_gdl + 1))
print("K en kN/mm:")
mostrar_matriz(K / 1000, todos)

K en kN/mm:
 gdl      1      2      3      4      5      6
   1   95.5   38.4  -66.7    0.0  -28.8  -38.4
   2   38.4   51.2    0.0    0.0  -38.4  -51.2
   3  -66.7    0.0   95.5  -38.4  -28.8   38.4
   4    0.0    0.0  -38.4   51.2   38.4  -51.2
   5  -28.8  -38.4  -28.8   38.4   57.6    0.0
   6  -38.4  -51.2   38.4  -51.2    0.0  102.4

7.5 El ensamblaje con `np.ix\_`

NumPy permite seleccionar de una vez un bloque de filas y columnas de una matriz. La función np.ix_(filas, columnas) indica que queremos todas las combinaciones de esas filas y esas columnas: con los cuatro índices de una barra, el bloque de $4\times4$ que le corresponde en la matriz global.

g = indices(gdl_de_barra(barras[2]))
print(g)
K[np.ix_(g, g)].shape

[2, 3, 4, 5]
(4, 4)

Con esta selección, las 16 sumas del doble bucle se escriben en una sola línea, y el ensamblaje completo cabe en una función corta. Es la versión que usará nuestra calculadora:

def ensamblar_rigidez(nodos, barras, rigideces):
    """Matriz de rigidez global de la armadura."""
    n_gdl = 2 * len(nodos)
    K = np.zeros((n_gdl, n_gdl))
    for numero, barra in barras.items():
        g = indices(gdl_de_barra(barra))
        K[np.ix_(g, g)] += rigideces[numero]
    return K


K_global = ensamblar_rigidez(nodos, barras, rigideces)

Ambas versiones deben dar el mismo resultado:

np.allclose(K_global, K)

True

7.6 Verificación

La matriz global hereda las propiedades de las matrices de las barras. Debe ser simétrica y, como la armadura todavía no tiene apoyos, singular: puede trasladarse en $x$ , en $y$ y girar sin deformarse. Esos tres movimientos de cuerpo rígido hacen que su rango sea $2n - 3$ :

simetrica = np.allclose(K_global, K_global.T)
rango = np.linalg.matrix_rank(K_global)
print(f"Simétrica: {simetrica}")
print(f"Rango: {rango} de {n_gdl}")

Simétrica: True
Rango: 3 de 6

En consecuencia, una traslación de toda la armadura no produce fuerzas:

# Toda la armadura se traslada 1 mm en x
traslacion_x = np.array([1.0, 0.0, 1.0, 0.0, 1.0, 0.0])
K_global @ traslacion_x

array([0., 0., 0., 0., 0., 0.])

Comprobemos también dos términos a mano. Al grado de libertad 5 (nodo 3, dirección $x$ ) llegan las barras 2 y 3, cada una con $EA/L\,c^2 = 80 \cdot 0.36 = 28.8\,\mathrm{kN/mm}$ : su suma, 57.6 kN/mm, es el término de la fila y columna 5. Al grado de libertad 3 (nodo 2, dirección $x$ ) llegan las barras 1 y 2: $66.7 + 28.8 = 95.5\,\mathrm{kN/mm}$ .

Problema resuelto 7.1 · La matriz global del ejemplo B

Ensamblar la matriz de rigidez global del ejemplo B y verificar sus propiedades.

Solución

Encadenamos las funciones de los capítulos anteriores. Como la matriz es de $8\times8$ , la mostramos sin decimales y con columnas más angostas:

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)

print("K del ejemplo B en kN/mm:")
mostrar_matriz(K_b / 1000, list(range(1, 9)),
               decimales=0, ancho=6)
print(f"Simétrica: {np.allclose(K_b, K_b.T)}")
print(f"Rango: {np.linalg.matrix_rank(K_b)} de 8")

K del ejemplo B en kN/mm:
 gdl     1     2     3     4     5     6     7     8
   1   162    38  -133     0     0     0   -29   -38
   2    38   151     0     0     0  -100   -38   -51
   3  -133     0   162   -38   -29    38     0     0
   4     0     0   -38   151    38   -51     0  -100
   5     0     0   -29    38   162   -38  -133     0
   6     0  -100    38   -51   -38   151     0     0
   7   -29   -38     0     0  -133     0   162    38
   8   -38   -51     0  -100     0     0    38   151
Simétrica: True
Rango: 5 de 8

Verificación. La matriz es simétrica y su rango es $5 = 2 \cdot 4 - 3$ , como corresponde a una armadura libre en el plano. En el término de la fila y columna 1 (nodo 1, dirección $x$ ) se suman la barra 1, horizontal, con $EA/L = 133.3\,\mathrm{kN/mm}$ , y la diagonal 5, con $80 \cdot 0.36 = 28.8\,\mathrm{kN/mm}$ : 162 kN/mm. La barra vertical 4 no aporta nada en $x$ .

Repaso y resumen del capítulo

La matriz de rigidez global es la suma de las matrices de las barras: el término $k_{ab}$ de una barra se suma en la fila $g_a$ y la columna $g_b$ de $\mK$ .
np.zeros crea una matriz de ceros; sus términos se leen y se modifican con matriz[fila, columna], y += acumula.
range(4) produce 0, 1, 2 y 3; los bucles pueden anidarse.
np.ix_(filas, columnas) selecciona un bloque de una matriz y permite ensamblar cada barra en una línea.
mostrar_matriz imprime una matriz con sus grados de libertad numerados desde 1.
La matriz global es simétrica y, sin apoyos, singular: su rango es $2n - 3$ .

Ejercicios

Ejercicio 7.1

Calcular a mano el término de la fila y columna 6 de la matriz global del ejemplo A, en kN/mm, sumando los aportes de las barras que llegan al nodo 3. Comprobarlo con el programa.

Ver la solución

Al nodo 3 llegan las barras 2 y 3. En ambas, el grado de libertad 6 es el cuarto ( $u_4$ ), y su término es $EA/L\,s^2 = 80 \cdot 0.64 = 51.2\,\mathrm{kN/mm}$ :

# Al gdl 6 (nodo 3, dirección y) llegan las barras 2 y 3
# Su gdl local es u4 (posición 3 de la matriz de 4 x 4)
aporte_2 = rigideces[2][3, 3] / 1000
aporte_3 = rigideces[3][3, 3] / 1000
fila, columna = indices([6, 6])
termino = K_global[fila, columna] / 1000
print(f"Barra 2: {aporte_2:.1f} kN/mm")
print(f"Barra 3: {aporte_3:.1f} kN/mm")
print(f"Suma:    {aporte_2 + aporte_3:.1f} kN/mm")
print(f"K:       {termino:.1f} kN/mm")

Barra 2: 51.2 kN/mm
Barra 3: 51.2 kN/mm
Suma:    102.4 kN/mm
K:       102.4 kN/mm

Ejercicio 7.2

Ensamblar una matriz global del ejemplo A que contenga solo la barra 2 y mostrarla. Comparar el resultado con la figura 7.1.

Ver la solución

solo_barra_2 = {2: barras[2]}
K_2 = ensamblar_rigidez(nodos, solo_barra_2, rigideces)
mostrar_matriz(K_2 / 1000, list(range(1, 7)))

 gdl      1      2      3      4      5      6
   1    0.0    0.0    0.0    0.0    0.0    0.0
   2    0.0    0.0    0.0    0.0    0.0    0.0
   3    0.0    0.0   28.8  -38.4  -28.8   38.4
   4    0.0    0.0  -38.4   51.2   38.4  -51.2
   5    0.0    0.0  -28.8   38.4   28.8  -38.4
   6    0.0    0.0   38.4  -51.2  -38.4   51.2

Solo el bloque de filas y columnas 3 a 6 es distinto de cero, como en la figura 7.1.

Ejercicio 7.3

Verificar que una traslación vertical de toda la armadura del ejemplo A tampoco produce fuerzas.

Ver la solución

traslacion_y = np.array([0.0, 1.0, 0.0, 1.0, 0.0, 1.0])
K_global @ traslacion_y

array([0., 0., 0., 0., 0., 0.])

Ejercicio 7.4

Ensamblar la matriz global del ejemplo A recorriendo las barras en orden inverso (3, 2, 1) y comprobar que el resultado no cambia. ¿Por qué?

Ver la solución

al_reves = {3: barras[3], 2: barras[2], 1: barras[1]}
K_al_reves = ensamblar_rigidez(nodos, al_reves, rigideces)
np.allclose(K_al_reves, K_global)

True

La suma es conmutativa: el orden en que las barras aportan sus términos no cambia el resultado.

Ensamblaje de la matriz global

7.1 La idea del ensamblaje

7.2 Matrices de ceros e índices

7.3 El ensamblaje, término a término

7.4 Mostrar una matriz con sus grados de libertad

7.5 El ensamblaje con np.ix\_

7.6 Verificación

Repaso y resumen del capítulo

Ejercicios

7.5 El ensamblaje con `np.ix\_`