ESCUELA POLITÉCNICA NACIONAL - Repositorio...
175
ESCUELA POLITÉCNICA NACIONAL ESCUELA DE INGENIERÍA DESARROLLO DE SOFTWARE PARA EL ANÁLISIS BIDIMENSIONAL DE PLACAS MEDIANTE EL MÉTODO DE LOS ELEMENTOS FINITOS PROYECTO PREVIO A LA OBTENCIÓN DEL TÍTULO DE INGENIERO MECÁNICO JORGE EDUARDO ROMERO ROMERO DIRECTOR: Dr. VÍCTOR H. GUERRERO Quito, mayo 2006
Transcript of ESCUELA POLITÉCNICA NACIONAL - Repositorio...
ESCUELA POLITÉCNICA NACIONAL
ESCUELA DE INGENIERÍA
DESARROLLO DE SOFTWARE PARA EL ANÁLISIS
BIDIMENSIONAL DE PLACAS MEDIANTE EL MÉTODO DE LOS
ELEMENTOS FINITOS
PROYECTO PREVIO A LA OBTENCIÓN DEL TÍTULO DE INGENI ERO
MECÁNICO
JORGE EDUARDO ROMERO ROMERO
DIRECTOR: Dr. VÍCTOR H. GUERRERO
Quito, mayo 2006
ii
DECLARACIÓN
Yo, Jorge Eduardo Romero Romero, declaro bajo juramento que el trabajo aquí
descrito es de mi autoría, que no ha sido previamente presentado para ningún
grado o calificación profesional, y, que he consultado las referencias bibliográficas
que se incluyen en este documento.
A través de la presente declaración cedo mis derechos de propiedad intelectual
correspondientes a este trabajo a la Escuela Politécnica Nacional, según lo
establecido por la ley de Propiedad Intelectual, por su Reglamento y por la
normatividad institucional vigente.
Jorge Eduardo Romero Romero
iii
CERTIFICACIÓN
Certifico que el presente trabajo fue desarrollado por Jorge Eduardo Romero
Romero, bajo mi supervisión.
Dr. Víctor H. Guerrero
DIRECTOR DEL PROYECTO
iv
AGRADECIMIENTOS
Quiero agradecer al Dr. Víctor Guerrero por ser la persona que me impulsó y guió
durante el desarrollo de este trabajo, comenzando inicialmente como un proyecto
final para el curso de Elementos Finitos dictado por él, y luego como proyecto de
tesis.
v
DEDICATORIA
A mis padres, por su constante apoyo, supieron darme el aliento que necesitaba
para concluir este trabajo.
vi
CONTENIDO
CAPÍTULO 1: INTRODUCCIÓN........................... ................................................. 1
1.1 INTRODUCCIÓN AL MÉTODO DE LOS ELEMENTOS FINITOS ........................................... 1
1.2 COMPARACIÓN DEL MEF CON OTROS MÉTODOS ............................................................ 2
1.3 FUNCIONAMIENTO DEL MEF ................................................................................................ 3
1.4 LIMITACIONES........................................................................................................................ 7
CAPÍTULO 2: ANÁLISIS DE PLACAS ..................... ............................................ 8
2.1 FLEXIÓN DE PLACAS ............................................................................................................. 8
2.1.1 TEORÍA CLÁSICA DE PLACAS DELGADAS.................................................................. 9
2.1.2 DESARROLLO DEL ELEMENTO RECTANGULAR PARA FLEXIÓN EN PLACAS ..... 10
2.1.3 ELEMENTO TRIANGULAR DKT ................................................................................... 20
2.2 PLACAS PLANAS .................................................................................................................. 26
2.2.1 ESFUERZO PLANO....................................................................................................... 27
2.2.2 ELEMENTO PLANO RECTANGULAR .......................................................................... 27
2.2.3 ELEMENTO TRIANGULAR CST ................................................................................... 32
2.3 CONVERGENCIA .................................................................................................................. 36
2.4 CRITERIO GENERAL DE CEDENCIA EN TRES DIMESIONES........................................... 36
2.4.1 POSTULADO 1: EXISTE UNA FUNCIÓN GENERAL DE CEDENCIA ......................... 37
2.4.2 POSTULADO 2: EL MATERIAL ES ISOTRÓPICO ....................................................... 37
2.4.3 POSTULADO 3: LA CEDENCIA ES INDEPENDIENTE DEL ESFUERZO
HIDROSTÁTICO...................................................................................................................... 38
2.4.4 POSTULADO 4: COMPORTAMIENTO IDÉNTICO A TENSIÓN Y COMPRESIÓN...... 38
2.5 HIPÓTESIS DE FALLA PARA MATERIALES DÚCTILES ..................................................... 38
2.5.1 HIPÓTESIS DEL ESFUERZO CORTANTE MÁXIMO (CRITERIO DE TRESCA)......... 39
2.5.2 HIPÓTESIS DE LA ENERGÍA DE DISTORSIÓN (CRITERIO DE VON MISES) .......... 39
CAPÍTULO 3: GENERACIÓN DE MALLAS NO ESTRUCTURADAS .. .............. 41
3.1 INTRODUCCIÓN ................................................................................................................... 41
3.2 DEFINICIÓN DE SUPERFICIES............................................................................................ 42
3.3 MÉTODOS PARA GENERAR MALLAS NO ESTRUCTURADAS ......................................... 42
3.3.1 MÉTODO DEL AVANCE FRONTAL .............................................................................. 43
3.3.2 MÉTODO DE TRIANGULACIÓN DE DELAUNAY......................................................... 45
3.4 POST-PROCESO .................................................................................................................. 47
3.4.1 ALISADO ........................................................................................................................ 47
3.4.2 LIMPIEZA ....................................................................................................................... 49
3.4.3 REFINAMIENTO............................................................................................................. 51
vii
CAPÍTULO 4: PROGRAMACIÓN........................... ............................................. 53
4.1 DESCRIPCIÓN DEL PROGRAMA......................................................................................... 53
4.1.1 SELECCIÓN DEL MÉTODO DE MALLADO.................................................................. 53
4.1.2 SELECCIÓN DE LOS ELEMENTOS ............................................................................. 54
4.1.3 OPTIMIZACIÓN DE LA MEMORIA NECESARIA .......................................................... 54
4.2 LENGUAJE ............................................................................................................................ 55
4.3 PROGRAMACIÓN EN MATLAB ............................................................................................ 56
4.3.1 SUBPROGRAMAS CON INTERFAZ GRÁFICA ............................................................ 57
4.3.2 SUBPROGRAMAS SIN INTERFAZ GRÁFICA .............................................................. 63
4.4 PROGRAMACIÓN EN JAVA.................................................................................................. 70
4.4.1 SUBPROGRAMAS CON INTERFAZ GRÁFICA ............................................................ 71
4.4.2 SUBPROGRAMAS SIN INTERFAZ GRÁFICA .............................................................. 79
4.4.3 EJECUCIÓN DEL PROGRAMA..................................................................................... 91
4.4.4 LIBRERÍAS EXTERNAS USADAS POR EL PROGRAMA ............................................ 91
CAPÍTULO 5: PRUEBAS DE VERIFICACIÓN................ .................................... 96
5.1 EJEMPLOS DE PRUEBA PARA ELEMENTOS TIPO PLACA............................................... 96
5.1.1 EJEMPLO 1: PLACA RECTANGULAR.......................................................................... 96
5.1.2 EJEMPLO 2: DISCO....................................................................................................... 98
5.1.3 EJEMPLO 3: VIGA ....................................................................................................... 100
5.2 EJEMPLOS DE PRUEBA PARA ELEMENTOS TIPO 2D .................................................... 102
5.2.1 EJEMPLO 1: HERRADURA ......................................................................................... 102
5.2.2 EJEMPLO 2: LÁMINA DELGADA ................................................................................ 103
5.2.3 EJEMPLO 3: PLACA .................................................................................................... 105
5.3 COMPARACIÓN DE LOS TRES PROGRAMAS ................................................................. 107
5.3.1 COMPARACIÓN DEL MALLADO ................................................................................ 107
5.3.2 COMPARACIÓN DE LA SOLUCIÓN ........................................................................... 108
5.3.3 COMPARACIÓN DE LA CONVERGENCIA................................................................. 110
CAPÍTULO 6: CONCLUSIONES Y TRABAJO FUTURO.......... ........................ 112
6.1 CONCLUSIONES................................................................................................................. 112
6.2 TRABAJO FUTURO............................................................................................................. 113
BIBLIOGRAFÍA ....................................... .......................................................... 115
ANEXOS ............................................................................................................ 118
viii
LISTA DE FIGURAS
Figura 2.1 Placa delgada ..................................................................................... 10
Figura 2.2 Corte diferencial a una placa............................................................... 11
Figura 2.3 Esfuerzos actuantes en la placa ......................................................... 13
Figura 2.4 Cargas actuantes en la placa.............................................................. 14
Figura 2.5 Elemento rectangular de 12 grados de libertad................................... 15
Figura 2.6 Características de convergencia del elemento rectangular [22].......... 19
Figura 2.7 Elemento DKT..................................................................................... 21
Figura 2.8 Geometría de un elemento triangular.................................................. 23
Figura 2.9 Condiciones de esfuerzo plano [11] .................................................... 27
Figura 2.10 Elemento plano rectangular .............................................................. 28
Figura 2.11 Notación para un elemento en estuerzo plano.................................. 32
Figura 2.12 Hipótesis de falla ............................................................................... 40
Figura 3.1 (a) Frente inicial, (b) frente actual, (c) frente actualizado con el nuevo
elemento............................................................................................................... 43
Figura 3.2 Triangulación de Delaunay (2D).......................................................... 45
Figura 3.3 Malla (a) antes del alisado Laplaciano, (b) después ........................... 48
Figura 3.4 Ejemplo de un refinamiento de Delaunay............................................ 52
Figura 3.5 Ejemplo de un refinamiento local triangular ........................................ 52
Figura 4.1 Ventana principal del programa .......................................................... 57
Figura 4.2 Ventana Acerca................................................................................... 58
Figura 4.3 Ventana Análisis.................................................................................. 59
Figura 4.4 Ventana Apoyos.................................................................................. 59
Figura 4.5 Ventana Cargas .................................................................................. 60
Figura 4.6 Ventana Geometría ............................................................................. 61
Figura 4.7 Ventana Malla ..................................................................................... 62
Figura 4.8 Ventana Opciones............................................................................... 62
Figura 4.9 Ventana de propiedades ..................................................................... 63
Figura 4.10 Ventana Tipo de Elemento................................................................ 63
Figura 4.11 Ventana Placa................................................................................... 71
Figura 4.12 Ventana Acerca................................................................................. 72
Figura 4.13 Ventana Análisis................................................................................ 73
ix
Figura 4.14 Ventana Apoyos................................................................................ 73
Figura 4.15 Ventana Barras de herramientas....................................................... 74
Figura 4.16 Ventana Cambios.............................................................................. 74
Figura 4.17 Ventana Cargas ................................................................................ 75
Figura 4.18 Ventana Ayuda.................................................................................. 75
Figura 4.19 Ventana GraficoPlano ....................................................................... 76
Figura 4.20 Ventana InfoNodo ............................................................................. 76
Figura 4.21 Ventana Librería de materiales ......................................................... 77
Figura 4.22 Ventana Malla ................................................................................... 77
Figura 4.23 Ventana Opciones............................................................................. 78
Figura 4.24 Ventana Propiedades........................................................................ 79
Figura 4.25 Ventana Unidades............................................................................. 79
Figura 4.26 Numeración de nodos y elementos ................................................... 90
Figura 5.1 Placa empotrada ................................................................................. 96
Figura 5.2 Resultados del ejemplo 1. (a) Matlab, (b) Java, (c) Algor.................... 97
Figura 5.3 Disco ................................................................................................... 98
Figura 5.4 Resultados del ejemplo 2. (a) Matlab, (b) Java, (c) Algor.................... 99
Figura 5.5 Viga ................................................................................................... 100
Figura 5.6 Resultados del ejemplo 3. (a) Matlab, (b) Java, (c) Algor.................. 101
Figura 5.7 Herradura .......................................................................................... 102
Figura 5.8 Resultados del ejemplo 1. (a) Matlab, (b) Java, (c) Algor.................. 103
Figura 5.9 Lámina delgada................................................................................. 104
Figura 5.10 Resultados del ejemplo 2. (a) Matlab, (b) Java, (c) Algor................ 105
Figura 5.11 Placa ............................................................................................... 105
Figura 5.12 Resultados del ejemplo 3. (a) Matlab, (b) Java, (c) Algor................ 106
Figura 5.13 Tiempos de mallado para los ejemplos de (a) sección 5.1.1 y (b)
sección 5.2.1 ...................................................................................................... 108
Figura 5.14 Tiempos de cálculo para los ejemplos de (a) sección 5.1.1 y (b)
sección 5.2.1 ...................................................................................................... 109
x
LISTA DE TABLAS
Tabla 4.1 Memoria ocupada por la matriz de rigidez del sistema......................... 54
Tabla 5.1 Resultados del ejemplo 1 ..................................................................... 97
Tabla 5.2 Resultados del ejemplo 2 ..................................................................... 99
Tabla 5.3 Resultados del ejemplo 3 ................................................................... 101
Tabla 5.4 Resultados del ejemplo 1 ................................................................... 103
Tabla 5.5 Resultados del ejemplo 2 ................................................................... 104
Tabla 5.6 Resultados del ejemplo 3 ................................................................... 106
Tabla 5.7 Convergencia del esfuerzo de von Mises (ejemplo sección 5.1.1)..... 110
Tabla 5.8 Convergencia del esfuerzo de von Mises (ejemplo sección 5.2.1)..... 111
xi
RESUMEN
En este trabajo se han desarrollado dos programas para el análisis bidimensional
de placas, uno realizado en Matlab 7.0 y el otro en Java empleando el J2SE
Development Kit 5.0, Java 3D 1.4.0 y NetBeans 5.0. Estos programas emplean el
método de los elementos finitos para realizar el mencionado análisis.
La versión del programa en Java, al ser la más completa, merece especial
atención. Posee una interfaz gráfica de uso intuitivo. Tanto la geometría como los
apoyos y cargas se ingresan de forma gráfica, el mallado se realiza de manera
automática y el análisis, que emplea dos librerías de libre distribución para el
manejo de matrices, puede realizarse en poco tiempo. Además los resultados se
pueden visualizar en forma gráfica y numérica.
Para poder implementar el método de los elementos finitos como técnica de
análisis computacional se ha investigado los elementos más comunes utilizados
para discretizar un sistema plano continuo, resultando en cuatro elementos
principales: el elemento rectangular, el elemento triangular DKT, el elemento
plano rectangular y el elemento triangular CST. Adicionalmente se estudió las dos
técnicas predominantes de generación de mallas no estructuradas: el método del
avance frontal y la triangulación de Delaunay. Estos estudios sirvieron como punto
de partida para el desarrollo de los algoritmos que hacen funcionar el programa.
Finalmente, se ha estudiado una serie de ejemplos de prueba para demostrar la
validez de los programas desarrollados al comparar los resultados con los
obtenidos por Algor v19. Dichos resultados demuestran que el software funciona
perfectamente siempre y cuando los problemas que se analicen puedan ser
representados por la teoría de placas. El estudio concluye con la comparación de
los tiempos de mallado y de solución, y la comparación de la convergencia de los
resultados si se aumenta la densidad de la malla en los tres programas.
xii
PRESENTACIÓN
Este trabajo tiene como objetivo el desarrollar un programa para el análisis
bidimensional de placas empleando el método de los elementos finitos. En
nuestro país, este tema ha sido objeto de un estudio muy superficial en algunos
de sus tópicos y casi nulo en otros, enfocándose principalmente a mostrar el uso
de software comercial mediante la solución de problemas relativamente simples.
Para poder abordar el desarrollo de tal programa es necesario comprender de qué
se trata el método de los elementos finitos, cómo discretizar un sistema continuo y
con qué tipo y forma de elementos se realiza tal discretización. Además, es
necesario determinar qué lenguaje de programación es el idóneo. Estos
requerimientos se exponen en sucesivos capítulos.
En el capítulo 1 se da una breve introducción al método de los elementos finitos y
el porqué de su reciente preferencia frente a otros métodos en el análisis de
problemas de ingeniería. El capítulo concluye con la enumeración de los pasos
necesarios para llegar a la solución de un problema por el método de los
elementos finitos y una descripción de sus limitaciones.
En el capítulo 2 se deducen las matrices de rigidez y de esfuerzo para elementos
tipo placa y tipo 2D. Adicionalmente, en anexos, se muestran las matrices ya
calculadas para estos elementos. También se da una breve explicación sobre las
dos teorías de falla de materiales dúctiles más importantes actualmente. Esta
información es usada por el programa para calcular desplazamientos, reacciones,
deformaciones y esfuerzos.
En el capítulo 3 se describe el proceso de mallado, haciendo hincapié en el
mallado no estructurado y en los métodos más populares para realizar la malla.
Además se describen los procesos posteriores al mallado como alisado, limpieza
y refinamiento.
xiii
En el capítulo 4 se describe el programa desarrollado. Primero se lo hace en
Matlab porque es un lenguaje fácil de aprender y especialmente apto para trabajar
con matrices. Luego se presenta la versión en Java, con la que se consigue que
el programa sea independiente de la plataforma de trabajo que se escoja. Para
ambas versiones del programa se presentan los algoritmos más significativos y
una breve explicación de todas las funciones (Matlab) y clases (Java). En anexos
se incluyen los diagramas de flujo de las subrutinas más importantes.
En el capítulo 5 se resuelven varios problemas empleando los programas
desarrollados en Matlab y en Java. Los resultados obtenidos con estos programas
se comparan con aquellos obtenidos por Algor v19. También se comparan los
tiempos que emplean para realizar el mallado no estructurado y para calcular la
solución.
En el capítulo 6 se exponen las conclusiones de este trabajo y se dan varias ideas
de posibles modificaciones y ampliaciones futuras al programa, que podrían ser
desarrolladas principalmente en nuevos proyectos de titulación.
1
CAPÍTULO 1.
INTRODUCCIÓN
1.1 INTRODUCCIÓN AL MÉTODO DE LOS ELEMENTOS FINITOS
El método de los elementos finitos MEF (o FEM por sus siglas en inglés) es una
técnica de análisis numérico que permite obtener soluciones aproximadas a una
extensa variedad de problemas de ingeniería. Estos problemas caen en el grupo
de los llamados problemas con valor en la frontera, que no son sino problemas
matemáticos en los cuales una o más variables dependientes deben satisfacer
una ecuación diferencial dentro del dominio de variables independientes y
satisfacer condiciones específicas en la frontera del dominio. A menudo a los
problemas de valor en la frontera también se los conoce como problemas de
campo.
El MEF se ha convertido rápidamente en la técnica de solución dominante en casi
todos los campos de la ingeniería y la ciencia aplicada. El alto grado de
aceptación del método se puede atribuir a los siguientes factores principales:
1. El MEF implantado computacionalmente permite modelar de manera casi
exacta la forma real del dominio de un problema, las restricciones y las
condiciones de carga.
2. La amplia disponibilidad y el incremento en la potencia de las
computadoras personales han facilitado el uso de métodos numéricos en el
análisis y diseño en ingeniería y han puesto el MEF al alcance de la mayor
parte de los ingenieros.
3. Una vez que un programa para el análisis estructural se escribe, este
puede ser fácilmente modificado para resolver diversos problemas
ingenieriles, ya que las operaciones matriciales básicas involucradas son
las mismas.
2
4. El MEF es versátil y flexible. El dominio puede tener forma, restricciones y
condiciones de carga arbitrarias. Aún en este caso, el MEF brinda buenas
aproximaciones a la solución exacta de un problema.
Las aplicaciones del MEF incluyen problemas estáticos lineales y no lineales,
dinámicos y de estabilidad, mecánica de fluidos, transferencia de calor,
biomecánica, geomecánica, en los campos estructural, mecánico, naval, minería e
ingeniería aeroespacial. Este amplio rango de aplicaciones del MEF se debe en
parte a su formulación basada en los principios variacionales.
1.2 COMPARACIÓN DEL MEF CON OTROS MÉTODOS
Cada vez más situaciones de ingeniería demandan la obtención de soluciones
numéricas aproximadas a los problemas que plantean, dada la dificultad para
obtener la solución exacta. Por ejemplo, hablemos del caso en que se quiere
encontrar la capacidad de carga de una placa que tiene varios orificios de formas
complejas. Al escribir las ecuaciones que gobiernan el sistema y las condiciones
de frontera para este problema se nota que es muy difícil encontrar una solución
analítica simple. La dificultad está en que la geometría es irregular y puede en
ciertos casos calificarse de “arbitraria”. Soluciones analíticas a problemas de este
tipo rara vez existen. Sin embargo, esta es la clase de problemas que los
ingenieros están llamados a resolver.
Varias alternativas han sido propuestas para superar este dilema. Una posibilidad
es hacer suposiciones que simplifiquen y reduzcan el problema a uno que pueda
ser abordado de manera más fácil. Algunas veces este procedimiento funciona,
pero a menudo conduce a respuestas erróneas en mayor o menor medida. Ahora
que las computadoras están ampliamente disponibles, una alternativa viable es
mantener la complejidad del problema y encontrar una solución numérica
aproximada.
Diversos métodos de análisis numérico aproximado han evolucionado a través de
los años. Un método comúnmente utilizado es el método de diferencias finitas. El
3
modelo de diferencias finitas da una aproximación de las ecuaciones gobernantes
de un problema. Este modelo se mejora mientras más puntos se empleen para
formar la malla que describe el dominio del problema. Con las técnicas de
diferencias finitas se pueden tratar algunos problemas complejos. Sin embargo,
cuando se trata con geometrías irregulares o condiciones de frontera inusuales, el
método de diferencias finitas es difícil de aplicar.
En contraste con el método de diferencias finitas, en el cual se modela la región
de la solución como un arreglo de puntos dentro de una malla, el método de los
elementos finitos modela la región de la solución como el conjunto de muchas
pequeñas subregiones o elementos interconectados. Un modelo de elementos
finitos de un problema da una aproximación del problema físico a resolver y no
solo de las ecuaciones gobernantes. La idea básica del método de elementos
finitos es que una región de solución puede ser modelada analíticamente o ser
aproximada reemplazándola con un ensamble de elementos discretos. Como
estos elementos pueden ser juntados de diversas maneras, pueden ser usados
para representar formas sumamente complejas.
1.3 FUNCIONAMIENTO DEL MEF
En un problema continuo de cualquier dimensión, la variable de campo posee un
infinito número de valores que corresponden a todos y cada uno de los puntos
que se encuentran en la región de solución. Consecuentemente, el problema tiene
un número infinito de incógnitas. El procedimiento de discretización por elementos
finitos reduce el problema a uno con un número finito de incógnitas al dividir la
región de solución en elementos y expresar la variable de campo desconocida en
términos de funciones de aproximación asumidas para cada elemento. Las
funciones de aproximación (llamadas también funciones de interpolación) están
definidas por los valores de las variables de campo en puntos específicos
llamados nodos o puntos nodales. Los nodos usualmente están en el contorno del
elemento, conectando los elementos adyacentes. Un elemento puede tener
también nodos interiores. Los valores de la variable de campo y las funciones de
4
interpolación para los elementos definen completamente el comportamiento de la
variable de campo dentro de los elementos.
Para la representación de un problema por elementos finitos los valores nodales
de la variable de campo se convierten en las incógnitas. Una vez que se
encuentran los valores de estas incógnitas, las funciones de interpolación
permiten calcular el valor de la variable de campo en cualquier punto del
ensamble de elementos.
Claramente, la naturaleza de la solución y el grado de aproximación dependen no
solamente del tamaño y el número de elementos usados, sino también de las
funciones de interpolación seleccionadas. Como se podría esperar, no se puede
escoger funciones arbitrariamente, porque ciertas condiciones de compatibilidad
deber ser satisfechas. A menudo las funciones se escogen de manera que la
variable de campo o sus derivadas sean continuas a través del contorno de
elementos adyacentes.
También hay que mencionar una importante característica del método de los
elementos finitos que lo aparta de otros métodos numéricos. Esta característica
es la habilidad de formular soluciones para elementos individuales antes de
ponerlos juntos para representar la totalidad del problema. Esto significa, por
ejemplo, que si se trata con un problema de análisis de esfuerzos, se hallan las
ecuaciones que describen el comportamiento de cada elemento individual y
entonces se ensamblan los elementos para describir el comportamiento de toda la
estructura. En esencia, un problema complejo se reduce a considerar una serie de
muchos problemas simplificados.
Sin tener en cuenta las aproximaciones empleadas para hallar las propiedades
del elemento, la solución a un problema continuo por el método de los elementos
finitos siempre sigue un proceso ordenado paso a paso. Estos pasos se listan en
términos generales a continuación:
5
1. Discretizar el sistema continuo. El primer paso es dividir el continuo o
región de solución en elementos. Varias formas de elementos pueden ser
usadas, y diferentes formas de elementos pueden ser empleadas en la
misma región de solución. Por ejemplo, cuando se analiza una estructura
elástica que tiene diferentes tipos de componentes tales como placas y
vigas, no solo es deseable, sino muy necesario usar diferentes elementos
en la misma solución.
2. Seleccionar las funciones de interpolación. El siguiente paso es asignar los
nodos a cada elemento y entonces escoger la función de interpolación para
representar la variación de la variable de campo sobre el elemento. La
variable de campo puede ser un escalar, un vector, o un tensor de orden
determinado. Frecuentemente se selecciona polinomios como funciones de
interpolación para la variable de campo debido a que son fáciles de integrar
y diferenciar. El grado del polinomio escogido depende del número de
nodos asignados al elemento, la naturaleza y número de incógnitas en
cada nodo, y ciertos requerimientos de continuidad impuestos en los nodos
y a lo largo de la frontera del elemento. La magnitud de la variable de
campo así como la magnitud de sus derivadas pueden ser las incógnitas
en los nodos.
3. Hallar las propiedades del elemento. Una vez que ha sido establecido el
modelo por elementos finitos, se debe determinar la ecuación matricial que
exprese las propiedades de los elementos individuales. Para esta tarea se
puede usar uno de los procedimientos siguientes: el enfoque directo, el
enfoque variacional, o el enfoque de residuos ponderados.
4. Ensamblar las propiedades de los elementos para obtener el sistema de
ecuaciones. Para encontrar las propiedades de todo el sistema modelado
por la red de elementos se debe “ensamblar” todas las propiedades de los
elementos. En otras palabras, se combinan las ecuaciones matriciales que
expresan el comportamiento de los elementos y forman la ecuación
matricial que expresa el comportamiento del sistema total. La ecuación
6
matricial del sistema incluye a todos los nodos del sistema. La base para el
procedimiento de ensamblaje parte de que en un nodo en donde los
elementos se interconectan, el valor de la variable de campo es el mismo
para cada elemento que comparte ese nodo. Una característica única del
método de los elementos finitos es que el sistema de ecuaciones que
describe el problema se genera ensamblando las ecuaciones de los
elementos individuales. En contraste, en el método de diferencias finitas el
sistema de ecuaciones se genera escribiendo las ecuaciones nodales.
5. Imponer las condiciones de frontera. Antes de que el sistema de
ecuaciones esté listo para ser resuelto, se lo debe modificar para tomar en
cuenta las condiciones de frontera del problema. En esta etapa se imponen
los valores nodales conocidos de las variables dependientes.
6. Resolver el sistema de ecuaciones. El proceso de ensamblado resulta en
un sistema de ecuaciones simultáneas que se resuelve para obtener los
valores nodales desconocidos del problema. Si el problema describe un
estado estable o en equilibrio, entonces se debe resolver un sistema de
ecuaciones algebraicas lineales o no lineales. Si el problema es inestable,
las incógnitas nodales son una función del tiempo, y se debe resolver un
sistema de ecuaciones diferenciales ordinarias lineales o no lineales.
7. Calcular variables adicionales. Muchas veces se usa la solución del
sistema de ecuaciones para calcular otros parámetros importantes. Por
ejemplo, en un problema estructural las incógnitas nodales son las
componentes del desplazamiento. A partir de estos desplazamientos se
calculan los esfuerzos y deformaciones de los elementos.
7
1.4 LIMITACIONES
Las grandes ventajas del método de los elementos finitos son la habilidad de
manejar geometrías arbitrarias y materiales no homogéneos. Estas dos
características permiten que se pueda modelar formas complejas constituidas por
regiones de varios materiales.
Sin embargo, el método está basado en la técnica de dividir la región o dominio
del problema en un número finito de elementos y entonces encontrar la mejor
solución posible que sea continua dentro de los elementos, pero que puede no ser
continua en la frontera entre elementos. La cantidad de “saltos” en los resultados
se usa a menudo para evaluar la precisión de la solución. Esta precisión depende
del número y tamaño de los elementos y del grado de la función de interpolación
usada dentro del elemento [14].
8
CAPÍTULO 2.
ANÁLISIS DE PLACAS
Una placa es una forma particular de un sólido tridimensional. Es un elemento
estructural plano con un espesor mucho menor que las otras dimensiones (largo y
ancho). Es uno de los más importantes elementos estructurales y es usado para
modelar y analizar estructuras tales como puertas, partes de automóviles, etc. [23]
Las placas se pueden clasificar en placas delgadas y gruesas. En forma rigurosa,
para que una placa sea considerada delgada, la razón entre el espesor y la
longitud más corta de la placa debe ser menor a 1/10. Las placas tratadas en este
trabajo están hechas de un material que es homogéneo e isotrópico. Si las
propiedades del material son las mismas en todas partes, se dice que el material
es homogéneo, mientras que si las propiedades del material son idénticas en
todas las direcciones, se dice que es isotrópico.
En este capítulo se presenta la formulación de cuatro elementos típicamente
empleados en el análisis de placas. Estos elementos se escogieron porque su
formulación es fácil de entender y son los elementos más simples para el análisis.
Por ejemplo, el número de grados de libertad en cada nodo es 3 para flexión en
placas y 2 para placas planas. Asimismo, el elemento rectangular está formado
por 4 nodos y el elemento triangular por 3. Otros tipos de elementos comúnmente
usados involucran elementos con un número mayor de grados de libertad por
nodo y/o un mayor número de nodos.
2.1 FLEXIÓN DE PLACAS
Una placa en flexión puede ser considerada la extensión a dos dimensiones de
una viga en flexión simple. Tanto placas como vigas soportan cargas
transversales o perpendiculares a su plano geométrico. Una viga tiene un único
momento de flexión, mientras que una placa resiste flexión sobre dos ejes y tiene
9
un momento de torsión. Se considerará la teoría clásica de placas delgadas o
teoría de placas de Kirchhoff [26].
2.1.1 TEORÍA CLÁSICA DE PLACAS DELGADAS
Un análisis matemático exacto de esfuerzos en una placa delgada, sujeta a
cargas actuantes normales a su superficie, requiere la solución de ecuaciones
diferenciales de elasticidad en tres dimensiones. En la mayoría de los casos, sin
embargo, tal enfoque implicaría dificultades matemáticas difíciles de superar.
Afortunadamente, para la mayoría de las aplicaciones técnicas, la teoría clásica
de Kirchhoff para placas delgadas proporciona resultados suficientemente exactos
sin la necesidad de llevar a cabo un análisis de esfuerzos tridimensional. Las
simplificaciones usadas en la derivación de la ecuación de placas están basadas
en las siguientes suposiciones:
1. El material es homogéneo, isotrópico y linealmente elástico, esto es, sigue
la ley de Hooke.
2. La placa es inicialmente plana.
3. La superficie media de la placa permanece sin deformarse durante la
flexión.
4. El espesor de la placa, h, es pequeño comparado con las otras
dimensiones. Esto es, la dimensión lateral más corta de la placa es al
menos 10 veces más larga que el espesor de esta.
5. Las deflexiones transversales w(x, y) son pequeñas comparadas con el
espesor de la placa. Una deflexión máxima de 1/10 del espesor se
considera el límite.
6. Las pendientes de la superficie media son pequeñas comparadas a la
unidad.
7. Las secciones tomadas normales a la superficie media antes de la
deformación permanecen planas y normales a la superficie media
flexionada. Consecuentemente, las deformaciones cortantes se
desprecian.
10
8. El esfuerzo normal σz en la dirección transversal a la superficie de la placa
puede ser despreciado.
Con la ayuda de estas suposiciones, el problema original en tres dimensiones se
reduce a un problema de placas en dos dimensiones [22].
2.1.2 DESARROLLO DEL ELEMENTO RECTANGULAR PARA FLEX IÓN EN
PLACAS
Se han desarrollado un gran número de formulaciones de elementos para placas
en flexión. En este trabajo se presenta la derivación de la matriz de rigidez para
uno de los elementos más comunes para flexión en placas [26].
2.1.2.1 Comportamiento básico de la geometría y la deformación
Considere la placa delgada mostrada en la figura 2.1. Esta placa está situada en
el plano xy, tiene un espesor h medido en la dirección z y está sometida a una
carga q.
Figura 2.1 Placa delgada
La superficie de la placa está entre z = ± h/2, y la superficie media está en z = 0.
La geometría básica de la placa es como sigue:
11
1. El espesor de la placa es mucho menor que sus dimensiones b y c (esto
es, h << b o c). Si h es mayor que aproximadamente 1/10 la extensión de la
placa, la deformación de corte transversal debe ser tomada en cuenta y se
dice que la placa es gruesa.
2. La deflexión w es mucho menor que el espesor h (esto es, w/h << 1).
2.1.2.2 Suposiciones de Kirchhoff
Considere el corte diferencial de la placa por planos perpendiculares al eje x,
como se muestra en la figura 2.2:
Figura 2.2 Corte diferencial a una placa
La carga q causa que la placa se deforme lateralmente o hacia arriba en la
dirección z. La deflexión w del punto P se asume que es una función solo de x e y,
esto es w = w(x, y) y la placa no se estira en la dirección z. La línea ab dibujada
perpendicular a la superficie de la placa antes de ser cargada permanece
perpendicular a la superficie después de aplicada la carga. Estas condiciones son
consistentes con las suposiciones de Kirchhoff:
1. Normales permanecen normales. Esto implica que las deformaciones
cortantes transversales γyz = 0 y γxz = 0. Sin embargo γxy no es igual a cero.
Ángulos rectos en el plano de la placa no pueden permanecer rectos
después de la carga. La placa debe torcerse en el plano.
2. Cambios en el espesor pueden ser despreciados. Esto significa que εz = 0.
12
3. El esfuerzo normal σz no tiene efecto en las deformaciones εx y εy en las
ecuaciones esfuerzo-deformación y es considerado despreciable.
4. Fuerzas dentro del plano son despreciables, y la resistencia de esfuerzo
plano puede ser impuesta más tarde. Por lo tanto, se asume que los
desplazamientos en el plano de la superficie media, h = 0, en las
direcciones x e y son cero; u(x, y, 0) = 0 y v(x, y, 0) = 0.
Basado en las suposiciones de Kirchhoff, en cualquier punto P el desplazamiento
en la dirección x debido a una pequeña rotación α es:
∂∂−=−=
x
wzzu α
En el mismo punto, el desplazamiento en la dirección y es:
∂∂−=−=
y
wzzv α
Las curvaturas de la placa están entonces dadas como una tasa de cambio de los
desplazamientos angulares de las normales y son definidas como:
2
2
x
wkx ∂
∂−= 2
2
y
wky ∂
∂−= yx
wkxy ∂∂
∂−=22
Usando las definiciones para los esfuerzos dentro del plano, junto con las
relaciones para la curvatura, las ecuaciones para la deformación/desplazamiento
en el plano son:
2
2
x
wzx ∂
∂−=ε 2
2
y
wzy ∂
∂−=ε yx
wzxy ∂∂
∂−=2
2γ
La primera de las ecuaciones anteriores es usada en la teoría de vigas. Las dos
ecuaciones restantes son nuevas para la teoría de placas.
2.1.2.3 Relaciones esfuerzo-deformación
Basadas en la tercera suposición de Kirchhoff, las ecuaciones de esfuerzo plano
que relacionan los esfuerzos en el plano con las deformaciones en el plano para
un material isotrópico son:
( )yxx
E νεεν
σ +−
=21
( )xyy
E νεεν
σ +−
=21
xyxy Gγτ =
(2.1a)
(2.1b)
(2.2)
(2.3)
(2.4)
13
Los esfuerzos normales en el plano y esfuerzos de corte se muestran actuando
en los bordes de la placa mostrada en la figura 2.3:
Figura 2.3 Esfuerzos actuantes en la placa
Similar a la variación de esfuerzos en una viga, los esfuerzos varían linealmente
en la dirección z desde la superficie media de la placa. Los esfuerzos cortantes
transversales τyz y τxz están también presentes, aún cuando la deformación
cortante transversal es despreciable. Estos esfuerzos varían cuadráticamente a lo
largo del espesor de la placa.
Los momentos de flexión actuando a lo largo del borde de la placa pueden ser
relacionados con los esfuerzos por:
∫−
=2
2
h
h
xx dzzM σ ∫−
=2
2
h
h
yy dzzM σ ∫−
=2
2
h
h
xyxy dzzM τ
Expresando los esfuerzos en términos de deformaciones se obtiene:
( )∫−
+−
=2
221
h
h
yxx dzE
zM νεεν
( )∫−
+−
=2
221
h
h
xyy dzE
zM νεεν
∫−
=2
2
h
h
xyxy dzzGM γ
Usando las relaciones deformación-curvatura, la expresión para el momento es:
( )yxx kkDM ν+= ( )xyy kkDM ν+= ( )
xyxy kD
M2
1 ν−=
Donde ( )[ ]23 112 ν−= EhD es la rigidez a la flexión de la placa.
(2.5a)
(2.5b)
(2.6)
14
La mayor magnitud del esfuerzo normal de cada borde de la placa se localiza en
la parte superior e inferior, en z = ± h/2.
2
6
h
M xx =σ
Las ecuaciones de equilibrio para placas en flexión son importantes en la
selección de los campos de desplazamiento del elemento. Las ecuaciones
diferenciales gobernantes son:
0=+∂
∂+
∂∂
qy
Q
x
Q yx
0=−∂
∂+
∂∂
xxyx Q
y
M
x
M
0=−∂
∂+
∂∂
yxyy Q
x
M
y
M
Donde q es la carga distribuida transversal y Qx y Qy son las cargas lineales de
corte transversales como se muestra en la figura 2.4:
Figura 2.4 Cargas actuantes en la placa
Sustituyendo las expresiones momento-curvatura en las dos últimas ecuaciones
diferenciales anteriores, resolviendo para Qx y Qy, y sustituyendo los resultados
en la primera ecuación listada arriba, la ecuación diferencial parcial gobernante
para placas delgadas en flexión isotrópicas puede ser derivada como:
(2.7)
(2.8)
15
qy
w
yx
w
x
wD =
∂∂+
∂∂∂+
∂∂
4
4
22
4
4
4 2
Donde la solución para flexión de placas es una función del desplazamiento
transversal w. Si se desprecia la diferenciación con respecto a la dirección y, la
ecuación anterior se simplifica a la ecuación para una viga y la rigidez de flexión D
de la placa se reduce al EI de la viga cuando el efecto de Poisson es fijado en
cero.
2.1.2.4 Energía potencial de una placa
La energía potencia total de una placa está dada por:
( )∫ ++=V
xyxyyyxx dVU γτεσεσ2
1
La energía potencial puede ser expresada en términos de momentos y curvaturas
como:
( )∫ ++=A
xyxyyyxx dAkMkMkMU2
1
2.1.2.5 Derivación de la matriz de rigidez del elemento
En esta sección se derivan las ecuaciones necesarias para el elemento
rectangular de 12 grados de libertad mostrado en la figura 2.5.
Figura 2.5 Elemento rectangular de 12 grados de libertad
(2.9)
(2.10a)
(2.10b)
16
Paso 1. Discretizar y seleccionar el tipo de elemento
Considere el elemento de 12 grados de libertad mostrado en la figura 2.5. Cada
nodo tiene 3 grados de libertad: un desplazamiento transversal w en la dirección
z, una rotación θx sobre el eje x, y una rotación θy sobre el eje y.
Los desplazamientos nodales en el nodo i son:
=
yi
xi
iw
d
θθ
Donde las rotaciones se relacionan con los desplazamientos transversales por:
y
wx ∂
∂=θ x
wy ∂
∂−=θ
El signo negativo en θy se debe a que se requiere un desplazamiento negativo de
w para producir una rotación positiva sobre el eje y. La matriz de desplazamiento
total del elemento es:
=
n
m
j
i
d
d
d
d
d
Paso 2. Seleccionar las funciones de desplazamiento
Como el elemento placa tiene 12 grados de libertad, se selecciona un polinomio
de 12 términos en x e y tal como:
( ) 312
311
310
29
28
37
265
24321, xyayxayaxyayxaxayaxyaxayaxaayxw +++++++++++=
La función dada arriba es un polinomio incompleto de cuarto orden. Sin embargo,
vale la pena notar que el polinomio está completo hasta el tercer orden (primeros
10 términos), y se debe seleccionar los dos últimos términos de entre los
restantes 5 términos del cuartillo completo. La selección de x3y y xy3 asegura que
habrá continuidad en el desplazamiento en los contornos entre elementos. Los
términos x4 y y4 ocasionarían discontinuidades a lo largo de los contornos entre
elementos. El término final x2y2 no puede ser emparejado con ningún término así
que también es despreciado. La aproximación de la función de desplazamiento
(2.11)
(2.13)
(2.14)
(2.12)
17
también satisface la ecuación diferencial básica sobre la parte no cargada de la
placa. Además, la función describe el movimiento de cuerpo rígido y deformación
constante en la placa. Sin embargo, la discontinuidad en la pendiente entre
elementos a lo largo de las fronteras comunes a elementos no está asegurada.
Las constantes a1 a la a12 pueden ser determinadas expresando las 12
ecuaciones simultáneas conectando los valores de w y su pendiente en los nodos
cuando las coordenadas toman los valores apropiados.
−−−−−−−=
∂∂−
∂∂
12
3
2
1
3222
2322
33322322
202202010332020100
1
a
aaa
yyxyxyxyxxyxyxyxyxxyyxyxyyxxyxyxyx
x
wy
w
w
⋮
O también expresado como:
Pa=Ψ
Donde P es la matriz de orden 3x12 del lado derecho de la ecuación anterior.
Entonces se evalúa la matriz en cada punto nodal:
−−−−−−−−
−−−−−−−−−=
12
4
3
2
1
322
33322322
3222
2322
33322322
202202010
1
202202010
332020100
1
a
a
a
a
a
yyxyyxxyx
yxyxyyxyxxyyxxyx
yyxyyxxyx
yxxyyxxyx
yxyxyyxyxxyyxxyx
yw
w
nnnnnnnnn
jjjjjjjjjjjjjj
iiiiiiiii
iiiiiiiii
iiiiiiiiiiiiiiii
yn
j
yi
xi
i
⋮⋮⋮
θ
θθ
En una forma compacta las ecuaciones anteriores vienen dadas por:
Cad =
Entonces, las constantes a pueden ser obtenidas por:
dCa 1−=
Sustituyendo la forma anterior en la forma general de la matriz se obtiene:
dPC 1−=Ψ o Nd=Ψ
Donde 1−= PCN es la función de forma de la matriz.
(2.15)
(2.16b)
(2.17)
(2.16a)
18
Paso 3. Definir las relaciones deformación (curvatura)/desplazamiento y esfuerzo
(momento)/curvatura
Recordando la forma general de las curvaturas:
2
2
x
wkx ∂
∂−= 2
2
y
wky ∂
∂−= yx
wkxy ∂∂
∂−=22
La matriz de curvatura puede ser escrita como:
−−−−−−−−−−−−−
=
212
211985
121096
11874
66442
6622
6262
yaxayaxaa
xyayaxaa
xyayaxaa
k
k
k
xy
y
x
O en una forma matricial como:
Qak =
Donde Q es la matriz de coeficientes multiplicada por a en la matriz de curvaturas.
−−−−−−−−−
−−−−=
12
3
2
1
22 660440020000606200200000060026002000
a
aaa
yxyxxyyx
xyyxQ
⋮
Por tanto:
Qak = dQCk 1−=⇒ o Bdk =
Donde: 1−= QCB
La matriz momento/curvatura para una placa está dada por:
DBd
k
k
k
D
M
M
M
M
xy
y
x
xy
y
x
=
=
=
Donde D para materiales isotrópicos está dada por:
( ) ( )
−−=
21000101
112 2
3
ν
ν
νv
EhD
Paso 4. Derivación de la matriz de rigidez del elemento y ecuaciones
La matriz de rigidez está dada por la forma usual de la matriz de rigidez que es:
∫∫= dxdyDBBk T
(2.18)
(2.19)
(2.20)
(2.21)
(2.22)
(2.23)
(2.24)
19
La matriz de rigidez para el elemento rectangular de 4 nodos es de orden 12x12.
La expresión explícita para la matriz de rigidez para el elemento rectangular
consta en el anexo A. Las buenas características de convergencia de este
elemento de flexión en placas se muestran en la figura 2.6, en donde se aprecia
que para un número relativamente pequeño de elementos por cuarto de placa el
error se vuelve prácticamente cero.
Figura 2.6 Características de convergencia del elemento rectangular [22].
2.1.2.6 Cálculo de esfuerzos
Generalmente los esfuerzos se calculan en los puntos nodales de los elementos
mediante la ecuación general:
( ) ( )( ) ( )( ) ( )( )Nee
Ne
NN dSEDNdE === εσ
También se pueden escribir los esfuerzos resultantes en los puntos nodales
como:
( )
( )N
xy
y
xN
mmm
m
=
(2.25)
20
Y para el elemento enésimo:
( ) ( )( )Nee
Ne dSm =
Donde Se(N) es la matriz de esfuerzo de la placa. Una versión más explícita de la
matriz de esfuerzo es:
( ) ( )( )
( )
( )Ne
N
i
T
T
T
N
ieN
i d
Nyx
Ny
Nx
EEDNdm ⋅
∂∂∂−
∂∂−
∂∂−
==2
2
2
2
2
2
Donde N es la matriz columna de la función de forma correspondiente y E
representa la matriz de elasticidad para placas isotrópicas y homogéneas.
( ) ( )
−−
=2100
01
01
112 2
3
νν
ν
νEh
E
Por lo tanto, la matriz de esfuerzos puede escribirse como:
( ) [ ]( )
( )
( )( )Nee
N
NNe dS
d
d
d
d
SSSSm =
=
4
3
2
1
4321
En el anexo B consta la expresión para Se(N) para un elemento rectangular.
Debido a que los esfuerzos se calculan en los nodos de cada elemento, esto
puede causar que en nodos comunes el valor varíe considerablemente. Estas
variaciones se pueden resolver promediando los resultados de los elementos que
compartan el mismo nodo [22].
2.1.3 ELEMENTO TRIANGULAR DKT
Este procedimiento para la formulación de elementos triangulares de flexión para
placas delgadas está basado en la teoría discreta de Kirchhoff. El elemento
resultante se conoce como el elemento triangular discreto de Kirchhoff (DKT por
sus siglas en inglés) [2] [22] [23].
(2.26)
(2.27)
(2.28)
(2.29)
21
Las características de convergencia del elemento DKT han sido extensivamente
examinadas. Debido a su excelencia probada, ha permanecido como uno de los
mejores elementos triangulares para flexión de placas. El enfoque DKT alternativo
usa la teoría de Reissner – Mindlin para placas moderadamente gruesas y asume
que las deflexiones nodales y rotaciones son independientes unas de otras.
De esta manera, las funciones de forma para los componentes del
desplazamiento pueden ser continuas en las fronteras entre elementos. Como la
placa a ser analizada no es moderadamente gruesa sino delgada, la deformación
transversal al cortante debe ser cero en puntos específicos. Las suposiciones de
la teoría de placas de Kirchhoff se introducen en puntos discretos a lo largo del
contorno del elemento. Las funciones de forma están diseñadas para mantener
compatibilidad, así el elemento es conforme. El elemento DKT se muestra en la
figura 2.7.
Figura 2.7 Elemento DKT
La energía de flexión se puede representar en la forma siguiente:
∫=A
bT
b dxdykEkU2
1
Donde:
( ) ( )
−−
=2100
01
01
112 2
3
ν
ν
νv
EhEb
(2.30)
(2.31)
22
Las curvaturas están dadas por:
+=
xyyx
yy
xx
k
,,
,
,
θθθθ
Los componentes del desplazamiento u, v y w en cualquier punto se pueden
representar como:
( )yxzu x ,θ= ( )yxzv y ,θ= ( )yxww ,=
Donde w es el desplazamiento transversal y θx y θy son las rotaciones en la
dirección normal a los planos xz y yz respectivamente.
El elemento DKT usa componentes de desplazamiento nodales w, θx, θy, para
cada uno de los tres nodos esquineros. Las rotaciones θx y θy se aproximan por:
∑=
=6
1ixiix N θθ ∑
=
=6
1iyiiy N θθ
Las funciones de forma Ni, formadas con coordenadas naturales como polinomios
cuadráticos, se expresan como:
( )( )ηξηξ −−−−= 21121N
( )122 −= ξξN
( )123 −= ηηN
ξη44 =N
( )ηξη −−= 145N
( )ηξξ −−= 146N
Donde ξ y η son coordenadas de área.
Para la obtención del elemento empleando la teoría discreta de Kirchhoff, se debe
imponer las suposiciones de Kirchhoff en puntos particulares. Para el elemento
mostrado en la figura 2.8, las suposiciones resultan en:
0,
, =
++
=yy
xx
ww
θθ
γ
En los nodos de las esquinas, y
0, =+ sktk wθ
(2.32)
(2.33)
(2.34)
(2.35)
(2.36)
(2.37)
23
Donde k = 4, 5, 6. Se asume que la rotación sobre la dirección s en la figura 2.8
en los nodos del medio es un valor promedio de las rotaciones de los nodos de
los extremos:
( )ajaiak θθθ +=2
1
Donde los nodos de las esquinas son ij = 23, 31 y 12.
Figura 2.8 Geometría de un elemento triangular
(a) Elemento triangular y nodos, (b) Cantidades de una frontera particular,
(c) Relación entre θx, θy y θa, θt
La variación de los desplazamientos transversales se representa por:
sjjij
siiij
sk wwl
wwl
w ,,, 4
1
2
3
4
1
2
3 ++−−=
Al sustituir las ecuaciones anteriores en la ecuación 2.34 se obtienen las
funciones de forma para θx y θy, lo que resulta en:
ixx vN=θ i
yy vN=θ
La matriz de rigidez del elemento DKT se formula a partir del principio de trabajo
virtual, dado por la ecuación:
( )∫=A
Ti dAEDuDuW δδ
( ) ( )∫∫ +=A V
T
A KBT
Ki dAuDEuDdADEDW γγδθθδδ
(2.39)
(2.40)
(2.38)
24
∫∫ +=A V
TT
A KBTT
Ki dAuDEDudADEDW γγδθδθδ
Como el efecto de la deformación al cortante se desprecia, las funciones de forma
para θx y θy se sustituyen en la primera parte de la ecuación 2.41, resultando:
iiT
A
iB
TiTi
A
iKB
TTK
iT
A KBTT
Ki
kvvdAvBEBvW
dAvNDEDNvdADEDW
δδδ
δθδθδ
==
==
∫
∫∫
Entonces, la matriz de rigidez está dada por:
∫ ∫−
=1
0
1
0
2η
ηξddBEBAk bT
En donde la matriz deformación – desplazamiento está dada por:
++−−−−
+=
Ty
Ty
Tx
Tx
Ty
Ty
Tx
Tx
NyNyNxNx
NxNx
NyNy
AB
ηξηξ
ηξ
ηξ
,12,31,12,31
,12,31
,12,31
2
1
Donde A es el área del elemento y 311212312 yxyxA −=
Las matrices Nx,ξ, Ny,η, Nx,ξ y Ny,η están dadas por:
( ) ( )( ) ( )
( ) ( ) ( )( ) ( )
( ) ( )( ) ( )( )
( )( )
−−−+−
−+−++−−−−++−−
+−−+++−+−−−+−
=
45
54
45
646
466
646
656
656
656
,
21622121
21642121
rrqqPP
rrrqqqPPP
rrrqqqPPP
Nx
ηηη
ηξξηξηξ
ηξηξηξηξ
ξ
( ) ( )( ) ( )
( ) ( )( ) ( )
( ) ( )( ) ( )
( )( )( )
−−−+−
−+−−−−+−
++−−+−−−+−−+−
−+−
=
54
54
54
646
646
646
656
656
656
,
21211
2121
21121
qqrrtt
qqqrrr
tttqqqrrr
ttt
N y
ηηη
ηξηξ
ηξηξ
ηξηξ
ξ
(2.42)
(2.43)
(2.44)
(2.45a)
(2.45b)
(2.41)
25
( ) ( )( ) ( )
( ) ( ) ( )( )( )
( )( ) ( )( ) ( )
( ) ( )
−+−++−−+−+−−
−−−+
+−−+++−+−−−−−−
=
545
545
545
46
64
64
655
655
565
,
21622121
21642121
rrrqqqPPP
rrqqPP
rrrqqqPPP
Nx
ξηηξηξη
ξξξ
ξηηξξηξη
η
( ) ( )( ) ( )
( ) ( )( )( )( )
( ) ( )( ) ( )
( ) ( )
−+−−−+−−
+−−−−
−+
++−−+−−+
−−−−
=
545
545
545
64
64
64
655
655
565
,
21211
21
2121121
qqqrrr
tttqq
rrtt
qqqrrr
ttt
N y
ξηξη
ξηξξξ
ξηξηξη
η
Donde: 26 ijijk lxP −= 23 ijijijk lyxq =
223 ijijk lyr = 26 ijijk lyt −=
k = 4, 5, 6 para ij = 23, 31, 12, respectivamente.
Para facilitar la manipulación de estas matrices y la integración [2], se introduce la
forma explícita para la matriz de rigidez:
αα EA
K Te 2
1=
Donde A representa el área del triángulo y E es la matriz del material elástico,
isotrópico y homogéneo, y está dada por:
( )
−
=
RvD
DRvDR
vDRDR
E
2
100
0
0
24
1
La matriz D es la rigidez flexional de la placa dada por:
( )2
3
112 v
EhD
−=
(2.45c)
(2.45d)
(2.46)
(2.47)
(2.48)
26
Además:
=
211121112
R
En el anexo C consta la expresión explícita para la matriz α.
2.1.3.1 Cálculo de esfuerzos
Conociendo los desplazamientos nodales de del elemento en su sistema de
coordenadas locales, los esfuerzos se pueden calcular a partir de:
( ) ( )( )
( )( )Neb
N
xy
y
xN
e dLEmmm
m αηξ24
1, =
=
Donde Eb es la matriz de elasticidad dada por:
( ) ( )
−−=
21000101
112 2
3
vv
v
v
EhEb
Además:
=l
ll
L00
0000
Con:
0000 = ξηξ −−= 1l
Dado que la matriz L depende de ξ y η, el cálculo de los esfuerzos en cualquier
punto del elemento es posible [22].
2.2 PLACAS PLANAS
Una placa plana es un elemento bidimensional de manera que dos coordenadas
definen cualquier posición del elemento. Existen dos tipos de análisis con placas
planas: esfuerzo plano y deformación plana. En este estudio solo se tratará el
primer caso.
(2.49)
(2.50)
(2.51)
(2.52)
27
2.2.1 ESFUERZO PLANO
Se define como un estado de esfuerzo en el cual el esfuerzo normal σz y los
esfuerzos al cortante τxz y τyz dirigidos en dirección perpendicular al plano se
asumen igual a cero. Se define por las siguientes suposiciones:
1. El cuerpo es pequeño en una dirección (z por convención) en comparación
con las otras dimensiones, además el espesor en z es uniforme y es menor
a 1/10 de la dimensión más pequeña en el plano xy.
2. El cuerpo está sujeto a cargas solo en el plano xy.
3. El material es linealmente elástico, isotrópico y homogéneo.
En la figura 2.9 se muestra una ilustración de las condiciones de esfuerzo plano.
Figura 2.9 Condiciones de esfuerzo plano [11]
2.2.2 ELEMENTO PLANO RECTANGULAR
Considérese el elemento rectangular con cuatro nodos y ocho grados de libertad
como se muestra en la figura 2.10. Para describir el campo de desplazamientos
en este elemento bidimensional, se asume que los desplazamientos están en el
plano xy, además se representa el desplazamiento u con dos componentes
rectilíneos ux y uy [23], así:
28
( )( )
= ηξηξ,,
y
x
uu
u
Se introduce una solución de prueba u con 8 incógnitas, con 8 grados de libertad
(GDL). Se asume que un polinomio puede representar adecuadamente cada uno
de los desplazamientos ux y uy. Para el desplazamiento en la dirección x, se
escoge un polinomio bilineal:
( ) [ ]
=+++==
4
3
2
1
4321
ˆˆˆˆ
1ˆˆˆˆˆ,
uuuu
uuuuuNu xuxx ηξηξηξηξηξ
Donde el origen de coordenadas ξ, η está en el nodo 1. Este polinomio contiene
las cuatro incógnitas 4321 ˆ,ˆ,ˆ,ˆ uuuu . De forma similar, los desplazamientos en la
dirección y se escogen para tener la forma del polinomio:
( ) ηξηξηξ 8765 ˆˆˆˆˆ, uuuuuNu yuyy +++==
De manera que los desplazamientos asumidos son:
=
8
1
ˆ
ˆ
1001
u
u
uu
y
x ⋮⋮⋮
ηξηξηξηξ
Figura 2.10 Elemento plano rectangular
Si se considera solamente ux, que en el nodo 1 toma el valor:
(2.53)
(2.54a)
(2.54c)
(2.54b)
29
[ ]
=
4
3
2
1
1
ˆˆˆˆ
0001
uuuu
ux
Donde ux1 es el desplazamiento del nodo 1 en la dirección x. Entonces la relación
para las cuatro esquinas es:
xux
x
x
x
x
x uN
uuuu
uuuu
v ˆˆ
ˆˆˆˆ
1001111100110001
4
3
2
1
4
3
2
1
=
=
=
Entonces:
xxxuxx vGvNu == −1ˆˆ
Y, por lo tanto, de la ecuación 2.54a:
( ) xxxxuxx vNvGNu ==ηξ ,
Con:
[ ] ( )( ) ( ) ( )[ ]ξηηξηξηξηξηξ −−−−=
−−−
− 1111
10011111
00110001
1
uxN ( )1ˆ −= uxx NG xN
De manera que en la ecuación 2.58 los desplazamientos supuestos ux se han
expresado en términos de los desplazamientos nodales desconocidos. De igual
forma se procede para uy. Si los desplazamientos ux y uy se colocan juntos:
( )( ) ( ) ( )( )( ) ( ) ( )
−−−−
−−−−=
4
3
2
1
4
3
2
1
1111001111
y
y
y
y
x
x
x
x
y
x
uuuuuuuu
uu
ξηξηηξηξξηξηηξηξ
Que también se expresa como:
ivNNNN
NNNNu
=4321
4321
Los desplazamientos asumidos contienen 8 GDL con 2 GDL por nodo.
(2.55)
(2.56)
(2.57)
(2.58)
(2.59)
(2.60b)
(2.60a)
30
2.2.2.1 Principio del trabajo virtual
La matriz de rigidez para un elemento se deriva al sustituir el campo de
desplazamientos asumidos de la ecuación 2.60 en la ecuación que representa el
principio del trabajo virtual. En notación matricial, el principio del trabajo virtual se
expresa como:
0=−−=− ∫∫∫���������������
III
S
T
II
V VT
I
V
T
P
dSpudVpudVW δδσδεδ
Se supone que solo existen cargas en la superficie, por lo que no es necesario
considerar la segunda integral. En términos de M elementos, la relación para el
trabajo virtual puede ser expresada como:
01
=
−=− ∑ ∫∫=
M
i
i
III
S
T
I
V
T
P
dSpudVW����������
δσδεδ
Las deformaciones se obtienen de la función de desplazamiento u empleando la
relación:
Du=ε
Donde el operador matricial es:
∂∂∂
∂==
xy
y
x
uDD 00
Donde, como x = aξ, y = bη:
ξξ
∂=∂∂=
∂∂=∂
aaxx
11
ηη
∂=∂∂=
∂∂=∂
bbyy
11
Las deformaciones en términos de las funciones de prueba son:
ii
y
x
xy
y
x
xy
y
x
BvDNvuu ==
∂∂∂
∂=
00
γεε
Donde B se obtiene al aplicar las derivadas en D a las funciones de prueba N.
(2.61)
(2.62)
(2.63)
(2.64)
(2.65)
31
( )( ) ( ) ( )( )( ) ( ) ( )
−−−−
−−−−
∂∂
∂
∂
ξηξηηξηξξηξηηξηξ
ξη
η
ξ
1111000000001111
11
10
01
ab
b
a
D N
−−−−−−−−
−−−−
−−−−
=
aaaabbbb
bbbb
aaaa
ηηηηξξξξ
ξξξξ
ηηηη
1111
110000
000011
B
2.2.2.2 Ley del material
La ecuación 2.65 es la relación deformación – desplazamiento necesaria en la
primera integral de la ecuación 2.62. Además, para expresar el esfuerzo en
términos de los desplazamientos, una relación esfuerzo – deformación, se
necesita la ley del material. Si se asume una condición de esfuerzo plano para
este caso, entonces la ley del material es:
( )
−−=
xy
y
x
xy
y
x
vv
v
v
E
γεε
τσσ
21000101
1 2
2.2.2.3 Formación de la matriz de rigidez del elemento
Si se expresa el principio del trabajo virtual en términos de los desplazamientos
nodales, para un solo elemento la integral de volumen I es:
( ) udVkuudVEDDuudVEDuDdVV
DTuV
Tu
TuV
TuV
T
∫∫∫∫ === δδδσδε
La matriz de rigidez del elemento puede ser obtenida a partir de esta expresión.
Ahora, la integral que contiene kD debe ser expresada en forma discreta. Para
este problema bidimensional, h es el espesor constante del elemento y dV =
hdxdy = habdξdη. Por esto la integral puede ser expresada como:
(2.66)
(2.67)
(2.68)
32
( ) i
k
TiT
V
iu
Tu
iT
V uT
uT vdEBdBabhvNdVvEDNDvudVEDDu
i
��� ���� ��∫ ∫∫∫ ==1
0
1
0
ηξδδδ
Donde la matriz de rigidez es igual a:
∫ ∫∫ ==1
0
1
0
ηξdEBdBabhEBdVBk T
V
Ti
En el anexo D está dada de forma explícita la matriz de rigidez para el caso de
elementos rectangulares con 8 GDL. La matriz de rigidez del anexo D está
ensamblada de forma tal que los primeros 4 desplazamientos corresponden a la
dirección x, mientras que los siguientes 4 están en la dirección y [23].
2.2.3 ELEMENTO TRIANGULAR CST
El elemento triangular CST (constant strain triangle) para esfuerzo plano es el
elemento más fácil para ser desarrollado matemáticamente. Está formado por 3
nodos con 2 desplazamientos en el plano por nodo, además la variación de los
desplazamientos dentro del elemento es lineal [11].
Figura 2.11 Notación para un elemento en estuerzo plano
En la figura 2.11 se describe un elemento triangular con tres nodos que
representa un cuerpo sujeto a esfuerzo plano. Los nodos del elemento se
numeran como se muestra en la figura, los desplazamientos nodales en la
dirección x son u1, u2 y u3, mientras que los desplazamientos en la dirección y son
v1, v2 y v3. Para el elemento triangular en esfuerzo plano los desplazamientos son.
(2.69)
(2.70)
33
( ) ( ) ( ) ( ) NuuyxNuyxNuyxNyxu =++= 332211 ,,,,
( ) ( ) ( ) ( ) NvvyxNvyxNvyxNyxv =++= 332211 ,,,,
Donde N1, N2 y N3 son las funciones de interpolación. Usando la representación
discretizada del campo de desplazamientos, los componentes de las
deformaciones del elemento son entonces:
33
22
11 u
x
Nu
x
Nu
x
N
x
ux ∂
∂+
∂∂+
∂∂=
∂∂=ε
33
22
11 v
y
Nv
y
Nv
y
N
y
vy ∂
∂+
∂∂
+∂
∂=
∂∂=ε
33
22
11
33
22
11 v
x
Nv
x
Nv
x
Nu
y
Nu
y
Nu
y
N
x
v
y
uxy ∂
∂+
∂∂
+∂
∂+
∂∂
+∂
∂+
∂∂
=∂∂+
∂∂=γ
Si se define el campo de desplazamientos como:
( )
=
3
2
1
3
2
1
vvvuuu
eδ
La matriz de deformación del elemento puede ser expresada como:
( )eB
vvvuuu
x
N
x
N
x
N
y
N
y
N
y
Ny
N
y
N
y
Nx
N
x
N
x
N
δε =
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
=
3
2
1
3
2
1
321321
321
321
000
000
Donde B es la matriz de derivadas parciales de las funciones de interpolación,
también conocida como la matriz de deformación – desplazamiento. Se puede
observar que las derivadas parciales son constantes ya que las funciones de
interpolación son lineales. Por lo tanto, las componentes de las deformaciones
son constantes a lo largo del volumen del elemento. Consecuentemente, el
elemento triangular de tres nodos para esfuerzo plano se conoce como un
triángulo con deformación constante CST (constant strain triangle). La energía
elástica de deformación del elemento es:
( ) ( )
( )
( ) ( ) ( ) ( )( ) ( )eTeTe
V
eTTe
V
eTee DBBVdVDBBdVDU
ee
δδδδεε2
1
2
1
2
1 === ∫∫∫∫∫∫
(2.71)
(2.72)
(2.73)
(2.74)
(2.75)
34
Donde V(e) es el volumen total del elemento.
Considerando las fuerzas como en la figura 2.11b, empleando la notación:
=
y
y
y
x
x
x
ffffff
f
3
2
1
3
2
1
De manera que se puede expresar el trabajo de las fuerzas externas como:
fW Tδ=
La energía potencial total para un elemento es entonces:
fDBBV
WU TTTe
e δδδ −=−=Π2
Considerando que el elemento se encuentra en equilibrio, se tiene:
0=∂
Π∂
iδ 6,,1⋯=i
Por lo que resulta la relación matricial:
fDBBV Te =δ
Esta ecuación matricial es de la forma:
fk =δ
Donde k es la matriz de rigidez del elemento dada por:
DBBVk Te=
Donde D es la matriz del material:
( )
−−=
21000101
1 2v
vv
v
ED
2.2.3.1 Evaluación de la matriz de rigidez
La matriz de rigidez dada por la ecuación 2.82 se puede evaluar conociendo las
funciones de interpolación, las cuales, para el elemento triangular, son:
(2.76)
(2.77)
(2.78)
(2.79)
(2.80)
(2.81)
(2.82)
(2.83)
35
( ) ( ) ( ) ( )[ ] ( )
( ) ( ) ( ) ( )[ ] ( )
( ) ( ) ( ) ( )[ ] ( )yxA
yxxxyyyxyxA
yxN
yxA
yxxxyyyxyxA
yxN
yxA
yxxxyyyxyxA
yxN
333122112213
222311331132
111233223321
2
1
2
1,
2
1
2
1,
2
1
2
1,
γβα
γβα
γβα
++=−+−+−=
++=−+−+−=
++=−+−+−=
Y las derivadas parciales requeridas son:
( )A
yyAx
N
22
1 132
1 β=−=∂
∂ ( )
Axx
Ay
N
22
1 123
1 γ=−=
∂∂
( )A
yyAx
N
22
1 213
2 β=−=∂
∂ ( )
Axx
Ay
N
22
1 231
2 γ=−=
∂∂
( )A
yyAx
N
22
1 321
3 β=−=
∂∂
( )A
xxAy
N
22
1 312
3 γ=−=
∂∂
Entonces B está dada por:
−−−−−−−−−
−−−=
211332123123
123123
211332
000000
2
1
yyyyyyxxxxxxxxxxxx
yyyyyy
AB
=
321321
321
321
000000
2
1
βββγγγγγγ
βββ
AB
Nótese que para espesor constante el volumen del elemento es hA, al sustituir en
la ecuación 2.82 resulta en:
( )
−
−=
321321
321
321
33
22
11
33
22
11
2000
000
2
100
01
01
0
0
0
0
0
0
14βββγγγγγγ
βββ
νν
ν
βγβγβγγβγβγβ
νA
Ehk
En el anexo E está dada de forma explícita la matriz de rigidez para el elemento
triangular CST en esfuerzo plano.
(2.86)
(2.84)
(2.85)
(2.87)
36
2.3 CONVERGENCIA
Una de las características que debe cumplir el tipo de elemento escogido es que
al realizar el análisis este converja a la solución exacta a medida que el tamaño
de los elementos tiende a cero. Es decir, la exactitud de la solución se incrementa
a medida que la malla de elementos finitos es continuamente refinada. Para
conseguir convergencia, el elemento debe ser completo y compatible (o
conforme). Esto significa que las funciones de desplazamiento deben ser capaces
de representar el desplazamiento de cuerpo rígido y los estados de deformación
constante. La compatibilidad asegura que no ocurran saltos dentro de los
elementos y entre los elementos cuando el sistema de elementos se ensambla y
se carga. Se considera las siguientes características en más detalle.
1. Las funciones de prueba deben ser capaces de representar
desplazamientos a los que los elementos están sometidos como un cuerpo
rígido sin desarrollar esfuerzos.
2. Las funciones de desplazamiento de un elemento deben ser tales que la
deformación en cada elemento se aproxime a un valor constante en el
límite cuando el elemento se acerque a un tamaño muy pequeño.
Los elementos que no satisfacen los requerimientos de compatibilidad son
llamados incompatibles o no conformes. Si la incompatibilidad desaparece al
incrementar el refinamiento de la malla, los elementos aún pueden ser aceptables
si convergen a la solución correcta [23].
2.4 CRITERIO GENERAL DE CEDENCIA EN TRES DIMESIONES
Una vez que se ha calculado el estado de esfuerzos en el que se encuentra un
sólido, es necesario determinar si el material del que está hecho este sólido
soportará los esfuerzos que en el se han inducido. Si esto no sucede se dice que
el material ha fallado y que es necesario utilizar otro material para garantizar que
el sólido pueda soportar las cargas a las que está sometido. En el caso de la
teoría clásica de la plasticidad se considera que la falla se produce en el momento
37
en que las cargas actuantes son suficientes para llevar al material al estado de
cedencia o fluencia. Para poder establecer un criterio de falla en el caso de
esfuerzos tridimensionales, es necesario plantear una serie de postulados en
relación al inicio de la cedencia [19].
2.4.1 POSTULADO 1: EXISTE UNA FUNCIÓN GENERAL DE CEDENCIA
Se supone que la cedencia depende solo del estado de esfuerzos y no de cómo
se ha alcanzado este estado. Esta ley define el límite de elasticidad y será
descrita como un criterio para la cedencia. Se asume la existencia de una función
de cedencia ( )ijf τ tal que:
Material de comportamiento elástico
si ( )ijf τ < 0
o si ( )ijf τ = 0 y ( )ijf τɺ < 0
Donde ( )ijf τ = 0 define la superficie de cedencia en el espacio de esfuerzos y
( )ijf τɺ < 0 indica que se está descargando el material.
Material de comportamiento plástico
si ( )ijf τ = 0 y ( )ijf τɺ ≥ 0
La función no positiva f, debido a la simetría del tensor de esfuerzos, puede ser
expresada como:
( )123123332211 ,,,,, ττττττff =
Donde se tienen 6 variables de esfuerzo independientes. Además, la función f
puede ser expresada como una función de los esfuerzos principales τ1, τ2, τ3 y 3
ángulos α1, α2, α3 que definen la dirección de los ejes principales, así:
( )321321 ,,,,, ααατττff =
2.4.2 POSTULADO 2: EL MATERIAL ES ISOTRÓPICO
Si el material es isotrópico, no puede haber direcciones preferidas y la función
debe tener la misma forma sin importar la orientación de los ejes. Esto significa
(2.88)
(2.89)
(2.90)
(2.91)
38
que la función f de la ecuación 2.91 no es una función de los ángulos α1, α2, α3
que definen la dirección de los ejes principales; es decir:
( )321 ,, τττff =
Además, la función f no debería cambiar si los ejes son intercambiados. Por lo
que se puede concluir:
( ) ( ) ( )231312321 ,,,,,, τττττττττ fff == , etc.
2.4.3 POSTULADO 3: LA CEDENCIA ES INDEPENDIENTE DEL ESFUERZO
HIDROSTÁTICO
Esto significa que la función f depende solo de los esfuerzos principales
deviatóricos s1, s2, s3. Esto es:
( )321 ,, sssff =
Además, se pueden usar los invariantes de esfuerzos deviatóricos J1 = Is, J2 = -IIs,
J3 = IIIs, en lugar de los esfuerzos principales deviatóricos, Nótese que por
definición J1 = 0, entonces:
( )32 , JJff =
2.4.4 POSTULADO 4: COMPORTAMIENTO IDÉNTICO A TENSIÓ N Y
COMPRESIÓN
Se requiere que el valor de la función de cedencia no cambie cuando los signos
de los esfuerzos cambien. Esto es:
( ) ( )ijij ff ττ −=
Estos son los postulados claves para el desarrollo de una función de cedencia
( )ijf τ en la teoría clásica de plasticidad.
2.5 HIPÓTESIS DE FALLA PARA MATERIALES DÚCTILES
Cuando un elemento se somete a una carga de manera que el estado de
esfuerzos es uniaxial, entonces el esfuerzo y la resistencia se pueden comparar
(2.92)
(2.93)
(2.94)
(2.95)
(2.96)
39
directamente para determinar si se producirá la falla. Sin embargo, no existe una
formulación teórica que represente exactamente el estado de esfuerzos para la
condición de cedencia para el caso tridimensional. A continuación se exponen dos
hipótesis para la predicción de fallas que han demostrado ser bastante acertadas
[20].
2.5.1 HIPÓTESIS DEL ESFUERZO CORTANTE MÁXIMO (CRITE RIO DE
TRESCA)
Esta hipótesis expresa que la fluencia comienza cuando el esfuerzo cortante
máximo de cualquier elemento iguala al esfuerzo cortante máximo en una probeta
de ensayo a la tensión del mismo material cuando esa probeta comienza a fluir.
Si se ordenan los esfuerzos normales principales como σ1 > σ2 > σ3, entonces la
hipótesis del esfuerzo cortante máximo predice que la fluencia ocurrirá cuando:
2y
máx
S≥τ o yS≥− 31 σσ
En donde la resistencia a la fluencia en cortante está dada por la ecuación:
ysy SS ⋅= 5.0
2.5.2 HIPÓTESIS DE LA ENERGÍA DE DISTORSIÓN (CRITER IO DE VON
MISES)
La hipótesis de la energía de distorsión predice que la fluencia ocurrirá cuando la
energía de distorsión en un volumen unitario iguale la energía de distorsión en el
mismo volumen cuando se someta a un esfuerzo uniaxial hasta la resistencia de
fluencia. Por lo tanto predice que la fluencia ocurrirá cuando:
( ) ( ) ( ) 21213
232
221
2'
−+−+−=
σσσσσσσ
yS≥'σ
Para el estado de esfuerzos biaxial, sean σ1 y σ2 los dos esfuerzos principales.
Entonces, de la ecuación 2.99, se obtiene:
(2.97)
(2.98)
(2.99)
40
( ) 212221
21' σσσσσ +−=
Esta ecuación representa una elipse de esfuerzos principales tal como se ilustra
en la figura 2.12.
Figura 2.12 Hipótesis de falla
Como se muestra en la figura 2.12 el criterio de Tresca es más conservador que
el criterio de von Mises ya que está dentro de la elipse. Adicionalmente, el criterio
de von Mises proporciona una estimación razonable de la falla por fatiga,
especialmente en casos de repetidas cargas de tensión y tensión – corte.
(2.100)
41
CAPÍTULO 3.
GENERACIÓN DE MALLAS NO ESTRUCTURADAS
3.1 INTRODUCCIÓN
El advenimiento de métodos numéricos basados en la subdivisión espacial de un
dominio en elementos, inmediatamente trae como consecuencia la tarea de la
generación automática de una malla. Esto principalmente para el caso de técnicas
numéricas basadas en los mallados no estructurados. Mallas típicas consisten de
elementos triangulares o tetraédricos que pueden tener un número arbitrario de
elementos vecinos. Otros elementos incluyen trapezoides, ladrillos (bricks) o
prismas. Con la disponibilidad de computadoras más potentes, se ha ansiado la
simulación de sistemas geométricos y físicos cada vez más complejos. Durante la
última década ha habido una considerable cantidad de esfuerzo dedicado a la
generación automática de mallas, resultando en dos técnicas muy poderosas: el
avance frontal y la triangulación de Delaunay.
Cualquier generador automático de mallas requiere como entrada la siguiente
información:
a) Definición de la superficie, es decir la descripción del contorno de las
superficies del dominio a ser mallado.
b) Tamaño y forma de la malla, es decir la descripción de cómo estarán en el
espacio el tamaño del elemento, forma y orientación.
c) Tipo de elemento, por ejemplo triangular, tetraédrico, etc.
d) Técnica de generación de malla, es decir la selección de un método
adecuado para conseguir la generación de la malla deseada.
Estrictamente hablando, una malla estructurada puede ser reconocida porque
todos los nodos interiores de la malla tienen un mismo número de elementos
adyacentes. La malla generada por un generador de mallas estructuradas
típicamente está compuesta de elementos rectangulares o hexagonales. Los
42
algoritmos empleados generalmente involucran complejas técnicas de suavizado
iterativas que intentan alinear los elementos con los contornos o dominios físicos.
En la generación de mallas no estructuradas, en cambio, está permitido cualquier
número de elementos conectados en un nodo. Las mallas triangulares y
tetraédricas son las más comunes. Las mallas no estructuradas son mucho más
flexibles ya que no están sometidas a ningún posicionamiento especial de sus
nodos.
El desarrollo de las técnicas de generación de mallas no estructuradas ha
avanzado de forma casi totalmente independiente al de las mallas estructuradas.
Esto se debe a que los métodos de resolución por diferencias finitas utilizan
mallas estructuradas mientras que las mallas no estructuradas se emplean
generalmente con el método de los elementos finitos.
3.2 DEFINICIÓN DE SUPERFICIES
El primer paso en la generación de una malla consiste en definir la forma de la
superficie del dominio. Para ello hay que modelar mediante algún tipo de curva la
geometría del objeto que se desea calcular.
La información disponible sobre la forma de un objeto sobre el cual se va a
generar una malla consiste en la determinación de una serie de puntos de
definición situados sobre su contorno. Por esta razón debe utilizarse algún tipo de
curva para interpolar la geometría del contorno entre esos puntos. Es
recomendable el utilizar tipos de curvas estándar, como las empleadas por
programas CAD para definir la geometría. Funciones tales como B-Splines y
NURBS se emplean para definir la superficie del dominio [1].
3.3 MÉTODOS PARA GENERAR MALLAS NO ESTRUCTURADAS
Si se considera la tarea de rellenar un dominio dado con elementos, solo existen
dos caminos para conseguirlo:
43
• Generación simultánea de nodos y elementos. Los algoritmos de este
procedimiento generan a la vez los nodos y los elementos, de forma que
cada vez que se genera un nuevo punto, éste se conecta inmediatamente
al resto de la malla. El método más utilizado dentro de este grupo es el
método del avance frontal AFM (advancing front method).
• Generación de nodos y elementos en fases distintas. En estas técnicas
primero se realiza la distribución de nodos sobre la malla y luego se hace la
triangulación de la malla para formar los elementos mediante la conexión
de los nodos ya existentes. La técnica más utilizada es el método de
triangulación de Delaunay.
Ambos grupos de procedimientos producen mallas de elementos triangulares [1].
3.3.1 MÉTODO DEL AVANCE FRONTAL
El método del avance frontal consiste en avanzar dentro del espacio no mallado
añadiendo un elemento cada vez. La región que separa la porción mallada de
espacio de la no mallada se llama frente.
Figura 3.1 (a) Frente inicial, (b) frente actual, (c) frente actualizado con el nuevo
elemento
El método del avance frontal recibe su nombre del hecho de que en cada etapa
del proceso existe un frente de generación formado por todos los posibles lados
44
sobre los cuales se va a generar un nuevo triángulo. De forma repetitiva el
proceso consta de una etapa de generación de un triángulo y de otra etapa de
actualización del frente, teniendo en cuenta la existencia del nuevo triángulo.
El algoritmo, que está bosquejado en la figura 3.1, puede ser resumido como
sigue [15]:
1. Definir los contornos del dominio a ser mallado. Sin entrar en detalle, se
asumirá alguna forma general consistente de zonas y líneas que rodean o
delimitan los contornos, y puntos en la intersección de líneas.
2. Definir la variación espacial del tamaño del elemento, direcciones de
estrechamiento para los elementos que se creen.
3. Empleando la información dada para la distribución del tamaño y forma del
elemento en el espacio y las definiciones de las líneas, se generan los
lados a lo largo de las líneas que conectan las zonas de superficie. Estos
lados forman el frente inicial de la triangulación de las zonas de superficie.
4. Se triangula las superficies usando la información dada para la distribución
del tamaño y forma del elemento en el espacio, los lados ya generados y la
definición de la superficie. Esto proporciona el frente inicial de caras.
5. Hallar los parámetros de generación (tamaño del elemento, direcciones de
estrechamiento) para estas caras.
6. Seleccionar la cara siguiente a ser borrada del frente, para evitar
elementos largos cruzando sobre regiones de elementos pequeños, la cara
que forme el elemento nuevo de menor tamaño se selecciona como la
siguiente cara a ser borrada de la lista de caras.
7. Para la cara que va a ser borrada: se selecciona la posición del “mejor
punto” para la introducción de un nuevo punto, luego se determina si existe
un punto en la malla recién generada que puede ser usado en lugar del
nuevo punto. Si existe tal punto, se coloca este punto en la lista y se
continúa buscando. Determinar si el elemento formado con el punto
seleccionado cruza alguna cara existente. En caso afirmativo, se
selecciona un nuevo punto y se vuelve a probar.
8. Añadir un nuevo elemento, punto y caras a sus respectivas listas.
45
9. Borrar las caras conocidas de la lista de caras.
10. Añadir las nuevas caras al frente.
11. Si hay alguna cara no seleccionada en el frente, ir a 6.
3.3.2 MÉTODO DE TRIANGULACIÓN DE DELAUNAY
El criterio de Delaunay no es un algoritmo para la generación de mallas, sino que
provee de un criterio para la conexión de un conjunto de puntos existentes en el
espacio. Por ello es necesario proporcionar adicionalmente un método para la
generación de las localizaciones de los nodos dentro de la geometría.
Figura 3.2 Triangulación de Delaunay (2D)
La técnica de triangulación de Delaunay tiene una larga historia en las
Matemáticas, Geofísica e Ingeniería. Dada una colección de puntos P = x1, x2,…,
xn, se puede definir un conjunto de regiones o volúmenes V = v1, v2,…, vn
asignados a cada uno de los puntos, que satisfacen la siguiente propiedad:
cualquier localización dentro de vi está más cerca de xi que cualquier otro de los
puntos:
jii xxxxv −<−Ρ= : ij ≠∀
Este conjunto de volúmenes V, que cubre el dominio completamente, se conoce
como la teselación de Dirichlet. Los volúmenes vi son poliedros convexos, y se
conocen como regiones de Voronoi. Juntando todos los pares de puntos xi, xj de
un lado a otro de los contornos de los poliedros resulta en una triangulación del
polígono convexo de P. Esta triangulación se conoce comúnmente como la
triangulación de Delaunay. El conjunto de triángulos (tetraedros) que forman la
46
triangulación de Delaunay satisfacen la propiedad que ningún otro punto está
contenido dentro del circuncírculo (circunesfera) formado por los nodos del
triángulo (tetraedro). Aunque esta y otras propiedades de las triangulaciones de
Delaunay han sido estudiadas por algún tiempo, procedimientos prácticos de
triangulación solo han aparecido en la última década.
El algoritmo de Delaunay es un proceso secuencial. Cada nuevo punto se
introduce en la estructura ya existente, la cual se rompe y se reconecta para
formar una nueva triangulación de Delaunay. Un procedimiento de triangulación
se resume aquí para el caso bidimensional [1]:
1. Definir un polígono convexo que contendrá todos los puntos. Especificar
cuatro puntos junto con su diagrama de Voronoi asociado.
2. Introducir un nuevo punto dentro del polígono xn+1.
3. Hallar todos los vértices del diagrama de Voronoi que serán eliminados.
Todos los vértices cuyo circuncírculo asociado contenga el nuevo punto
añadido serán eliminados.
4. Hallar los puntos de formación de los vértices de Voronoi a eliminar. Estos
puntos serán los que se unirán al punto recién añadido para formar nuevos
triángulos y rehacer la triangulación. La unión de cada uno de estos puntos
con el nuevo punto formará cada una de las nuevas caras de los nuevos
triángulos.
5. Determinar los vértices de Voronoi inmediatos a los vértices que van a ser
eliminados que a su vez no sean eliminables.
6. Determinar los puntos de formación de los nuevos vértices de Voronoi.
Estarán formados por el nuevo punto añadido junto con dos puntos más
que serán contiguos entre sí y formarán parte de alguno de los triángulos
ya existentes.
7. Determinar los vértices de Voronoi inmediatos a los nuevos vértices. Para
cada nuevo vértice se hace una búsqueda sobre los puntos de formación
de todos los vértices inmediatos hallados en el paso 5. En caso de que un
par de puntos de formación coincida para un nuevo vértice y para otro de
47
los hallados en 5 entonces estos dos vértices pasan a tener una conexión
inmediata.
8. Reordenar la estructura de datos eliminando las filas correspondientes a
los vértices eliminados.
9. Repetir los pasos 2 al 8 para cada nuevo punto.
La complejidad de la técnica de triangulación de Delaunay es ( )( )NNO log , donde
N expresa el número de elementos.
3.4 POST-PROCESO
Sería muy raro que algún algoritmo de generación de mallas sea capaz de definir
una malla que sea óptima sin alguna forma de post-proceso para mejorar la
calidad de los elementos. Las dos categorías principales de mejora de mallas
incluyen alisado y limpieza. Las técnicas de alisado incluyen cualquier método
que ajuste la localización de los nodos mientras se mantiene la conectividad de
los elementos. La limpieza generalmente se refiere a cualquier proceso que
cambie la conectividad de los elementos.
3.4.1 ALISADO
La mayoría de los procesos de alisado involucran alguna forma de proceso
iterativo que reposiciona los nodos individuales para mejorar la calidad local de
los elementos. Una gran variedad de técnicas de alisado han sido propuestas,
una de las más empleadas es el alisado Laplaciano [25].
3.4.1.1 Alisado Laplaciano
Geométricamente, el alisado Laplaciano es relativamente fácil de visualizar en el
espacio bidimensional. Como se muestra en la figura 3.3, un generador
automático de mallas puede generar elementos no uniformes. Algunos elementos
muestran tasas de deformación que causan imprecisión numérica en el análisis.
48
Figura 3.3 Malla (a) antes del alisado Laplaciano, (b) después
El objetivo del alisado Laplaciano es reubicar los nodos internos de la malla de
manera que los elementos sean lo más equiláteros posible.
El alisado Laplaciano se hace simplemente al promediar los vértices vecinos
iterativamente. Un paso del proceso se describe por la siguiente fórmula:
∑
−+=
n
xxxx
oldoldioldnew 0
00 λ 10 << λ
El alisado Laplaciano elimina el ruido de alta frecuencia rápidamente, pero
requiere muchas iteraciones para bajas frecuencias. Otro problema es que el
alisado Laplaciano no garantiza un mejoramiento en la calidad de los elementos,
especialmente cuando los ángulos de los elementos son muy grandes o muy
pequeños
3.4.1.2 Alisado de Taubin
La idea del procedimiento de Taubin está basada en la analogía con el análisis de
Fourier para una y dos señales regularmente muestreadas dimensionalmente.
Recordando que cualquier función f(x) se puede escribir como una combinación
de funciones seno y coseno. Por simplicidad se asume que la función es simétrica
y solo funciones coseno son empleadas:
( ) ( ) ( )∫∞
∞−
= ωωω dxfxf cos⌢
(3.1)
(3.2)
49
Donde ( )ωf⌢
se denomina espectro. La observación de Taubin fue que si se
considera las funciones básicas del espectro como funciones propias de un
continuo Laplaciano, se puede también intentar construir una representación
espectral similar usando vectores propios del Laplaciano discreto. Se puede
reescribir la ecuación del alisado Laplaciano en la forma:
( ) oldoldoldnew xKIKxxx λλ +=+=
Donde K es una matriz; esta es la matriz Laplaciana discreta.
Si u1,…uN son los vectores propios de K, donde N es el número total de vértices
en la malla. Entonces se puede expandir el vector de la posición inicial del vértice
como:
∑=
=N
iii
init uax1
Los valores propios de los vectores propios corresponden a las frecuencias en el
caso continuo. El alisado corresponde a eliminar los términos en la suma para los
cuales los valores propios λ, son grandes.
Taubin muestra que aplicando dos pasos del alisado Laplaciano, uno con
coeficiente positivo y otro con coeficiente negativo, repetidamente, se tiene el
efecto deseado. Los coeficientes µ y λ para pasos secuenciales deben ser
reemplazados por:
frecuenciadeCorte=+µλ11
Donde λ es positivo y µ es negativo.
El número de pasos se escoge empíricamente como en el caso del alisado
Laplaciano.
3.4.2 LIMPIEZA
Al igual que el suavizado, existe una gran variedad de métodos actualmente
empleados para mejorar la calidad de la malla al realizar cambios locales a las
(3.3)
(3.4)
(3.5)
50
conectividades de los elementos. Los métodos de limpieza generalmente aplican
algún criterio que debe cumplirse para realizar una operación local. Los criterios
en general pueden ser definidos como:
1. Mejoramiento de la forma, y
2. Mejoramiento topológico
Además, las operaciones de limpieza generalmente no son hechas
individualmente, sino que son usadas junto con el suavizado [16].
3.4.2.1 Mejoramiento de la forma
Para mallados triangulares, a menudo se realizan simples cambios de las
diagonales. Para cada lado en la triangulación se hace un chequeo para
determinar en que posición el lado mejoraría la forma de los dos triángulos
adyacentes. El criterio de Delaunay puede ser también usado para determinar la
posición de un lado.
3.4.2.2 Mejoramiento topológico
Un método común para mejorar el mallado es intentar optimizar el número de
bordes compartidos por un único nodo. Esto se conoce como la valencia del nodo
o el grado. Se asume que la forma de los elementos locales será mejorada. Para
una malla triangular debería haber óptimamente 6 lados junto a cada nodo. Donde
quiera que haya un nodo que no tenga una valencia ideal, la calidad de los
elementos del alrededor también será menor que la óptima. Realizando
transformaciones locales a los elementos se puede mejorar la topología y de esa
forma la calidad de los elementos. Varios métodos se han propuesto para mejorar
la valencia de los nodos para mallados triangulares y cuadriláteros.
51
3.4.3 REFINAMIENTO
Los procedimientos de refinamiento de elementos son numerosos. El refinamiento
se define como cualquier operación realizada a la malla que efectivamente
reduzca el tamaño de los elementos localmente. La reducción en el tamaño puede
ser requerido para captar un fenómeno físico local, o puede ser hecho
simplemente para mejorar la calidad de los elementos localmente [16].
Algunos métodos de refinamiento pueden ser considerados como algoritmos de
generación de mallas. Un procedimiento de refinamiento se aplica hasta que la
densidad nodal se haya alcanzado. Frecuentemente, los algoritmos de
refinamiento se usan como parte de un proceso de solución adaptativo, donde los
resultados de una solución anterior proveen un criterio para el refinamiento de la
malla. Dos de los principales métodos para refinamiento triangular y tetraédrico
incluyen:
1. Inserción de puntos
2. Plantillas (templates)
3.4.3.1 Inserción de puntos
Un enfoque simple de refinamiento consiste en insertar un único nodo en el
centroide de un elemento existente, dividiendo el triángulo en tres o un tetraedro
en cuatro. Este método generalmente no proporciona elementos de buena
calidad, particularmente después de varias iteraciones. Para mejorar el
procedimiento, una triangulación de Delaunay puede ser usada para borrar los
triángulos localmente y conectar el nodo a la triangulación manteniendo el criterio
de Delaunay.
52
Figura 3.4 Ejemplo de un refinamiento de Delaunay
3.4.3.2 Plantillas
Una plantilla se refiere a una descomposición específica del triángulo. Un ejemplo
es descomponer un único triángulo en cuatro triángulos similares insertando un
nuevo nodo en cada uno de sus lados como se muestra en la figura 3.5. Para
mantener una malla conforme, también se puede definir plantillas adicionales
basadas en el número de lados que han sido divididos.
Figura 3.5 Ejemplo de un refinamiento local triangular
53
CAPÍTULO 4.
PROGRAMACIÓN
Una vez que se ha completado el estudio de los elementos para el análisis de
placas y las técnicas de generación de mallas, se aplicará esta teoría en el
desarrollo de un programa que facilite la resolución rápida de problemas.
4.1 DESCRIPCIÓN DEL PROGRAMA
El programa deberá cumplir ciertos requisitos que garanticen que será de utilidad
en las áreas educativa y de investigación:
• Interfaz gráfica de uso intuitivo
• Capacidad de resolver problemas formados por elementos tipo placa y
elementos tipo 2D (placas planas)
• Posibilidad de mallado estructurado y no estructurado en dos dimensiones
• Presentación de resultados en forma gráfica y numérica
4.1.1 SELECCIÓN DEL MÉTODO DE MALLADO
Un tema muy importante a decidir es el método de mallado que se empleará.
Garantizar una buena calidad de la malla significa asegurar que los resultados no
difieran demasiado de los valores reales. También debe tomarse en cuenta el
tiempo requerido por el mallador automático para generar la malla.
En el capítulo 3 se ha visto que los métodos de mallado no estructurado más
populares son el método de triangulación de Delaunay y el método del avance
frontal. Este último es el que se usará para el mallado en el programa. La razón
fundamental consiste en que el método de triangulación de Delaunay no calcula la
distribución de los nodos sobre una superficie, solo su conexión, mientras que el
método del avance frontal realiza ambos pasos. Si se escogiera el primer método
54
sería necesario idear un algoritmo adicional para la colocación de los nodos, lo
que incrementaría la complejidad del programa.
Para mejorar la calidad de la malla se ha optado por aplicar el alisado Laplaciano,
mismo que se ha escogido por ser un método bastante simple de implementar.
Adicionalmente se utilizarán las optimizaciones topológicas descritas en [1].
Cuando la geometría del problema sea rectangular se podrá generar una malla
estructurada, tanto de elementos rectangulares como triangulares. En este caso el
tiempo requerido se reduce considerablemente.
4.1.2 SELECCIÓN DE LOS ELEMENTOS
En el capítulo 2 se han deducido las matrices de rigidez y esfuerzo para 4
elementos. Como el programa podrá resolver problemas con elementos tipo placa
y tipo 2D, se usarán los 4 elementos. Los elementos rectangulares solo se podrán
utilizar cuando la malla generada sea estructurada y los elementos triangulares se
podrán emplear siempre.
4.1.3 OPTIMIZACIÓN DE LA MEMORIA NECESARIA
Los requerimientos de memoria, sobre todo para almacenar la matriz de rigidez
del sistema a medida que se ensambla, pueden crecer desmesuradamente
dependiendo del número de elementos generados por el mallador.
Tabla 4.1 Memoria ocupada por la matriz de rigidez del sistema
Número de elementos
Número de nodos
Orden de la matriz de rigidez del sistema (3 GDL por nodo)
Memoria ocupada (8 Bytes por elemento)
100 121 M363x363 1.01 MB 400 441 M1323x1323 13.35 MB 900 961 M2883x2883 63.41 MB 1600 1681 M5043x5043 194.03 MB
55
Como se ve en la tabla 4.1, la memoria necesaria que debe tener el equipo de
computación donde se ejecute el programa se eleva muy rápidamente. Por esto
es necesario emplear algoritmos de almacenamiento que permitan reducir la
cantidad de memoria demandada durante el ensamblaje de la matriz del sistema
principalmente. Para ello se parte de la característica de esta matriz de ser muy
dispersa, es decir que la mayor parte de sus elementos tendrán el valor de cero.
Las matrices dispersas pueden ser almacenadas en memoria de diferente manera
que las matrices llenas: solamente se guardan los elementos distintos de cero,
junto con la posición que ocupan en la matriz. Además se debe recurrir a un
algoritmo de numeración nodal eficiente, esto determina el ancho de banda de la
matriz de rigidez global. Para el caso de placas con 3 grados de libertad (GDL)
por nodo y conociendo que la matriz de rigidez es simétrica, el ancho de banda,
BW, puede ser calculado así:
( )13 += dBW
Donde d es la mayor diferencia en la numeración nodal para todos los elementos
del ensamblado. La ecuación anterior indica que d debe ser pequeño para
obtener el ancho de banda óptimo. De manera que el ancho de banda más
pequeño se puede obtener simplemente al numerar los nodos a lo largo de la
dimensión más corta de la placa. Al reducir el ancho de banda, los requerimientos
de almacenamiento en la computadora y el tiempo de solución de la matriz del
sistema son igualmente reducidos.
4.2 LENGUAJE
Para facilitar la elaboración del programa, primero se desarrollará una versión en
Matlab, y posteriormente se la desarrollará en Java. Comenzar por Matlab tiene la
ventaja de que su lenguaje de programación es sencillo y de alto nivel, que
incorpora una gran cantidad de funciones de cálculo predefinidas y está orientado
a la manipulación de matrices. Además Matlab permite crear una relativamente
amplia variedad de gráficos sin mucha dificultad. Esto significa que se puede
completar el programa en Matlab más rápido que en otros lenguajes
convencionales. Estas ventajas se aprovecharán para, en una primera etapa,
(4.1)
56
concentrarse en la definición de algoritmos, subrutinas y formas de interacción
que resultarían en un código óptimo. Todo esto dejando de lado las dificultades
que existen en el estudio de la sintaxis de un lenguaje de programación de más
bajo nivel, el desarrollo de algoritmos para la resolución de sistemas de
ecuaciones, esquemas de presentación gráfica de resultados, etc.
Después se pasarán los algoritmos desarrollados a Java. Java tiene la ventaja de
que un programa codificado en este lenguaje puede ejecutarse en cualquier
plataforma que disponga de la máquina virtual Java, como son Windows, Unix -
Linux y Mac OS, entre los principales. Además, el uso de Java permite codificar
programas que pueden ser distribuidos y ejecutados sin tener que recurrir a la
adquisición de licencias de uso de software comercial costoso tal como es el caso
de Matlab.
El equipo de computación que se empleará para el desarrollo es un clon de las
siguientes características: procesador AMD Athlon 64 3000+ de 1800 MHz, 1 GB
de memoria RAM, disco duro SATA de 160 GB, monitor de pantalla plana de 17
pulg. y una unidad óptica DVD±R/RW. Los programas que se utilizarán son:
plataforma Microsoft Windows XP Profesional, Matlab 7.0, J2SE Development Kit
5.0, Java 3D 1.4.0, NetBeans 5.0 y Algor 19.
4.3 PROGRAMACIÓN EN MATLAB
El programa se lo desarrollará en la versión disponible más reciente de Matlab,
versión 7.0. Esta versión de Matlab cuenta con nuevos componentes que
permiten crear interfaces más elaboradas, además de mejoras en el lenguaje de
programación y nuevas funciones. Lamentablemente esto no permite que el
programa pueda ejecutarse en versiones anteriores de Matlab.
El programa consta de 36 subprogramas de los cuales 10 tienen interfaz gráfica.
Todos los subprogramas se explican a continuación.
57
4.3.1 SUBPROGRAMAS CON INTERFAZ GRÁFICA
Placa.- Es la aplicación principal. Al ejecutarse el programa se muestra la ventana
gráfica desde donde se puede llamar a otros subprogramas. Al inicio las opciones
disponibles son escoger el tipo de elemento (por defecto es placa) y editar la
geometría. Una vez editada la geometría se activan las opciones para configurar
propiedades y el mallado. Al realizar la malla se activan las opciones para
ingresar los apoyos y las cargas. Para que se active la opción Analizar se deben
ingresar al menos los apoyos. Al ejecutar Analizar se calculan desplazamientos,
cargas, deformaciones, esfuerzos y reacciones y se muestra la gráfica de la
geometría deformada del problema. Existe la opción de guardar en disco el
archivo del problema o solo los resultados en un archivo de texto. En la figura 4.1
se muestra la ventana una vez terminada.
Figura 4.1 Ventana principal del programa
58
Acerca.- Muestra una ventana informativa sobre la versión del programa
(actualmente versión 1.0), en la figura 4.2 se muestra esta ventana.
Figura 4.2 Ventana Acerca
Analisis.- Despliega una ventana gráfica para iniciar el análisis del problema. A
medida que transcurre el análisis muestra las etapas de cálculo y el tiempo
empleado por cada una.
Dependiendo del tipo de elemento (placa o 2D) se llama a las subrutinas
calcMatRigSis o calcMatRigSis2D que ensamblan la matriz de rigidez del
sistema a partir de los datos recolectados por la aplicación principal. Luego se
aplican las condiciones de frontera, se calculan los desplazamientos y las cargas,
y a partir de estos resultados se calculan deformaciones, esfuerzos y reacciones
con las subrutinas calcMatDefEsf o calcMatDefEsf2D . En la figura 4.3 se
muestra esta ventana.
59
Figura 4.3 Ventana Análisis
Apoyos.- Permite imponer las restricciones de manera gráfica mediante la
selección de uno o más nodos. Se pueden restringir los desplazamientos en los
ejes X, Y y Z y las rotaciones alrededor de los ejes X y Y, dependiendo del tipo de
elemento. En la figura 4.4 se muestra esta ventana.
Figura 4.4 Ventana Apoyos
60
Cargas.- Al igual que Apoyos , muestra una ventana para imponer cargas sobre la
geometría del problema. Para elementos tipo placa, fuerzas actuantes en la
dirección del eje Z y momentos alrededor de X y Y. Para elementos tipo 2D,
fuerzas actuantes en la dirección de los ejes X y Y. En la figura 4.5 se muestra
esta ventana.
Figura 4.5 Ventana Cargas
Geometria.- Ventana para la edición de la geometría del problema, se pueden
dibujar líneas, rectángulos, arcos, círculos y elipses. Para el caso de líneas se
guardan sus coordenadas inicial y final, además se considera al rectángulo como
el dibujo de 4 líneas sucesivas. Para el caso de arcos, círculos y elipses se
guarda su radio mayor y menor, su centro y los ángulos inicial y final.
En esta misma ventana se pueden crear puntos de refinamiento de la malla.
Además, mediante la creación de varias partes de geometría se puede tener
regiones que posean diferentes propiedades y espesor. En la figura 4.6 se
muestra esta ventana.
61
Figura 4.6 Ventana Geometría
Malla.- Permite configurar los parámetros del mallado como forma de los
elementos, tamaño de los elementos y el factor de cercanía. Si se seleccionó la
opción Estructurado para modelos rectangulares, y la geometría es rectangular, el
mallado será estructurado, caso contrario el mallado será no estructurado.
Siempre que sea posible debería optarse por el mallado estructurado,
especialmente si el número de elementos a generar es alto (más de 1000), ya que
realizar el mallado no estructurado puede consumir un tiempo considerable. En la
figura 4.7 se muestra esta ventana.
El factor de cercanía es un parámetro que determina como se forman los
elementos de la malla. Un valor de 1 o menor favorece la creación de elementos
triangulares equiláteros, pero puede conducir a errores en geometrías muy
complejas. Un valor mayor que 1 hará que los elementos tiendan a formarse
empleando los nodos de elementos ya creados, esto en cambio puede crear
mallas de pobre calidad.
62
Figura 4.7 Ventana Malla
Opciones.- Ventana donde se pueden configurar dos parámetros: si se mostrarán
los detalles del análisis y el factor de escala usado en la representación de los
resultados. Se muestra en la figura 4.8.
Figura 4.8 Ventana Opciones
Propiedades.- Permite establecer las propiedades de cada parte de la geometría.
Las propiedades son el módulo de elasticidad, módulo de Poisson y el espesor.
Existe la posibilidad de seleccionar un material de la lista desplegable que
contiene algunos materiales comúnmente usados en ingeniería. En la figura 4.9
se muestra esta ventana.
63
Figura 4.9 Ventana de propiedades
TipoElemento.- Permite escoger el tipo de elemento del que está formado la
geometría del problema, existen dos tipos de elementos: placa y 2D. Si se escoge
placa, las cargas actuantes deberán ser normales al plano de la geometría
(dirección Z) y el desplazamiento será igualmente en esa dirección. Si se escoge
2D las cargas deberán estar en el plano de la geometría y el desplazamiento solo
se podrá dar en este mismo plano. Se puede revisar la teoría de placas empleada
en el capítulo 2. Esta ventana se muestra en la figura 4.10.
Figura 4.10 Ventana Tipo de Elemento
4.3.2 SUBPROGRAMAS SIN INTERFAZ GRÁFICA
areaPoligono.- Calcula el área de un polígono, si las coordenadas del polígono
están dadas en sentido antihorario el área es positiva, caso contrario el área es
negativa.
64
El área se calcula dividiendo el polígono de n lados en n - 2 triángulos, entonces
el área es la suma de las áreas triangulares definidas por:
( ) 231121231 yxyxA ⋅−⋅=
Donde x1y1, x2y2, x3y3 son las coordenadas del triángulo. El código es el siguiente:
for i = 2:size(P,1)-1 A0 = ((P(i,3)-P11)*(P12-P(i,2)) - (P11-P(i,1))* (P(i,4)-P12)) / 2; % Suma de las areas triangulares A = A + A0; end
buscarVecMat.- Busca un vector dentro de una matriz y devuelve la posición del
vector, si no se encuentra devuelve cero. El código es:
% Para que el vector sea igual deben coincidir los dos elementos del % vector con los elementos de la fila de la matriz MLogica = all([Mat(:,1)==Vec(1) Mat(:,2)==Vec(2)],2 ); Ind = find(MLogica); % Un indice 0 indica que no existe ese vector dentr o de la matriz if isempty(Ind), Ind = 0; end
calcMatDefEsf.- Calcula las deformaciones y esfuerzos del elemento tipo placa.
Para elementos rectangulares llama a la subrutina matDefEsfRec y para
elementos triangulares llama a matDefEsfTriDKT repetidamente para cada
elemento. Esto se puede ver en el siguiente algoritmo ejecutado iterativamente
para cada elemento rectangular:
% Para los elementos rectangulares es suficiente co n conocer su % ancho y su longitud a = Nxy(Elem(i,2),1)-Nxy(Elem(i,1),1); b = Nxy(Elem(i,3),2)-Nxy(Elem(i,2),2); % Desplazamientos del elemento, son 4 nodos con 3 d esplazamientos % cada uno: Uz, Rx y Ry for j = 1:4 d(3*j-2) = Despl(3*Elem(i,j)-2); d(3*j-1) = Despl(3*Elem(i,j)-1); d(3*j ) = Despl(3*Elem(i,j) ); end % Calculo de las deformaciones, esfuerzos y momento s flectores Se = matDefEsfRec(E(Elem(i,5)),v(Elem(i,5)),h(Elem( i,5)),a,b,d); % V contiene las deformaciones, esfuerzos y momento s flectores, % ademas un contador para poder calcular el promedi o de los valores
(4.2)
65
for j = 1:4 V(3*Elem(i,j)-2,:) = V(3*Elem(i,j)-2,:)+[Se(3*j -2,:) 1]; V(3*Elem(i,j)-1,:) = V(3*Elem(i,j)-1,:)+[Se(3*j -1,:) 1]; V(3*Elem(i,j) ,:) = V(3*Elem(i,j) ,:)+[Se(3*j ,:) 1]; end
Este código promedia los valores de esfuerzos, deformaciones y momentos
flectores en cada nodo, según lo recomendado en la sección 2.1.2.6.
calcMatDefEsf2D.- Calcula las deformaciones y esfuerzos del elemento tipo 2D.
Para elementos rectangulares llama a la subrutina matDefEsfRec2D y para
elementos triangulares llama a matDefEsfTriCST repetidamente para cada
elemento. El algoritmo es similar que para calcMatDefEsf .
calcMatRigSis.- Ensambla la matriz de rigidez del sistema para elementos tipo
placa. Iterativamente, para cada elemento, determina su conectividad, calcula su
matriz de rigidez y la agrega a la matriz de rigidez del sistema. El cálculo de la
matriz de rigidez de cada elemento se hace con las subrutinas matRigRec y
matRigTriDKT , dependiendo de si los elementos son rectangulares o
triangulares respectivamente. El siguiente fragmento de código se ejecuta
repetidamente para cada elemento rectangular:
% Conectividad del elemento, 4 nodos con 3 GDL por nodo Ce = [3*Elem(i,1)-2 3*Elem(i,1)-1 3*Elem(i,1) ... 3*Elem(i,2)-2 3*Elem(i,2)-1 3*Elem(i,2) ... 3*Elem(i,3)-2 3*Elem(i,3)-1 3*Elem(i,3) ... 3*Elem(i,4)-2 3*Elem(i,4)-1 3*Elem(i,4)]; % Para los elementos rectangulares es suficiente co n conocer su % ancho y su longitud a = Nxy(Elem(i,2),1)-Nxy(Elem(i,1),1); b = Nxy(Elem(i,3),2)-Nxy(Elem(i,2),2); % Calculo de la matriz de rigidez del elemento, de orden 12x12 Ke = matRigRec(E(Elem(i,5)),v(Elem(i,5)),h(Elem(i,5 )),a,b); % Se agrega la matriz de rigidez del elemento a la del sistema % Por ahora solo se agrega la diagonal superior, la otra parte es % simetrica y se agrega al final for j = 1:12 for k = j:12 K(Ce(j),Ce(k)) = K(Ce(j),Ce(k))+Ke(j,k); end end
66
calcMatRigSis2D.- Igual que calcMatRigSis para elementos tipo 2D, las
subrutinas empleadas para el cálculo de la matriz de rigidez de cada elemento
son matRigRec2D y matRigTriCST , dependiendo de si los elementos son
rectangulares o triangulares respectivamente.
calcNodosFrente.- Subrutina llamada por calcNodosFrenteArea para calcular
las coordenadas de los nodos del contorno y el frente de generación inicial para
cada parte de geometría.
calcNodosFrenteArea.- Subrutina llamada por genMalla2D para calcular las
coordenadas de los nodos del contorno, el frente de generación inicial y el área de
la geometría.
genMalla2D.- Realiza el mallado no estructurado del problema. El algoritmo de
mallado emplea el método del avance frontal, tal como se describe en [1]. A partir
de la geometría del problema y de los datos proporcionados por Malla , se van
generando elementos triangulares sobre todo el contorno hasta “llenar”
completamente la superficie del problema. El algoritmo es el siguiente:
% Generacion de los elementos triangulares while size(Fre,1) % Generar mient ras existan segmentos % 1. Eleccion del segmento del frente de genera cion xyg = Nxy(Fre(1,1:2),:); [xyr Fi] = rotarCoordXYFi(xyg); % S egemento en SCL % 2. Obtencion del valor de espaciamiento d d = xyr(2,1)-xyr(1,1); % El tamaño se ajusta de manera que el elemento resultante no sea % demasiado grande ni demasiado pequeño if Tam < 0.55*d d = 0.55*d; elseif Tam>2*d && (Fre(1,1)>NxyC && Fre(1,2)>Nx yC) d = 1.8*d; end % 3. Determinacion del punto C1 h = 0.866025404*d; % A ltura del triangulo C = [(xyr(2,1)+xyr(1,1))/2 xyr(1,2)+h]; % C 1 esta en SCL % 4. Determinacion de los puntos C2, C3, C4, C5 y C6 % Estos puntos estan equiespaciados en el segme nto que une C1 y M for i = 1:5 C(end+1,:) = [C(1,1) C(1,2)-i*h/6];
67
end C = rotarCoordXYFi(C,-Fi); % P aso de C de SCL a SCG % 5. Nodos pertenecientes al frente de generaci on r = 4*Tam; % Radio en el q ue se buscaran nodos d = sqrt((Nxy(:,1)-C(1,1)).^2+(Nxy(:,2)-C(1,2)) .^2); % Distancia % De estos nodos se escogen los que esten situa dos sobre el semiplano % AB y que contiene a C1 Nf = Nxy(d<r,:); d = d(d<r); NfAB = rotarCoordXY([xyg ; Nf]); Nf = [d Nf]; NfAB = Nf(NfAB(3:end,2)>0.1*tol,:); % 6. Colocacion de C1...C6 en la lista d1 = P(5)*h; % Penalizacion de una distancia d1 Nf = [NfAB ; [[1 7/6 8/6 9/6 10/6 11/6]'.*d1 C] ]; % Colocacion de acuerdo a su distancia a C1 Nf = sortrows(Nf); Nf(:,1) = []; NfAB(:,1) = [] ; % 7. Determinacion del punto de conexion aj % Para ahorrar tiempo se comprobara solo con lo s frentes de la % parte actual f = Fre(Fre(:,3)==Fre(1,3),1:2); % Se toma el primer nodo de la lista Nf que cum ple con las % condiciones 7a y 7b for i = 1:size(Nf,1) % 7a. Sin nodos dentro del triangulo ABaj e xcepto Ci Nfs = NfAB(NfAB(:,1)~=Nf(i,1) | NfAB(:,2)~= Nf(i,2),:); xy1 = rotarCoordXY([xyg ; Nfs]); % Para el lado 12 xy2 = rotarCoordXY([xyg(2,:) ; Nf(i,:) ; Nf s]); % Para el lado 23 xy3 = rotarCoordXY([Nf(i,:) ; xyg(1,:) ; Nf s]); % Para el lado 31 % Como xy1, xy2 y xy3 estan en SCL, entonce s para determinar si % los nodos estan dentro del triangulo bast a con comprobar si sus % valores son mayores a cero if all(xy1(3:end,2)<-tol | xy2(3:end,2)<-to l | xy3(3:end,2)<-tol) % 7b. El segmento ajM no corta ninguna cara existente en el % frente de generacion % Por facilidad se rotan las coordenada s de manera que ajM % coincide con el eje de las abscisas m = rotarCoordXY([sum(xyg)/2 ; Nf(i,:) ; Nxy(f(:,1),:) ; Nxy(f (:,2),:)]); s = size(f,1); % No se incluyen los segmentos ajM y AB en el chequeo p = m(4:s+2,:); p(:,3:4) = m(s+4:end,:) ; % Se eliminan las caras que incluyen a aj % Si no queda ninguna cara se encontro aj p = p((p(:,1)~=m(2,1) | p(:,2)~=m(2,2)) & (p(:,3)~=m(2,1) | p(:,4)~=m(2,2)) ,:); if isempty(p), aj = Nf(i,:); break, end % Se eliminan las caras que pasan por ' debajo' o por 'arriba' % del segmento ajM p = p(~((p(:,2)<-tol & p(:,4)<-tol) |
68
(p(:,2)>tol & p(:,4)>tol)),:); if isempty(p), aj = Nf(i,:); break, end % Se eliminan las caras que pasan por l a 'izquierda' o por la % 'derecha' del segmento ajM p = p(~((p(:,1)<-tol & p(:,3)<-tol) | (p(:,1)>m(2,1)+tol & p(:,3)>m(2 ,1)+tol)),:); if isempty(p), aj = Nf(i,:); break, end % Se emplea la ecuacion b = x1-(x2-x1)/ (y2-y1)*y1 % Las caras no cortan ajM si b < 0 y b > ajM p = p(:,1)-(p(:,3)-p(:,1))./(p(:,4)-p(: ,2)).*p(:,2); if all(p<-tol | p>m(2,1)+tol), aj = Nf( i,:); break, end end end % 8. Formacion y almacenamiento del nuevo eleme nto j = buscarVecMat(aj,Nxy); % Ya ex iste el nodo? if j == 0 Nxy(end+1,:) = aj; j = size(Nxy,1); % Creac ion de un nuevo nodo end Elem(end+1,:) = [Fre(1,1) Fre(1,2) j Fre(1,3)]; % Nuevo elemento % 9. Actualizacion del frente de generacion j0 = buscarVecMat([j Fre(1,1)],Fre(:,1:2)); % Y a existe la cara 31? if j0, Fre(j0,:) = []; else Fre(end+1,:) = [Fre (1,1) j Fre(1,3)]; end j0 = buscarVecMat([Fre(1,2) j],Fre(:,1:2)); % Y a existe la cara 23? if j0, Fre(j0,:) = []; else Fre(end+1,:) = [j F re(1,2) Fre(1,3)]; end Fre(1,:) = []; aj = []; % E limina la cara 12 end
Posteriormente se refina la malla en los puntos de refinamiento ingresados en
Geometria . Para esto se emplea una plantilla que divide cada triángulo que
comparte el punto de refinamiento en 4 triángulos internos [16]. Finalmente se
realizan algunas operaciones de cosmética como son la optimización topológica y
el alisado de la malla. Para este último se emplea el alisado Laplaciano [25] [28],
que corresponde a la siguiente porción de código:
% 3. Alisado de la malla for m = 1:3 for i = NxyNC E = Elem(any(Elem(:,1:3)==i,2),1:3); V = unique(E(:)); V = V(V~=i); % Alisado Laplaciano % La formula es: x0new = x0old+h*sum((xiold -x0old)/n), 0 < h < 1 % Se escoge h = 0.8 Nxy(i,1) = Nxy(i,1)+0.8*sum((Nxy(V,1)-Nxy(i ,1))./length(V)); Nxy(i,2) = Nxy(i,2)+0.8*sum((Nxy(V,2)-Nxy(i ,2))./length(V)); end end
69
Este algoritmo se ejecuta 3 veces. Este valor se escogió luego de realizar varias
mallas con diferente número de iteraciones. Un valor mayor aumenta el tiempo de
generación de la malla excesivamente sin una mejora notable de su calidad.
graficar2D.- Grafica los desplazamientos para elementos tipo 2D, es decir, en el
mismo plano de la geometría.
graficar3D.- Grafica los resultados en el espacio tridimensional, esto hace que
Matlab automáticamente asigne una coloración de acuerdo al valor en Z del
gráfico, haciendo más fácil la interpretación del resultado del problema.
graficarGeometria.- Grafica la geometría del problema, es decir se grafican las
líneas, arcos, elipses y puntos de refinamiento que definen la geometría.
graficarModelo.- Grafica la geometría ya mallada del problema, también es capaz
de representar los apoyos y cargas que se hayan ingresado.
matDefEsfRec, matDefEsfRec2D.- Calculan las deformaciones y esfuerzos para el
elemento rectangular tipo placa y 2D. Los cálculos se realizan a partir de las
propiedades del elemento y de las coordenadas que definen su forma. Las
matrices de esfuerzos para implementar los algoritmos se obtuvieron de varias
fuentes [22] [23].
matDefEsfTriCST, matDefEsfTriDKT.- Calculan las deformaciones y esfuerzos para
el elemento triangular tipo placa y 2D respectivamente [2] [11] [22] [23].
matRigRec, matRigRec2D.- Calculan la matriz de rigidez para el elemento
rectangular tipo placa y 2D. Los cálculos se realizan a partir de las propiedades
del elemento y de las coordenadas que definen su forma. Las matrices de rigidez
de las subrutinas y otras matrices usadas en el cálculo se obtuvieron de varias
fuentes [22] [23].
70
matRigTriDKT, matRigTriCST.- Calculan la matriz de rigidez para el elemento
triangular tipo placa y 2D respectivamente [11] [22] [23].
modeloRec.- Comprueba si la geometría es rectangular, esto puede hacerse para
decidir si es factible realizar el mallado estructurado o si necesariamente el
mallado será no estructurado.
polarRec.- Pasa de coordenadas polares a rectangulares. El algoritmo es muy
simple como se ve a continuación:
X = Rx.*cos(Ang); Y = Ry.*sin(Ang);
rotarCoordXY.- Rota la matriz de coordenadas dada al sistema local definido por
las dos primeras coordenadas de la matriz y ajusta el origen a la primera fila de la
matriz.
rotarCoordXYFi.- Rota la matriz de coordenadas dada al sistema local definido por
un ángulo Fi.
verificar.- Comprueba si el valor ingresado es un número válido. Este método es
llamado por todas las subrutinas con interfaz gráfica en donde se ingrese algún
valor numérico.
4.4 PROGRAMACIÓN EN JAVA
La versión del programa en Java se la desarrolla empleando J2SE Development
Kit 5.0, Java 3D 1.4.0 y el entorno de programación NetBeans 5.0.
El programa consta de 51 ficheros fuente, 20 ficheros de ayuda de extensión htm
y 132 imágenes de extensión png distribuidos en 6 packages. Además emplea las
librerías de distribución gratuita llamadas colt versión 1.2.0 y mtj versión 0.9.6.
Ambas librerías permiten operar con matrices llenas y dispersas. Finalmente el
programa gratuito IzPack versión 3.7.2 crea el instalador para el programa.
71
4.4.1 SUBPROGRAMAS CON INTERFAZ GRÁFICA
Todas las clases aquí listadas pertenecen al package edu.epn.gui . NetBeans
genera ficheros de extensión form, lo que facilita la construcción de interfaces
gráficas con el uso del ratón.
Placa.- Clase principal, muestra la ventana inicial desde donde se puede invocar
otras ventanas. A diferencia de su contraparte en Matlab, esta versión en Java
ofrece muchas más opciones entre las que se destacan la capacidad de poder
dibujar la geometría con el puntero del ratón, asignar apoyos y cargas
directamente a los nodos seleccionados en el gráfico en 3 dimensiones y la
capacidad de poder seleccionar el sistema de unidades. Para más información se
puede consultar la ayuda integrada en el programa. Esta ventana se muestra en
la figura 4.11.
Figura 4.11 Ventana Placa
72
Acerca.- Esta clase simplemente muestra la versión del programa y los nombres
del autor y director del programa. En la figura 4.12 se muestra esta ventana.
Figura 4.12 Ventana Acerca
Analisis.- Clase que muestra las etapas del proceso de resolución y el tiempo
transcurrido. El código que ejecuta el análisis es el siguiente:
private class InicioAnalisis extends Thread { private Indicador ind; public InicioAnalisis(Indicador ind) { super("análisis"); this.ind = ind; } public void run() { try { new AnalisisModelo(ind, parent.config.r esolutor). iniciarAnalisis(parent.modelo); } catch (IterativeSolverNotConvergedExcepti on e) { btnAceptar.setEnabled(true); lblAnalizando.setText("Terminado:"); Utilidades.mostrarErrorInfo(Analisis.th is, "El resolutor no converge"); return; } tiempoTrans.detener(); lblAnalizando.setText("Terminado:"); btnAceptar.setEnabled(true); } }
73
En el código se ve que esta subclase deriva de Thread y, por lo tanto, es un hilo
de ejecución independiente. Además si el resolutor es iterativo puede no
converger a la solución del problema, en este caso se mostrará un mensaje de
error. En la figura 4.13 se muestra la ventana de análisis.
Figura 4.13 Ventana Análisis
Apoyos.- Clase que permite asignar apoyos a uno o más nodos seleccionados, las
restricciones disponibles dependen del tipo de elemento. En la figura 4.14 se
muestra esta ventana.
Figura 4.14 Ventana Apoyos
74
BarrasHerramientas.- Clase que permite escoger las barras de herramientas
visibles y el tamaño de los íconos. Esta ventana se muestra en la figura 4.15.
Figura 4.15 Ventana Barras de herramientas
CambiosAyuda.- Clase que muestra los cambios sufridos por el programa desde la
última versión. En la figura 4.16 se muestra esta ventana.
Figura 4.16 Ventana Cambios
75
Cargas.- Clase que permite asignar fuerzas y momentos a uno o más nodos
seleccionados, los campos disponibles dependen del tipo de elemento. En la
figura 4.17 se muestra esta ventana.
Figura 4.17 Ventana Cargas
ContenidoAyuda.- Clase que muestra la ayuda del programa, se divide en 4
categorías: Introducción, General, Resolución del modelo y Características
adicionales. En la figura 4.18 se muestra esta ventana.
Figura 4.18 Ventana Ayuda
76
GraficoPlano.- Clase que muestra el gráfico de un plano de corte hecho a los
resultados. En la figura 4.19 se muestra esta ventana.
Figura 4.19 Ventana GraficoPlano
InfoNodo.- Clase que muestra el número y coordenadas del nodo seleccionado.
Esta ventana se muestra en la figura 4.20.
Figura 4.20 Ventana InfoNodo
LibreriaMateriales.- Clase que permite ver, crear y borrar materiales, que luego
pueden ser empleados para definir las propiedades de la geometría. Esta ventana
se muestra en la figura 4.21.
77
Figura 4.21 Ventana Librería de materiales
Malla.- Clase que permite configurar los parámetros de la malla. Estos parámetros
son los mismos que para la versión en Matlab, salvo la opción de poder definir el
número de iteraciones del alisado Laplaciano. Esta opción, como se indicó al
exponer la subrutina genMalla2D para la versión hecha en Matlab, no contribuye
a mejorar la calidad de la malla de manera evidente para un valor mayor a 10.
Esta ventana se muestra en la figura 4.22.
Figura 4.22 Ventana Malla
78
Opciones.- Clase que permite configurar varias opciones del programa que se
guardarán por medio de la clase Configuracion que se describe más adelante.
Las opciones están agrupadas en 4 pestañas: General, Gráficos, Análisis y
Resultados. Para más información consultar la ayuda del programa. En la figura
4.23 se muestra esta ventana.
Figura 4.23 Ventana Opciones
Propiedades.- Clase que permite definir las propiedades de una o más partes de
geometría. Las propiedades son el módulo de elasticidad y de Poisson del
material, y el espesor. También está la posibilidad de escoger un material de una
de las dos librerías de materiales disponibles. En la figura 4.24 se muestra esta
ventana.
79
Figura 4.24 Ventana Propiedades
Unidades.- Clase que permite definir las unidades del problema. Por defecto las
unidades son las correspondientes al SI. En la figura 4.25 se muestra esta
ventana.
Figura 4.25 Ventana Unidades
4.4.2 SUBPROGRAMAS SIN INTERFAZ GRÁFICA
En el anexo F se puede consultar los diagramas de flujo de los métodos más
importantes.
80
ArbolContenido.- Clase requerida por ContenidoAyuda (ver más adelante) para
mostrar el árbol de ayuda del programa. El método más útil es mostrarPagina ,
que recibe como parámetros un objeto de la clase JEditorPane y una cadena
de texto que indica la dirección de la página en formato htm a mostrar. Este
método se muestra a continuación:
public void mostrarPagina(JEditorPane editor, Strin g pagina) { try { editor.setPage(getClass().getResource(pagin a)); } catch (java.io.IOException e1) { JOptionPane.showMessageDialog(editor, "Error al cargar la página de ayuda ", "Error", JOptionPane.ERROR_MESSAGE) ; } }
ArbolModelo.- Clase requerida por ArbolGrafico (ver más adelante) para
mostrar el árbol del modelo, conteniendo las categorías Unidades, Geometría,
Malla, Apoyos y Cargas.
Configuracion.- Clase que guarda varios parámetros del funcionamiento del
programa, por ejemplo: posición de la ventana principal, estilo visual, color del
modelo, apoyos, cargas y grilla, tamaño de íconos, etc. El método
guardarOpciones guarda la instancia de esta clase al archivo placa.opc
ubicado en el mismo directorio del programa.
Materiales.- Clase abstracta de la que derivan MaterialesLista y
MaterialesUsuario y que define dos métodos adicionales que deben ser
redefinidos: getModuloElasticidad y getModuloPoisson .
MaterialesLista.- Clase que deriva de Materiales y que contiene una lista de
nombres de materiales y otra lista de sus propiedades. Redefine los métodos de
Materiales . Por ejemplo, el código del método que permite conocer el módulo
de Poisson es:
public double getModuloPoisson(int idx) { return PROP_MATERIALES[idx][1]; }
81
MaterialesUsuario.- Clase que deriva de Materiales y que permite crear una lista
de materiales particular, para ello dispone de varios métodos. Por ejemplo, para
agregar un nuevo material, dispone del método insertarMaterial , su código
es:
public void insertarMaterial(int idx, String nombre , double elast, double poiss) { Material m = new Material(nombre, elast, poiss) ; if (! LIST_MATERIALES.contains(m)) { LIST_MATERIALES.add(idx, m); } }
Modelo.- Clase que guarda los datos del problema que se analizará, como son su
geometría, propiedades, tipo y forma de elemento, coordenadas nodales, apoyos,
cargas, etc. Un objeto de esta clase se puede leer y escribir a un archivo de
extensión p2d.
ModeloFilter.- Clase que deriva de FileFilter y que es usada por
fchAbrirGuardar de la clase Placa para filtrar el tipo de archivos que se
pueden abrir y guardar a solo aquellos con extensión p2d.
NodoModelo.- Clase que es empleada por ArbolModelo para guardar el nombre,
posición y tipo de cada uno de los nodos agregados al árbol.
RenderAtrib.- Clase que deriva de RenderingAttributes y que es usada por
algunas clases del package edu.epn.graficos3D para habilitar la escritura de
la visibilidad. Su único constructor muestra como se realiza este proceso:
public RenderAtrib() { setCapability(RenderingAttributes.ALLOW_VISIBLE _WRITE); }
RendererModelo.- Clase que deriva de DefaultTreeCellRenderer y que es
usada por ArbolGrafico para mostrar diferentes íconos dependientes del tipo
de nodo del árbol.
82
RendererTema.- Igual que la clase anterior, usado por ContenidoAyuda .
SimpleMouse.- Clase que controla la respuesta de los gráficos en 3 dimensiones al
moverse el puntero del ratón manteniendo presionado un botón del mismo.
Existen 3 tipos de comportamientos: traslación, rotación y acercamiento o
alejamiento. Por ejemplo, si está activado el comportamiento de traslación, el
siguiente fragmento de código modifica la posición del gráfico:
if ((Math.abs(dy) > 50) || (Math.abs(dx) > 50)) bre ak; traslacion.x = dx * X_FACTOR_TRASLACION; traslacion.y = -dy * Y_FACTOR_TRASLACION; transformX.set(traslacion); currXform.mul(transformX, currXform);
Apoyo3D.- Clase que crea el gráfico de un apoyo tridimensional. Su constructor
toma como parámetros las posiciones en X, Y, y Z del nuevo apoyo.
BarraColor3D.- Clase que crea el gráfico de una barra de colores. Esta clase es
usada por PanelResultados (ver más adelante).
Carga3D.- Clase que crea el gráfico de una carga tridimensional. Su constructor
toma como parámetros las posiciones en X, Y y Z de la nueva carga.
Ejes3D.- Clase que crea el gráfico de ejes en 3 dimensiones. Esta clase es usada
por PanelGrafico3D (ver más adelante).
Info3D.- Clase que muestra una cadena de texto dentro de un cuadro en el
espacio tridimensional. Usada por PanelResultados .
Linea3D.- Permite crear una línea en el espacio tridimensional. Esta clase es
usada por PanelResultados para definir un plano de corte.
Modelo3D.- Clase que crea una geometría en 3 dimensiones, dispone de métodos
para cambiar su geometría, apariencia, color, escala, rotación y traslación. Por
ejemplo el código para cambiar la escala es el siguiente:
83
public void cambiarEscala(double escala) { mouseModelo3D.escala = escala; Transform3D transform = new Transform3D(); tgRotacion.getTransform(transform); transform.setScale(escala); tgRotacion.setTransform(transform); }
Seleccion3D.- Clase que crea un rectángulo o un círculo de selección en el espacio
tridimensional. Esta clase es usada por PanelModelo que más adelante se
explica.
Texto3D.- Clase que muestra dos cadenas de texto en el espacio tridimensional.
Usada por PanelResultados .
ArbolGrafico.- Clase derivada de JTree y que crea un árbol usado por Placa en
donde se muestra información relativa al problema actual. Los datos mostrados
pueden ser: partes de geometría, figuras geométricas, apoyos, cargas, unidades y
malla.
PanelGeometria.- Clase que deriva de PanelGrafico2D , es el panel en donde se
dibuja la geometría, redefine el método paint de manera que cada vez que es
llamado este método se redibuja la geometría.
Cada vez que se ingresa una figura geométrica se llama al método
agregarGeometria de la clase Placa. Este método guarda los datos de la
figura geométrica dibujada. Para el caso de líneas se guardan sus coordenadas
inicial y final. Para el caso de rectángulos se almacena como si se tratase de 4
líneas sucesivas. Los arcos guardan las coordenadas de su centro, su radio y sus
ángulos inicial y final. Para círculos y elipses se guardan las coordenadas de su
centro y sus radios mayor y menor (iguales, si es un círculo). Adicionalmente a
estos datos, para todas las figuras geométricas se almacena la parte a la que
pertenece y el valor de un parámetro que indica si la geometría dibujada es
externa o interna.
84
PanelGrafico2D.- Clase que deriva de JPanel , es una clase base para el
tratamiento de gráficos en 2 dimensiones.
PanelGrafico3D.- Clase que deriva de PanelGrafico2D y agrega métodos
básicos para poder trabajar con gráficos en 3 dimensiones.
PanelModelo.- Clase que deriva de PanelGrafico3D , es el panel en donde se
dibuja la geometría ya mallada, implementa la selección de nodos lo que permite
asignar apoyos y cargas.
PanelResultados.- Clase que deriva de PanelGrafico3D , en este panel se
grafican los resultados. Sus dos métodos más importantes son crearModelo y
mostrarResultados , el primero crea el gráfico de la geometría desplazada
mientras que el segundo cambia la coloración de la geometría para reflejar el tipo
de resultado.
SplashImagen.- Clase que muestra una gráfica mientras se carga la aplicación
principal. Esta clase es llamada en el método main de Placa , este fragmento de
código se muestra a continuación:
// Mientras se carga la interfaz principal se muest ra el // logo del programa SplashImagen.mostrar(Placa.class.getResource( "/edu/epn/imagenes/Splash.png")); java.awt.EventQueue.invokeLater(new Runnable() { public void run() { new Placa().setVisible(true); } }); // Una vez cargada la interfaz, el logo desaparece SplashImagen.cerrar();
Utilidades.- Clase que dispone de varios métodos utilizados por las clases con
interfaz gráfica.
85
AnalisisModelo.- Clase que realiza el análisis del modelo, su constructor se detalla
a continuación:
public AnalisisModelo(Indicador ind, byte tipo, boo lean precondicionar) { this.ind = ind; this.tipo = tipo; this.precondicionar = precondicionar; }
Este constructor toma como parámetros un objeto de la clase Indicador
(descrito más adelante), un byte que indica el tipo de resolutor a emplear, que
puede ser utilizando la descomposición LU o el método iterativo de los gradientes
biconjugados, y un boolean que indica si se precondiciona el método iterativo.
Una vez creada la instancia de esta clase se debe llamar al método
iniciarAnalisis que es el que efectivamente resuelve el problema. El
siguiente fragmento de código es ejecutado para elementos tipo placa:
// Ensamblaje de la matriz de rigidez del sistema SparseDoubleMatrix2D kSist = calcMatRigSis(mod); indicar(); // Número de restricciones int r = 0; for (int i = 0; i < mod.apoyos.length; i++) { if (mod.apoyos[i][2]) r++; // Restricc ión Uz if (mod.apoyos[i][3]) r++; // Restricc ión Rx if (mod.apoyos[i][4]) r++; // Restricc ión Ry } // Aplicación de las restricciones a las cargas // Generación del vector de restricciones int n = 3 * mod.apoyos.length - r; int[] restr = new int[n]; double[] Fr = new double[n]; for (int i = 0, k = 0; i < mod.apoyos.length; i++) { for (int j = 0; j < 3; j++) { if (! mod.apoyos[i][j + 2]) { Fr[k] = mod.cargas[i][j + 2]; restr[k] = 3 * i + j; k++; } } } DenseDoubleMatrix1D despl = new DenseDoubleMatrix1D ( 3 * mod.apoyos.length); // Emplea el resolutor definido en opciones // Existen dos tipo de resolutor: usando la descomp osición LU de la
86
// matriz y empleando el método iterativo de los gr adientes // biconjugados, aunque este método puede no conver ger a la solución y // generar error switch (tipo) { case Configuracion.RESOLUTOR_LU: // Aplicación de las restricciones a la mat riz de rigidez SparseDoubleMatrix2D KrLU = new SparseDoubl eMatrix2D(n, n); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { val = kSist.getQuick(restr[i], rest r[j]); if (val != 0) KrLU.setQuick(i, j, v al); } } KrLU.trimToSize(); indicar(); // Cálculo de los desplazamientos LUDecompositionQuick LUD = new LUDecomposit ionQuick(); DenseDoubleMatrix1D UrLU = new DenseDoubleM atrix1D(Fr); LUD.decompose(KrLU); LUD.solve(UrLU); for (int i = 0; i < n; i++) { despl.setQuick(restr[i], UrLU.getQuick( i)); } break; default: // Aplicación de las restricciones a la mat riz de rigidez FlexCompRowMatrix KrIt0 = new FlexCompRowMa trix(n, n); for (int i = 0, k = 1; i < n; i++) { for (int j = 0; j < n; j++) { val = kSist.getQuick(restr[i], rest r[j]); if (val != 0) KrIt0.set(i, j, val); } if (i > k * 500) { // Esto evita que la matriz consuma demasiada memoria KrIt0.compact(); k++; } } KrIt0.compact(); DenseVector FrD = new DenseVector(Fr); DenseVector UrIt = new DenseVector(n); CompRowMatrix KrIt = new CompRowMatrix(KrIt 0); KrIt0 = null; System.gc(); indicar(); // Cálculo de los desplazamientos BiCG sol = new BiCG(FrD); if (precondicionar) { Preconditioner pre = new DiagonalPrecon ditioner(n); pre.setMatrix(KrIt); sol.setPreconditioner(pre); }
87
sol.solve(KrIt, FrD, UrIt); for (int i = 0; i < n; i++) { despl.setQuick(restr[i], UrIt.get(i)); } } // Ensamblaje de la matriz de desplazamientos mod.desplazamientos = new double[mod.apoyos.length] [3]; for (int i = 0; i < mod.desplazamientos.length; i++ ) { for (int j = 0; j < 3; j++) { mod.desplazamientos[i][j] = despl.getQuick( 3 * i + j); } } indicar(); // Cálculo de las cargas DoubleMatrix1D cargas = new Algebra().mult(kSist, d espl); // Ensamblaje de la matriz de cargas mod.reaccionesCargas = new double[mod.apoyos.length ][3]; for (int i = 0; i < mod.reaccionesCargas.length; i+ +) { for (int j = 0; j < 3; j++) { mod.reaccionesCargas[i][j] = cargas.getQuic k(3 * i + j); } } indicar(); // Cálculo de las deformaciones y esfuerzos calcMatDefEsf(mod); indicar();
Se puede ver que este algoritmo llama a su vez a varios métodos que realizan
tareas específicas: calcMatRigSis calcula la matriz de rigidez del sistema,
indicar es usado por la clase Analisis para mostrar el tiempo consumido en
cada etapa del análisis y calcMatDefEsf calcula las deformaciones y esfuerzos.
Por ejemplo, el algoritmo empleado por calcMatRigSis para calcular la matriz
de rigidez del sistema para elementos rectangulares es:
for (int[] elm : mod.elem) { // Conectividad del elemento, 4 nodos con 3 GDL por nodo Ce[0] = 3 * elm[0]; Ce[1] = 3 * elm[0] + 1; Ce[2] = 3 * elm[0] + 2; Ce[3] = 3 * elm[1]; Ce[4] = 3 * elm[1] + 1; Ce[5] = 3 * elm[1] + 2; Ce[6] = 3 * elm[2]; Ce[7] = 3 * elm[2] + 1; Ce[8] = 3 * elm[2] + 2; Ce[9] = 3 * elm[3];
88
Ce[10]= 3 * elm[3] + 1; Ce[11]= 3 * elm[3] + 2; // Para los elementos rectangulares es suficien te con conocer su // ancho y su longitud a = mod.nxy.get(elm[1])[0] - mod.nxy.get(elm[0] )[0]; b = mod.nxy.get(elm[2])[1] - mod.nxy.get(elm[1] )[1]; // Cálculo de la matriz de rigidez del elemento , de orden 12x12 prop = mod.propiedad.get(elm[4]); Ke = ElementoPlaca.matrizRigidezRec( prop[0], prop[1], prop[2], a, b); // Se agrega la matriz de rigidez del elemento a la del sistema for (int i = 0; i < 12; i++) { for (int j = 0; j < 12; j++) { kSist.setQuick(Ce[i], Ce[j], kSist.getQuick(Ce[i], Ce[j]) + Ke[i ][j]); } } }
Este algoritmo es bastante sencillo: para cada elemento se calcula el vector de
conectividades, su largo y ancho, se calcula su matriz de rigidez y esta se agrega
a la matriz de rigidez del sistema.
Elemento2D.- Clase que provee de métodos para el cálculo de las matrices de
rigidez y esfuerzo para el elemento tipo 2D, los métodos disponibles para ello son:
matrizRigidezRec , matrizRigidezTriDKT , matrizEsfuerzoRec y
matrizEsfuerzoTriDKT .
ElementoPlaca.- Esta clase, al igual que la anterior, proporciona métodos para el
cálculo de las matrices de rigidez y esfuerzo para el elemento tipo placa, los
métodos disponibles son: matrizRigidezRec2D , matrizRigidezTriCST ,
matrizEsfuerzoRec2D y matrizEsfuerzoTriCST .
Indicador.- Clase que permite que se pueda conocer el progreso del análisis,
dispone de dos métodos sincronizados indicar y obtener . Su uso se puede
mostrar con la siguiente fracción de código del constructor de la clase Analisis :
ind = new Indicador(); inicio = new InicioAnalisis(ind); detalles = new DetallesAnalisis(ind);
89
Primero se crea una nueva instancia de la clase Indicador llamada ind y luego
se crean dos instancias de InicioAnalisis y DetallesAnalisis . Ambas
clases derivan de Thread y, por lo tanto, son hilos de ejecución independientes.
Al pasarles ind a los constructores de ambas clases estas pueden interactuar:
inicio ejecuta ind.indicar() y luego detalles llama a ind.obtener().
De esta forma, secuencialmente, se conoce el tiempo transcurrido y este se
muestra en la ventana de análisis.
Malla2D.- Clase que realiza el mallado estructurado y no estructurado. Sus dos
métodos public son mallaEstructurada y mallaNoEstructurada , que
reciben como único parámetro una instancia de la clase Modelo (explicada
posteriormente). El código de mallaEstructurada no es muy largo y, por lo
tanto, se lista completo a continuación:
public static void mallaEstructurada(Modelo mod) { int nx, ny; double[] n1 = mod.linea.get(0); double[] n2 = mod.linea.get(1); // Inicializa los nodos y elementos mod.nxy = new ArrayList<double[]>(); mod.elem = new ArrayList<int[]>(); // Número de elementos en X y en Y if (mod.densidadTamano) { double A0 = (n1[2] - n1[0]) * (n2[3] - n2[1 ]) / mod.densidad * (mod.formaElemen to + 1); nx = (int) Math.ceil((n1[2] - n1[0]) / Math .sqrt(A0)); ny = (int) Math.ceil((n2[3] - n2[1]) / Math .sqrt(A0)); } else { nx = (int) Math.ceil((n1[2] - n1[0]) / mod. tamano); ny = (int) Math.ceil((n2[3] - n2[1]) / mod. tamano); } // Generación de las coordenadas de los nodos mod.nxy.ensureCapacity((nx + 1) * (ny + 1)); for (int i = 0; i < nx + 1; i++) { for (int j = 0; j < ny + 1; j++) { mod.nxy.add(new double[]{n1[0] + i * (n 1[2] - n1[0]) / nx, n2[1] + j * (n2[3] - n2[1]) / ny}); } } // Numeración de los nodos: Considerando el pla no XY, los nodos se // numeran de abajo hacia arriba y de izquierda a derecha if (mod.formaElemento == Modelo.ELEMENTO_RECTAN GULAR) { mod.elem.ensureCapacity(nx * ny); for (int i = 0; i < nx; i++) {
90
for (int j = 0; j < ny; j++) { // Para cada elemento la numeración es en sentido // antihorario mod.elem.add(new int[]{j + i * (ny + 1), j + (i + 1) * (ny + 1), j + (i + 1) * (ny + 1) + 1, j + i * (ny + 1) + 1, 0}); } } } else { mod.elem.ensureCapacity(2 * nx * ny); for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) { // Para cada elemento la numeración es en sentido // antihorario mod.elem.add(new int[]{j + i * (ny + 1), j + (i + 1) * (ny + 1), j + i * (ny + 1) + 1, 0}); mod.elem.add(new int[]{j + (i + 1) * (ny + 1), j + (i + 1) * (ny + 1) + 1, j + i * (ny + 1) + 1, 0}); } } } }
Este código inicia calculando el número de elementos que tendrá la malla en X y
en Y. Luego crea los nodos sobre la superficie del modelo partiendo de la
coordenada inferior izquierda de la geometría hasta llegar a la coordenada
superior derecha. El último paso es crear los elementos, los cuales están
formados por elementos numerados en sentido antihorario. Por ejemplo, en la
figura 4.26 se muestra como quedarían numerados los nodos y elementos de una
geometría con una malla estructurada formada por elementos rectangulares.
Figura 4.26 Numeración de nodos y elementos
91
MatrizOp.- Clase que dispone de varios métodos generales para operar con
vectores y matrices. Entre las operaciones que realizan está conocer los valores
máximo y mínimo de vectores, transpuesta de matrices, suma, resta y
multiplicación de vectores y matrices y varias operaciones para obtener
submatrices.
4.4.3 EJECUCIÓN DEL PROGRAMA
Para la ejecución del programa se requiere que esté instalada la máquina virtual
Java versión 5.0 o posterior y el módulo java3D versión 1.3.1 o posterior en el
equipo a utilizar. Para la instalación es suficiente con dar doble clic al fichero
InstalarPlaca.jar (solo válido en Windows) o escribir en una terminal de línea de
comandos:
PATH_JAVA/java -jar PATH_INSTALADOR/InstalarPlaca.j ar
Donde PATH_JAVA es la ruta hasta el ejecutable java y PATH_INSTALADOR es la ruta
que indica donde está InstalarPlaca.jar. Una vez instalado, para iniciar el
programa se puede dar doble clic al acceso directo del programa, o ejecutar
EjecutarPlaca.bat en Windows y EjecutarPlaca.sh en Linux. Para otras
plataformas se debe escribir en una terminal de línea de comandos:
PATH_JAVA/java -Xmx256m -jar PATH_INSTALADOR/Placa. jar
4.4.4 LIBRERÍAS EXTERNAS USADAS POR EL PROGRAMA
Como se indicó anteriormente, el programa emplea 3 librerías de distribución
gratuita, completamente desarrolladas en Java. A continuación se da una breve
descripción de cada una.
colt.- Es una librería que provee una infraestructura para la computación técnica y
científica en java. Contiene algoritmos para el análisis de datos, álgebra lineal,
estadística, programación paralela y concurrente, entre otros. Para más
información se puede consultar su página Web: http://dsd.lbl.gov/~hoschek/colt
92
La parte empleada por el programa es la referente al álgebra lineal. Por ejemplo,
el cálculo de los desplazamientos se realiza con la siguiente porción de código:
// Cálculo de los desplazamientos LUDecompositionQuick LUD = new LUDecompositionQuick (); DenseDoubleMatrix1D UrLU = new DenseDoubleMatrix1D( Fr); LUD.decompose(KrLU); LUD.solve(UrLU);
mtj.- Es una colección detallada de estructuras de datos de matrices, resolutores
lineales, métodos de mínimos cuadrados y descomposiciones de vectores y
valores propios. Está diseñada para ser usada como una librería para el
desarrollo de aplicaciones numéricas. Su página Web es: http://rs.cipr.uib.no/mtj
Su capacidad para operar con matrices dispersas es empleada para el
almacenaje de la matriz del sistema y para el cálculo de los desplazamientos
empleando un método iterativo, tal como se muestra a continuación:
// Aplicación de las restricciones a la matriz de r igidez FlexCompRowMatrix KrIt0 = new FlexCompRowMatrix(n, n); for (int i = 0, k = 1; i < n; i++) { for (int j = 0; j < n; j++) { val = kSist.getQuick(restr[i], restr[j]); if (val != 0) KrIt0.set(i, j, val); } if (i > k * 500) { // Esto evita que la matriz consuma demasia da memoria KrIt0.compact(); k++; } } KrIt0.compact(); DenseVector FrD = new DenseVector(Fr); DenseVector UrIt = new DenseVector(n); CompRowMatrix KrIt = new CompRowMatrix(KrIt0); KrIt0 = null; System.gc(); indicar(); // Cálculo de los desplazamientos BiCG sol = new BiCG(FrD); if (precondicionar) { Preconditioner pre = new DiagonalPreconditioner (n); pre.setMatrix(KrIt); sol.setPreconditioner(pre); } sol.solve(KrIt, FrD, UrIt);
93
El uso del método iterativo reduce considerablemente el tiempo de cálculo si se
compara con el método de la descomposición LU empleado por la librería colt,
pero en cambio, al tratarse de un proceso iterativo, puede no converger a la
solución del problema.
IzPack.- Es una herramienta que permite construir instaladores que se ejecuten en
cualquier sistema operativo. Su diseño es modular, permitiendo seleccionar las
aplicaciones que se instalarán para una plataforma específica. Su página Web es:
http://www.izforge.com/izpack
IzPack emplea archivos de extensión xml para describir la instalación. Por
ejemplo, el siguiente fragmento de código del archivo InstalarPlaca.xml genera el
instalador para el programa:
<?xml version="1.0" encoding="iso-8859-1" standalon e="yes" ?> <!-- Instalador del programa Placa Creado con IzPack versión 3.7.2 --> <installation version="1.0"> <!-- Sección de información --> <info> <appname>Placa</appname> <appversion>1.0</appversion> <authors> <author name="Jorge Eduardo Romero" email="[email protected]"/> </authors> <url>http://www.epn.edu.ec</url> <javaversion>1.5.0</javaversion> </info> <!-- Preferencias de la instalación Instalador de dimensiones 640x480. Sin opción de cambio --> <guiprefs width="640" height="480" resizable="n o"> <laf name="metouia"> <os family="unix" /> </laf> </guiprefs> <!--
94
Sección de localización Idiomas disponibles: español e inglés --> <locale> <langpack iso3="spa"/> </locale> <!-- Sección de recursos --> <resources> <res id="HTMLInfoPanel.info" src="Leame.htm "/> <res id="shortcutSpec.xml" src="AccesoWindows.xm l"/> <res id="Unix_shortcutSpec.xml" src="AccesoUnix. xml"/> </resources> <!-- Sección de paneles Se usarán los paneles en el orden en que es tán indicados a continuación --> <panels> <panel classname="HelloPanel"/> <panel classname="HTMLInfoPanel"/> <panel classname="TargetPanel"/> <panel classname="PacksPanel"/> <panel classname="InstallPanel"/> <panel classname="ShortcutPanel"/> <panel classname="SimpleFinishPanel"/> </panels> <!-- Sección de paquetes --> <packs> <pack name="Placa" required="yes"> <description>Programa compilado en NetB eans IDE 5.0 </description> <file src="Placa.jar" targetdir="$INSTA LL_PATH"/> <file src="Leame.htm" targetdir="$INSTA LL_PATH"/> <file src="MaterialesUsuario.msr" targetdir="$INS TALL_PATH"/> <file src="bin" targetdir="$INSTALL_PATH"/> <file src="iconos" targetdir="$INSTALL_PATH"/> <file src="lib" targetdir="$INSTALL_PATH"/> <parsable targetfile="$INSTALL_PATH/bin/EjecutarP laca.sh"/> <parsable targetfile= "$INSTALL_PATH/bin/DesinstalarPlaca.sh "/> <executable targetfile= "$INSTALL_PATH/bin/EjecutarPlaca.sh" s tage="never"/> <executable targetfile= "$INSTALL_PATH/bin/DesinstalarPlaca.sh " stage="never"/> </pack> <pack name="Ejemplos" required="no" presele cted="yes"> <description>Ejemplos realizados con el programa Placa </description> <file src="ejemplos" targetdir="$INSTAL L_PATH"/> </pack> </packs> <!--
95
Sección de archivos nativos Permite crear accesos directos en Windows --> <native type="izpack" name="ShellLink.dll"/> </installation>
Para crear el instalador considerando que IzPack está instalado en D:\Archivos de
programa\ se debe escribir en una terminal de línea de comandos:
D:\Archiv~1\IzPack\bin\compile InstalarPlaca.xml -b .
Esta línea de código crea el archivo InstalarPlaca.jar. Para mayor información se
puede consultar la ayuda incluida con IzPack.
96
CAPÍTULO 5.
PRUEBAS DE VERIFICACIÓN
Los ejemplos que a continuación se detallan sirven para probar la capacidad de
cálculo del programa. Además, para verificar que todos los algoritmos
desarrollados trabajan correctamente, se comparará las respuestas con las
obtenidas por el programa Algor versión 19.
5.1 EJEMPLOS DE PRUEBA PARA ELEMENTOS TIPO PLACA
5.1.1 EJEMPLO 1: PLACA RECTANGULAR
Este ejemplo consiste en una placa rectangular de dimensiones 1 x 2 m y espesor
5 cm. empotrada en una pared. El material del que está hecha la placa tiene un
módulo de elasticidad de 200x109 N/m2 y un módulo de Poisson de 0.3. Sobre la
placa actúa una carga perpendicular de 1000 N tal como se muestra en la figura
5.1. Se pide calcular la deflexión en el punto de aplicación de la carga y los
esfuerzos de von Mises y Tresca máximos.
Figura 5.1 Placa empotrada
97
Comparación de resultados
Los resultados obtenidos se muestran en la tabla 5.1 y corresponden a mallas con
una densidad aproximada de 400 elementos rectangulares o cuadriláteros. Como
se puede observar, para el caso de la deflexión, la diferencia en los resultados es
despreciable.
Tabla 5.1 Resultados del ejemplo 1
Matlab Java Algor
Forma elemento Rectangular Rectangular Cuadrilateral
Número elementos 435 435 406
Deflexión 0.0013658 m 0.00136575 m 0.00136508 m
Esfuerzo von Mises 5.5913e6 Pa 5.59126e6 Pa 5.32326e6 Pa
Esfuerzo Tresca 2.8401e6 Pa 2.84013e6 Pa 2.99774e6 Pa
En la figura 5.2 se muestra la geometría deformada tal como es presentada por
los tres programas.
Figura 5.2 Resultados del ejemplo 1. (a) Matlab, (b) Java, (c) Algor
(a) (b)
(c)
98
Nótese que no es posible hacer coincidir el número de elementos utilizado
durante el análisis en el programa desarrollado en Matlab o Java y en Algor v19.
Esto se debe a que en ambos casos el usuario solo está en capacidad de indicar
al programa la densidad aproximada de la malla que desea utilizar durante el
análisis. Las subrutinas de discretización de cada programa son las que
determinan la densidad final de la malla en función de los algoritmos que
incorporan.
5.1.2 EJEMPLO 2: DISCO
Un disco como el de la figura 5.3 tiene un espesor de 2.5 mm y está hecho de
acero inoxidable AISI 202 (E = 206.84x109 N/m2, ν = 0.3). Se busca calcular los
valores máximos para la deflexión, las deformaciones y esfuerzos en las
direcciones X y Y.
Figura 5.3 Disco
Comparación de resultados
Los resultados corresponden a una malla con una densidad aproximada de 3000
elementos. En este caso la relativamente alta densidad de las mallas utilizadas
ayudó para que los resultados obtenidos por los diferentes programas sean
similares, como se muestra en la tabla 5.2.
99
Tabla 5.2 Resultados del ejemplo 2
Matlab Java Algor
Forma elemento Triangular Triangular Triangular
Número elementos 2970 2970 3056
Desplazamiento Z 1.2323e-4 m 1.23229e-4 m 1.24118e-4 m
Deformación X 7.2291e-4 7.22908e-4 7.6983e-4
Deformación Y 5.7295e-4 5.72975e-4 5.32813e-4
Esfuerzo X 1.5596e8 Pa 1.55959e8 Pa 1.55862e8 Pa
Esfuerzo Y 1.3304e8 Pa 1.33036e8 Pa 1.21026e8 Pa
Adicionalmente, en la figura 5.4 se puede observar que los tres programas
presentan gráficas muy parecidas para la geometría deformada.
Figura 5.4 Resultados del ejemplo 2. (a) Matlab, (b) Java, (c) Algor
(a) (b)
(c)
100
5.1.3 EJEMPLO 3: VIGA
La viga de la figura 5.5 es de acero ASTM-A572 (E = 199.95x109 N/m2, ν = 0.29) y
tiene un espesor constante de 7 cm. La carga distribuida sobre la viga ejerce una
fuerza total de 400 N. Calcular los valores máximos del desplazamiento en la
dirección Z, las deformaciones y esfuerzos en la dirección de los ejes X y Y y los
esfuerzos de von Mises y Tresca máximos.
Figura 5.5 Viga
Comparación de resultados
En la tabla 5.3 se puede ver que hay una pequeña variación en los valores
obtenidos aún entre las versiones de Matlab y Java del programa desarrollado.
Esto se debe a que el programa en Java emplea un resolutor iterativo. Sin
embargo, si se refina la malla las diferencias en los resultados se vuelven cada
vez más pequeñas, tal como se demuestra en la sección 5.3.3.
101
Tabla 5.3 Resultados del ejemplo 3
Matlab Java Algor
Forma elemento Triangular Triangular Triangular
Número elementos 3054 3054 2996
Desplazamiento Z 1.6667e-3 m 1.66667e-3 m 1.66585e-3 m
Deformación X 2.8706e-5 2.88726e-5 2.78834e-5
Deformación Y 8.5595e-6 8.57438e-6 8.51359e-6
Esfuerzo X 5.8006e6 Pa 5.83038e6 Pa 5.54861e6 Pa
Esfuerzo Y 1.0244e6 Pa 1.02442e6 Pa 9.74914e5 Pa
Esfuerzo von Mises 5.6992e6 Pa 5.73461e6 Pa 5.60106e6 Pa
Esfuerzo Tresca 2.7958e6 Pa 2.81671e6 Pa 2.8438e6 Pa
En la figura 5.6 se muestran las gráficas para la deflexión generadas por los tres
programas.
Figura 5.6 Resultados del ejemplo 3. (a) Matlab, (b) Java, (c) Algor
(a) (b)
(c)
102
5.2 EJEMPLOS DE PRUEBA PARA ELEMENTOS TIPO 2D
5.2.1 EJEMPLO 1: HERRADURA
La herradura que se muestra en la figura 5.7 es de acero AISI 1015 (E =
204.77x109 N/m2 y ν = 0.28), tiene un espesor de 1 mm. Calcular el
desplazamiento horizontal máximo y los esfuerzos de von Mises y Tresca
máximos cuando se aplica una fuerza de 70 N.
Figura 5.7 Herradura
Comparación de resultados
En la tabla 5.4 se muestran los resultados obtenidos utilizando los programas
desarrollados y el software comercial. Como se puede ver, en este caso existe
una pequeña diferencia en los resultados. Debe recordarse que la teoría de
placas es válida si el espesor es menor a 1/10 de las dimensiones de la placa y
que los desplazamientos resultantes deben ser menores a 1/10 del espesor.
Cuanto más difieran estos valores, los resultados serán más inexactos. Las
gráficas presentadas por los tres programas se muestran en la figura 5.8.
103
Tabla 5.4 Resultados del ejemplo 1
Matlab Java Algor
Forma elemento Triangular Triangular Triangular
Número elementos 1008 1008 970
Despl. horizontal 0.00018275 m 0.000182759 m 0.00017954 m
Esfuerzo von Mises 1.7015e8 Pa 1.70197e8 Pa 1.86317e8 Pa
Esfuerzo Tresca 8.2713e7 Pa 8.27425e7 Pa 9.392e7 Pa
Figura 5.8 Resultados del ejemplo 1. (a) Matlab, (b) Java, (c) Algor
5.2.2 EJEMPLO 2: LÁMINA DELGADA
Una lámina de espesor uniforme está sujeta a la acción de las cargas que se
muestran en la figura 5.9. La lámina esta hecha de acero ASTM-A36 (E =
199.95x109 N/m2 y ν = 0.29) y tiene un espesor de 0.5 mm. Calcular los valores
máximos para el desplazamiento, deformación y esfuerzo en las direcciones X y Y
y los esfuerzos de von Mises y Tresca máximos.
(a) (b) (c)
104
Figura 5.9 Lámina delgada
Comparación de resultados
Para este caso, al ser la densidad de la malla bastante alta, los resultados difieren
muy poco. Además, como se dijo anteriormente, si se refina la malla los valores
tienden a converger. Los resultados se muestran en la tabla 5.5 y en la figura
5.10.
Tabla 5.5 Resultados del ejemplo 2
Matlab Java Algor
Forma elemento Triangular Triangular Triangular
Número elementos 3064 3064 3000
Desplazamiento X 4.9998-6 m 4.99977e-6 m 4.99472e-6 m
Desplazamiento Y 5.0263e-7 m 5.02733e-7 m 5.0283e-7 m
Deformación X 1.0435e-4 1.04351e-4 1.02118e-4
Deformación Y 0 0 -1.43346e-36
Esfuerzo X 2.1173e7 Pa 2.1173e7 Pa 2.06995e7 Pa
Esfuerzo Y 4.4692e6 Pa 4.46928e6 Pa 4.50737e6 Pa
Esfuerzo von Mises 2.0926e7 Pa 2.09262e7 Pa 2.05068e7 Pa
Esfuerzo Tresca 1.0236e7 Pa 1.02357e7 Pa 1.044635e7 Pa
105
Figura 5.10 Resultados del ejemplo 2. (a) Matlab, (b) Java, (c) Algor
5.2.3 EJEMPLO 3: PLACA
La placa de la figura 5.11 tiene un espesor de 2 cm. y está hecha de aluminio
2024-T3 (E = 73.1x109 N/m2, ν = 0.33). Se busca calcular los desplazamientos,
deformaciones y esfuerzos máximos en las direcciones X y Y y los esfuerzos de
von Mises y Tresca máximos.
Figura 5.11 Placa
(a) (b)
(c)
106
Comparación de resultados
En la tabla 5.6 y en la figura 5.12 se resumen los resultados obtenidos. Se puede
notar que la mayor disparidad se presenta en el esfuerzo de Tresca.
Tabla 5.6 Resultados del ejemplo 3
Matlab Java Algor
Forma elemento Triangular Triangular Triangular
Número elementos 9872 9872 10018
Desplazamiento X 2.7825e-4 m 2.78246e-4 m 2.76509e-4 m
Desplazamiento Y 5.8059e-4 m 5.80589e-4 m 5.77171e-4 m
Deformación X 4.403e-5 4.40304e-5 4.30953e-5
Deformación Y 3.0571e-5 3.05693e-5 2.93243e-5
Esfuerzo X 3.3756e6 Pa 3.3756e6 Pa 3.33089e6 Pa
Esfuerzo Y 2.5355e6 Pa 2.53543e6 Pa 2.48163e6 Pa
Esfuerzo von Mises 9.3996e6 Pa 9.39959e6 Pa 9,40915e6 Pa
Esfuerzo Tresca 4.4826e6 Pa 4.48263e6 Pa 5.01215e6 Pa
Figura 5.12 Resultados del ejemplo 3. (a) Matlab, (b) Java, (c) Algor
(a) (b)
(c)
107
5.3 COMPARACIÓN DE LOS TRES PROGRAMAS
En las secciones anteriores se ha verificado mediante ejemplos que el programa,
en sus versiones de Matlab y Java, funciona correctamente. Ahora se va a
comparar los tiempos requeridos para realizar las dos operaciones que consumen
más tiempo: mallado y solución del problema. También se comparará si al refinar
la malla los resultados de los tres programas tienden a converger a una solución
común.
5.3.1 COMPARACIÓN DEL MALLADO
El programa desarrollado en este trabajo emplea el método del avance frontal
para realizar el mallado triangular no estructurado, además puede mallar
estructuradamente placas rectangulares. Algor solo puede realizar mallas no
estructuradas, pero en cambio estas pueden estar formadas por elementos
triangules, cuadriláteros o una combinación de ambos.
En la figura 5.13 se muestran dos gráficas de curvas de tiempos empleados por
los tres programas para realizar el mallado triangular no estructurado como
función del número de elementos utilizados para el análisis. Las geometrías
usadas para la comparación corresponden a las de los ejemplos de las secciones
5.1.1 y 5.2.1.
En ambas gráficas se puede notar que Algor es muy eficiente para mallar
geometrías, mientras que la versión del programa en Matlab es el más lento.
108
0
10
20
30
40
50
60
0 2000 4000 6000 8000 10000 12000
Número de elementos
Tie
mpo
(s)
Matlab Java Algor
0
10
20
30
40
50
60
70
0 2000 4000 6000 8000 10000 12000
Número de elementos
Tie
mpo
(s)
Matlab Java Algor
Figura 5.13 Tiempos de mallado para los ejemplos de (a) sección 5.1.1 y (b)
sección 5.2.1
5.3.2 COMPARACIÓN DE LA SOLUCIÓN
Procediendo en igual forma que la sección anterior, se resolverá los ejemplos de
las secciones 5.1.1 y 5.2.1 pero variando la densidad de la malla para poder
registrar el tiempo requerido por cada programa para encontrar la solución de los
(a)
(b)
109
problemas. Estos resultados se muestran en la figura 5.14. Para la versión en
Java se empleó el método iterativo precondicionado como resolutor.
0
20
40
60
80
100
120
140
0 2000 4000 6000 8000 10000 12000
Número de elementos
Tie
mpo
(s)
Matlab Java Algor
0
5
10
15
20
25
30
35
0 2000 4000 6000 8000 10000 12000
Número de elementos
Tie
mpo
(s)
Matlab Java Algor
Figura 5.14 Tiempos de cálculo para los ejemplos de (a) sección 5.1.1 y (b)
sección 5.2.1
(a)
(b)
110
En este caso se puede ver que para el ejemplo de la sección 5.1.1 formado por
elementos tipo placa (3 GDL por nodo), la versión del programa en Java
incrementa el tiempo requerido para la solución rápidamente. En cambio, para el
ejemplo de la sección 5.2.1 formado por elementos tipo 2D (2 GDL por nodo),
este tiempo es mucho menor, inclusive menor que el tiempo de Algor. En ambos
casos la versión del programa en Matlab es la más rápida.
5.3.3 COMPARACIÓN DE LA CONVERGENCIA
Para comparar la convergencia se volverá a resolver los ejemplos de las
secciones 5.1.1 y 5.2.1 pero variando la densidad de las mallas. Solo se compara
el esfuerzo de von Mises máximo tal como se muestra en las tablas 5.7 y 5.8.
Como se observa en las tablas 5.7 y 5.8, a medida que aumenta el número de
elementos de la malla, los resultados del programa desarrollado, en sus versiones
en Matlab y Java presentan una menor diferencia al ser comparados con los
resultados producidos por Algor. En la tabla 5.8 se nota que no constan los
valores para 100 y 250 elementos usando Algor, esto es debido a que el mallador
no permitió densidades tan bajas.
Tabla 5.7 Convergencia del esfuerzo de von Mises (ejemplo sección 5.1.1)
# elementos (aprox.) Matlab Java Algor
100 5.5462e6 Pa 5.54619e6 Pa 5.21803e6 Pa
250 5.6248e6 Pa 5.62482e6 Pa 5.30214e6 Pa
500 5.5502e6 Pa 5.55016e6 Pa 5.33605e6 Pa
1000 5.4576e6 Pa 5.4576e6 Pa 5.36716e6 Pa
2000 5.4599e6 Pa 5.45995e6 Pa 5.38113e6 Pa
3000 5.46e6 Pa 5.46004e6 Pa 5.41223e6 Pa
4000 5.4599e6 Pa 5.4599e6 Pa 5.42447e6 Pa
5000 5.4604e6 Pa 5.46035e6 Pa 5.4166e6 Pa
111
Tabla 5.8 Convergencia del esfuerzo de von Mises (ejemplo sección 5.2.1)
# elementos (aprox.) Matlab Java Algor
100 5.8363e7 Pa 5.83632e7 Pa -
250 1.1611e8 Pa 1.16112e8 Pa -
500 1.5148e8 Pa 1.51494e8 Pa 8.7488e7 Pa
1000 1.7015e8 Pa 1.70197e8 Pa 1.86317e8 Pa
2000 1.8367e8 Pa 1.83671e8 Pa 2.02481e8 Pa
3000 1.8928e8 Pa 1.89275e8 Pa 2.10061e8 Pa
4000 1.9214e8 Pa 1.9214e8 Pa 1.94843e8 Pa
5000 1.9415e8 Pa 1.94147e8 Pa 1.98387e8 Pa
Cabe aclarar que dependiendo del estado de cargas al que está sujeto una
geometría, el número de elementos de la malla necesario para que el error en los
resultados sea despreciable puede variar ampliamente. En todo caso, cuando se
esperen esfuerzos locales muy altos es recomendable realizar mallas más
densas.
112
CAPÍTULO 6.
CONCLUSIONES Y TRABAJO FUTURO
6.1 CONCLUSIONES
El presente trabajo ha demostrado que en el Ecuador se puede usar de forma
exitosa el lenguaje de programación Java para el desarrollo de aplicaciones
científicas y técnicas. En concreto, se ha demostrado que Java puede ser
utilizado de manera eficiente y efectiva para desarrollar software para resolver
problemas de análisis y diseño de estructuras mediante el método de los
elementos finitos. A pesar de que Java es un lenguaje interpretado (se ejecuta
algo más lento que C++, por ejemplo), sus ventajas (lenguaje seguro, fácil de
aprender, multiplataforma, gran cantidad de APIs disponibles) superaron con
mucho a esta limitación.
Los ejemplos de prueba estudiados en el capítulo 5 de este trabajo permiten
afirmar que el software desarrollado funciona correctamente, siempre y cuando el
problema siendo analizado cumpla las condicionantes de la teoría de placas
empleada por el programa. Aún así, se puede notar que mientras la malla sea
más densa los resultados serán más satisfactorios, aunque podría bastar con
hacer un refinamiento local a la malla en las regiones donde se espere obtener
valores críticos de desplazamientos, deformaciones o esfuerzos, para conseguir
buenos resultados.
Al comparar los tiempos empleados para discretizar un dominio determinado
utilizando Algor y las versiones del programa desarrolladas en Matlab y en Java,
se comprobó que Algor produjo mallas de mejor calidad en un tiempo mucho
menor. Sin embargo, el hecho de que mallar en la versión desarrollada en Java
aproximadamente 10000 elementos tome alrededor de 30 segundos indica que sí
se puede emplear este programa para generar mallas bastante densas de manera
eficiente.
113
Por otro lado, el tiempo requerido en Java para la solución de los sistemas de
ecuaciones derivados resultó ser, en general, el peor de entre los tres programas
estudiados. Si bien este hecho constituye una limitante para la densidad del
mallado que se puede emplear, se debe reconocer que el tener que esperar
menos de 2 minutos para analizar un problema cuya geometría está formada por
10000 elementos es insignificante comparado con el tiempo que requeriría un
método analítico, o un proceso manual, sobre todo en geometrías complejas.
Finalmente, el hecho de que la versión desarrollada en Matlab sea la más rápida
para la resolución de sistemas de ecuaciones, solo refleja que Matlab emplea
algoritmos de cálculo altamente optimizados, sobre todo para el cálculo matricial.
Esta es precisamente la ventaja que ha permitido a Matlab colocarse a la cabeza
de las ventas en el segmento de software comercial empleado para realizar
cálculos matemáticos.
6.2 TRABAJO FUTURO
Este programa ha sido diseñado de forma que pueda ser modificado, mejorado y
ampliado posteriormente por cualquier persona con un conocimiento intermedio
del lenguaje de programación Java. Durante este desarrollo del programa existen
varios puntos que podrían ser considerados. A continuación se describen los
principales:
Mejorar el mallador.- El mallador utilizado por el programa emplea el método del
avance frontal. Este no es el único método disponible y sería muy beneficioso
para el programa si se dispusiera de, por ejemplo, un mallador que utilice el
método de la triangulación de Delaunay. También se puede buscar un método de
alisado que produzca mallas de mejor calidad que el alisado Laplaciano
actualmente utilizado.
Uso de geometrías complejas.- Las B-Splines y NURBS (nonuniform racional b-
spline) son representaciones de curvas que permiten construir complejas
geometrías a partir de un conjunto de puntos en el espacio. El uso de estos tipos
114
de curvas ampliarían el tipo de problemas que pueden ser analizados por el
programa.
Resolutor propio.- La actual implementación del programa desarrollado emplea dos
librerías generales y de libre distribución para la resolución de matrices dispersas.
Si se creasen subrutinas de cálculo específicas para el programa desarrollado, de
seguro se ganaría en velocidad de análisis y se reduciría el tamaño del programa.
Tipos de elementos.- Este es un punto que puede significativamente mejorar el
programa. Actualmente el programa solo es capaz de resolver problemas de
mecánica de sólidos que involucren el uso de elementos tipo placa y 2D. Sin
embargo, el método de los elementos finitos permite resolver problemas que van
desde simples estructuras definidas en una dimensión hasta el análisis dinámico
no lineal de complejas geometrías 3D.
Sería bastante sencillo aumentar los tipos de elementos que se puedan utilizar en
la resolución de problemas con este programa, empezando con elementos tipo
viga columna en 2 y 3 dimensiones, elementos tipo membrana (shell) y placas
gruesas. Otros tipos de elementos, tales como los tipo prisma o bricks pudieran
requerir mucho más tiempo.
115
BIBLIOGRAFÍA
[1] XII CURSO DE MÁSTER EN MÉTODOS NUMÉRICOS PARA CÁLCULO Y
DISEÑO EN INGENIERÍA; Técnicas de Pre y Postproceso Gráfico; 2001.
[2] BATOZ, J. L.; An Explicit Formulation for an Efficient Triangular Plate-
Bending Element; 1982.
[3] BOUVIER, Dennis; Getting Started with the Java 3D API; Sun
Microsystems; 2001.
[4] CAMPIONE, Mary; WALRATH, Kathy; HUML, Alison; The Java Tutorial;
Third Edition; Pearson Education.
[5] CERVANTES, Guillermo; MEJÍA, Carlos; Precondicionamiento de Métodos
Iterativos; Revista de la Academia Colombiana de Ciencias; Volumen 28;
2004.
[6] ECKEL, Bruce; Thinking in Java; Second Edition; 2000.
[7] GARCÍA, Javier; Aprenda Java Como si Estuviera en Primero; 2000.
[8] GARCÍA, Javier; Aprenda Matlab 6.5 Como si Estuviera en Primero;
Madrid; 2004.
[9] GOSLING, James; The Java Language Specification; Second Edition;
Addison Wesley; 2000.
[10] HOMMEL, Scott; Convenciones de Código para el Lenguaje de
Programación JAVA; Sun Microsystems Inc.; 1999.
[11] HUTTON, David V.; Fundamentals of Finite Element Analysis; McGraw-Hill;
2004.
116
[12] HYER, M. W.; Laminated Plate and Shell Theory; 2000.
[13] KANSARA, Kaushalkumar; TESIS; Development of Membrane, Plate and
Flat Shell Elements in Java; Virginia Polytechnic Institute; 2004.
[14] KUNWOO, Lee; Principles of CAD/CAM/CAE Systems; ADDISON
WESLEY; 1999.
[15] LÖHNER, Rainald; Automatic Unstructured Grid Generators; 1997.
[16] OWEN, Steven J.; A Survey of Unstructured Mesh Generation Technology;
Proceedings 7th International Meshing Roundtable, Dearborn, MI, October
1998.
[17] PADRÓN, Miguel A.; SUÁREZ, José P.; PLAZA, Angel; Adaptive
Techniques for Unstructured Nested Meshes; 2004.
[18] S. H. Lo; Finite Element Mesh Generation and Adaptive Meshing; 2002.
[19] SHAMES, Irving H.; COZARELLI, Francis A.; Elastic and Inelastic Stress
Analysis; Revised Printing; Taylor & Francis; 1997.
[20] SHIGLEY, Joseph E.; MISCHKE, Charles R.; Diseño en Ingeniería
Mecánica; Sexta edición; McGraw Hill; 2002.
[21] SUN MICROSYSTEMS; The Java 3D API Specification; 2002.
[22] SZILARD, Rudolph; Theories and Applications of Plate Analysis; John
Wiley & Sons; 2004.
[23] WUNDERLICH, Walter; Mechanics of Structures Variational and
Computational Methods; Second edition; CRC Press; 2003.
117
[24] ZIENKIEWICZ, O. C.; Finite Element Procedures in the Solution of Plate
and Shell Problems; 1965.
[25] ZORÍN, Denis; Curvature and Geodesics, Discrete Laplacian and Related
Smoothing Methods; 2002.
[26] http://www.ce.memphis.edu/7117/notes.htm
[27] http://www.inf.ethz.ch/personal/cetin/thesis/thesis/node18.html
[28] http://www.me.cmu.edu/faculty1/shimada/gm98/project/judy/project/
118
ANEXOS
Anexo A. Matriz de rigidez para el elemento rectangular ................................... 119
Anexo B. Matriz de esfuerzo para el elemento rectangular................................ 120
Anexo C. Matriz α para el elemento DKT .......................................................... 121
Anexo D. Matriz de rigidez para el elemento plano rectangular ......................... 122
Anexo E. Matriz de rigidez para el elemento triangular CST.............................. 123
Anexo F. Diagrama de bloques del programa.................................................... 124
Anexo G. Diagramas de flujo.............................................................................. 125
119
Anexo A. Matriz de rigidez para el elemento rectangular
( )
( ) ⋅−
=21180 v
EhK N
e
−−−
−−
−−
−−−−−
−−−−−−
V
ZR
HGF
XNV
TMZR
NMLHGF
YQWKV
UPSJZR
QPOKJIHGF
WKYQXNV
SJUPTMZR
KJIQPONMLHGF
φφ
φφφφ
φφφφφφ
Donde los coeficientes de rigidez individuales son:
( )ab
hvF
222 60601242 −++−= ρρ
( )b
hvG
211 12330 −− ++= ρρρ
( )a
hvH
21 12330 ρρρ ++=−
( )
ab
hvI
222 30601242 −+−+−= ρρ
( )[ ]b
hvJ
211330 −−+= ρρ
( )a
hvK
21 12315 ρρρ −−=−
( )ab
hvL
222 30601242 ρρ +−+−=−
( )
b
hvM
211 12315 −− ++−= ρρρ
( )[ ]a
hvN
21 1330 ρρ −+=−
( )
ab
hvO
222 30301242 −−−−= ρρ
( )[ ]b
hvP
211315 −−+−= ρρ
( )[ ]a
hvQ
21 1315 ρρ −−=−
( )[ ] 211420 hvR −−+= ρρ ( )[ ] 21110 hvS −−−= ρρ
( )[ ] 211410 hvT −−−= ρρ ( )[ ] 2115 hvU −−+= ρρ
( )[ ] 21 1420 hvV ρρ −+= − ( )[ ] 21 1410 hvW ρρ −−= −
( )[ ] 21 110 hvX ρρ −−= − ( )[ ] 21 15 hvY ρρ −+= −
( )215vhZ = 0=φ
b
a=ρ a
b=−1ρ
Simétrica
120
Anexo B. Matriz de esfuerzo para el elemento rectangular
( )
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
−−−−−−−−−−−−−−+−−−−−−+−−−−
−−−−−−−−−−−−−−+−−−−+−−
−−−−−−−−−−−−−+−−+−
−−−−−−−−−−−−−−−−−+−−−−+
=
−−
−−
−−
−−
−−
−−
−−
−−
avbvvbvvvavvbvavbvavvbavvbva
bvvavbvvavvvbbvavvavvbvvbava
vavvavbvvbvvvavbvavbavbavvbv
avvvvvavbvvvavbbvavavbvvbav
ab
DS N
e
111011001101446620600002644662060000260111111010012064466026000206446602600000110111101100002644662060000264466206
10100101111102600020644660260002064466
11
11
11
11
11
11
11
11
ρρρρρρρρ
ρρρρρρρρ
ρρρρρρρρ
ρρρρρρρρ
Donde:
( )2
3
112 v
EhD
−=
b
a=ρ
121
Anexo C. Matriz αααα para el elemento DKT
( ) ( ) ( ) ( ) ( )( )( )
( ) ( )( )( )
( ) ( ) ( )( ) ( ) ( ) ( )
−++−++−−−++−++−−−−+−−−−−−
−−+−+−+−−−+−+−−+−−−+−−−−
−−−−+−−−+−
−−+−−−−−
−−−
=
3452354354434343533553
35236352343523343443235523523
24424242342342623233623
255252336352233525263
5234324352343523434343523523523
424242424234223
525252352522352
543543543434343535353
363363
363363
41212
242224
11100100
2200040200002040
yqqxyrryttqyqxtyqyyrtyrxrxqxqxpxpxxryrpxxrqxpx
xrqxpxrxxqxypxpxxypxxrqxpxxypxrxxyqxpxpxqxqxxrxrxtxtxqxrxtxqxrxtx
qxrxtxqxrxxtxxqxrxtxxqxrxxtx
rryqqyppyryqypyryqypyypyypyypyypy
α
Con:
jiij xxx −= jiij yyy −= 222ijijij yxl +=
223234 6 lxp −= 2
3135 6 lxp −= 212126 6 lxp −=
223234 6 lyt −= 2
3135 6 lyt −= 22323234 3 lyxq =
231335 3 lyxq = 2
232234 3 lyr = 2
312315 3 lyr =
NOTA: Las filas 9 y 10 forman parte de la misma expresión.
122
Anexo D. Matriz de rigidez para el elemento plano rectangular
( ) ⋅−
=2124 v
Ehvk ii
−−
−−−−
−−−−−−
4
3
2
1
4
3
2
1
2323
3232
2323
3232
22
22
y
y
y
y
x
x
x
x
uuuuuuuu
ABAACACABA
vvvvAvvvvCA
vvvvABAvvvvBACA
βα
βαβα
βαβαβα
βαβαβαβα
αβ
αβαβ
αβαβαβ
αβαβαβαβ
Donde:
ab=α ba=β vv −= 11 ( )vv −= 132 ( )vv 3133 −=
148 vA βααβ += 144 vB βααβ −= 128 vC βααβ +−=
148 vA αββα += 144 vB αββα −= 128 vC αββα +−=
Simétrica
123
Anexo E. Matriz de rigidez para el elemento triangular CST
( ) ⋅−
=214 vA
Ehk
++++++
++++
+++++
++++++
23
23
323222
22
313121212
121
333223311323
23
2332222112323222
22
13311221113131212121
21
2
12
12
1
βγββγγβγββγγββγγβγ
γβγβγβγβγβγβ
γβγβγβγβγβγγββγβ
γβγβγβγβγβγγββγγββγβ
CCCCCC
vCvCvC
Cvv
CvCC
CvCvv
CCC
Donde:
321 yy −=β 132 yy −=β 213 yy −=β
231 xx −=γ 312 xx −=γ 123 xx −=γ
2
1 vC
−=
Simétrica
124
Anexo F. Diagrama de bloques del programa
Menú principal
Archivo Modelo Análisis Resultados Ver
Contenido Nuevo modelo
Abrir
Guardar como
Guardar resultados
Salir
Tipo de elemento
Propiedades
Malla
Apoyos
Cargas
Analizar
Tipo de análisis
Desplaza-mientos
Reacciones y cargas
Deformacio-nes
Esfuerzos
Orientación
Visualización
Ejes
Apoyos y cargas
Mover
Expandir
Acerca de Placa
Momentos flectores
Acercar / alejar
Rotación 3D
Guardar
Edición
Deshacer
Rehacer
Cortar
Copiar
Pegar
Borrar
Cambios Editor de
geometrías
Selección
Consulta
Herramientas
Librería de materiales
Geometría
Modelo
Resultados
Unidades
Barras de herramientas
Opciones
Ayuda
125
Anexo G. Diagramas de flujo
Los diagramas de flujo que a continuación se presentan son los correspondientes
a los algoritmos desarrollados en Java. Los diagramas expuestos son los más
importantes, porque de ellos depende la realización del mallado y del análisis. Los
diagramas de flujo presentados corresponden a los métodos:
1. mallaEstructurada
2. mallaNoEstructurada
3. iniciarAnalisis
mallaEstructurada es un método de la clase Malla2D dentro del package
edu.epn.operaciones encargado de realizar el mallado estructurado.
mallaNoEstructurada es un método de la clase Malla2D que realiza el
mallado triangular no estructurado. Este método llama a varias subrutinas como
son calcNodosFrenteArea , getListEsc , buscar , buscarNodo y unico
dentro de la misma clase y a anguloFi , mult , matTransf y buscar de la clase
MatrizOp .
iniciarAnalisis es un método de la clase AnalisisModelo que se encarga
de efectuar el análisis. Durante su operación llama a varias subrutinas como
indicar , calcMatRigSis y calcMatDefEsf para elementos tipo placa y
calcMatDefEsf2D y calcMatRigSis2D para elementos tipo 2D, estas dos
últimas subrutinas no se detallarán por ser muy similares a sus contrapartes para
elementos tipo placa. Adicionalmente se ocupan varios métodos de las librerías
colt y mtj, como se indicó en la sección 4.4.4.
126
Inicio de mallaEstructurada
Se conoce el # de elementos
Si No
Calcula el ancho y alto de la geometría
rectangular
Calcula el número de elementos en X y Y a
partir del área
Calcula el área de la geometría rectangular
Calcula el número de elementos en X y Y
i < # elementos en X + 1
i = 0
No
j = 0
j < # elementos en Y + 1
Si
Agrega las coordenadas del
nuevo nodo
j = j + 1
No
i = i + 1
1A
Gen
erac
ión
de la
s co
orde
nada
s de
los
nodo
s
Si
1. Subrutina mallaEstructurada
127
Elementos de forma rectangular
Si No
Fin de mallaEstructurada
1A
i < # elementos en X
i = 0
No
j = 0
j < # elementos en Y
Si
Agrega nuevo elemento
rectangular
j = j + 1
No
Si
i = i + 1
i < # elementos en X
i = 0
No
j = 0
j < # elementos en Y
Si
Agrega nuevo elemento triangular
inferior
j = j + 1
No
Si
i = i + 1
Agrega nuevo elemento triangular
superior Gen
erac
ión
de lo
s el
emen
tos
128
2. Subrutina mallaNoEstructurada
Inicio de mallaNoEstructurada
Calcular la posición de los nodos del contorno y el
frente de generación inicial calcNodosFrenteArea
Si
2A
Seleccionar segmento del frente de generación
Calcula el área de la geometría
calcNodosFrenteArea
Se conoce la densidad de la
malla
Si No
Calcula el tamaño de los elementos a
partir del área
Lee el tamaño de los elementos
Se asigna la memoria suficiente para
guardar nodos y elementos
Queda algún frente de
generación
No 2B
Calcula el ángulo del segmento
MatrizOp.anguloFi
Pasa el segmento al sistema de
coordenadas locales
2C
Calcula el número inicial de nodos → nxyC
129
Determinación del punto de conexión C1
Paso de C al sistema de coordenadas
global
2C
i = 1
i < 6
Agrega nuevo punto Ci
i = i + 1
No
Si
Det
erm
inac
ión
de lo
s pu
ntos
C
2, C
3, C
4, C
5 y
C6
Calcula el radio en el que se
buscarán nodos
Queda algún nodo
Si
Calcula la distancia desde
C1 hasta el nodo
Distancia < radio
Si No
Calcula el área entre el nodo
actual y el segmento
Area positiva
Si No
Agrega el nodo y su distancia a nf
No
Nod
os p
erte
neci
ente
s al
fr
ente
de
gene
raci
ón
Agrega C1 … C6 a nf → nfC
Ordena los nodos de nfC de acuerdo a su
distancia
Extrae el frente perteneciente a la parte actual getListEsc
Remueve el frente actual de la lista → freP
2D
Obtención del valor de
espaciamiento d
130
2D
Quedan nodos en nfC
Si
No
Calcula el área entre el nodo y el segmento → A0
Si No
Quedan nodos en nf
Si
Es el mismo nodo de nfC
Calcula las áreas entre el nodo de nf y un lado del triángulo
nodo nfC - segmento
Suma de áreas triangulares < A0
No
Calcula el ángulo y longitud del segmento ajM
Si A
Com
prob
ar q
ue n
o ha
yan
nodo
s de
ntro
del
triá
ngul
o A
Baj
Queda algún frente en freP
Si
Obtiene las coordenadas del
frente
El frente no corta el segmento ajM
Si
A
No
No
aj = nodo de nfC
ajM
no
cort
a ni
ngun
a ca
ra
del f
rent
e de
gen
erac
ión
A
Se encontró aj
Si No
Error Busca aj en la lista de nodos
buscar
Se encontró aj
Si No
Agrega aj a la lista de nodos
Crea el nuevo elemento
2E
131
2E
Busca la cara 31 buscar
Se encontró la cara 31
Si No
Borra el frente 31 Agrega el frente 31
Busca la cara 23 buscar
Se encontró la cara 23
Si No
Borra el frente 23 Agrega el frente 23
Borra la cara actual
2A
132
2B
i = nxyC + 1
i < número de nodos
Si
No
Busca los elementos que contiene al nodo i
buscarNodo
idxn = parte de los elementos
# elementos = 3
Si
Halla el vector de nodos unico
No
# elementos = 4
Si No
Añade el nuevo elemento
Agrega posición del nodo a nod
Halla el vector de nodos → V
unico
Añade los nuevos elementos
Agrega posición del nodo a nod
i = i + 1
Opt
imiz
ació
n to
poló
gica
2F
133
Tamaño de nod > 0
Si No
i = tamaño de nod + 1
2F
i > -1
Si
No
Borra los elmentos que contiene a i
Actualiza la numeración de
los nodos
i = i - 1
Act
ualiz
ació
n de
nod
os y
ele
men
tos
i = nxyC + 1 - # puntos de
refinamiento
i < nxyC - 1
Si
No
Busca los elementos que contiene al nodo i
→ elm buscarNodo
idxn = parte de los elementos
Halla el vector de nodos → V
unico
Quedan nodos en V
Si
No
Agrega nuevo nodo entre el punto de refinamiento y un
vértice de las caras
2H
2G
2I
Ref
inam
ient
o de
la m
alla
en
los
punt
os d
e re
finam
ient
o
134
2G
j = 0
j < tamaño de elm
Si
No
Agrega nodo en el centro de la cara
El lado es un borde
Si No
Agrega nuevo nodo del contorno
Agrega 2 nuevos elementos externos
Agrega nuevo nodo del contorno
Agrega 4 nuevos elementos internos
Borra los elementos divididos
2H
j = j + 1
135
2I
Extrae los nodos que no pertenecen al
contorno → nxyNC
i = 0
i < # iteraciones del alisado
Si
No
Queda algún nodo en nxyNC
Si
Busca los elementos que comparten el
nodo → elm buscarNodo
Obtiene el vector de nodos → V
Obtiene el vector de nodos → V
Centra el nodo en el polígono formado por V aplicando el alisado
Laplaciano
No
i = i + 1
Alis
ado
de la
mal
la
Fin de mallaNoEstructurada
136
2.1. Subrutina calcNodosFrenteArea
Queda algún arco
Inicio de calcNodosFrenteArea
Calcula el ángulo inicial, final y paso del arco
Si
No
Agrega las líneas de geometría a la
lista lineas
tam0 <> 0 Si No
t = menor(t, tam0)
Calcula el tamaño del
elemento → t
Agrega el vector de tamaño inicial
tam0 a ta m
Convierte de coordenadas polares a rectangulares → xy
i = 1
i < # filas de matriz xy
Si
Agrega nuevo segmento de
arco a lineas
i = i + 1
No
Con
vier
te lo
s ar
cos
en p
olíg
onos
21A
Agrega t a tam
137
21A
Queda alguna elipse
Calcula el ángulo de paso de la
elipse
Si
No
tam0 <> 0 Si No
t = menor(t, tam0)
Calcula el tamaño del
elemento → t
Convierte de coordenadas polares a rectangulares → xy
i = 1
i < # filas de matriz xy
Si
Agrega nuevo segmento de
elipse a lineas
i = i + 1
No
Con
vier
te la
s el
ipse
s en
pol
ígon
os
21B
Calcula el radio promedio
Agrega t a tam
Extrae el primer vector de lineas y lo agrega a pol
138
21B
Quedan vectores en lineas
Si
No
P = último vector de pol
i = 0
i < longitud de lineas
Agrega vector a pol y lo
remueve de linea
Siguiente vector = p
Si No
Inverso de vector = p
Si No
Agrega inverso de vector a pol y
lo remueve de linea
No
No se agregó ningún vector o
long. linea = 0
Si
No
Calcula el área de pol
Si
Calcula las coordenadas nodales y el frente de
generación inicial de pol
Longitud de lineas > 0
Si
Borra pol y le agrega el primer vector de
lineas
No
21C
Ord
ena
el p
olíg
ono
actu
al
A
A
139
Fin de calcNodosFrenteArea
21C
Tam0 <> 0 Si No
Queda algún punto de refinamiento
Si
Agrega el punto de refinamiento a la lista
de coordenadas nodales
No
Devuelve el área total
Agr
ega
los
punt
os
de r
efin
amie
nto
140
2.2. Subrutina getListEsc
Queda algún vector en la lista
Inicio de getListEsc
Devuelve la selección
Si
No
Fin de getListEsc
Lee lista de vectores, idx , val
Elemento de posición idx del
vector = val
Si No
Agregar vector a la selección
141
2.3. Subrutina buscar
Queda algún vector en la lista
Inicio de buscar
Devuelve -1
Si
No
Fin de buscar
Lee lista de vectores, vec
Vector = vec Si No
Devuelve posición del
vector
A
A
142
2.4. Subrutina buscarNodo
Queda algún vector en la lista
Inicio de buscarNodo
Devuelve la selección
Si
No
Fin de buscarNodo
Lee lista de vectores, nodo
Algún elemento del vector = nodo
Si No
Agrega vector a la selección
143
2.5. Subrutina unico
Queda algún vector en la lista
Inicio de unico
Lee lista de vectores, idx
Si
No
Fin de unico No
Copia todos menos el último elemento del
vector en elems
i = 0, j = 0
i < elems
Si
elems[i] = idx
Si
k = 0
k < j
Si
No elems[i ] = selec[k]
Si
selec[j ] = elems[i]
j = j + 1
k = k + 1
A
No
i = i + 1
A
A
No
i = 0
i < j
Si
No
unicos[i] = selec[i]
i = i + 1
Devuelve unicos
144
2.6. Subrutina MatrizOp.anguloFi
Inicio de MatrizOp.anguloFi
Fin de MatrizOp.anguloFi
Lee matriz de coordenadas
Si No
Devuelve PI/2
Calcula diferencia en X y en Y de las 2
primeras filas → Lx, Ly
Lx = 0 y Ly > 0
Lx < 0 y Ly = 0
No Si
Devuelve PI
Lx = 0 y Ly < 0
No Si
Devuelve 3*PI/2
Lx < 0 No Si
Devuelve atan(Ly/Lx)+PI
Devuelve atan(Ly/Lx)
145
2.7. Subrutina MatrizOp.mult
Inicio de MatrizOp.mult
Devuelve M
Si
No
Fin de MatrizOp.mult
Lee matrices mat1, mat2
sum = 0, i = 0
i < # filas de mat1
j = 0
j < # columnas de mat2
No
Si
k = 0
k < # columnas de mat1
Si
No
sum = sum + mat1[i][k] *
mat2[k][j]
k = k + 1
M[i][j] = sum
sum = 0
j = j + 1
i = i + 1
146
2.8. Subrutina MatrizOp.matTransf
Inicio de MatrizOp.matTransf
Devuelve M
Fin de MatrizOp.matTransf
Lee ángulo Fi
M = {{cos(Fi),-sin(Fi)}, {sin(Fi), cos(Fi)}}
147
2.9. Subrutina MatrizOp.buscar
i < longitud del vector
Inicio de MatrizOp.buscar
Si
No
Fin de MatrizOp.buscar
Lee vector, escalar
Elemento i del vector = escalar
Si No
Devuelve i
i = 0
A i = i + 1
Devuelve -1
A
148
3. Subrutina iniciarAnalisis
Inicio de iniciarAnalisis
Elementos tipo placa
Si No
Ensamblaje de la matriz de rigidez del sistema calcMatRigSis
Liberar memoria para el análisis
Calcula el tiempo de inicio
Indica el tiempo que demoró
indicar
# restricciones r = 0, i = 0
i < # apoyos
Si
Apoyo está restringido
Si No
r = r + 1
No
Ensamblaje de la matriz de rigidez del sistema
calcMatRigSis2D
Indica el tiempo que demoró
indicar
3A
Núm
ero
de r
estr
icci
ones
i = i + 1
# restricciones r = 0, i = 0
i < # apoyos
Si
Apoyo está restringido
Si No
r = r + 1
No
3B
Núm
ero
de r
estr
icci
ones
i = i + 1
149
3A
i < # apoyos
i = 0
No
j = 0
J < 3
Si
j = j + 1
No
Si
i = i + 1
Apoyo no restringido
Si
Agrega la carga al vector de
cargas reducido
Almacena la posición de la
restricción
No
3C
Apl
icac
ión
de la
s re
stric
cion
es a
l ve
ctor
de
carg
as
150
3C
No Si
i < orden matriz reducida
i = 0
No
j = 0
j < orden matriz reducida
Si
Obtiene siguiente elemento de la matriz
del sistema
j = j + 1
No
Si
i = i + 1
Elemento = 0
No Si
Almacena elemento en la matriz reducida
i < orden matriz reducida
i = 0
No
j = 0
j < orden matriz reducida
Si
Obtiene siguiente elemento de la matriz
del sistema
j = j + 1
No
Si
i = i + 1
Elemento = 0
No Si
Almacena elemento en la matriz reducida
i es múltiplo de 500
Si No
Compactar matriz reducida
3E 3D
Apl
icac
ión
de la
s re
stric
cion
es a
la
mat
riz d
e rig
idez
del
sis
tem
a
Resolutor LU
151
Cálculo de los desplazamientos usando colt 1.2.0
Reconstruir el vector de
desplazamientos
Cálculo de los desplazamientos usando mtj 0.9.6
Reconstruir el vector de
desplazamientos
Indica el tiempo que demoró
indicar
Indica el tiempo que demoró
indicar
3D 3E
Precondicionar Si No
Aplicar precondicionador
i < # apoyos
i = 0
No
j = 0
j < 3
Si
Almacena siguiente
desplazamiento
j = j + 1
No
Si
i = i + 1
Indica el tiempo que demoró
indicar
Ens
ambl
aje
de la
mat
riz
de d
espl
azam
ient
os
3F
152
Cálculo de las cargas
3F
i < # apoyos
i = 0
No
j = 0
j < 3
Si
Almacena siguiente carga
j = j + 1
No
Si
i = i + 1
Indica el tiempo que demoró
indicar
Ens
ambl
aje
de la
mat
riz
de c
arga
s
Cálculo de las deformaciones y
esfuerzos calcMatDefEsf
Fin de iniciarAnalisis
Liberar la memoria usada
3G
153
3B
i < # apoyos
i = 0
No
j = 0
J < 2
Si
j = j + 1
No
Si
i = i + 1
Apoyo no restringido
Si
Agrega la carga al vector de
cargas reducido
Almacena la posición de la
restricción
No
3H
Apl
icac
ión
de la
s re
stric
cion
es a
l ve
ctor
de
carg
as
154
3H
Resolutor LU
No Si
i < orden matriz reducida
i = 0
No
j = 0
j < orden matriz reducida
Si
Obtiene siguiente elemento de la matriz
del sistema
j = j + 1
No
Si
i = i + 1
Elemento = 0
No Si
Almacena elemento en la matriz reducida
i < orden matriz reducida
i = 0
No
j = 0
j < orden matriz reducida
Si
Obtiene siguiente elemento de la matriz
del sistema
j = j + 1
No
Si
i = i + 1
Elemento = 0
No Si
Almacena elemento en la matriz reducida
i es múltiplo de 500
Si No
Compactar matriz reducida
3J 3I
Apl
icac
ión
de la
s re
stric
cion
es a
la
mat
riz d
e rig
idez
del
sis
tem
a
155
Cálculo de los desplazamientos usando colt 1.2.0
Reconstruir el vector de
desplazamientos
Cálculo de los desplazamientos usando mtj 0.9.6
Reconstruir el vector de
desplazamientos
Indica el tiempo que demoró
indicar
Indica el tiempo que demoró
indicar
3I 3J
Precondicionar Si No
Aplicar precondicionador
i < # apoyos
i = 0
No
j = 0
j < 2
Si
Almacena siguiente
desplazamiento
j = j + 1
No
Si
i = i + 1
Indica el tiempo que demoró
indicar
Ens
ambl
aje
de la
mat
riz
de d
espl
azam
ient
os
3K
156
Cálculo de las cargas
3K
i < # apoyos
i = 0
No
j = 0
j < 2
Si
Almacena siguiente carga
j = j + 1
No
Si
i = i + 1
Indica el tiempo que demoró
indicar
Ens
ambl
aje
de la
mat
riz
de c
arga
s
Cálculo de las deformaciones y
esfuerzos calcMatDefEsf2D
3G
157
3.1. Subrutina indicar
Indicar tiempo transcurrido
Inicio de indicar
Resta: tiempo actual – tiempo
anterior
Si No
Indica resultado
Calcula nuevo tiempo inicial
Fin de indicar
158
3.2. Subrutina calcMatRigSis
Inicio de calcMatRigSis
Elementos de forma rectangular
Si No
Quedan elementos
Si
No
Obtiene sus propiedades
Calcula su matriz de rigidez
matrizRigidezRec
Lee siguiente elemento
No
Devuelve la matriz de rigidez
del sistema
Fin de calcMatRigSis Calcula su vector de conectividad
Calcula el factor de conversión de
las unidades
Crea la matriz dispersa vacía
Calcula su ancho y longitud
Calcula el tamaño de la matriz de rigidez
del sistema
i = 0
i < 12
Si
j = 0
j < 12
Si
Suma siguiente elemento de la matriz
de rigidez del elemento a la matriz
de rigidez del sistema
j = j + 1
i = i + 1
No
Agr
ega
la m
atriz
de
rigid
ez d
el
elem
ento
a la
del
sis
tem
a
32B
32A
159
Quedan elementos
Si
No
Obtiene sus propiedades
Calcula su matriz de rigidez matrizRigidezTriDKT
Lee siguiente elemento
No
Calcula su vector de conectividad
i = 0
i < 9
Si
j = 0
j < 9
Si
Suma siguiente elemento de la matriz
de rigidez del elemento a la matriz
de rigidez del sistema
j = j + 1
i = i + 1
No
Agr
ega
la m
atriz
de
rigid
ez d
el
elem
ento
a la
del
sis
tem
a
32B
32A
160
3.3. Subrutina calcMatDefEsf
Inicio de calcMatDefEsf
Elementos de forma rectangular
Si No
Quedan elementos
Si
No
Calcula su ancho y longitud
Calcula sus deformaciones, esfuerzos y momentos
flectores matrizEsfuerzoRec
Lee siguiente elemento
Calcula el factor de conversión de
las unidades
Obtiene sus propiedades
i = 0
i < 4
Agrega siguiente desplazamiento
al vector de desplazamientos
i = i + 1
Si
No
Ens
ambl
a el
vec
tor
de d
espl
azam
ient
os
33A
33B
i = 0
i < 4
j = 0
j < 3
Suma siguiente elemento a la matriz de deformaciones,
esfuerzos y momentos flectores
Aumenta en 1 el contador de
deformaciones, esfuerzos y momentos
flectores
j = j + 1
i = i + 1
Si
No
No
Si
Sum
a la
mat
riz d
e de
form
acio
nes,
es
fuer
zos
y m
omen
tos
flect
ores
161
Quedan elementos
Si
No
Calcula sus deformaciones, esfuerzos y momentos
flectores matrizEsfuerzoTriDKT
Lee siguiente elemento
Obtiene sus propiedades
i = 0
i < 3
Agrega siguiente desplazamiento
al vector de desplazamientos
i = i + 1
Si
No
Ens
ambl
a el
vec
tor
de d
espl
azam
ient
os
33B
i = 0
i < 3
j = 0
j < 3
Suma siguiente elemento a la matriz de deformaciones,
esfuerzos y momentos flectores
Aumenta en 1 el contador de
deformaciones, esfuerzos y momentos
flectores
j = j + 1
i = i + 1
Si
No
No
Si
Sum
a la
mat
riz d
e de
form
acio
nes,
es
fuer
zos
y m
omen
tos
flect
ores
33A
162
Fin de calcMatDefEsf
Si
No
i = 0
33B
i < # apoyos
j = 0
j < 3
Promedia las deformaciones
del nodo
Promedia los esfuerzos del
nodo
Promedia los momentos flectores
del nodo
j = j + 1
i = i + 1
Si
No P
rom
edio
de
defo
rmac
ione
s, e
sfue
rzos
y
mom
ento
s fle
ctor
es
