Notebook

Contenido bajo licencia Creative Commons BY 4.0 y código bajo licencia MIT. © Juan Gómez y Nicolás Guarín-Zapata 2020. Este material es parte del curso Modelación Computacional en el programa de Ingeniería Civil de la Universidad EAFIT.

Ensamblaje de la matriz global¶

Este notebook explica cómo construir el sistema de ecuaciones para un sistema de masas y resortes, proceso que suele denominarse ensamblaje de ecuaciones. Para ello, se presenta un ejemplo paso a paso en el que se calculan las matrices de rigidez de cada elemento y luego se suman (ensamblan) a la matriz global.

Adicionalmente, el notebook auxiliar presenta de manera automática las matrices de rigidez locales para cada elemento y la matriz de rigidez global luego de ensamblar el elemento correspondiente.

Nota: Este notebook utiliza los mismos archivos de texto del sistema de resortes del notebook principal y disponibles en la carpeta files del repositorio. Para poder continuar con el resto del notebook es necesario que estos archivos estén disponibles en memoria.

In [1]:

import numpy as np

In [2]:

def readin():
    nodes = np.loadtxt('files/sprnodes.txt', ndmin=2)
    mats = np.loadtxt('files/sprmater.txt', ndmin=2)
    elements = np.loadtxt('files/spreles.txt', ndmin=2, dtype=np.int)
    loads = np.loadtxt('files/sprloads.txt', ndmin=2)
    return nodes, mats, elements, loads

Introducción¶

En este notebook se describen los detalles del proceso de ensamblaje de la matriz de rigidez global en el cual se suma la contribución de cada uno de los resortes. Mientras mas elementos se sumen, mas rígido será el sistema resultante. El proceso simplemente obtiene la "dirección" en la matriz de rigidez global en la que se deben ubicar los diferentes valores de las matrices locales de cada uno de los resortes. Estas direcciones se encuentran almacenadas, en forma de numero de ecuación, en un arreglo que tiene tantas filas como resortes tenga el sistema.

Al completar este notebook usted debería estar en la capacidad de:

Entender el papel que desempeñan en el programa los arreglos de condiciones de frontera, de conectividades y de ensamblaje de ecuaciones, IBC, IELCON y DME_mat respectivamente.
Entender los detalles del algoritmo de ensamblaje de matrices de rigidez globales.

Sistema estructural = sistema de ecuaciones¶

Todo programa de análisis estructural basado en el método de rigidez resuelve un sistema de ecuaciones de la forma general

$[K_G] \{U_G\} = \{F_G\}\, ,$

en el que $[K_G]$ es la matriz de coeficientes del sistema y representa la rigidez de toda la estructura, la cual resulta del aporte de todos los elementos del sistema; $\{U_G\}$ es el vector de desplazamientos (o grados de libertad) del sistema; y $\{F_G\}$ es el vector de fuerzas externas aplicadas al sistema.

Una vez resuelto el sistema de ecuaciones y se determinen los desplazamiento es posible usar relaciones fuerza-desplazamiento para cada elemento y determinar así las fuerzas internas. En las secciones que siguen se explican los pasos necesarios para armar el sistema de ecuaciones en un programa típico de análisis estructural. A esta operación se le conoce como ensamblaje del sistema.

Estructura de un programa y arreglos fundamentales¶

Considere el sistema masa-resortes discutido en el Notebook principal y mostrado nuevamente en la figura. En esta sección discutiremos solamente el proceso de ensamblaje de la de matriz rigidez global mediante la suma de las contribuciones de las rigideces de cada uno de los resortes.

En la figura los números encerrados en hexágonos azules representan masas o puntos de conexión de resortes, también denominados nodos. En el sistema del ejemplo el nodo $0$ se encuentra "amarrado" y por ende no puede desplazarse, impidiendo a la vez que bajo la acción de cargas el sistema se desplace como un cuerpo rígido.

Además, hemos nombrado los resortes con letras para facilitar la explicación. Para referirnos a los elementos de la matriz local de cada resorte usaremos subíndices para denotar la posición del elemento en la matriz local y un superíndice para indicar el nombre del resorte al que pertenece el elemento. Por ejemplo la matriz de rigidez local del resorte $A$ la representaremos como:

$k^A=k\begin{bmatrix} 1 &-1\\ -1 &1 \end{bmatrix} \equiv\begin{bmatrix} k_{00}^A &k_{01}^A\\ k_{10}^A &k_{11}^A \end{bmatrix}$

de manera que $k_{01}^A$ hace referencia a un coeficiente de rigidez que pertenece al resorte $A$ y que en su matriz local esta ubicado en la fila $0$ columna $1$ . Nótese, además, que hemos iniciado el conteo de elementos desde $0$ para ser consistentes con lo usado en Python.

El objetivo del ejemplo es indicar en que posiciones de la matriz global o del sistema de ecuaciones debemos colocar los elementos de las matrices de cada resorte.

Para realizar el ensamblaje necesitamos usar varios arreglos que se explican a continuación.

Arreglo de identificadores de ecuaciones: `IBC`¶

El arreglo de ecuaciones IBC contiene los identificadores de las ecuaciones asignadas a cada nodo o masa. Los nodos que se encuentren restringidos no tendrán ecuación asignada (pues su valor, nulo, es conocido) mientras que a cada nodo cuyo desplazamiento sea desconocido se le asignará un numero de ecuación. En este caso esta información se especifica al programa mediante un código que indica si el nodo se puede desplazar (valor de $0$ ) o si se encuentra restringido (valor de $-1$ ). Esta información hace parte del modelo y se encuentra en el archivo de nodos. El programa extrae estos códigos del archivo de nodos y los pasa a una primera versión del arreglo IBC. Este tendrá tantas filas como nodos y tantas columnas como grados de libertad tenga el problema. En este ejemplo cada nodo podrá tener asignado un solo grado de libertad.

En el archivo de nodos del presente ejemplo identifique donde se especifica la información de ecuaciones asignadas a cada nodo. Compare este información contra el sistema mostrado en la figura.

Posteriormente, la función eqcounter() utiliza la primera versión del arreglo IBC y asigna números de ecuaciones a los grados de libertad activos (indicados por código 0) y los almacena en el mismo arreglo IBC. En otras palabras, esta función cambia los ceros por números de ecuaciones.

Por ejemplo, el nodo $0$ se encuentra restringido como se indica con el valor de $-1$ . Por lo tanto, a este nodo no se le asignará una ecuación. Por su parte, el nodo 1 tiene asignado un valor de $0$ y a este se le asignará el primer número de ecuación correspondiente también a $0$ . Siguiendo con el conteo, el nodo 2 también tiene asignado un valor de $0$ y dado que es el siguiente nodo en la lista al mismo se le asignará la ecuación $1$ .

La rutina es la siguiente.

In [3]:

def eqcounter(nodes):
    nn = nodes.shape[0]
    IBC = np.zeros((nn, 1), dtype=np.integer)
    neq = 0
    for cont in range(nn):
        IBC[cont] = int(nodes[cont, 2])
        if IBC[cont] == 0:
            IBC[cont] = neq
            neq = neq + 1
    return neq, IBC

In [4]:

nodes, mats, elements, loads = readin()
neq, IBC = eqcounter(nodes)
IBC

Out[4]:

array([[-1],
       [ 0],
       [ 1],
       [ 2]])

Arreglo de conectividades de los elementos: `IELCON`¶

El siguiente arreglo necesario en el proceso es el arreglo que almacena los nodos a los que está conectado cada resorte. En el lenguaje comúnmente usado en el método de rigideces, a los nodos que definen un resorte se le conocen como las conectividades del elemento. En este programa nos referiremos al arreglo de conectividades como IELCON. Nuevamente, note que cada fila del arreglo almacena la información de un elemento. Por ejemplo, según la información almacenada en IELCON el resorte B está definido por los nodos 1 y 2.

En el archivo de elementos del presente ejemplo identifique donde se especifica la información de conectividades de cada elemento. Compare este información contra el sistema del ejemplo.

Arreglo de ensamblaje: `DME_mat`¶

El siguiente paso en el proceso es la creación de la matriz ensambladora denominada como DME_mat. Este arreglo es una combinación del arreglo de conectividades IELCON y del arreglo de ecuaciones nodales IBC.

Cada fila de esta matriz contiene los identificadores de ecuaciones asociados con cada uno de los nodos del elemento y será la que indique como distribuir los diferentes elementos de las matrices locales en la matriz global.

Para formar esta matriz se recorre un elemento a la vez el arreglo de conectividades IELCON; se identifican los nodos correspondientes a cada elemento y posteriormente se extraen las ecuaciones correspondientes a cada uno de estos nodos del arreglo IBC.

El siguiente bloque de código ensambla la matriz DME_mat.

In [5]:

def DME(nodes, elements):
    nels = elements.shape[0]
    DME_mat = np.zeros((nels, 2), dtype=np.integer)
    neq, IBC = eqcounter(nodes)
    ndof = 2
    nnodes = 2
    for ele in range(nels):
        for node in range(nnodes):
            pos = elements[ele, node + 3]
            DME_mat[ele, node] = IBC[pos]
    return DME_mat, IBC, neq

In [6]:

DME_mat, IBC, neq = DME(nodes, elements)
print(DME_mat)

[[-1  0]
 [ 0  1]
 [ 0  1]
 [ 1  2]]

Localización de la matriz de un elemento en la matriz global¶

En el proceso de ensamble cada elemento de la matriz local es colocado en la dirección de la matriz global extraída del arreglo DME_mat.

La matriz resultante para el problema del ejemplo se muestra en la figura y en lo que sigue se muestra el proceso paso a paso.

Resorte A¶

El resorte A esta conectado a los nodos 0 y 1. Sin embargo, en el nodo 0 no hay ecuación asignada como se indica por el valor de -1 en el arreglo DME_mat, mientras que al nodo 1 se le asigna la ecuación 0. Por lo tanto, como se ilustra en la figura el único valor de la matriz local del resorte A que es necesario ensamblar es el $k_{11}^A$ . En este proceso, en el que las componentes de la matriz asociadas a grados de libertad restringidos no se ensamblan, evita tener que eliminar dichos grados de libertad de la matriz global en un paso posterior.

Resorte B¶

El resorte B esta conectado a los nodos 1 y 2 los cuales tienen asignadas las ecuaciones 0 y 1. Por lo tanto esta matriz local se localizará en las direcciones $(0,0)$ , $(0,1)$ y $(1,0)$ .

Resorte C¶

Nótese que los resortes B y C se encuentran dispuestos en paralelo, lo que hace el sistema más rígido y por ende sus valores se localizan en las misma posiciones de la matriz global.

Resorte D¶

El resorte D está conectado a los nodos 2 y 3 los cuales tienen asignadas las ecuaciones 1 y 2. Por lo tanto, esta matriz local se localizará en las direcciones $(1,1)$ , $(1,2)$ y $(2,2)$ . Nótese ahora que cada que adicionamos al sistema un resorte en serie este se conecta a través de una sola posición con la matriz de rigidez global la cual a la vez aumenta su tamaño.

La función assembly() que se se presenta a continuación realiza la operación de cálculo de las matrices de rigidez locales y al mismo tiempo las ensambla en la matriz de rigidez global. Note que el algoritmo está conformado por un ciclo principal que recorre el sistema elemento por elemento. Una vez dentro del ciclo la rutina obtiene de los arreglos que almacenan el modelo la información del elemento actual, calcula su matriz local (ver uelspring()) y posteriormente obtiene su dirección en la matriz global. La última parte del ciclo asigna cada elemento a su posición en la matriz.

In [7]:

def assembly(elements, mats, nodes, neq, DME_mat):
    IELCON = np.zeros([2], dtype=np.integer)
    KG = np.zeros((neq, neq))
    nels = elements.shape[0]
    nnodes = 2
    ndof = 2
    for el in range(nels):
        elcoor = np.zeros((nnodes))
        im = np.int(elements[el, 2])
        par = mats[im]
        for j in range(nnodes):
            IELCON[j] = elements[el, j+3]
            elcoor[j] = nodes[IELCON[j], 1]
        kloc = uelspring(par)
        dme = DME_mat[el, :ndof]
        for row in range(ndof):
            glob_row = dme[row]
            if glob_row != -1:
                for col in range(ndof):
                    glob_col = dme[col]
                    if glob_col != -1:
                        KG[glob_row, glob_col] = KG[glob_row, glob_col] +\
                                                 kloc[row, col]
    return KG

La función uelspring calcula para cada elemento su matriz de rigidez local para que esta sea ensamblada posteriormente usando las direcciones extraídas del arreglo DME_mat.

In [8]:

def uelspring(kcof):
    """
    Esta rutina calcula la matriz de rigidez para
    un elemento tipo resorte conformado por 2 nodos.

    Parámetros
    ----------
    kcof : float
      Coeficiente de rigidez del resorte (>0).

    Retorna
    -------
    kloc : ndarray
      Matriz de coeficientes de rigidez local para
      el elemento tipo resorte (2, 2).


    """
    kloc = np.array([
        [kcof, -kcof],
        [-kcof, kcof]])
    return kloc

Similarmente, la función loadassem ensambla las cargas aplicadas sobre cada masa. Para esto la función utiliza el arreglo IBC que almacena el identificador de ecuación correspondiente a cada masa.

In [9]:

def loadasem(loads, IBC, neq, nl):
    """
    Ensambla el vector de cargas o de valores conocidos
    en el sistema global de ecuaciones.

    Parámetros
    ----------
    loads : ndarray
      Arreglo con las cargas impuestas a las partículas.
    IBC : ndarray (int)
      Arreglo que almacena el identificador de grado de libertad
      a cada partícula.
    neq : int
      Numero de ecuaciones en el sistema.
    nl : int
      Numero de cargas en el sistema.

    Retorna
    -------
    rhs_glob : ndarray
      Arreglo con las cargas impuestas a las partículas,
      rhs se refiere a lado derecho (del inglés right-hand-side).

    """
    rhs_glob = np.zeros((neq))
    for cont in range(nl):
        il = int(loads[cont, 0])
        ilx = IBC[il]
        if ilx != -1:
            rhs_glob[ilx] = loads[cont, 1]
    return rhs_glob

Programa principal¶

El programa principal consiste entonces en un ensamblador de la matriz global, o lo que es lo mismo, en un ensamblador de la matriz de coeficientes del sistema global de ecuaciones y en un ensamblador de un vector de valores conocidos, correspondiente en este problema a el vector de cargas.

In [10]:

KG = assembly(elements, mats, nodes, neq, DME_mat)
RHSG = loadasem(loads, IBC, neq, 3)

Una vez ensamblado el sistema de ecuaciones este se resuelve para determinar el vector de desplazamientos conocidos.

In [11]:

UG = np.linalg.solve(KG, RHSG)
print(UG)

[0.002  0.0025 0.0045]

Un aspecto importante a tener en cuenta es que en muchos programas de análisis estructural se ensambla la contribución de todos los elementos sin tener en cuenta los grados de libertad restringidos. Si en este punto se intenta resolver el sistema de ecuaciones resultante, la matriz de coeficientes será singular y no será posible invertirla reflejando el hecho de que el sistema no se encuentra en equilibrio ya que se puede desplazar como un cuerpo rígido: en términos matemáticos esto equivale a que el sistema de ecuaciones tiene infinitas soluciones. Para poder resolver el sistema es necesario, entonces, modificar la matriz de coeficientes y el vector de cargas de manera que se tengan en cuenta los grados de libertad restringidos (valores iguales a 0).

Alternativamente, se puede proceder como se ha discutido en este notebook, en donde en el proceso de ensamblaje ya se han tenido en cuenta las condiciones de apoyo o restricciones de cuerpo rígido del sistema. En este procedimiento el sistema de ecuaciones es soluble sin tener que realizar modificaciones adicionales.

Actividades para la clase¶

En las figuras se muestra una parte de un sistema de masas conectadas por resortes conjuntamente con su bloque en el arreglo DME_mat almacenando las direcciones de ensamblaje. Usando la notación adoptada en el notebook para referirse a los diferentes elementos de las matrices locales se pide realizar el ensamblaje de cada subsistema.

Subsistema 1

Subsistema 2

Glosario de términos¶

Ensamblaje: Operación mediante la cual se considera el aporte de un sistema de resortes a la rigidez de todo el sistema.

Conectividades: Identificadores de nodos que definen un elemento y se almacenan en un arreglo denominado IELCON. Su representación en términos de grados de libertad indica las direcciones de la matriz de rigidez global en las que se ensambla el elemento.

Operador ensamblador DME_mat: Arreglo que tiene tantas filas como elementos en el sistema y en cada fila contiene las direcciones en las que se ensambla la matriz local correspondiente al elemento en la global.

Formato del notebook¶

La siguiente celda cambia el formato del Notebook.

In [12]:

from IPython.core.display import HTML
def css_styling():
    styles = open('./nb_style.css', 'r').read()
    return HTML(styles)
css_styling()

Out[12]:

In [ ]: