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.

Cerchas planas¶

Introducción¶

En este Notebook se parte del programa de resortes simple discutido anteriormente y tras realizar cambios menores en algunos parámetros este se modifica para que resuelva estructuras correspondientes a cerchas planas. Además se introducen los conceptos de sistemas de referencia locales a los elementos y global para la estructura así como la relación entre ambos.

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

Entender el algoritmo de solución de estructuras como un proceso general de ensamblaje de rigideces en un sistema global de ecuaciones.
Reconocer la diferencia entre los sistemas de referencia local (propio de cada elemento) y global (para toda la estructura) y la necesidad de expresar las rigideces en único sistema de referencia.
Expresar relaciones fuerza-desplazamiento en los sistemas de referencia local y global.
Reconocer las modificaciones necesarias para convertir un programa fundamental de ensamblaje de resortes en uno para estructuras mas complejas.

Ensamblaje de elementos cercha¶

En lo que sigue se usarán las siguientes condiciones.

Se asumirán estructuras correspondientes a ensamblajes de elementos barra conectados por articulaciones en sus uniones.
Por sus condiciones geométricas se asume que los elementos solo tienen rigidez axial y que por ende estos pueden entenderse como equivalentes a resortes de rigidez $k=\frac{AE}{l}$ .

La siguiente figura muestra un elemento típico. El eje de referencia $x$ dispuesto en la dirección longitudinal del elemento representa el sistema local de referencia. Nótese que si se estudia el elemento en sus sistema local este es completamente equivalente al resorte discutido anteriormente.

Desplazamientos en el sistema de referencia local

Usando $k = \frac{AE}{l}$ se tiene que la relación fuerza-desplazamiento en el sistema local de referencia es de la forma:

$\begin{Bmatrix}f_1\\f_2\end{Bmatrix} = \frac{AE}{l} \begin{bmatrix} 1&-1\\ -1&1\end{bmatrix} \begin{Bmatrix}u_1\\u_2\end{Bmatrix}\, .$

Considere ahora la siguiente estructura o ensamblaje conformado por 2 elementos tipo barra. Asumiendo que se conocen las relaciones fuerza-desplazamiento para cada una de las barras de la cercha se requiere determinar el desplazamiento del vértice superior de la cercha y las fuerzas internas en los elementos.

Nótese que además del sistema de referencia local de cada elemento ahora se tiene un sistema de referencia único denotado como $X-Y$ y denominado el sistema de referencia Global.

Para obtener la rigidez total de la estructura es necesario considerar la contribución de todos los elementos en el sistema de referencia de la cercha. Será necesario entonces transformar la matriz de rigidez de cada elemento del sistema local al sistema común o global. Para esto es conveniente definir:

$U, F$ : Desplazamientos (o grados de libertad) y fuerzas en el sistema de referencia global.
$u, f$ : Desplazamientos (o grados de libertad) y fuerzas en el sistema de referencia local.

Estas variables se relacionan mediante la matriz de rotación $\lambda$ de acuerdo con;

$u=\lambda U\, .$

Suponiendo que aplicamos ahora un desplazamiento $\delta u$ y determinamos el trabajo de las fuerzas a lo largo de este desplazamiento de acuerdo con:

$W = \delta u^T f\, .$

Ahora, dado que el trabajo es una cantidad escalar y por lo tanto independiente del sistema de referencia es válido escribir:

$W = \delta U^T F=\delta u^T f\, .$

Usando la ecuación de transformación bajo rotación en esta última se obtiene:

$\delta U^T F=\delta U^T \lambda^T f$

de donde se concluye que

$F =\lambda^T f\, .$

Si ahora se consideran las relaciones fuerza-desplazamiento

$f = ku\, ,$

donde $k$ es la matriz de rigidez local se tiene que podemos escribir:

$\begin{align} &\lambda^Tf = \lambda^T k u\\ &\lambda^Tf = \lambda^T k \lambda U\\ &F = KU \end{align}$

de donde se obtiene la ecuación de transformación bajo rotación de las relaciones fuerza-desplazamiento en un sistema de referencia local y un sistema global

$K=\lambda^T k\lambda\, .$

Es importante observar que en el sistema global el elemento parece tener 2 grados de libertad por nodo mientras que en el sistema local solo existe el desplazamiento axial.

La información requerida para calcular $K$ es entregada como parámetros de entrada a la función uel como se describe a continuación.

Encuentre la matriz de transformación bajo rotación $\lambda$ requerida para la formulación de la matriz de rigidez en el sistema de referencia global.

Ejemplo¶

Las modificaciones que se deben aplicar al programa de resortes son aquellas relacionadas con el hecho de tener 2 grados de libertad por nodo.

(Los archivos de texto con los datos de entrada para las partículas, elementos, coeficientes de rigidez y cargas de este problema están almacenados en la carpeta files de este repositorio).

In [1]:

import numpy as np

Lea los archivos de entrada de la carpeta files.

In [2]:

def readin():
    nodes = np.loadtxt('files/Cnodes.txt', ndmin=2)
    mats = np.loadtxt('files/Cmater.txt', ndmin=2)
    elements = np.loadtxt('files/Celes.txt', ndmin=2)
    loads = np.loadtxt('files/Cloads.txt', ndmin=2)
    return nodes, mats, elements, loads

La subrutina eqcounter cuenta ecuaciones y genera el arreglo de condiciones de frontera.

In [3]:

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

Ahora la función DME calcula la matriz de ensamblaje de ecuaciones.

In [4]:

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

La función assembly usa el modelo y la matriz DME para calcular la matriz de rigidez global.

In [5]:

def assembly(elements, mats, nodes, neq, DME):
    IELCON = np.zeros((2), dtype=np.integer)
    KG = np.zeros((neq, neq))
    nels = elements.shape[0]
    nnodes = 2
    ndof = 4
    for el in range(nels):
        elcoor = np.zeros((nnodes, 2))
        im = np.int(elements[el , 2])
        par0 = mats[im, 0]
        par1 = mats[im, 1]
        for j in range(nnodes):
            IELCON[j] = elements[el , j+3]
            elcoor[j , 0] = nodes[IELCON[j], 1]
            elcoor[j , 1] = nodes[IELCON[j], 2]
        kloc = ueltruss2D(elcoor, par0, par1)
        dme = DME[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 rutina ueltruss2D usa las coordenadas de los nudos y los parámetros de material para calcular la matriz de rigidez local ya transformada al sistema de referencia global.

Agregue un comentario a cada línea relevante de los códigos de las siguientes rutinas y úselos para escribir los pseudocódigos correspondientes. En particular identifique el calculo de la matriz de transformación bajo rotación $\lambda$ .

In [6]:

def ueltruss2D(coord, A, Emod):
    vec = coord[1, :] - coord[0, :]
    length = np.linalg.norm(vec)
    nx = vec[0]/length
    ny = vec[1]/length
    Q = np.array([
        [nx, ny , 0 , 0],
        [0,  0, nx , ny]])
    kl = (A*Emod/length) * np.array([
        [1, -1],
        [-1, 1]])
    kG = Q.T @ kl @ Q
    return kG

La rutina loadassem forma el vector de cargas en los nudos.

In [7]:

def loadasem(loads, IBC, neq, nl):
    RHSG = np.zeros([neq])
    for cont in range(nl):
        il = int(loads[cont, 0])
        ilx = IBC[il , 0]
        ily = IBC[il , 1]
        if ilx != -1:
            RHSG[ilx] = loads[cont, 1]
        if ily != -1:
            RHSG[ily] = loads[cont, 2]
    return RHSG

El programa principal mantiene la misma estructura que el programa de resortes, es decir, se tienen siguientes pasos:

Lee el modelo.
Determina la matriz de ensamblaje DME.
Ensambla el sistema global de ecuaciones.
Determina los desplazamientos globales UG tras resolver el sistema global.

In [8]:

nodes, mats, elements, loads = readin()
DME, IBC, neq = DME(nodes, elements)
KG = assembly(elements, mats, nodes, neq, DME)
RHSG = loadasem(loads, IBC, neq, 1)
UG = np.linalg.solve(KG, RHSG)
print(UG)

[ 0.         -0.15018948]

Problemas propuestos¶

Problema 1¶

Implemente una función que calcule las fuerzas nodales en cada elemento y que verifique el equilibrio del sistema.

Problema 2¶

Determine la rigidez lateral de la estructura usando la relación:

$k = \frac{P}{\delta}\, .$

Problema 3¶

Introduzca las modificaciones necesarias a los archivos de entrada para adicionar una barra adicional al modelo y que conecte los nodos $0$ y $1$ . Luego use este nuevo modelo para determinar los desplazamientos. Comente los resultados.

Problema 4¶

Repare la cercha mostrada en la figura adicionando elementos o imponiendo restricciones apropiadas a los desplazamientos. (Cree un nuevo paquete de archivos de datos).

Referencias¶

Klaus-Jürgen Bathe (2006). Finite element procedures. Klaus-Jurgen Bathe. Prentice Hall International.
Juan Gómez, Nicolás Guarín-Zapata (2018). SolidsPy: 2D-Finite Element Analysis with Python, https://github.com/AppliedMechanics-EAFIT/SolidsPy.

Formato del notebook¶

La siguiente celda cambia el formato del Notebook.

In [9]:

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

Out[9]:

In [ ]: