Capítulo 5

La rigidez de una barra

Con la geometría de todas las barras resuelta, entramos al corazón del método: la matriz de rigidez de cada barra, que relaciona las fuerzas en sus extremos con los desplazamientos de sus nodos. Para escribirla necesitamos dos herramientas nuevas. Las funciones, que nos permiten dar nombre a un cálculo y reutilizarlo sin copiarlo, y NumPy, la biblioteca que convierte a Python en una herramienta de cálculo matricial. Al final del capítulo, nuestra calculadora tendrá la matriz de rigidez de cada barra de cualquier armadura.

5.1 La matriz de rigidez de una barra

Una barra de armadura solo trabaja a fuerza axial. En su propio eje $x'$ , los desplazamientos de sus extremos $\vul = [u'_1,\ u'_2]\transp$ y las fuerzas axiales en ellos $\vf' = [f'_1,\ f'_2]\transp$ se relacionan por la ley de Hooke, con la matriz de rigidez local (local stiffness matrix):

\vf' = \mkl\,\vul, \qquad \mkl = \EAL \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}. \tag{5.1}

Los nodos, en cambio, se desplazan en las direcciones $x$ e $y$ globales: cada barra tiene cuatro grados de libertad (degrees of freedom), $\vu = [u_1,\ u_2,\ u_3,\ u_4]\transp$ (figura 5.1). El desplazamiento de cada extremo en la dirección de la barra es la proyección de su desplazamiento global sobre el eje $x'$ , lo que se escribe con la matriz de transformación (transformation matrix):

\vul = \mT\,\vu, \qquad \mT = \begin{bmatrix} c & s & 0 & 0 \\ 0 & 0 & c & s \end{bmatrix}. \tag{5.2}

Figura 5.1:

Grados de libertad de una barra en ejes globales y su eje local $x'$ .

Como las fuerzas se transforman con la transpuesta, $\vf = \mT\transp\vf'$ , llegamos a la matriz de rigidez de la barra en ejes globales:

\begin{aligned} \mk = \mT\transp\,\mkl\,\mT = \EAL \begin{bmatrix} c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2 \end{bmatrix} \end{aligned}

Podríamos programar directamente la forma desarrollada, término a término. Preferimos, sin embargo, programar el producto $\mT\transp\mkl\mT$ : es más corto, se parece más a la teoría y, como veremos en los próximos capítulos, las mismas matrices $\mT$ y $\mkl$ nos servirán para calcular las fuerzas internas.

5.2 Funciones

Una función (function) es un cálculo con nombre. Se define con def, seguido del nombre, de los parámetros (parameters) entre paréntesis y de dos puntos. El cuerpo va sangrado, como en un bucle, y la instrucción return indica el resultado:

def rigidez_axial(E, A, L):
    """Rigidez axial EA/L de una barra, en kN/m."""
    return E * A / L


rigidez_axial(2.0e8, 0.0020, 5.0)

80000.0

Definir la función no calcula nada: solo le enseña a Python cómo hacerlo. El cálculo ocurre cuando la llamamos, escribiendo su nombre con los valores de los parámetros, los argumentos (arguments), entre paréntesis. El texto entre triples comillas, debajo de la primera línea, es la documentación (docstring) de la función: explica qué hace y en qué unidades.

Con una función resolvemos el problema que dejamos pendiente en el capítulo 4: el bucle de la geometría estaba escrito para los diccionarios nodos y barras, y para el ejemplo B tuvimos que copiarlo. Convertido en función, recibe los diccionarios como parámetros:

def calcular_geometria(nodos, barras):
    """Longitud y cosenos directores de cada barra.

    Devuelve un diccionario {número de barra:
    {"longitud": L, "c": c, "s": s}}.
    """
    geometria = {}
    for numero, barra in barras.items():
        inicial = nodos[barra["nodo_inicial"]]
        final = nodos[barra["nodo_final"]]
        dx = final["x"] - inicial["x"]
        dy = final["y"] - inicial["y"]
        L = math.sqrt(dx**2 + dy**2)
        geometria[numero] = {"longitud": L,
                             "c": dx / L, "s": dy / L}
    return geometria

La misma función sirve ahora para cualquier armadura:

geometria = calcular_geometria(nodos, barras)
geometria_b = calcular_geometria(nodos_b, barras_b)
geometria_b[6]

{'longitud': 5.0, 'c': -0.6, 's': 0.8}

5.3 NumPy

Las listas de Python no saben multiplicar matrices. Para eso usamos NumPy, la biblioteca de cálculo numérico sobre la que se apoya casi todo el Python científico. Se importa con el nombre corto np, y con np.set_printoptions pedimos que muestre los resultados con un decimal y sin notación científica:

import numpy as np

np.set_printoptions(precision=1, suppress=True)

El objeto central de NumPy es el arreglo (array), que creamos con np.array a partir de una lista de filas. Su atributo shape indica sus dimensiones:

k_local = np.array([[1, -1],
                    [-1, 1]])
k_local.shape

(2, 2)

Las operaciones con un número se aplican a cada elemento. Así obtenemos $\mkl$ para una barra con $EA/L = 80\,000\,\mathrm{kN/m}$ :

80000 * k_local

array([[ 80000, -80000],
       [-80000,  80000]])

El producto matricial se escribe con el operador @, y la transpuesta de un arreglo T, con T.T. Comprobemos la ecuación (5.2) con la barra 3 del ejemplo A, cuyo nodo final se desplaza 3 mm en $x$ y 4 mm en $y$ :

# Barra 3 del ejemplo A: c = 0.6, s = 0.8
T = np.array([[0.6, 0.8, 0.0, 0.0],
              [0.0, 0.0, 0.6, 0.8]])
# El nodo final se mueve 3 mm en x y 4 mm en y
u = np.array([0.0, 0.0, 3.0, 4.0])
print(T.shape, u.shape, T.T.shape)
T @ u

(2, 4) (4,) (4, 2)
array([0., 5.])

El nodo final se aleja 5 mm del inicial en la dirección de la barra: el desplazamiento $(3, 4)$ es paralelo a ella, y su módulo es 5.

5.4 La matriz de rigidez en Python

Traducimos la teoría del apartado 5.1 a tres funciones, una por cada matriz. Cada una es corta y hace una sola cosa, lo que facilita leerlas y verificarlas por separado:

def rigidez_local(E, A, L):
    """Matriz de rigidez de la barra en sus ejes (2x2)."""
    return E * A / L * np.array([[1.0, -1.0],
                                 [-1.0, 1.0]])


def matriz_transformacion(c, s):
    """Matriz T (2 x 4) tal que u' = T u."""
    return np.array([[c, s, 0.0, 0.0],
                     [0.0, 0.0, c, s]])


def rigidez_barra(E, A, L, c, s):
    """Matriz de rigidez de la barra en ejes globales."""
    T = matriz_transformacion(c, s)
    return T.T @ rigidez_local(E, A, L) @ T

La última línea es casi idéntica a la ecuación clave: $\mT\transp\mkl\mT$ . Probemos con la barra 1 del ejemplo A, que es horizontal:

g = geometria[1]
k1 = rigidez_barra(2.0e8, 0.0020, g["longitud"],
                   g["c"], g["s"])
k1

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

Con $c = 1$ y $s = 0$ , solo sobreviven los términos en $c^2$ , que valen $EA/L = 2.0\times10^8 \cdot 0.0020/6 = 66666.7\,\mathrm{kN/m}$ : una barra horizontal no ofrece rigidez a los desplazamientos verticales de sus nodos.

Para obtener las matrices de todas las barras, envolvemos el bucle en una función más, que guarda cada matriz en un diccionario con el número de barra como clave:

def calcular_rigideces(barras, geometria):
    """Matriz de rigidez global (4 x 4) de cada barra."""
    rigideces = {}
    for numero, barra in barras.items():
        g = geometria[numero]
        rigideces[numero] = rigidez_barra(
            barra["elasticidad"], barra["area"],
            g["longitud"], g["c"], g["s"])
    return rigideces


geometria = calcular_geometria(nodos, barras)
rigideces = calcular_rigideces(barras, geometria)

Las funciones calcular_geometria y calcular_rigideces son las dos primeras etapas de nuestra calculadora. A partir de ahora, los datos fluyen de una a otra: los diccionarios del modelo entran a la primera, y su resultado, junto con las barras, entra a la segunda.

5.5 Verificación

Toda matriz de rigidez de una barra debe cumplir tres propiedades, que podemos comprobar con NumPy sin hacer ningún cálculo a mano:

Simetría. $\mk = \mk\transp$ . La función np.allclose compara dos arreglos admitiendo las pequeñas diferencias del redondeo.
Singularidad. Una barra suelta puede moverse sin deformarse, así que $\mk$ no se puede invertir. Su rango (rank), la cantidad de filas independientes, es 1: todas las filas son múltiplos de la misma.
Cuerpo rígido. Si la barra se traslada sin deformarse, no aparecen fuerzas: $\mk\,\vu = \vect{0}$ .

for numero, k in rigideces.items():
    simetrica = np.allclose(k, k.T)
    rango = np.linalg.matrix_rank(k)
    print(f"Barra {numero}: simétrica = {simetrica}, "
          f"rango = {rango}")

Barra 1: simétrica = True, rango = 1
Barra 2: simétrica = True, rango = 1
Barra 3: simétrica = True, rango = 1

# Traslación de la barra 2 como cuerpo rígido: 1 mm en x
traslacion = np.array([0.001, 0.0, 0.001, 0.0])
rigideces[2] @ traslacion

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

La singularidad no es un defecto: es la expresión matemática de que una barra sola no es una estructura. En el capítulo 8 veremos cómo los apoyos la eliminan.

Problema resuelto 5.1 · La barra 2 del ejemplo A

Calcular a mano los términos de la matriz de rigidez global de la barra 2 del ejemplo A y compararlos con el resultado de rigidez_barra.

Solución

La barra 2 tiene $L = 5\,\mathrm{m}$ , $c = -0.6$ y $s = 0.8$ , de modo que $EA/L = 2.0\times10^{8} \cdot 0.0020 / 5 = 80\,000\,\mathrm{kN/m}$ . Los tres términos distintos de la ecuación clave son:

g = geometria[2]
EA_L = 2.0e8 * 0.0020 / g["longitud"]
c, s = g["c"], g["s"]
print(f"EA/L = {EA_L:.0f} kN/m, c = {c:.1f}, s = {s:.1f}")
print(f"EA/L·c²  = {EA_L * c**2:9.1f}")
print(f"EA/L·cs  = {EA_L * c * s:9.1f}")
print(f"EA/L·s²  = {EA_L * s**2:9.1f}")
rigideces[2]

EA/L = 80000 kN/m, c = -0.6, s = 0.8
EA/L·c²  =   28800.0
EA/L·cs  =  -38400.0
EA/L·s²  =   51200.0
array([[ 28800., -38400., -28800.,  38400.],
       [-38400.,  51200.,  38400., -51200.],
       [-28800.,  38400.,  28800., -38400.],
       [ 38400., -51200., -38400.,  51200.]])

Verificación. Cada término de la matriz coincide, con su signo, con el patrón de la ecuación clave: $c^2$ y $s^2$ en la diagonal, y $cs$ negativo porque la barra sube hacia la izquierda.

Repaso y resumen del capítulo

La matriz de rigidez de una barra en ejes globales es $\mk = \mT\transp\mkl\mT$ , con $\mkl$ de $2\times2$ y $\mT$ de $2\times4$ .
Una función se define con def, recibe parámetros y devuelve un resultado con return. Sus variables internas son locales.
Con funciones, el mismo cálculo sirve para cualquier armadura: basta con llamar a calcular_geometria con otros diccionarios.
NumPy aporta los arreglos: np.array, shape, el producto matricial @ y la transpuesta .T. El operador * multiplica elemento por elemento.
La matriz de una barra es simétrica, de rango 1 y no produce fuerzas ante un movimiento de cuerpo rígido.

Ejercicios

Ejercicio 5.1

Escribir una función peso_barra(A, L, gamma=78.5) que devuelva el peso propio de una barra, en kN, y usarla para calcular el peso del ejemplo A. El valor después del signo igual es un valor por omisión: se usa si al llamar la función no se indica otro.

Ver la solución

def peso_barra(A, L, gamma=78.5):
    """Peso propio de una barra, en kN."""
    return gamma * A * L


peso_total = 0.0
for numero, barra in barras.items():
    L = geometria[numero]["longitud"]
    peso_total += peso_barra(barra["area"], L)
print(f"Peso propio: {peso_total:.3f} kN")

Peso propio: 2.512 kN

El resultado coincide con el del ejercicio 4.2, ahora con una función reutilizable.

Ejercicio 5.2

Mover el nodo final de la barra 3 del ejemplo A 1 mm en la dirección perpendicular a la barra, es decir, en la dirección $(-s, c)$ , y calcular las fuerzas que aparecen con $\mk\,\vu$ . Explicar el resultado.

Ver la solución

g = geometria[3]
c, s = g["c"], g["s"]
# El nodo final se mueve 1 mm perpendicular a la barra
perpendicular = 0.001 * np.array([0.0, 0.0, -s, c])
rigideces[3] @ perpendicular

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

Un desplazamiento perpendicular a la barra no cambia su longitud, de modo que no la deforma axialmente y no produce fuerzas. Los valores -0. son ceros: el signo es un detalle de la aritmética de punto flotante. Este movimiento corresponde a una pequeña rotación de la barra alrededor de su nodo inicial, otro movimiento de cuerpo rígido.

Ejercicio 5.3

Calcular la matriz de rigidez global de una barra vertical de 3 m, de acero y con 10 cm² de área, que va del nodo $(0, 0)$ al nodo $(0, 3)$ . ¿Qué términos son distintos de cero y por qué?

Ver la solución

# Barra vertical de 3 m: de (0, 0) a (0, 3)
k_vertical = rigidez_barra(2.0e8, 0.0010, 3.0, 0.0, 1.0)
k_vertical

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

Con $c = 0$ y $s = 1$ , solo sobreviven los términos en $s^2$ , asociados a los desplazamientos verticales: una barra vertical solo resiste desplazamientos en su propia dirección. $EA/L = 2.0\times10^{8} \cdot 0.0010 / 3 = 66666.7\,\mathrm{kN/m}$ .

Ejercicio 5.4

Calcular las matrices de rigidez de las seis barras del ejemplo B con las funciones del capítulo y verificar que todas son simétricas y de rango 1.

Ver la solución

geometria_b = calcular_geometria(nodos_b, barras_b)
rigideces_b = calcular_rigideces(barras_b, geometria_b)
for numero, k in rigideces_b.items():
    simetrica = np.allclose(k, k.T)
    rango = np.linalg.matrix_rank(k)
    print(f"Barra {numero}: simétrica = {simetrica}, "
          f"rango = {rango}")

Barra 1: simétrica = True, rango = 1
Barra 2: simétrica = True, rango = 1
Barra 3: simétrica = True, rango = 1
Barra 4: simétrica = True, rango = 1
Barra 5: simétrica = True, rango = 1
Barra 6: simétrica = True, rango = 1