{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introducción al Análisis Estructural"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introducción"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"En este Notebook se introduce de manera intuitiva el concepto de rigidez mediante el problema simple de un resorte sometido a una carga constante. Posteriormente, la idea se extiende al caso de un sólido (o medio continuo). Luego, se sume que el el sólido está conformado por pedazos triangulares, cada uno de ellos definido por 3 \"partículas\", conectadas por \"resortes\". Esta idealización permite representar el sólido como un sistema de partículas conectadas por resortes. Una vez el problema se reduce a un sistema de masas-resortes se hace evidente que la solución de ambos sistemas usa el mismo tipo de algoritmo. En la parte final del Notebook se implementa un programa general para resolver sistemas representados en términos de relaciones fuerza-desplazamiento.\n",
"\n",
"\n",
"**Al completar este notebook usted debería estar en la capacidad de:**\n",
"\n",
"* Entender el concepto generalizado de rigidez tanto en sistemas continuos como discretos.\n",
"\n",
"* Usar razonamientos mecánicos para resolver problemas de equilibrio (estático) conformados por partículas conectadas por resortes.\n",
"\n",
"* Reconocer la equivalencia entre el problema de un sólido sometido a cargas externas y el de un sistema de masas-resortes.\n",
"\n",
"* Entender los pasos para resolver problemas de sistemas masa-resorte o su equivalente.\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Metodología para usar el notebook\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"El presente notebook esta diseñado para ser usado bajo una metodología de aprendizaje activo y mediante un esquema de clase invertida. El NB está conformado por una parte teórica y una parte computacional. Para poder alcanzar de manera satisfactoria los objetivos de aprendizaje es necesario que antes de la clase presencial el estudiante lea con detalle los contenidos teóricos del NB.\n",
"\n",
"<div class=\"alert alert-warning\">\n",
"Los bloques resaltados en amarillo corresponden a actividades para la clase, las cuales en algunos casos deben ser realizadas con el apoyo del profesor.\n",
"</div>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Concepto intuitivo de rigidez\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Considere los 2 modelos mostrados en la siguiente figura. El primero corresponde a un resorte de constante $k$ el cual tras ser sometido a la acción de una fuerza externa $F$ experimenta un cambio en su longitud igual a $\\delta$. El mismo problema se ilustra en el modelo abstracto de un sólido representado por la \"papa\" y en el cual la línea continua describe la forma de la \"papa\" antes de aplicar la fuerza externa $F$, mientras que la línea punteada describe la forma de la \"papa\" tras aplicar la fuerza. Los apoyos esquemáticos (de triángulo y rodillo) que se muestran en la figura indican que la \"papa\" no se podrá desplazar como un cuerpo rígido.\n",
"\n",
"<center>\n",
" <img src=\"img/rigidez.svg\"\n",
" alt=\"Ejemplo de rigidez\"\n",
" style=\"width:800px\">\n",
"</center>\n",
"\n",
"\n",
"En el resorte se satisface la relación\n",
"\n",
"\\begin{equation}\n",
"F = k \\delta\n",
"\\end{equation}\n",
"\n",
"y en la cual el coeficiente de rigidez $k$ representa la fuerza necesaria para producir un desplazamiento unitario $\\delta = 1.0$.\n",
"\n",
"En el caso de la \"papa\", el problema es mas complejo ya que existen muchas direcciones posibles y puntos en los cuales aplicar las fuerzas y en los cuales medir los desplazamientos. ¿Cómo podemos cuantificar la rigidez para este caso, entonces?\n",
"\n",
"Sea $\\delta _i$ el desplazamiento en un punto cualquiera medido en una dirección $i$ y sea $F_j$ la fuerza aplicada en la dirección $j$ también en un punto arbitrario de la \"papa\". En este caso siguiendo la idea del resorte es posible escribir:\n",
"\n",
"\\begin{equation}\n",
"F_j = k \\delta_i\\, .\n",
"\\end{equation}\n",
"\n",
"Comparando esta relación con la correspondiente al problema del resorte se puede concluir que el coeficiente $k$ representa la fuerza necesaria a lo largo de la dirección $j$ para producir un desplazamiento unitario en la dirección $i.$\n",
"\n",
"\n",
"Nótese que si se cambia la dirección en la que se aplica la fuerza y la dirección en la que se mide el desplazamiento es esperable que el coeficiente de rigidez también cambie. Por lo tanto un modelo mas general sería:\n",
"\n",
"\\begin{equation}\n",
"F_j = k_{ji} \\delta_i\\, .\n",
"\\end{equation}\n",
"\n",
"y en el cual los subíndices $ij$ dan cuenta de que el coeficiente de rigidez $k_{ij}$ representa la fuerza necesaria a lo largo de la dirección $j$ para producir un desplazamiento unitario en la dirección $i.$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Como describir el sólido?\n",
"\n",
"El modelo propuesto para el caso del sólido es físicamente correcto pero operativamente existen infinitas direcciones posibles $i$ y $j$ y surge entonces la pregunta: **¿cómo cuantificar la rigidez para el caso de la \"papa\"?**\n",
"\n",
"Consideremos algunas posibles soluciones:\n",
"\n",
"1. Solución discreta: Aproximar el sólido mediante un sistema finito de partículas conectadas por resortes.\n",
"\n",
"2. Solución del continuo: Estudiar el comportamiento de un elemento diferencial generando un sistema de ecuaciones diferenciales (en las variables $u_i$, $\\epsilon_{ij}$ y $\\sigma_{ij}$).\n",
"\n",
"3. Solución numérica: Combinar **1** y **2** para hacer una aproximación discreta del continuo (métodos por elementos finitos). La comparación entre las soluciones **1** y **2** se describe en la siguiente tabla:\n",
"\n",
"<div class=\"alert alert-warning\">\n",
"\n",
"Discutir con sus compañeros de grupo las posibles ventajas y desventajas de las soluciones **1** a la **3**.\n",
"\n",
"</div>\n",
"\n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"| Planteamiento del continuo | Planteamiento discreto |\n",
"|:--------------------------:|:----------------------:|\n",
"| Número infinito de elementos diferenciales. | Número finito de partículas conectadas por resortes. |\n",
"| Desplazamientos, deformaciones y tensiones. | Desplazamientos y fuerzas. |"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Elementos finitos (aproximando el continuo).\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"La figura muestra nuevamente el problema esquemático de la \"papa\". Esta ha sido dividida en un numero **finito** de triángulos conectados en las esquinas. Cada triángulo esta definido por 3 puntos correspondientes a sus vértices. Denominemos un triangulo típico como $\\Omega_e$. Asuma que cada vértice del triángulo típico $\\Omega_e$ se comporta como una partícula y que esta puede experimentar desplazamientos (y fuerzas) en la dirección horizontal y vertical. Asuma además que para el triángulo típico las 3 partículas que lo definen están conectadas por resortes.\n",
"\n",
"<div class=\"alert alert-warning\">\n",
"\n",
"¿Cuántos coeficientes del tipo $k_{ij}$ existen para cada triángulo?\n",
"\n",
"</div>\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<center>\n",
" <img src=\"img/dominio_discreto.svg\"\n",
" alt=\"Dominio discretizado.\"\n",
" style=\"width:400px\">\n",
"</center>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Relación fuerza-desplazamiento para todo el sistema\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Se tiene entonces que para un triangulo típico $\\Omega_e$ es posible relacionar las 6 fuerzas con los 6 desplazamientos mediante una ecuación general como:\n",
"\n",
"\\begin{equation}\n",
"F_j = k_ {ji} \\delta_i\n",
"\\end{equation}\n",
"\n",
"y en la cual el término $k_{ji}$ representa 36 coeficientes de rigidez. Si se considera que el índice $j$ puede tomar 6 valores correspondientes a 6 fuerzas y que el índice $i$ puede tomar 6 valores correspondientes a 6 desplazamientos entonces se puede reconocer que se tiene un sistema de 6 ecuaciones el cual podemos escribir matricialmente como:\n",
"\n",
"\\begin{equation}\n",
"\\{ F \\}=[ K ] \\{ u \\}\\, .\n",
"\\end{equation}\n",
"\n",
"De manera similar, es posible anticipar que es posible escribir una relación que considere todos los puntos del dominio. Por ejemplo, si se tienen $N$ puntos en total habrán $2N$ fuerzas y $2N$ desplazamientos. Escribamos la relación entre fuerzas y desplazamientos para todo el sólido (representado por $N$ puntos) como:\n",
"\n",
"\\begin{equation}\n",
"F_j^G = k_ {ji}^G \\delta_i\\, .\n",
"\\end{equation}\n",
"\n",
"En esta relación hemos usado el superíndice $G$ para enfatizar que las fuerzas, desplazamientos y coeficientes de rigidez corresponden a toda el sólido y no a un elemento típico. En forma matricial se tiene:\n",
"\n",
"\\begin{equation}\n",
"\\{ F^G \\} = [ K^G ] \\{ U^G \\}\\, .\n",
"\\end{equation}\n",
"\n",
"<div class=\"alert alert-warning\">\n",
"\n",
"¿Cuántos coeficientes del tipo $k_{ij}$ existen para todo el sólido si este esta representado por $N$ puntos?\n",
"\n",
"</div>\n",
"\n",
"El superínidce $G$ también indica que las ecuaciones son válidas para todo el sistema global representado por todos los triángulos en los que ha sido dividido el sólido.\n",
"\n",
"<div class=\"alert alert-warning\">\n",
" \n",
"* Identifique el argumento de tipo físico que permite formar el sistema global de ecuaciones:\n",
"\n",
"\\begin{equation}\n",
"F_j ^ G = k_{ji} ^ G \\delta_i\\, .\n",
"\\end{equation}\n",
"\n",
"* Si el sólido se divide en $M$ triángulos unidos a través de un total de $N$ vértices, ¿cuántos desplazamientos del tipo $\\delta_i$ hay en el sistema?\n",
"\n",
"* Si el sólido se divide en $M$ triángulos unidos a través de un total de $N$ vértices, ¿cuántas fuerzas del tipo $F_j ^ G$ hay en el sistema?\n",
"</div>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"De acuerdo con la discusión anterior, tanto el equilibrio de un triángulo típico $\\Omega_e$ como el de un resorte pueden formularse mediante una relación fuerza-desplazamiento de la forma:\n",
"\n",
"\\begin{equation}\n",
"\\{ F \\} = [ K ] \\{ u \\}\\, .\n",
"\\end{equation}\n",
"\n",
"En el caso del resorte es sencillo determinar la matriz de rigidez $[K]$.\n",
"\n",
"En el caso del sólido el problema es un poco más elaborado y por el momento nos conformaremos con suponer que esta matriz se encuentra disponible. En cualquier caso, sea que se trata de un sistema de partículas conectadas por resortes, o un sólido representado por triángulos, el problema se resuelve ensamblando la matriz de rigidez global $[K^G]$ y resolviendo el sistema de ecuaciones:\n",
"\n",
"\\begin{equation}\n",
"\\{ F^G \\} = [ K^G ] \\{ U^G \\}\\, .\n",
"\\end{equation}\n",
"\n",
"\n",
"En la siguiente sección se discute con más detalle el problema usando un sistema de partículas unidas por resortes y se implementa su solución en el computador."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div class=\"alert alert-warning\">\n",
" \n",
"### Sistema particula-resorte\n",
"\n",
" \n",
"La siguiente actividad permite entender desde un punto de vista físico el proceso de solución del problema de varias partículas interactuando a través de un sistema de resortes como el que se muestra en la figura.\n",
"\n",
"<center>\n",
"<img src=\"img/sistema_resortes.svg\"\n",
" alt=\"Sistema de masas y resortes\"\n",
" style=\"width:400px\">\n",
"</center>\n",
"\n",
"El sistema de varias partículas se encuentra conectado por resortes con coeficiente de rigidez $k$. El sistema se encuentra sometido a diferentes fuerzas externas (representadas en múltiplos de $P$). Por efectos de visualización las partículas se representan como carritos rectangulares.\n",
"\n",
"Se requiere:\n",
"\n",
"1. Plantear las ecuaciones de equilibrio para una partícula típica mostrada en la figura. Asuma que la partícula esta conectada a un resorte $i$ por la izquierda y a un resorte $i+1$ por la derecha de manera que el diagrama de cuerpo libre para la partícula es como el que se muestra en la figura.\n",
"\n",
"<center>\n",
"<img src=\"img/dcl_masa.svg\"\n",
" alt=\"Diagrama de cuerpo libre para una partícula\"\n",
" style=\"width:200px\">\n",
"</center>\n",
"\n",
"\n",
"\n",
"2. Plantear las relaciones fuerza-desplazamiento para el resorte típico mostrado en la figura. En esta los subíndices 1 y 2 hacen referencia a las fuerzas y desplazamientos de los puntos 1 y 2 que definen el resorte.\n",
"\n",
"<center>\n",
" <img src=\"img/elemento_resorte.svg\"\n",
" alt=\"Balance de fuerzas para un resorte\"\n",
" style=\"width:300px\">\n",
"</center>\n",
"\n",
"3. Enumerar las partículas y resortes del sistema. (Asumir el empotramiento como una partícula).\n",
"\n",
"4. Escribir las ecuaciones de equilibrio para cada una de las 4 partículas.\n",
"\n",
"5. Usando las relaciones fuerza-desplazamiento de los diferentes resortes re-escribir las ecuaciones del numeral **4**.\n",
"\n",
"</div>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Algoritmo de solución"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"En la solución del problema anterior se usaron planteamientos mecánicos para formular las ecuaciones de equlibrio del sistema de partículas. Consideremos ahora su solución de una manera sistemática de tal forma que el problema se pueda codificar en un programa general de análisis estructural.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Equilbrio de los resortes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Considere el siguiente sistema de 3 partículas, las cuales están conectadas por resortes $i$ e $i+1$ con coeficientes de rigidez $K^i$ y $K^{i+1}$. De manera similar, las partículas tienen asignados los nombres $j-1$, $j$ y $j+1$ y sus respectivos desplazamientos se denotan por $u_{j-1}$, $u_{j}$ y $u_{j+1}$.\n",
"\n",
"<center>\n",
" <img src=\"img/ibc.svg\"\n",
" alt=\"files\"\n",
" style=\"width:300px\">\n",
"</center>\n",
"\n",
"Escribamos las ecuaciones de equilibrio para los resortes $i$ e $i+1$ en términos de los desplazamientos $u_{j - 1}$, $u_j$ y $u_{j + 1}$ como:\n",
"\n",
"$$\\begin{Bmatrix} f_1^i\\\\ f_2^i \\end{Bmatrix}\n",
"= \\begin{bmatrix}\n",
" k_{11}^i & k_{12}^i\\\\\n",
" k_{21}^i & k_{22}^i\n",
"\\end{bmatrix}\n",
"\\begin{Bmatrix} u_{j - 1}\\\\ u_j \\end{Bmatrix}\\, ,$$\n",
"\n",
"y\n",
"\n",
"$$\\begin{Bmatrix} f_1^{i + 1}\\\\ f_2^{i + 1}\\end{Bmatrix}\n",
"= \\begin{bmatrix}\n",
" k_{11}^{i + 1} & k_{12}^{i + 1}\\\\\n",
" k_{21}^{i + 1} & k_{22}^{i + 1} \\end{bmatrix}\n",
"\\begin{Bmatrix} u_j\\\\ u_{j + 1} \\end{Bmatrix}\\, .$$\n",
"\n",
"Note que hemos usado notación fila-columna en los índices para los coeficientes de rigidez para simplificar su implementación en el computador."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Equilibrio de una partícula"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Asumiendo que la partícula $j$ se encuentra sometida a una carga externa $P$, se tiene, de su diagrama de cuerpo libre, que:\n",
"\n",
"$$\n",
"k_{21}^i u_{j - 1} + (k_{22}^i + k_{11}^{i + 1}) u_j + k_{12}^{i + 1} u_{j + 1} = P_j.\n",
"$$\n",
"\n",
"Procediendo de manera similar para las partículas ${j-1}$ y ${j+1}$ y considerando la contribución de los resortes $K^i$ y $K^{i+1}$ obtenemos el siguiente bloque del sistema completo de ecuaciones:\n",
"\n",
"$$\\begin{bmatrix}\n",
" k_{11}^i &k_{12}^i &0\\\\\n",
" k_{21}^i &k_{22}^i + k_{11}^{i + 1} &k_{12}^{i + 1}\\\\\n",
" 0 &k_{21}^{i + 1} &k_{22}^{i + 1}\n",
"\\end{bmatrix}\\, . $$\n",
"\n",
"Al considerar el sistema completo de partículas y resortes se obtiene un sistema de ecuaciones lineales de la forma general:\n",
"\n",
"$$\n",
"[K_G] \\{U_G\\} = \\{F_G\\}\\, .\n",
"$$\n",
"\n",
"En este sistema cada ecuación representa el equilibrio de una partícula.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Ensamblaje"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"La construcción de las matrices globales que describen el equilibrio de cada masa (o partícula) en el sistema puede programarse en el computador mediante un algoritmo que considere el aporte (en términos de fuerzas) de cada elemento (o resorte).\n",
"\n",
"En Mecánica Computacional, a este proceso se le denomina como **Ensamblaje** de las ecuaciones globales. La operación de ensamblaje puede realizarse después de identificar la conexión entre los (desplazamientos o) grados de libertad globales y los (desplazamientos o) grados de libertad locales mediante una matriz que almacena en cada fila los identificadores de los grados de libertad globales asociados a cada elemento.\n",
"\n",
"Por ejemplo, en el sistema de 3-partículas los resortes $i$ y $i+1$ tienen como desplazamientos de sus extremos los denotados como $u_{j-1}$ y $u_j$ y $u_j$ y $u_{j+1}$. Dichos índices contienen toda la información necesaria para realizar el ensamblaje. La matriz que almacena los indices globales para todos los elementos en el modelo es la matriz `DME` dada por:\n",
"\n",
"$$ DME = \\begin{bmatrix}\n",
" j - 1 &j\\\\\n",
" j &j + 1\\\\\n",
"\\end{bmatrix}\\, .$$\n",
"\n",
"Nótese que en esta matriz los datos de la primera fila mostrada son $j - 1$ y $j$, los cuales corresponden a los identificadores de los grados de libertad asociados con el resorte $i$. De la misma manera los datos de la segunda fila mostrada son $j$ y $j+1$ los cuales corresponden a los identificadores de los grados de libertad asociados con el resorte $i+1$.\n",
"\n",
"Una vez se tiene disponible la matriz `DME`, el proceso de ensamblaje se realiza como se describe a continuación:\n",
"\n",
"\\begin{align}\n",
"&K_{j - 1,j - 1} \\leftarrow K_{j - 1,j - 1} + k_{11}^i\\\\\n",
"&K_{j - 1,j} \\leftarrow K_{j - 1,j} + k_{12}^i\\\\\n",
"&K_{j,j - 1} \\leftarrow K_{j,j - 1} + k_{21}^i\\\\\n",
"&K_{j,j} \\leftarrow K_{j,j} + k_{22}^i\n",
"\\end{align}\n",
"\n",
"y\n",
"\n",
"\\begin{align}\n",
"&K_{j,j} \\leftarrow K_{j,j} + k_{11}^{i + 1}\\\\\n",
"&K_{j,j + 1} \\leftarrow K_{j,j + 1} + k_{12}^{i + 1}\\\\\n",
"&K_{j + 1,j} \\leftarrow K_{j + 1,j} + k_{21}^{i + 1}\\\\\n",
"&K_{j + 1,j + 1} \\leftarrow K_{j + 1,j + 1} + k_{22}^{i + 1}\n",
"\\end{align}\n",
"\n",
"\n",
"Note la conexión entre los índices locales, correspondientes a $1$ y $2$, y las posiciones en la matriz global, correspondientes a $j-1$, $j$ y $j+1$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Programa simple de análisis estructural"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"De las secciones anteriores se concluye que el problema de equilibrio se resuelve tras ensamblar la contribución de cada resorte (o triángulo) en una matriz de rigidez global. En esta sección discutimos su solución en el computador. Para definir el modelo usaremos archivos de texto en los cuales indicaremos al programa las partículas que lo conforman, los resortes del sistema y su conexión con las partículas, la localización de las cargas externas y las propiedades de los resortes.\n",
"\n",
"Inicialmente se describen las rutinas o funciones que conforman el programa y en la parte final estas son llamadas desde el programa principal.\n",
"\n",
"El programa principal ejecuta los siguientes pasos:\n",
"\n",
"* Lee el modelo.\n",
"\n",
"* Calcula la matriz `DME`.\n",
"\n",
"* Ensambla el sistema global de ecuaciones.\n",
"\n",
"* Resuelve el sistema para determinar los desplazamientos globales $UG$.\n",
"\n",
"(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)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" La rutina `readin()` lee los archivos de entrada que se encuentran en la carpeta `files` del respositorio y que deben ser almacenados localmente bajo esta misma estructura."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def readin():\n",
" nodes = np.loadtxt('files/sprnodes.txt', ndmin=2)\n",
" mats = np.loadtxt('files/sprmater.txt', ndmin=2)\n",
" elements = np.loadtxt('files/spreles.txt' , ndmin=2, dtype=np.int)\n",
" loads = np.loadtxt('files/sprloads.txt', ndmin=2)\n",
" return nodes, mats, elements, loads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"En el archivo de noos, en el cual se almacena la información de cada partícula, hay un código con valor $-1$ o $0$ indicando si la partícula se encuentra restringida (como en el empotramiento) o libre. Con esa información la rutina `eqcounter()` asigna identificadores de grados de libertad o números de ecuaciones a cada una de las partículas declaradas como libres.\n",
"\n",
"<div class=\"alert alert-warning\">\n",
"Agregue un comentario a cada línea relevante del código y use los comentarios para escribir el pseudocódigo de la función. \n",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def eqcounter(nodes):\n",
" nn = nodes.shape[0]\n",
" IBC = np.zeros((nn, 1), dtype=np.integer)\n",
" neq = 0\n",
" for cont in range(nn):\n",
" IBC[cont] = int(nodes[cont, 2])\n",
" if IBC[cont] == 0:\n",
" IBC[cont] = neq\n",
" neq = neq + 1\n",
" return neq, IBC"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"El numero de ecuación asignado a cada partícula es usado ahora para formar la matriz `DME`. Cada fila de esta matriz contiene los identificadores de los desplazamientos de los extremos del resorte.\n",
"\n",
"<div class=\"alert alert-warning\">\n",
"Agregue un comentario a cada línea relevante del código y use los comentarios para escribir el pseudocódigo de la función. \n",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def DME(nodes, elements):\n",
" nels = elements.shape[0]\n",
" DME_mat = np.zeros((nels, 2), dtype=np.integer)\n",
" neq, IBC = eqcounter(nodes)\n",
" ndof = 2\n",
" nnodes = 2\n",
" for ele in range(nels):\n",
" for node in range(nnodes):\n",
" pos = elements[ele, node + 3]\n",
" DME_mat[ele, node] = IBC[pos]\n",
" return DME_mat, IBC, neq"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Usando la matriz `DME` es posible ensamblar la matriz global de coeficientes de rigidez en términos de ecuaciones como:\n",
"\n",
"$$\n",
"K_{j - 1,j - 1} \\leftarrow K_{j - 1,j - 1} + k_{11}^i\\, .\\\\\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def assembly(elements, mats, nodes, neq, DME_mat):\n",
" IELCON = np.zeros([2], dtype=np.integer)\n",
" KG = np.zeros((neq, neq))\n",
" nels = elements.shape[0]\n",
" nnodes = 2\n",
" ndof = 2\n",
" for el in range(nels):\n",
" elcoor = np.zeros((nnodes))\n",
" im = np.int(elements[el, 2])\n",
" par = mats[im]\n",
" for j in range(nnodes):\n",
" IELCON[j] = elements[el, j+3]\n",
" elcoor[j] = nodes[IELCON[j], 1]\n",
" kloc = uelspring(par[0])\n",
" dme = DME_mat[el, :ndof]\n",
" for row in range(ndof):\n",
" glob_row = dme[row]\n",
" if glob_row != -1:\n",
" for col in range(ndof):\n",
" glob_col = dme[col]\n",
" if glob_col != -1:\n",
" KG[glob_row, glob_col] = KG[glob_row, glob_col] +\\\n",
" kloc[row, col]\n",
"\n",
" return KG"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"La función `uelspring` calcula la matriz de rigidez de un resorte típico, mientras que la función `loadassem` ensambla el vector de cargas externas aplicadas sobre las partículas.\n",
"\n",
"<div class=\"alert alert-warning\">\n",
"Agregue un comentario a cada línea relevante de los códigos de las siguientes función y úselos para escribir los pseudocódigos correspondientes. \n",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def uelspring(kcof):\n",
" \"\"\"\n",
" Esta rutina calcula la matriz de rigidez para\n",
" un elemento tipo resorte conformado por 2 nodos.\n",
"\n",
" Parámetros\n",
" ----------\n",
" kcof : float\n",
" Coeficiente de rigidez del resorte (>0).\n",
"\n",
" Retorna\n",
" -------\n",
" kloc : ndarray\n",
" Matriz de coeficientes de rigidez local para\n",
" el elemento tipo resorte (2, 2).\n",
"\n",
"\n",
" \"\"\"\n",
" kloc = np.array([\n",
" [kcof, -kcof],\n",
" [-kcof, kcof]])\n",
" return kloc"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def loadasem(loads, IBC, neq, nl):\n",
" \"\"\"\n",
" Ensambla el vector de cargas o de valores conocidos\n",
" en el sistema global de ecuaciones.\n",
"\n",
" Parámetros\n",
" ----------\n",
" loads : ndarray\n",
" Arreglo con las cargas impuestas a las partículas.\n",
" IBC : ndarray (int)\n",
" Arreglo que almacena el identificador de grado de libertad\n",
" a cada partícula.\n",
" neq : int\n",
" Numero de ecuaciones en el sistema.\n",
" nl : int\n",
" Numero de cargas en el sistema.\n",
"\n",
" Retorna\n",
" -------\n",
" rhs_glob : ndarray\n",
" Arreglo con las cargas impuestas a las partículas,\n",
" rhs se refiere a lado derecho (del inglés right-hand-side).\n",
"\n",
" \"\"\"\n",
" rhs_glob = np.zeros((neq))\n",
" for cont in range(nl):\n",
" il = int(loads[cont, 0])\n",
" ilx = IBC[il]\n",
" if ilx != -1:\n",
" rhs_glob[ilx] = loads[cont, 1]\n",
" return rhs_glob"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Programa principal"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.002 0.0025 0.0045]\n"
]
}
],
"source": [
"nodes, mats, elements, loads = readin()\n",
"DME, IBC, neq = DME(nodes, elements)\n",
"KG = assembly(elements, mats, nodes, neq, DME)\n",
"RHSG = loadasem(loads, IBC, neq, 3)\n",
"UG = np.linalg.solve(KG, RHSG)\n",
"print(UG)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div class=\"alert alert-warning\">\n",
"\n",
"### Actividades para resolver con el computador\n",
"\n",
"* Usando los desplazamientos encontrados en la solución del sistema determine las fuerzas resultantes en cada resorte e indique que elementos están a compresión y que elementos están a tensión.\n",
"\n",
"* Usando las fuerzas encontradas en el paso anterior verifique el equilibrio del sistema.\n",
"\n",
"* Libere el punto correspondiente al empotramiento y calcule nuevamente la solución. Discuta sus resultados.\n",
"\n",
"* Resuelva nuevamente el sistema haciendo la rigidez del resorte mas a la derecha igual a $5k$. Discuta como varían las fuerzas internas en comparación con las del sistema original.\n",
"\n",
"</div>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Glosario de términos"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Coeficiente de rigidez:** Fuerza experimentada por un resorte empotrado en su extremo y sometido en su extremo opuesto a un desplazamiento unitario.\n",
"\n",
"**Matriz de rigidez:** Matriz cuyos elementos son coeficientes de rigidez relacionando las fuerzas y los desplazamientos de un sistema de partículas.\n",
"\n",
"**Sistema discreto:** Sistema mecánico cuyas ecuaciones de equilibrio son representables de manera directa mediante relaciones fuerza-desplazamiento por medio de una matriz de rigidez.\n",
"\n",
"**Sistema continuo:** Sistema mecánico cuyas ecuaciones de equilibrio corresponden a ecuaciones diferenciales parciales que se satisfacen sobre el dominio del sistema.\n",
"\n",
"**Ensamblaje:** Operación mediante la cual se considera el aporte de un sistema de resortes a la rigidez de todo el sistema."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Formato del notebook"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"La siguiente celda cambia el formato del Notebook."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
"<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
"<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
"\n",
"<style>\n",
"\n",
"/*\n",
"Template for Notebooks for Modelación computacional.\n",
"\n",
"Based on Lorena Barba template available at:\n",
"\n",
" https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css\n",
"*/\n",
"\n",
"/* Fonts */\n",
"@font-face {\n",
"font-family: \"Computer Modern\";\n",
"src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
"}\n",
"\n",
"/* Text */\n",
"div.cell{\n",
"width:800px;\n",
"margin-left:16% !important;\n",
"margin-right:auto;\n",
"}\n",
"h1 {\n",
"font-family: 'Alegreya Sans', sans-serif;\n",
"}\n",
"h2 {\n",
"font-family: 'Fenix', serif;\n",
"}\n",
"h3{\n",
"font-family: 'Fenix', serif;\n",
"margin-top:12px;\n",
"margin-bottom: 3px;\n",
"}\n",
"h4{\n",
"font-family: 'Fenix', serif;\n",
"}\n",
"h5 {\n",
"font-family: 'Alegreya Sans', sans-serif;\n",
"}\t\n",
"div.text_cell_render{\n",
"font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
"line-height: 135%;\n",
"font-size: 120%;\n",
"width:600px;\n",
"margin-left:auto;\n",
"margin-right:auto;\n",
"}\n",
".CodeMirror{\n",
"font-family: \"Source Code Pro\";\n",
"font-size: 90%;\n",
"}\n",
"/* .prompt{\n",
"display: None;\n",
"}*/\n",
".text_cell_render h1 {\n",
"font-weight: 200;\n",
"font-size: 50pt;\n",
"line-height: 100%;\n",
"color:#CD2305;\n",
"margin-bottom: 0.5em;\n",
"margin-top: 0.5em;\n",
"display: block;\n",
"}\t\n",
".text_cell_render h5 {\n",
"font-weight: 300;\n",
"font-size: 16pt;\n",
"color: #CD2305;\n",
"font-style: italic;\n",
"margin-bottom: .5em;\n",
"margin-top: 0.5em;\n",
"display: block;\n",
"}\n",
".warning{\n",
"color: rgb( 240, 20, 20 )\n",
"}\n",
"</style>\n",
"\n",
"<script>\n",
"/* Equations */\n",
"\n",
"MathJax.Hub.Config({\n",
"TeX: {\n",
"extensions: [\"AMSmath.js\"]\n",
"},\n",
"tex2jax: {\n",
"inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
"displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
"},\n",
"displayAlign: 'center', // Change this to 'center' to center equations.\n",
"\"HTML-CSS\": {\n",
"styles: {'.MathJax_Display': {\"margin\": 4}}\n",
"}\n",
"});\n",
"</script>\n",
"\n",
"\n"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = open('./nb_style.css', 'r').read()\n",
" return HTML(styles)\n",
"css_styling()"
]
}
],
"metadata": {
"celltoolbar": "Raw Cell Format",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.1"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}