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.

Pórticos planos¶

Introducción¶

En este Notebook se agrega un elemento adicional al programa de análisis estructural, iniciado con el problema simple de resortes unidimesionales y extendido para considerar armaduras planas. En particular, ahora se adicionan elementos tipo viga para abordar problemas en los cuales los elementos funcionen principalmente a flexión. Adicionalmente, en una implementación independiente, es posible acoplar el comportamiento a carga axial de la cercha con el comportamiento a flexión de la viga.

El modelo de viga que se discute en este notebook corresponde a un modelo de Euler-Bernoulli en el cual se desprecian las deformaciones por cortante y por lo tanto válido para elementos con secciones transversales de baja esbeltez.

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

Reconocer las modificaciones necesarias para convertir un programa fundamental de ensamblaje de resortes en uno para estructuras conformadas por vigas.
Resolver problemas simples de estructuras conformadas por vigas sometidas a cargas puntuales.

Ensamblaje de elementos tipo viga¶

Consideremos un elemento tipo viga en su sistema local de referencia como el mostrado en la figura. En el sistema local el elemento tiene 2 grados de libertad por nodo, correspondientes al desplazamiento transversal y una rotación con respecto a un eje perpendicular al plano de la imagen.

El vector de grados de libertad (o desplazamientos generalizados) del elemento es:

$u^T=\begin{bmatrix}v_1 &\theta_1 &v_2 &\theta_2\end{bmatrix}$

mientras que el vector de fuerza (momentos y cortantes) esta dado por:

$f^T=\begin{bmatrix}f_1 &m_1 &f_2 &m_2\end{bmatrix}$

En el sistema de referencia local la matriz de rigidez relacionando fuerzas con desplazamientos es:

$\begin{Bmatrix} f_1\\ m_1\\ f_2\\ m_2\end{Bmatrix} = \begin{bmatrix} 12\frac{EI}{\mathcal l^3} &6\frac{EI}{\mathcal l^3} &-12\frac{EI}{\mathcal l^3} &6\frac{EI}{\mathcal l^2}\\ 6\frac{EI}{\mathcal l^3}&4\frac{EI}{\mathcal l}&-6\frac{EI}{\mathcal l^2}&2\frac{EI}{\mathcal l}\\ 12\frac{EI}{\mathcal l^3}&-6\frac{EI}{\mathcal l^2}&12\frac{EI}{\mathcal l^3}&-6\frac{EI}{\mathcal l^2}\\ 6\frac{EI}{\mathcal l^2}&2\frac{EI}{\mathcal l}&-6\frac{EI}{\mathcal l^2}&4\frac{EI}{\mathcal l} \end{bmatrix} \begin{Bmatrix} v_1\\ \theta_1\\ v_2\\ \theta_2\end{Bmatrix}$

Nota: La matriz de rigidez puede formularse por diferentes métodos que serán cubiertos en el curso de Análisis de Estructuras.

Para obtener la matriz de rigidez global de la estructura es necesario considerar nuevamente la contribución de todos los elementos en el sistema global de referencia. Con ese propósito procedemos como con el problema de la cercha. Usando la matriz de transformación bajo rotación $\lambda$ se tiene que:

$K = \lambda^T k\lambda$

donde $K$ es la matriz de rigidez para el elemento tipo viga en el sistema global de referencia. Note que en este sistema el elemento tiene ahora 3 grados de libertad por nodo.

Estructura aporticada simple¶

Considere el siguiente modelo simple conformado por un ensamblaje de 3 elementos. (Los archivos de datos de entrada están disponibles en la carpeta files).

Se requiere determinar el desplazamiento lateral de la estructura cuando es sometida a la carga lateral $P$ .

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.

Las modificacions que se deben aplicar al código de resortes (o al de cerchas) solo estan relacionadas con el hecho de que ahora hay 3 grados de libertad en cada nodo.

In [1]:

import matplotlib.pyplot as plt
import numpy as np
import sympy as sym

In [2]:

%matplotlib inline

Lea los archivos de entrada de la carpeta files.

In [3]:

def readin():
    nodes = np.loadtxt('files/Fnodes.txt', ndmin=2)
    mats = np.loadtxt('files/Fmater.txt', ndmin=2)
    elements = np.loadtxt('files/Feles.txt', ndmin=2)
    loads = np.loadtxt('files/Floads.txt', ndmin=2)
    return nodes, mats, elements, loads

La función eqcounter cuenta ecuaciones y genera el arreglo de condiciones de frontera.

In [4]:

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

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

In [5]:

def DME(nodes, elements):

    nels = elements.shape[0]
    IELCON = np.zeros([nels, 2], dtype=np.integer)
    DME_mat = np.zeros([nels, 6], dtype=np.integer)
    neq, IBC = eqcounter(nodes)
    ndof = 6
    nnodes = 2
    for i in range(nels):
        for j in range(nnodes):
            IELCON[i, j] = elements[i, j+3]
            kk = IELCON[i, j]
            for l in range(3):
                DME_mat[i, 3*j+l] = IBC[kk, l]
    return DME_mat, IBC, neq

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

In [6]:

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 = 6
    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 = uelbeam2DU(elcoor, par0, par1)
        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 uelbeam2D usa las coordenadas de los nodos 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 funciones y úselos para escribir los pseudocódigos correspondientes. En particular identifique el calculo de la matriz de transformación bajo rotación $\lambda$ .

In [7]:

def uelbeam2DU(coord, I, Emod):
    """2D-2-noded beam element
       without axial deformation

    Parametros
    ----------
    coord : ndarray
      Cordenadas nodales (2, 2).
    A : float
      Area de la seccion transversal.
    Emod : float
      Modulo de elasticidad (>0).

    Returna
    -------
    kl : ndarray
      Matriz de rigidez local para el elemento (4, 4).

    """
    vec = coord[1, :] - coord[0, :]
    nx = vec[0]/np.linalg.norm(vec)
    ny = vec[1]/np.linalg.norm(vec)
    L = np.linalg.norm(vec)
    Q = np.array([
        [-ny, nx, 0,  0,  0, 0],
        [0,  0, 1.0,  0,  0, 0],
        [0,  0, 0, -ny, nx, 0],
        [0,  0, 0,  0,  0, 1.0]])
    kl = (I*Emod/(L*L*L)) * np.array([
        [12.0, 6, -12.0, 6*L],
        [6,  4*L*L, -6*L, 2*L*L],
        [-12.0,  -6*L, 12.0, -6*L],
        [6*L,  2*L*L, -6*L, 4*L*L]])
    kG = np.dot(np.dot(Q.T, kl), Q)
    return kG

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

In [8]:

def loadasem(loads, IBC, neq, nl):
    RHSG = np.zeros([neq])
    for i in range(nl):
        il = int(loads[i, 0])
        ilx = IBC[il, 0]
        ily = IBC[il, 1]
        ilT = IBC[il, 2]
        if ilx != -1:
            RHSG[ilx] = loads[i, 1]
        if ily != -1:
            RHSG[ily] = loads[i, 2]
        if ilT != -1:
            RHSG[ilT] = loads[i, 3]
    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_mat.
Ensambla el sistema global de ecuaciones.
Determina los desplazamientos globales UG tras resolver el sistema global.

In [9]:

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

[ 2.25000000e+01  2.00000000e+01  1.27557539e-16 -1.25918779e-15
  2.00000000e+01  4.88970566e-16]

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 siguiente expresión:

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

Problema 3¶

Repotencie la estructura de tal forma que la rigidez lateral se incremente por un factor de 2.0

Problema 4¶

Repare el pórtico mostrado en la figura adicionando elementos y/o imponiendo restricciones apropiadas a los desplazamientos. (Cree un nuevo paquete de archivos de datos).

Problema 5¶

Para el portico de la figura asuma que la conexión del nudo 5 es articulada e indique como esta condición cambia los resultados.

Esquema del pórtico para el ejercicio propuesto.

Referencias¶

Bathe, Klaus-Jürgen. (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 [10]:

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

Out[10]: