Anexo D

Validación con OpenSees

A lo largo del libro verificamos la calculadora con el cálculo a mano (ejemplo A), con propiedades que deben cumplirse (simetría, equilibrio) y con el método de las secciones (ejemplo C). Este anexo agrega una verificación independiente: analizar los ejemplos B y C con otro programa, escrito por otras personas, y comparar los resultados.

D.1 OpenSees

OpenSees es un programa de análisis estructural de código abierto, desarrollado en la Universidad de California en Berkeley y muy usado en la investigación en ingeniería sísmica; su diseño se describe en McKenna (2011). Es gratuito y se usa desde Python con la biblioteca openseespy, de modo que cualquiera puede repetir esta comparación en Colab.

Para validar la calculadora, OpenSees tiene una ventaja sobre un programa comercial de pórticos: su elemento Truss es exactamente la barra de armadura del libro, articulada en sus extremos y con solo fuerza axial. No hay que liberar momentos ni anular el peso propio; basta con describir la misma armadura.

D.2 El modelo en OpenSees

La función analizar_con_opensees recibe los mismos diccionarios que analizar_armadura (capítulo 12) y devuelve los desplazamientos, las reacciones y las fuerzas en el mismo formato, para compararlos directamente:

import openseespy.opensees as ops


def analizar_con_opensees(nodos, barras, apoyos, cargas):
    """Analiza la armadura con OpenSees. Recibe los mismos
    datos que analizar_armadura y devuelve U, R y N en
    el mismo formato, para compararlos."""
    ops.wipe()  # borra el modelo anterior
    ops.model("basic", "-ndm", 2, "-ndf", 2)
    for nodo, punto in nodos.items():
        ops.node(nodo, punto["x"], punto["y"])
    for nodo, apoyo in apoyos.items():
        ops.fix(nodo, int(apoyo["x"]), int(apoyo["y"]))
    for numero, barra in barras.items():
        # Un material elástico por barra, con su E
        ops.uniaxialMaterial("Elastic", numero,
                             barra["elasticidad"])
        ops.element("Truss", numero,
                    barra["nodo_inicial"],
                    barra["nodo_final"], barra["area"],
                    numero)
    ops.timeSeries("Linear", 1)
    ops.pattern("Plain", 1, 1)
    for nodo, carga in cargas.items():
        ops.load(nodo, carga["x"], carga["y"])

    # Análisis estático y lineal, en un solo paso
    ops.system("BandSPD")
    ops.numberer("RCM")
    ops.constraints("Plain")
    ops.integrator("LoadControl", 1.0)
    ops.algorithm("Linear")
    ops.analysis("Static")
    ops.analyze(1)
    ops.reactions()

    # Los resultados, ordenados como los de la calculadora
    U = np.zeros(2 * len(nodos))
    for nodo in nodos:
        ix, iy = indices(gdl_del_nodo(nodo))
        U[ix], U[iy] = ops.nodeDisp(nodo)
    _, restringidos = clasificar_gdl(nodos, apoyos)
    R = []
    for g in restringidos:
        nodo = (g + 1) // 2
        direccion = 1 if g % 2 == 1 else 2  # x o y
        R.append(ops.nodeReaction(nodo, direccion))
    N = {}
    for numero in barras:
        fuerza = ops.eleResponse(numero, "axialForce")
        N[numero] = fuerza[0]
    return {"U": U, "R": np.array(R), "N": N,
            "restringidos": restringidos}

Cada paso tiene su equivalente en la calculadora. ops.node y ops.fix definen los nodos y los apoyos (con 1 para un grado de libertad restringido); ops.element("Truss", ...) crea cada barra con su área y un material elástico con su módulo $E$ ; ops.load aplica las cargas. Las seis líneas siguientes eligen cómo resolver el sistema (un análisis estático y lineal en un solo paso, que equivale a resolver $\mK_{ll}\,\vU_l = \vF_l$ una vez), ops.analyze(1) lo resuelve y ops.reactions() calcula las reacciones. Al final, los resultados se ordenan como los de la calculadora: clasificar_gdl solo se usa para listar las reacciones en el mismo orden.

La función comparar_resultados muestra los dos análisis lado a lado:

def comparar_resultados(nodos, calculadora, opensees,
                        tolerancia=1e-9):
    """Muestra los dos análisis lado a lado e indica si
    coinciden."""
    filas = []
    for nodo in nodos:
        for direccion, g in zip("xy", gdl_del_nodo(nodo)):
            i = indices([g])[0]
            filas.append([f"u{direccion} {nodo} (mm)",
                          calculadora["U"][i] * 1000,
                          opensees["U"][i] * 1000])
    restringidos = calculadora["restringidos"]
    for posicion, g in enumerate(restringidos):
        filas.append([f"R{g} (kN)",
                      calculadora["R"][posicion],
                      opensees["R"][posicion]])
    for numero in calculadora["N"]:
        filas.append([f"N{numero} (kN)",
                      calculadora["N"][numero],
                      opensees["N"][numero]])
    mostrar_tabla(["Resultado", "Calculadora",
                   "OpenSees"], filas, decimales=3,
                  ancho=13)
    diferencia = 0.0
    for _, valor_calculadora, valor_opensees in filas:
        diferencia = max(diferencia, abs(
            valor_calculadora - valor_opensees))
    if diferencia < tolerancia:
        print("\nLos dos análisis coinciden.")
    else:
        print(f"\nDifieren hasta en {diferencia:.3g}.")

resultados_b = analizar_armadura(nodos_b, barras_b,
                                 apoyos_b, cargas_b)
opensees_b = analizar_con_opensees(nodos_b, barras_b,
                                   apoyos_b, cargas_b)
comparar_resultados(nodos_b, resultados_b, opensees_b)

    Resultado  Calculadora     OpenSees
    ux 1 (mm)        0.000        0.000
    uy 1 (mm)        0.000        0.000
    ux 2 (mm)        0.179        0.179
    uy 2 (mm)        0.000        0.000
    ux 3 (mm)        0.900        0.900
    uy 3 (mm)       -0.081       -0.081
    ux 4 (mm)        0.854        0.854
    uy 4 (mm)       -0.481       -0.481
      R1 (kN)      -30.000      -30.000
      R2 (kN)        0.000        0.000
      R4 (kN)       80.000       80.000
      N1 (kN)       23.889       23.889
      N2 (kN)      -48.148      -48.148
      N3 (kN)       -6.111       -6.111
      N4 (kN)       -8.148       -8.148
      N5 (kN)       10.185       10.185
      N6 (kN)      -39.815      -39.815

Los dos análisis coinciden.

D.3 Ejemplo B

El ejemplo B (figura 3.2) es hiperestático, de modo que sus fuerzas dependen de la rigidez de las barras: es la prueba más exigente para la calculadora. La tabla D.1 compara los desplazamientos, y la tabla D.2, las reacciones y las fuerzas axiales.

Tabla D.1: Ejemplo B: desplazamientos de los nodos, en milímetros.

Nodo	$u_x$ calculadora	$u_x$ OpenSees	$u_y$ calculadora	$u_y$ OpenSees
1	0.000	0.000	0.000	0.000
2	0.179	0.179	0.000	0.000
3	0.900	0.900	−0.081	−0.081
4	0.854	0.854	−0.481	−0.481

Tabla D.2: Ejemplo B: reacciones y fuerzas axiales, en kilonewtons.

Resultado	Calculadora	OpenSees
Reacción $x$ del nodo 1	−30.00	−30.00
Reacción $y$ del nodo 1	0.00	0.00
Reacción $y$ del nodo 2	80.00	80.00
Barra 1	23.89	23.89
Barra 2	−48.15	−48.15
Barra 3	−6.11	−6.11
Barra 4	−8.15	−8.15
Barra 5	10.19	10.19
Barra 6	−39.81	−39.81

Los resultados coinciden en todas las cifras que muestran las tablas. La figura D.1 lo muestra también a la vista: la deformada de OpenSees se dibuja con un trazo ancho y translúcido, y la de la calculadora, con uno fino encima. Si fueran distintas, el trazo fino se saldría del ancho; en cambio, corre por su centro en todas las barras. La figura D.2 es el diagrama de fuerzas axiales de OpenSees, con el grosor de cada barra proporcional a su fuerza: es el mismo que dibuja nuestra calculadora (figura 10.5).

Figura D.1:

Ejemplo B: la deformada de OpenSees (trazo ancho) y la de la calculadora (trazo fino) coinciden.

Figura D.2:

Ejemplo B: fuerzas axiales y reacciones calculadas con OpenSees.

D.4 Ejemplo C

El ejemplo C (figura 12.2) es isostático: sus fuerzas dependen solo del equilibrio, y ya las verificamos con el método de las secciones en el capítulo 12. Compararlas con OpenSees confirma que ambos modelos tienen la misma geometría y las mismas cargas; la flecha, en cambio, depende de las áreas y del módulo de elasticidad, y verifica las rigideces.

Tabla D.3: Ejemplo C: flecha y reacciones.

Resultado	Calculadora	OpenSees
Desplazamiento $u_y$ del nodo 4 (mm)	−6.736	−6.736
Reacción $y$ del nodo 1 (kN)	62.50	62.50
Reacción $y$ del nodo 7 (kN)	62.50	62.50

Tabla D.4: Ejemplo C: fuerzas axiales, en kilonewtons.

Barra	Calculadora	OpenSees	Barra	Calculadora	OpenSees
1	62.50	62.50	12	−88.39	−88.39
2	62.50	62.50	13	25.00	25.00
3	100.00	100.00	14	−12.50	−12.50
4	100.00	100.00	15	0.00	0.00
5	62.50	62.50	16	−12.50	−12.50
6	62.50	62.50	17	25.00	25.00
7	−100.00	−100.00	18	53.03	53.03
8	−112.50	−112.50	19	17.68	17.68
9	−112.50	−112.50	20	17.68	17.68
10	−100.00	−100.00	21	53.03	53.03
11	−88.39	−88.39

También aquí los resultados coinciden. En la deformada (figura D.3), el trazo fino de la calculadora recorre el centro del de OpenSees en las 21 barras, y el diagrama de fuerzas axiales (figura D.4) muestra el patrón típico de una armadura Pratt: el cordón superior comprimido, el inferior traccionado, y el montante central sin fuerza.

Figura D.3:

Ejemplo C: las deformadas de OpenSees y de la calculadora coinciden.

Figura D.4:

Ejemplo C: fuerzas axiales y reacciones calculadas con OpenSees.

D.5 Todas las comparaciones de un vistazo

Las tablas muestran dos o tres decimales, pero la comparación puede ir mucho más lejos. La figura D.5 reúne todos los resultados de los dos ejemplos: los desplazamientos de todos los nodos, las reacciones y las fuerzas de todas las barras. Para cada uno dibuja la diferencia relativa entre los dos programas: la diferencia dividida por el mayor valor de su grupo, en escala logarítmica.

Figura D.5:

Diferencia relativa entre la calculadora y OpenSees en cada resultado de los ejemplos B y C (entre paréntesis, cuántos hay de cada tipo). Los círculos vacíos son resultados exactamente iguales.

Una validación razonable aceptaría diferencias de hasta el 1 % (la franja sombreada); las de nuestra calculadora están unos catorce órdenes de magnitud por debajo, alrededor de $10^{-16}$ . Ese número no es casual: es la precisión de la máquina, el error de redondeo de los números decimales de 64 bits con los que calculan Python y OpenSees. Dos programas que resuelven el mismo sistema por caminos distintos solo pueden diferir en ese redondeo, y eso es lo que encontramos.

La figura D.6 muestra lo mismo de otra forma: cada resultado es un punto, con el valor de la calculadora en el eje horizontal y el de OpenSees en el vertical. Si coinciden, el punto cae sobre la recta $y = x$ ; los 65 resultados están sobre ella.

Figura D.6:

Cada resultado de los ejemplos B y C, calculado con la calculadora y con OpenSees: todos caen sobre la recta $y = x$ .

D.6 Conclusión

Para una armadura hiperestática y para una isostática de 21 barras, la calculadora que construimos en este libro reproduce los desplazamientos, las reacciones y las fuerzas axiales de OpenSees hasta el redondeo de la máquina. No es casualidad: ambos resuelven el mismo sistema $\mK\,\vU = \vF$ con el mismo método. La diferencia es que ahora sabemos qué hace cada paso.