Source: http://docplayer.es/2792013-Universidad-de-chile-facultad-de-ciencias-fisicas-y-matematicas-departamento-de-computacion.html
Timestamp: 2018-02-24 01:27:27+00:00

Document:
UNIVERSIDAD DE CHILE FACULTAD DE CIENCIAS FÍSICAS Y MATEMÁTICAS DEPARTAMENTO DE COMPUTACIÓN - PDF
UNIVERSIDAD DE CHILE FACULTAD DE CIENCIAS FÍSICAS Y MATEMÁTICAS DEPARTAMENTO DE COMPUTACIÓN
Download "UNIVERSIDAD DE CHILE FACULTAD DE CIENCIAS FÍSICAS Y MATEMÁTICAS DEPARTAMENTO DE COMPUTACIÓN"
Santiago Contreras Zúñiga
1 UNIVERSIDAD DE CHILE FACULTAD DE CIENCIAS FÍSICAS Y MATEMÁTICAS DEPARTAMENTO DE COMPUTACIÓN HERRAMIENTA PARA LA SOLUCIÓN DE EDP SOBRE DOMINIOS GENERALES EN 2-D MEDIANTE MÉTODOS ADAPTIVOS MEMORIA PARA OPTAR AL TÍTULO DE INGENIERO CIVIL EN COMPUTACIÓN EDUARDO MERCADER ORTA PROFESOR GUÍA: MARÍA CECILIA RIVARA ZUÑIGA MIEMBROS DE LA COMISIÓN : NANCY HITSCHFELD KAHLER LUIS MATEU BRULÉ Este trabajo ha sido parcialmente nanciado por FONDECYT ( ) SANTIAGO DE CHILE ABRIL, 2012
2 Resumen En este trabajo se discute un sistema experimental para resolver ecuaciones diferenciales parciales (EDP) elípticas o parabólicas sobre dominios en 2-D mediante el método de elementos nitos (MEF). La herramienta combina algoritmos de renamiento y desrenamiento de triangulaciones conformes sobre dominios generales con bordes curvos en las fronteras e interfaces de medios; el MEF, usando elementos lineales; una estructura de datos adecuada; el uso de estimadores del error cometido en la resolución numérica; una interfaz gráca sobre XWindows y un lenguaje declarativo para la denición de los problemas. La herramienta fue desarrollada en forma modular, para permitir la incorporación de nuevas opciones, como elementos de grado mayor en el uso del MEF, y utiliza librerías externas, como son, Sparse y SuperLU para la resolución de los sistemas de ecuaciones poco densos, pdraw para visualización de resultados en 3D y GLADE para construcción de la interfaz gráca; todas estas herramientas corresponden a software de código abierto. La herramienta tiene un uso potencial en una amplia gama de aplicaciones, por ejemplo en cálculo estructural, y mecánica de uidos y permite al usuario manejar en forma exible la adaptabilidad, pudiendo denir o modicar a través de la interfaz gráca, las triangulaciones, nodos, conexiones, moléculas, condiciones de borde y lados curvos. Esta interfaz, también permite al usuario denir o modicar tanto el problema que se desea resolver, como las regiones en las que el usuario desea dirigir ya sea el renamiento como el desrenamiento de la triangulación. Los estimadores de error utilizados permiten crear indicadores que dirigen el renamiento y desrenamiento en forma adaptativa, para mejorar la solución con la menor interacción del usuario. Con ello solo se debe denir una triangulación inicial conforme que representa una malla gruesa y luego por medio de los mecanismos de renamiento explícito o los procesos adaptativo, obtener una triangulación que provea de una malla mas na, que permitirá obtener una solución de mejor calidad. Se ilustra el uso del sistema con problemas de prueba, de solución conocida; se muestra la imagen de la malla inicial del dominio y de algunas iteraciones, la malla y solución asociada. Se concluye que esta herramienta constituye un software general, exible y sencillo de usar para resolver problemas de EDP sobre dominios en 2-D generales. 2
3 Dedicatoria Dedico este trabajo a Patricia, mi esposa, por su amor, paciencia y apoyo incondicional a Eduardo Antonio, mi hijo, por su conanza y fe en mí a Elena, mi madre, por su ejemplo y sus esfuerzos y a María Cecilia, mi profesora guía.
4 Agradecimientos Agradezco a María Cecilia Rivara, a Edgardo Krell, a Patricio Poblete y a todos los que me han apoyaron tanto en NIC Chile, en la Universidad de Chile y en mi familia para poder completar este proyecto. Agradezco también a Fondecyt (Proyecto N ).
5 Índice general Índice general Índice de cuadros Índice de guras i iv v 1. Introducción Antecedentes generales Justicación Objetivo Contenido por capítulos Diseño de la aplicación Objetivos de la aplicación Componentes o módulos Librería de rutinas de renamiento/desrenamiento Resolución de ecuaciones Control de la aplicación Estimadores de error Presentación de resultados Renamiento/Desrenamiento de mallas de triángulos Estructura Molecular Compacta Deniciones Algoritmos de Renamiento Algoritmo de Desrenamiento Lados curvos Bisección de un triángulo por su lado curvo Expresión para hallar el punto medio del lado curvo (primer caso) Expresión para hallar el punto medio del lado curvo (segundo caso) Expresión para hallar el punto medio del lado curvo (tercer caso) Método de Newton-Raphson i
6 4. Método de Elementos Finitos Descripción Problemas elípticos Planteamiento Discretización del problema Elementos triangulares lineales Funciones de forma para Triángulos lineales C Evaluación numérica de las integrales Expresión para el ujo Resolución numérica Problemas parabólicos Planteamiento Discretización del problema Elementos triangulares lineales Funciones de forma para Triángulos lineales C Evaluación numérica de las integrales Expresión para el ujo Resolución numérica Integración de elementos triangulares curvos Estimadores del Error Descripción Indicador del error Denición Aplicación Estimadores del Error Estimador basado en la solución del problema Estimador basado en la solución obtenida Estimador basado en la suavidad de la solución obtenida Problemas de prueba Problema Descripción Resultados Problema Descripción Resultados Problema Descripción Resultados Problema Descripción Aplicación Introducción Ventana principal Ventana de consulta/creación/modicación del problema ii
7 7.4. Ventana de consulta/creación/modicación de la malla Ventana de presentación de soluciones (3DPlot) Otras ventanas Ejemplo de denición de un problema Creación de la malla inicial Denición del problema Conclusiones 67 Bibliografía 68 A. Estructura de datos 69 A.1. Nodo A.2. Lado A.3. Triángulo A.4. Conexión A.5. Molécula A.6. Triangulación A.7. Región A.8. Expresión B. Formato de los archivos de datos 74 B.1. Estructura del archivo de descripción del problema a resolver B.2. Estructura del archivo de descripción de la malla a utilizar C. Expresiones aritméticas 82 iii
8 Índice de cuadros 6.1. Estadísticas ejecución problema Estadísticas ejecución problema Estadísticas ejecución problema B.1. Archivo de descripción del problema (1/3) B.2. Archivo de descripción del problema (2/3) B.3. Archivo de descripción del problema (3/3) B.4. Archivo de descripción de la malla (1/3) B.5. Archivo de descripción de la malla (2/3) B.6. Archivo de descripción de la malla (3/3) C.1. Descripción de los elementos de la gramática de expresiones iv
9 Índice de guras Esquema funcional de la aplicación Ejemplo de una Triangulación Representación de la triangulación conforme de la gura previa Algoritmos de renamiento Propagación del renamiento Triángulo a eliminar Desrenamiento Renamiento de un triángulo por el lado curvo Elemento triangular lineal C 0 y su solución de prueba Flujo a lo largo de la frontera de dos elementos con un lado interno común Flujo concentrado en un punto de la frontera Flujo para un lado interno entre elementos Ejemplo de ujo nodal Malla inicial, problema Iteración 4, problema Iteración 12, problema Iteración 23, problema Malla inicial, problema Iteración 5, problema Iteración 12, problema Malla inicial, problema Iteración 2, problema Iteración 9, problema Dominio D Problema Ventana principal aplicación gxmalla Ventana de visualización de un problema Ventana de visualización de una malla Presentación de la solución Abrir archivo Acerca de Creación de nodos de la malla inicial Creación de la triangulación, las moléculas y los vecinos Creación de las condiciones de borde v
10 7.7.4.Modicación de los nodos de la malla Denición del problema vi
11 Capítulo 1 Introducción 1.1. Antecedentes generales La resolución de ecuaciones diferenciales parciales (EDP) es una necesidad de las más diversas áreas de la Ciencia y la Ingeniería, pues permite entre otros modelar problemas de resistencia de materiales, estructuras, dinámica de uidos o termodinámica. El Método de los Elementos Finitos (MEF) es una técnica que permite resolver EDP o sistemas de EDP que modelan problemas físicos o de ingeniería. Para la utilización del MEF, el usuario debe proveer: una EDP que modele su problema, un dominio acotado sobre el cual se encontrará la solución a la ecuación que modela el sistema, condiciones de borde (CB) sobre el dominio, sean ˆ sobre la función, es decir, CB esenciales (CBE), o ˆ sobre la primera derivada de ésta, es decir, CB Naturales (CBN). En la utilización del método de los elementos nitos clásico, el proceso se divide en las siguientes etapas: En primer lugar se genera una discretización o malla del dominio dado, que corresponde a una colección nita de elementos preseleccionados. Para esto se dene el tipo de elementos, nodos y propiedades geométricas. Por ejemplo, triángulos y/o rectángulos en los dominios de dos dimensiones (2-D) o tetraedros y/o hexaedros en los dominios de tres dimensiones (3-D). En forma más precisa diremos que, en el caso 2-D, los elementos serán triángulos lineales, isoparamétricos, de segundo orden, etc. Los nodos del problema corresponderán a los vértices de la malla y otros puntos auxiliares que permiten calcular funciones de forma sobre el elemento. En segundo lugar se construye el sistema lineal de ecuaciones sobre la discretización. Esta etapa incluye la derivación de las ecuaciones sobre los elementos geométricos, el ensamble de las ecuaciones elementales y la imposición de las condiciones de borde. 1
12 En tercer lugar se encuentra la solución numérica, resolviendo el sistema de ecuaciones ensamblado. Como la precisión de la solución numérica depende de la discretización utilizada, la construcción de la malla es un aspecto clave en la utilización correcta del método de elementos nitos. De acuerdo al error estimado en la solución, se procede a repetir la primera etapa preparando una malla que sea mas na en los zonas en donde éste estime que se requiera mejorar la precisión del resultado. El proceso de preparar una malla mas na se llama renamiento. Los procesos de modicación de las mallas, resultan complejos, dicultando la resolución de los problemas. La complejidad en el uso del MEF, motiva el desarrollo de herramientas y métodos que automaticen o simpliquen su operación. Una solución a este problema es el uso de métodos adaptivos. Al usar software basado en métodos adaptivos, la modicación de la malla se realiza en forma dinámica durante el proceso de resolución del problema. Esto depende en forma crítica de la disponibilidad de algoritmos exibles que permitan modicar la malla, ya sea por renamiento y/o por desrenamiento, de acuerdo al error estimado sobre la discretización. Este trabajo realizado en forma automática. obtendrá resultados que dependerán de que el software incluya estimadores conables del error asociados a la resolución de la ecuación diferencial en consideración. En forma más especíca, el diseño de software adaptivo para análisis de elementos nitos, sigue el siguiente esquema general: En primer lugar se construye una malla inicial gruesa que representa adecuadamente la geometría del problema. Luego se calcula la solución de elementos nitos asociada a ésta malla. Se procede entonces a evaluar esta solución (discretización). Mientras la solución numérica calculada no es aceptable, se repite el siguiente ciclo: ˆ En base a la evaluación de la solución se dene una submalla de renamiento, es decir, se identica el conjunto de elementos de la malla que deberían ser renados para mejorar la solución. ˆ Se utiliza entonces esta información para renar la malla. ˆ Se calcula una nueva solución de elementos nitos, que a su vez se evalúa para decidir si es sucientemente precisa para concluir el proceso. En el caso de los problemas parabólicos, una vez obtenida una solución para un tiempo t i se procede a evaluar el problema para el tiempo t i+1, entonces se repite el siguiente ciclo: ˆ En base a la evaluación de la solución se dene una submalla de desrenamiento, es decir, se identica el conjunto de elementos de la malla que deberían ser desrenados pues su pequeño tamaño no ayuda a mejorar la solución. ˆ Se utiliza entonces esta información para desrenar la malla. 2
13 ˆ Se calcula una nueva solución de elementos nitos, que a su vez se evalúa para decidir si ya no es lo sucientemente precisa para concluir el proceso. Se procede a repetir el ciclo desde 1.1 en la página anterior para el tiempo t i+1. Como la solución numérica obtenida en cada paso depende de la calidad de la malla construida para discretizar el problema continuo, la evaluación de la malla se realiza en la práctica analizando dicha solución numérica. En el contexto adaptivo, el proceso completo se realiza en forma totalmente automática usando un estimador conable del error para evaluar la solución y un criterio que permita decidir qué elementos deberían renarse para obtener una mejor solución (submalla de renamiento). En el sistema descrito en este trabajo, la generación y modicación de las mallas, se hace mediante algoritmos de renación y desrenación de mallas triangulares conformes basados en los algoritmos diseñados y desarrollados por Rivara y Spencer [2][4][3], una primera referencia de este trabajo fue presentada en el Congreso Internacional sobre Métodos Numéricos en Ingeniería y Ciencias Aplicadas celebrado en la Universidad de Concepción, Chile del 16 al 20 de Noviembre de 1992[5] Justicación En el momento del inicio de este trabajo (año 1988), no habían herramientas disponibles para la resolución de EDP sobre dominios generales en 2-D mediante mecanismos adaptativos. El desarrollo inicial evolucionó junto con el desarrollo de las tecnologías, incorporándose nuevos elementos y módulos, como GUI para los usuarios. Si bien, actualmente existen varias herramientas disponibles para la resolución de EDP sobre dominios generales en 2-D 1, la aplicación puede usarse como: prototipo para apreciar el funcionamiento del MEF, plataforma de evaluación de indicadores de error y su impacto en el control de la adaptabilidad, como plataforma para probar nuevos desarrollos con métodos multimalla o como base para el soporte de resolución de problemas en 3-D Objetivo El objetivo del trabajo es desarrollar una herramienta para la solución de EDP sobre dominios generales en 2-D mediante métodos adaptivos, que utilizase los desarrollos realizados en el renamiento y desrenamiento de triangulaciones conformes desarrollado por María Cecilia Rivara para representar los dominios y apoyar el proceso de solución de los problemas elípticos y parabólicos en forma adaptiva. Esta herramienta debería utilizar métodos adaptativos y permitir la evaluación de estimadores de error para estos métodos y ser una aplicación útil para su uso en Ingeniería, que dispusiera de interfaces amigables para el usuario. otras 1 FreeFem++,http://www.freefem.org/ - Matlab, entre 3
14 1.4. Contenido por capítulos En el capítulo 2 se describe la forma en que se diseñó e implementó el sistema. En el capítulo 3 se describen los algoritmos de renación y desrenación. En el capítulo 4 se hace referencia a como se implementa el M.E.F. En el capítulo 5 se presentan los estimadores del error. En el capítulo 6 se presentan problemas test y los resultados obtenidos al resolverlos con el sistema desarrollado. En el capítulo 7se presenta la aplicación, mostrando la interfaz gráca y un ejemplo de como utilizarla para resolver un problema. Para terminar en el capítulo 8 se discute el trabajo realizado. 4
15 Capítulo 2 Diseño de la aplicación 2.1. Objetivos de la aplicación Se denió en primer lugar el rango de aplicaciones que el sistema debería apoyar. Luego se denieron las siguientes estrategias: se estableció que resolvería ecuaciones o sistemas de ecuaciones diferenciales parciales(en adelante E.D.P.) en 2-D mediante elementos nitos triangulares. Esta decisión permite trabajar con dominios generales, obtener una buena precisión en la aproximación y simplicidad en el tratamiento. se estableció que debería utilizar métodos adaptativos y permitir la evaluación de estimadores de error para estos métodos. se estableció que se debería desarrollar una aplicación útil para su uso en Ingeniería, que dispusiera de interfaces amigables para el usuario. se estableció que debería soportar en una segunda etapa la inclusión de métodos de resolución multimallas. se estableció que la aplicación debía ser modular, para facilitar la incorporación de nuevas funcionalidades. se decidió realizar el desarrollo de la aplicación en el lenguaje C, desarrollándolo para sistemas UNIX 1, y que las interfaces grácas se desarrollarían en el sistema X- Window Componentes o módulos A continuación se procedió a la denición de los módulos que permitirían la construcción de la aplicación. Se denieron para ello los siguientes módulos: 1. Librería de rutinas para el renamiento/desrenamiento 2. Resolución de ecuaciones 1 UNIX es una marca comercial registrada de The Open Group en los Estados Unidos y otros paises. 2 X-Window systems, 5
16 Figura 2.2.1: Esquema funcional de la aplicación Definicion Definicion Calculo Evaluacion Presentacion del del de Estimadores de Terminar Dominio Problema Solucion de Error Resultados Desrefina/ Refina 3. Control de la aplicación 4. Estimadores de error 5. Presentación de resultados Basado en lo anterior se denió el esquema de operación de la aplicación. La gura muestra un diagrama lógico que describe la relación entre las diferentes funcionalidades de la aplicación. Finalmente se denieron las estructuras de datos, que se describe en el apéndice A, y el formato de los archivos, que se describe en el apéndice B. Con posterioridad, dada la disponibilidad de nuevas tecnologías, se procedió a agregar una interfaz gráca, que permitiera al usuario un mayor control sobre la operación Librería de rutinas de renamiento/desrenamiento Este módulo corresponde al esqueleto del sistema que permite el renamiento y desre- namiento de triangulaciones sobre dominios generales, usando los algoritmos descritos en el capítulo 3. Para lograr una mejor aproximación en los problemas reales, los algoritmos se generalizaron para manejar triángulos de bordes curvos, pudiendo ser estos actualmente arcos de circunferencia o parábolas, aunque el diseño permite que en el futuro se agreguen nuevas opciones. Los triángulos de bordes curvos se usan solo en las fronteras de los dominios o en las interfaces entre materiales diferentes dentro del dominio. La librerías de renamiento y desrenamiento, corresponden a la implementación de las rutinas y métodos desarrollados por Rivara y Spencer[3]. Estas librerías contienen las rutinas requeridas para el renamiento y desrenamiento de mallas de triángulos conformes, denidas por medio de una estructura de datos denominada Estructura Molecular Compacta (EMC). Todas las rutinas contenidas en las librerías fueron escritas en el lenguaje C. Las rutinas se agrupan de acuerdo a su función, en dos librerías: 1. libref.a agrupa las rutinas usadas para el renamiento de las mallas, 2. libdes.a agrupa las rutinas usadas para el desrenamiento de las mallas. 6
17 Resolución de ecuaciones Sobre el esqueleto que representa la EMC y las librerías de renamiento y desrenamiento se desarrollaron los programas que implementan el MEF, usando en esta etapa funciones de forma lineales, aunque el diseño permite incorporar posteriormente funciones de forma de mayor grado. El módulo del MEF se ha diseñado de modo que la denición del problema sea realizada en forma externa por el usuario, sin necesidad de que éste sepa programar, o requiera recompilar parte de la herramienta para cada problema. Así el usuario puede denir externamente, como parámetros del sistema, todas las funciones que serán utilizadas, usando para ello el lenguaje que se describe en el apéndice C. En base a lo anterior el sistema construye la matriz de rigidez asociada al problema y al lado derecho. Para la resolución del sistema de ecuaciones se han utilizado librerías disponibles en el dominio público para la operación con matrices poco densas, que utilizan un método directo basado en la descomposición LU, que aprovecha las características de baja densidad de elementos de la matriz de rigidez, estas son la librería Sparse 3 y la librería SuperLU 4. Se dispone en forma opcional de un algoritmo de renumeración de nodos para optimizar la resolución del problema. En el futuro podría incorporar el uso opcional de un método multimalla para la resolución del MEF. Las rutinas desarrolladas para la resolución de las ecuaciones se agrupan de acuerdo a su función, en dos librerías: 1. libmef.a agrupa las rutinas usadas para construir la matriz de rigidez, el lado derecho del problema numérico que será resuelto y los indicadores del error. 2. libmat.a agrupa las rutinas que contienen herramientas de cálculo que son usadas en el resto de los módulos Control de la aplicación La aplicación desarrollada originalmente para ser usada en modo terminal a través de la línea de comandos del shell, evolucionó, transformándose en una aplicación gráca que permitió mejorar la usabilidad. El control de esta aplicación se realiza por medio de una interfaz gráca. Esta interfaz fue desarrollada inicialmente usando las librerías de dominio público XForms(The Forms Library) 5 que residen en la librería libforms.a del sistema. Estas librerías proveen funciones que permiten construir interfaces en Xwindows. Dado que dichas librerías dejaron de mantenerse, se reconstruyó la aplicación usando las librerías GTK+ 6 (GIMP Toolkit) y la herramienta para desarrollo rápido multiplataforma GLADE 7 Este módulo permite las siguientes funcionalidades: 3 Sparse 1.3 A Sparse Linear Equation Solver, Kenneth S. Kundert y Alberto Sangiovanni-Vincentelli, University of California, Berkeley. Estas librerías están disponibles a través de netlib ( 4 SuperLU 3.0[9], General purpose library for the direct solution of large, sparse, nonsymmetric systems of linear equations on high performance machines, 5 XForms es Copyright (c) by T.C. Zhao and Mark Overmars. (http://www.nongnu.org/xforms/) 6 GTK+ is free software and part of the GNU Project. The licensing terms for GTK+, are GNU LGPL (http://www.gtk.org/) 7 Glade is Free Software released under the GNU GPL License. (http://glade.gnome.org/) 7
18 Denición de los dominios de los problemas y de la malla inicial. Denición de los problemas Control de la ejecución, permitiendo una operación semiautomática basada en procedimientos adaptativos de control de la calidad de la solución o procedimientos interactivos de renamiento/desrenamiento utilizando regiones previamente conguradas Estimadores de error Una vez resuelto el sistema, se calculan los errores locales de aproximación, sobre cada elemento y se usa esta información tanto para decidir qué elementos de la malla actual deben modicarse, por renación y/o desrenación, como para calcular un estimador global del error que permita evaluar la precisión global de la solución actual. Si la solución actual no es aceptable se modica la malla, para entonces proceder a reconstruir la matriz de rigidez y volver a empezar el ciclo adaptivo. El sistema provee algunos estimadores de error que el usuario de la aplicación puede seleccionar. Se pueden denir nuevos estimadores para el control de la adaptatividad, pero para esto se requiere modicar una de las librerías del sistema Presentación de resultados La presentación de resultados se realiza a través de: 1. la aplicación, por medio de paneles en los que se despliega una imagen gráca de la malla de triángulos y que genera un archivo en postscript una versión para imprimir. 2. un programa de dominio público, 3DPlot 8 que permite mostrar en 3 dimensiones las soluciones obtenidas, entregándole al usuario la posibilidad de realizar acciones sobre la imagen resultante, como son rotaciones o la impresión postscript de la imagen desplegada. 8 Autor: Kenny K.H. Toh 8
19 Capítulo 3 Renamiento/Desrenamiento de mallas de triángulos 3.1. Estructura Molecular Compacta La estructura molecular compacta (EMC), es la estructura de datos usada para representar en forma compacta una secuencia de triangulaciones conformes anidadas. Esta estructura es capaz de soportar los algoritmos de renación y desrenación de las mismas. La estructura de datos se componen de los siguientes objetos geométricos: Nodo, elemento básico de la triangulación. Contiene información sobre: ˆ posición del nodo físico, ˆ posición dentro de la malla, esto es, si está ubicado en el borde o al interior, ˆ condiciones de borde que lo afectan, ˆ tipo de material(para el caso de dominios no homogéneos), ˆ edad relativa al resto de los triángulos, ˆ pertenencia a lado curvo Molécula, asociada a un nodo está formada por el conjunto de sus nodos vecinos Triángulo, grupo de tres nodos vecinos que se denen en base a las moléculas Triangulación, grupo de triángulos generados en una operación de renación o desrenación y se obtienen por bisección de triángulos antiguos, creando nuevos nodos que se agregan a la antigua malla. Lado Curvo, polinomio que permite representar bordes o interfaces en dominios generales. La EMC guarda la información de todas las mallas generadas, almacenando solo la modicación local de la malla actual con respecto a la anterior. En la gura 3.1.2, de la página 10 se ilustra la estructura molecular asociada a la triangulación de la gura 3.1.1, de la página 10. Así, por ejemplo, la molécula asociada al nodo 2 es (3,5,7,1). Observe que la formación de las moléculas se hace en sentido positivo o antihorario. 9
20 Figura 3.1.1: Ejemplo de una Triangulación Figura 3.1.2: Representación de la triangulación conforme de la gura previa T
21 3.2. Deniciones Para usar las triangulaciones con el método de los elementos nitos se exigen las siguientes cualidades: 1. Condición de Conformidad, la intersección de cualquier par de triángulos adyacentes es un lado común o bien un vértice común. 2. Condición de No Degeneración,los ángulos internos de todos los triángulos están acotados inferiormente (por un valor α>0). 3. Condición de Suavidad, para cualquier par de triángulos no disjuntos t 1 y t 2 (cuyos respectivos diámetros son: h 1 y h 2 ) se cumple: min(h 1, h 2 ) max(h 1, h 2 ) > δ (3.2.1) esta condición asegura las propiedades de convergencia del método de elementos nitos y permite el uso de discretizaciones irregulares para manejar ecuaciones parciales con soluciones singulares. Las condiciones anteriores se pueden formalizar con la siguiente denición: Denición 1 : T se dirá una familia (α 0, δ) admisible de triangulaciones irregulares, conformes y no degeneradas, si se cumplen las siguientes condiciones: 1. Para cualquier triangulación τ T y para cualquier par de triángulos adyacentes t 1, t 2 τ, se cumpla la condición de suavidad para la misma constante δ. 2. Si α t es el ángulo interior mínimo del triángulo t, entonces min α t α 0 > 0, sobre todos los triángulos de todas las triangulaciones τ T. Para dar soporte a métodos multimalla se requiere que la secuencia de triangulaciones sea anidada, esto es: Denición 2 :Sean T i y T j dos triangulaciones pertenecientes a la secuencia de mallas, en que (i > j) (T i es mas gruesa que T j ). Sean t i, t j triángulos tales que t i T i y t j T j. Se dice que la triangulación T i y T j son anidadas si: t j T j! t i T i / t j t i Para crear la secuencia de triangulaciones a partir de una triangulación inicial gruesa se requiere proceder a modicar la malla actual de modo que conserve las propiedades antes indicadas. Estos procedimientos se llamaran renamiento y desrenamiento. Denición 3 :Se llama renamiento a la operación que aplicada a una triangulación T i genera una triangulación T i+1 que es más na, esto es, si t i y t i+1 son triángulos de las triangulaciones antes indicadas, entonces: 11
22 considerando a como subconjunto disjunto t i T i no existe t i+1 tal que t i t i+1, Denición 4 : Se llama desrenamiento a la operación que aplicada a una triangulación T i genera una triangulación T i+1 que es más gruesa, esto es, si t i y t i+1 son triángulos de las triangulaciones antes indicadas, entonces: considerando a como subconjunto disjunto t i+1 T i+1 no existe t i tal que t i+1 t i, Denición 5 : Se llama segmento curvo S al conjunto de puntos denido por: S = {(x, y) / (x, y) = ( a x t 3 + b x t 2 + c x t + d x, a y t 3 + b y t 2 + c y t + d y ), t [0, 1] ) Denición 6 : Se llamará triángulo curvo a aquel que posea un lado contenido en un segmento curvo S. El lado contenido se denominará lado curvo. Esto es: ABC curvo lado s AB ABC /s AB S, S segmento curvo 3.3. Algoritmos de Renamiento Para realizar el renamiento de una triangulación T se requiere: 1. Seleccionar los triángulos t a renar. 2. Para cada triangulo t, realizar su división (o bisección). 3. Proceder a dejar la triangulación T conforme. Para realizar la bisección se dispone de dos algoritmos, que fueron bautizados algoritmo de 2 triángulos y algoritmo de 4 triángulos. Estos algoritmos tienen la cualidad de que, en el caso de que la triangulación inicial cumpla con las condiciones indicadas en la denición 1 de la página 11, siempre producirán triangulaciones anidadas que pertenecen a la misma familia (α 0, δ) admisible de triangulaciones irregulares T, donde las constantes α 0 y δ que denen la familia solo dependen de la triangulación inicial T 0. El algoritmo de 2 triángulos, se basa en efectuar bisecciones por el punto medio del lado mas largo, a estas bisecciones las llamaremos bisecciones simples. El algoritmo de 4 triángulos consiste en que además de realizar la bisección por el punto medio del lado mas largo, se procede a realizar la bisección de los dos triángulos generados a partir de nodo creado en la bisección previa, a esta bisección la llamaremos bisección generalizada. En la gura de la página 13 se muestran los dos algoritmos. En el caso de ambos algoritmos, luego de realizar las bisecciones se procede a hacer conforme la malla por los nodos creados. La gura de la página 14 muestra un ejemplo de la propagación del renamiento en una malla, para hacer la malla conforme. El desarrollo de estos algoritmos como la demostración de sus propiedades se puede ver en la referencia[3]. 12
23 Algoritmo de 2 triángulos Triángulo Inicial Figura 3.3.1: Algoritmos de renamiento 3.4. Algoritmo de Desrenamiento Algoritmo de 4 triángulos Consideremos una secuencia anidada de triangulaciones conformes {T i, 0 i r}. Supongamos que se desea desrenar un triángulo t i, i > 0, que llamaremos i-nuevo dado que fue creado por renamiento en la triangulación T i. El desrenamiento de t i deberá revisar todas las triangulaciones T k, i k r, en donde r es la triangulación mas na del dominio. El desrenamiento sobre la secuencia T i, T i+1,..., T r del triángulo t i-nuevo se desarrolla a continuación. Sea τ = { triángulos a ser eliminados por el desrenamiento} 1. Determinar el vértice n del triángulo t a eliminar de la malla(n es un nodo i-nuevo). 2. Determinar el triángulo t τ(t i 1 ) que incluye a t (los casos posibles se pueden observar en la gura de la página 15, en donde t = ABC ). 3. Según el origen del nodo tenemos dos situaciones: a) Si es nodo creado por bisección simple de t gura 3.4.1(a) denir τ = {t } b) Si es nodo creado por bisección generalizada de t gura 3.4.1(b), denir τ = {t, /t, fue formado por la bisección generalizada de t que incluye a n} 4. Denir τ = τ {t, τ(t i ) : t, es no conforme} 5. Para cada t, τ, eliminar todos los triángulos contenidos en t, de todas las mallas posteriores (t, T k, i k r). 13
24 Refinamiento por algoritmo 2 Triángulos 1. Triangulación inicial 2. Bisección 3. Dejo conforme 4. Bisección Refinamiento por algoritmo 4 Triángulos 1. Triangulación inicial 2. Bisección 3. Dejo conforme 4. Bisección 5. Dejo conforme 6. Bisección Figura 3.3.2: Propagación del renamiento 14
25 C Figura 3.4.1: Triángulo a eliminar t* = ABC C t* = ABC t A tt n B A n B (a) (b) La gura de la página 16 muestra una secuencia de desrenamiento de una triangulación, en la que se procede a eliminar el nodo n, y posteriormente los nodos más nuevos n 1 y n 2, para poder así conservar la conformidad de la malla y evitar la existencia de elemento trapezoidales. El algoritmo de desrenamiento, considera la propagación del mismo para asegurar la conformidad de la triangulación. El algoritmo de desrenamiento no corresponde a la aplicación de la operación inversa de la operación de renamiento, en particular esto casi siempre es cierto en el caso del algoritmo de 4 triángulos Lados curvos Para poder resolver problemas sobre dominios irregulares, se ha escogido usar polinomios de grado 3 (cúbicos) para representar los lados curvos, ver deniciones 5 y 6 en la página 12. Los segmentos curvos estarán asociados a uno o más triángulos, de acuerdo a su denición inicial y a los procesos de renamiento y desrenamiento que se realizan durante la búsqueda de una solución al problema. La representación escogida permite aproximar por medio de polinomios de interpolación cualquier forma compleja o conjunto de punto, usando por ejemplo, splines cúbicos Bisección de un triángulo por su lado curvo Para renar un triángulo curvo por su lado curvo se ha denido el siguiente procedimiento (ver gura en la página 17): Sea s AB el lado curvo de un triángulo ABC, para encontrar el punto de bisección: 1. se encuentra el punto medio del segmento recto denido por los puntos A y B extremos del lado S AB 2. se traza la recta perpendicular al punto medio del segmento AB 3. se dene como punto de bisección la intersección entre la recta perpendicular al punto medio del segmento AB y el lado curvo s AB 15
26 Figura 3.4.2: Desrenamiento n 1 n 1 n 2 n n 2 1. Triangulación Inicial 2. Elimino nodo n n 2 3. Elimino nodo n1 4. Elimino nodo n2 16
27 Figura 3.5.1: Renamiento de un triángulo por el lado curvo A A (x, y ) 0 0 (x, y ) m m B B (x, y ) 1 1 C C Expresión para hallar el punto medio del lado curvo (primer caso) Sean A = (x 0, y 0 ) y B = (x 1, y 1 ) los extremos del lado curvo s AB. Consideremos primero que x 0 x 1 y que y 0 y 1. El punto medio del segmento AB es ( x 0 +x1, y 0+y 1 ) 2 2, la pendiente del segmento AB sería m = y 1 y 0 x 1 x 0 por lo que la pendiente de la recta perpendicular será m = x 1 x 0 y 1 y 0. Al conocer la pendiente de la recta perpendicular y un punto de esta (el punto medio de AB) tenemos la ecuación de la recta perpendicular, esta es: y = x 1 x 0 x + (y2 1 y0) 2 + (x 2 1 x 2 0) y 1 y 0 2(y 1 y 0 ) Si usamos la denición 5, tenemos que: a y t 3 + b y t 2 + c y t + d y = x 1 x 0 (a x t 3 + b x t 2 + c x t + d x ) + (y2 1 y0) 2 + (x 2 1 x 2 0) y 1 y 0 2(y 1 y 0 ) de donde obtenemos: 2 (a y (y 1 y 0 ) + a x (x 1 x 0 )) t (b y (y 1 y 0 ) + b x (x 1 x 0 )) t (c y (y 1 y 0 ) + c x (x 1 x 0 )) t+ + 2 (d y (y 1 y 0 ) + d x (x 1 x 0 )) ( (y 2 1 y 2 0) + (x 2 1 x 2 0) ) = 0 (3.5.1) La ecuación 3.5.1es de tercer grado en t, y en nuestro caso es resuelta por un método numérico iterativo, el método de Newton-Raphson. 17
28 Expresión para hallar el punto medio del lado curvo (segundo caso) Consideremos el caso mostrado en la sección , pero en que tenemos que x 0 = x 1, entonces tenemos que el punto medio de AB es ( x 0, y 0+y 1 ) 2, la ecuación del segmento AB sería x = x 0, la recta perpendicular será y = y 0+y 1. Luego tenemos que para encontrar el 2 punto medio del lado curvo debemos resolver: a y t 3 + b y t 2 + c y t + d y = y = y 0 + y 1 2 2a y t 3 + 2b y t 2 + 2c y t + 2d y (y 0 + y 1 ) = 0 (3.5.2) Se puede observar que la ecuación es el caso particular de la ecuación cuando x 0 = x Expresión para hallar el punto medio del lado curvo (tercer caso) Consideremos el caso mostrado en la sección , pero en que tenemos que y 0 = y 1, entonces tenemos que el punto medio de AB es ( ) x 0 +x 1, y 2 0, la ecuación del segmento AB sería y = y 0, la recta perpendicular será x = x 0+x 1. Luego tenemos que para encontrar el 2 punto medio del lado curvo debemos resolver: a x t 3 + b x t 2 + c x t + d x = x = x 0 + x 1 2 2a x t 3 + 2b x t 2 + 2c x t + 2d x (x 0 + x 1 ) = 0 (3.5.3) Se puede observar que la ecuación es el caso particular de la ecuación cuando y 0 = y 1. Es necesario destacar que en los desarrollos anteriores no se consideró el caso x 0 = x 1 e y 0 = y 1, pues no estaríamos en el caso de un triángulo Método de Newton-Raphson Sea x 0 un valor aproximado a la raíz de f(x) = 0, sea h la diferencia entre el valor real de la raíz α y x 0, esto es: α = x 0 + h La expansión de Taylor de segundo grado en x 0 es: f(α) = f(x 0 + h) = f(x 0 ) + hf (x 0 ) + h2 2! f (ξ) en donde ξ = x 0 + θh, 0 < θ < 1,toma un valor entre α y x 0. Si ignoramos el término de segundo grado, y evaluamos f(α) = 0, tenemos: de donde f(x 0 ) + hf (x 0 ) 0 h f(x 0) f (x 0 ) 18
29 de ahí que una mejor aproximación a la raíz de f(x) = 0, sea: f(x 1 ) = x 0 f(x 0) f (x 0 ) lo cual se puede repetir sucesivamente hasta obtener una aproximación sucientemente buena de la raíz de f(x) = 0. Obtenido el valor de t, se procede a obtener por medio de las ecuaciones del lado s AB el punto de bisección. 19
30 Capítulo 4 Método de Elementos Finitos 4.1. Descripción En lo que sigue se desarrollará el método de elementos nitos, para la aplicación considerada, destacando: Descripción de los problemas a resolver. Discretización del dominio, Sistema de ecuaciones para el caso de elementos triangulares lineales El método se desarrollará en forma separada para los problemas elípticos y parabólicos. Es importante destacar que este trabajo no intenta realizar un completo tratamiento del método, más información con respecto del mismo se puede encontrar consultando [6], etc Problemas elípticos Planteamiento Sea Ω un dominio cerrado y acotado en R 2 y Γ = Ω su frontera. Sean Γ E y Γ N subconjuntos de Γ tales que: dos Γ = Γ E ΓN y Γ E ΓN = { 0 o, un conjunto finito de puntos El software esta diseñado para resolver problemas elípticos de la forma : en donde: x (α x u x ) (α y y u y u, f, α x, α y y β son funciones de Ω R. ) + βu = f en Ω u = g E en Γ E u n = g N en Γ N (4.2.1) g E (condición de borde esencial, o de Dirichlet) es una función de Γ E R 20
31 g N (condición de borde natural o de Neumann) es una función de Γ N R n es el vector unitario normal a la frontera del dominio Discretización del problema Consideremos el espacio de Hilbert H0(Ω) 1 con el producto interno usual: (u, v) H H < u, v >= (u(x, y)v(x, y) + u(x, y) v(x, y))dxdy Ω Si multiplicamos la ecuación diferencial por v H, obtenemos: ( x (α u x x ) y (α u y ) + βu)v = fv (4.2.2) y Podemos luego a integrar la ecuación sobre el dominio Ω, obtenemos: ( Ω x (α u x x ) y (α u y ) + βu)v dω = fv dω (4.2.3) y Ω Dado que: x (α u x x )v = x (α u x x v) α u v x x x Obtenemos al introducir en : (4.2.4) Ω ( x (α u x x v) + α u v x x x y (α u y y v) + α u v y + βuv) dω = y y Ω fv dω (4.2.5) Si tomamos el primer y tercer termino de la integral del lado izquierdo tenemos que: [ Ω x (α u x x v) + ( u y (α u αx y v)] dω = v ) x u y α Ω y v dω (4.2.6) y El teorema de Green en el plano dice que si F : Ω Γ R 2 R 2, F H 1 (Ω), entonces: F dω = F n dσ (4.2.7) en donde Γ es la frontera de Ω, y n = Ω ( nx n y Γ ) es la normal exterior unitaria a Γ. Aplicando el teorema de Green en la ecuación obtenemos: Ω u v (α x x x + α u v y + βuv) dω = y y Ω u fv dω + (α x x vn u x + α y y vn y)dσ De lo anterior podemos concluir que si u es una solución de 4.2.1, entonces u también es solución de la formulación débil (variacional): Hallar u H tal que v H a(u, v) = l(v) (4.2.8) 21 Γ
32 en donde: u v a(u, v) = (α x Ω x x + α u v y + βuv) dω y y u l(v) = fv dω + (α x x vn u x + α y y vn y)dσ Ω Γ a(u, v) es una forma bilineal y l(v) es lineal. También se cumple que a y l son continuas, y el problema está bien planteado si α x, α y, β, f L 2 (Ω). La formulación variacional del problema que acabamos de presentar es el punto de partida del método de elementos nitos. El método consiste en encontrar una solución aproximada en un subespacio de dimensión nita. Supongamos que Ω R es discretizado en un conjunto de n t elementos triangulares t i, para i = 1,..., n t, que forma una triangulación T h con n v vértices y n l lados L j, para j = 1,..., n l, que forman la discretización de la frontera Γ. Así: Ω = n t n l t i = t i, Γ = L j t i T h i=1 j=1 Luego, si K es un triángulo t i T h cualquiera, podemos plantear el problema variacional { como: Encontrar uh V h (α Ω x u v + α x x y u v + βuv) dω = fv dω + (α y y Ω Γ x uvn u x x + α y vn y y)dσ v V h en donde: V h = { v; v es continua en Ω, v K P 1 (K) y v = g E en Γ E } (4.2.9) con P 1 (K) el espacio de los polinomios en K de grado menor o igual a 1; que es el espacio de las funciones lineales en K. V h es un espacio de dimensiones nitas. Consideremos el espacio: Ṽ h = {v; v es continua en Ω, v K P 1 (K) V h, ) se tiene que dim (Ṽh = n v y que {φ i } nv i=1denidas por: φ i P 1 (K) y φ i (x j ) = δ ij i, j = 1,..., n v (4.2.10) es una base de Ṽh. De ahí podemos escribir v Ṽh n v v(x) = η i φ i (x), en donde η i = v(x i ) para x Ω. i=1 Hay que destacar que i = 1,..., n v, el soporte de la función φ i se reduce al conjunto de triángulos de los cuales x i es un vértice (el soporte de una función es el conjunto de puntos en donde la función es diferente de 0). Se puede notar que: n v j=1 Sea Ũ K una solución sobre un elemento K, φ j (x) = 1 x Ω 22
33 n v Ũ K = a j φ K j (x, y) (4.2.11) j=1 considerando las funciones de forma denidas en base para el espacio V h, tenemos Ũ K φ K i (α x K x x + α Ũ K φ K i y y y + βũ K φ K i ) dω = fφ K Ũ K Ũ K i dω + (α x K Γ x φk i n x + α y y φk i n y )dσ i = 1,... n v (4.2.12) podemos simplicar la expresión si consideramos las componentes para los ujos en x y en y: τ K x = α x Ũ K x (4.2.13) luego como: τ K y = α y Ũ K y (4.2.14) τ K n = τ K x n K x + τ K y n K y = τ K n en que el subscriptor n indica que el ujo es hacia el interior de la frontera del elemento. Resumiendo tenemos: luego: n v j=1 ( K Γ Ũ K Ũ K u (α x x φk i n x + α y y φk i n y )dσ = φ K j α x x K φ K i x dω + φ K j α y K y fφ K i dω + Matricialmente, estas ecuaciones para el elemento K serían: A K 11 A K 1n v a A K n v1 A K n vn v a = nv en donde: A K ij = Ω K φ K j φ K i α x x x dω + φ K j α y K y L K i = fφ K i dω + K K τ K nφ K i dσ φ K ) i y dω + βφ K j φ K i dω a j = K τ K nφ K i dσ i = 1,..., n v (4.2.15) 23 K L 1. L nv (4.2.16) φ K i y dω + βφ K j φ K i dω (4.2.17) K τ K nφ K i dσ (4.2.18)
34 ~ U y a 2 ~ K U a 1 (x,y ) a 3 3 (x,y ) 3 3 K 1 (x,y ) 1 1 x Figura 4.2.1: Elemento triangular lineal C 0 y su solución de prueba la ecuación matricial representa la discretización del problema variacional para un elemento de la triangulación T h. Para resolver el problema es necesario reunir las ecuaciones de cada elemento, en un gran sistema de ecuaciones, o lo que es lo mismo, una gran ecuación matricial. Esto lo haremos para un tipo de elemento especíco Elementos triangulares lineales Corresponden al más simple de todos los elementos que satisfacen las condiciones de convergencia. A continuación procederemos a: Desarrollar las funciones de forma Preparar las integrales para su evaluación numérica Preparar las expresiones para evaluar el ujo Funciones de forma para Triángulos lineales C 0 La función solución de prueba del elemento K para el problema es: Ũ K = a + bx + cy (4.2.19) polinomio que podemos ver grácamente en la gura de la página
35 necesitamos expresar Ũ K en función de los parámetros del elemento K. Dado que sabemos que: Ũ K (x i, y i ) = a i (4.2.20) si aplicamos en tenemos el siguiente sistema de ecuaciones: a + bx 1 + cy 1 = a 1 a + bx 2 + cy 2 = a 2 a + bx 3 + cy 3 = a 3 si despejamos los términos a, b y c en función de a 1, a 2 y a 3 obtenemos: a = a 1(x 2 y 3 x 3 y 2 ) + a 2 (x 1 y 1 x 1 y 3 + x 3 y 1 ) + a 3 (2x 1 y 1 x 1 y 2 x 2 y 1 ) 2 En donde: = x 1 y 1 1 x 2 y 2 1 x 3 y 3 b = a 1(y2 y3) + a 2 (y 3 y 1 ) + a 3 (y2 y1) 2 c = a 1(x 3 x 2 ) + a 2 (x 3 x 1 ) + a 3 (x 2 x 1 ) 2 = 1 2 [(x 2y 3 x 3 y 2 ) (x 1 y 3 x 3 y 1 ) + (x 1 y 2 x 2 y 1 )] es decir, es el área de elemento K. Al introducir estos resultados en y si reagrupamos términos para obtener la forma de interpolación presentada en , obtenemos que: φ K j (x, y) = a j + b j x + c j y 2 con j = 1, 2, 3 (4.2.21) a j = x k y l x l y k b j = y k y l c j = x l x k En los subscriptores j, k, l toman los valores de 1,2,3 respectivamente para φ K 1 (x, y) y son permutados cíclicamente para φ K 2 (x, y) y φ K 3 (x, y). Finalmente podemos expresar la función de prueba mostrada en , como: Ũ K (x, y) = 3 a i φ K i (x, y) (4.2.22) i=1 25
36 Evaluación numérica de las integrales En lo que sigue desarrollaré la evaluación de cada una de los distintos tipos de integrales que componen la discretización Estas son: 1. Integrales de rigidez y de carga interior (página 26) 2. Integrales de ujo a través de los lados del elemento (página 27) 3. Integrales de ujo en el caso de cargas puntuales (página 28) 4. Integrales de ujo en el caso de lados internos (página 29) Evaluación numérica de integrales de rigidez y de carga interior Para la evaluación numérica de las integrales presentes en necesitamos aproximar las funciones α x, α y, β y f por una constante, válida para el elemento K. Para calcular dichas constantes aproximaremos los valores de las funciones respectivas por el valor que presenten en el centroide del elemento, así: α K x = α x (x c, y c ) α K y = α y (x c, y c ) β K = β(x c, y c ) con: f K = f(x c, y c ) (4.2.23) x c = 1 3 (x 1 + x 2 + x 3 ) y c = 1 3 (y 1 + y 2 + y 3 ) y x i, y i, i = 1, 2, 3, son las coordenadas de los nodos del elemento K. Si consideramos primero las integrales de rigidez que contienen a α x y a α y en y la función de forma descrita en , tenemos: K α x φ K j x φ K i x dω + φ K j α y K y φ K i y dω = αk x b j b i 2 2 K dω + α K y c j c i 2 2 K dω = αk x 4 b jb i + αk y 4 c jc i (4.2.24) Si consideramos la integral de rigidez que contiene a β, dado que las funciones de forma φ K j (x, y) corresponden a las coordenadas de área ζ j de un triángulo y como: ζ1ζ l 2 m ζ n l!m!n! 3 dω = K (l + m + n + 2)! 2 entonces: 26
37 Figura 4.2.2: Flujo a lo largo de la frontera de dos elementos con un lado interno común 4 2 τ n t 2 t ^ s K βφ K j φ K i dω = β K K φ K j φ K i dω = β K K ζ j ζ i dω = { β K 6 i = j β K 12 i j (4.2.25) Si consideramos la integral de carga, utilizando la aproximación de la función f a la constante f K, y aprovechando la fórmula de integración basada en las coordenadas de área, tenemos que: fφ K i dω = f K φ K i dω = f K (4.2.26) 3 K Evaluación numérica de integrales de ujo a través de los lados del elemento Consideremos el caso de dos elementos t 1 y t 2 que compartan un lado común y tengan lados en la frontera de Ω, esto es, Li t j /L i lado de t j, L i Γ; i = 1, 2, 3; j = 1, 2 Sea τ n (ŝ) la componente normal interna del ujo aplicado en la frontera (sea ŝ la coordenada a lo largo de la frontera). Tenemos que τ t 1 n = τ n (ŝ) t1 y que τ t 2 n = τ n (ŝ) t2. Para mayor claridad la gura de la página 27 gráca la situación indicada. La integral del ujo para el elemento t 1 se puede descomponer como: τ t 1 nφ t 1 i dσ = τ t 1 nφ t 1 i dσ + τ t 1 nφ t 1 i dσ + τ t 1 nφ t 1 i dσ t en que los límites 1,2 y 3 corresponden a los nodos del elemento t 1 así, el vector del ujo lo podemos escribir como: K 27
38 2 τ t 1 1 nφ t 1 1 dσ + 3 τ t 1 2 nφ t 1 1 dσ + 1 τ t 1 3 nφ t 1 1 dσ 2 τ t 1 1 nφ t 1 2 dσ + 3 τ t 1 2 nφ t 1 2 dσ + 1 τ t 1 3 nφ t 1 2 dσ 2 τ t 1 1 nφ t 1 3 dσ + 3 τ t 1 2 nφ t 1 3 dσ + 1 τ t 1 3 nφ t 1 3 dσ dado que las funciones de forma toman el valor 0 en el lado opuesto a su nodo la matriz anterior queda como: 2 τ t 1 1 nφ t 1 1 dσ + 1 τ t 1 3 nφ t 1 1 dσ 2 τ t 1 1 nφ t 1 2 dσ + 3 τ t 1 2 nφ t 1 2 dσ 3 τ t 1 2 nφ t 1 3 dσ + (4.2.27) 1 τ t 1 3 nφ t 1 3 dσ En que la primera columna corresponde a integrales de ujo externo y las restantes a integrales de ujos internos. Para evaluar el ujo en los lados de la frontera supondremos que podemos aproximar τ t 1 n como una constante τ t 1 que resulta de evaluar la primera función en el punto medio del lado, esto es: τ t 1 = τ t 1 n(x m, y m ), con x m = x 1 + x 2 e y m = y 1 + y Si aprovechamos la siguiente fórmula de integración para lados de un triángulo: resumiendo: j i ζ l iζ m j dσ = l!m! j (l + m + 1)! L ij, con L ij = dσ i τ t 1 nφ t 1 1 dσ = τ t 1 L 12 2 τ t 1 nφ t 1 2 dσ = τ t 1 L 12 2 (4.2.28) (4.2.29) Lo que en la práctica resulta en que, el ujo sobre el lado, se reparte en partes iguales sobre los nodos del mismo. Evaluación numérica de integrales de ujo en el caso de cargas puntuales En el caso en que se presenten situaciones de ujos concentrados en algún punto de la frontera, como el que muestra la gura de la página 29, entonces podemos representar, en el caso de que sea el nodo 1 el que concentre el ujo, la función del ujo como: τ n (ŝ) = τ 1 δ(ŝ ŝ 1 ) (4.2.30) Si reemplazamos la expresión en las integrales de la primera columna de la matriz obtendremos: que en la práctica aplica todo el ujo al nodo 1. 1 τ t 1 nφ t 1 1 dσ = τ 1 (4.2.31) τ t 1 nφ t 1 2 dσ = 0 (4.2.32) 28
39 2 t 1 1 τ 1 3 ^ s Figura 4.2.3: Flujo concentrado en un punto de la frontera 4 2 t 2 n t 2 n t 1 t Figura 4.2.4: Flujo para un lado interno entre elementos Evaluación numérica de integrales de ujo en el caso de lados internos Para evaluar el caso del ujo interno, consideremos la situación de dos elementos t 1 y t 2 vecinos, esto es que comparten un lado común, como se representa la gura de la página 29. El ujo en estos elementos queda representado por la ecuación para el elemento t 1 y para el elemento t 2, si consideramos a σ, a la coordenada a lo largo del lado del elemento, tenemos que el vector del ujo queda presentado por: 4 τ t 2 2 n t 2 φt 2 2 dσ, + 2 τ t 2 3 n t 2 φt 2 2 dσ, 3 τ t 2 4 n t 2 φt 2 3 dσ, + 2 τ t 2 3 n t 2 φt 2 3 dσ, 4 τ t 2 2 n t 2 φt 2 4 dσ, + (4.2.33) 3 τ t 2 4 n t 2 φt 2 4 dσ, Si ensamblamos las ecuaciones y las ecuaciones obtenemos: 29
40 τ t 1 n t 1 φt 1 1 dσ + 1 τ t 1 3 n t 1 φt 1 1 dσ τ t 1 1 n t 1 φt 1 2 dσ + 3 τ t 1 2 n t 1 φt 1 2 dσ + 2 τ t 2 3 τ t 1 3 n t 1 φt 1 3 dσ + 3 τ t 1 2 n t 1 φt 1 4 τ t 2 2 n t 2 φt 2 4 dσ, + 3 τ t 2 4 n t 2 φt 2 4 dσ, 3 dσ + 2 τ t 2 3 n t 2 φt 2 n t 2 φt 2 2 dσ, + 4 τ t dσ, + 3 τ t 2 4 n t 2 φt 2 2 dσ, n t 2 φt 2 3 dσ, (4.2.34) Pero dada la notación mostrada en la gura de la página 29 tenemos que para el lado común L 23 se cumple : dσ = dσ, n t 1 = n t 2 `τ t 2 n t 2 = `τ t 2 n t 1 Ademas dada las condiciones de continuidad de las funciones de forma, tenemos que: y φ t 1 2 L23 = φ t 2 2 L23 φ t 1 3 L23 = φ t 2 3 L23 Considerando lo anterior en la ecuación para aquellas partes que corresponde al lado común L 23, tenemos para las ecuaciones asociadas al nodo 2 que: 3 τ t 1 2 n t 1 φt 1 2 dσ + 2 τ t 2 3 n t 2 φt 2 2 dσ, = 3 τ t 1 2 n t 1 φt 1 2 dσ + 2 ( τ t 2 φ t 1 3 n t 1 ) 2 ( dσ) = 3 τ t 1 2 n t 1 φt 1 2 dσ + 2 τ t 2 3 n t 1 φt 1 2 dσ = 3 ( τ t 1 τ t 2 n t 2 1 n t 1 )φt 1 2 dσ y para las ecuaciones asociadas al nodo 3 que: (4.2.35) 3 τ t 1 2 n t 1 φt 1 3 dσ + 2 τ t 2 3 n t 2 φt 2 3 dσ, = 3 τ t 1 2 n t 1 φt 1 3 dσ + 2 ( τ t 2 φ t 1 3 n t 1 ) 3 (dσ) = 3 τ t 1 2 n t 1 φt 1 3 dσ + 2 τ t 2 3 n t 1 φt 1 3 dσ = 3 ( τ t 1 τ t 2 n t 2 1 n t 1 )φt 1 3 dσ (4.2.36) Si la carga interior f no está concentrada a lo largo del lado L 23, o es continua a lo largo del lado, tendremos que: τ t 1 n t 1 τ t 2 n t 1 = 0 (4.2.37) si la carga interior f se concentra a lo largo del lado L 23 entonces tendremos que: τ t 1 n t 1 τ t 2 n t 1 = τ 23 (4.2.38) en que τ 23 corresponderá a una condición de borde del problema. 30
41 Expresión para el ujo De las ecuaciones y prepararé las expresiones del ujo, así si t 1 es un elemento tal que t 1 Ω, tenemos que: τ t 1 x (x, y) = α x (x, y) Ũ t 1 (x, y) x τ t 1 y (x, y) = α y (x, y) Ũ t 1 (x, y) y Usando la ecuación tenemos que: (4.2.39) Ũ t 1 (x, y) x = 3 i=1 a i φ t 1 i (x, y) x = 3 i=1 a i b i 2 luego Ũ t 1 (x, y) y = 3 i=1 a i φ t 1 i (x, y) y = 3 i=1 a i c i 2 (4.2.40) τ t 1 x (x, y) = α x (x, y) τ t 1 y (x, y) = α y (x, y) 3 i=1 3 i=1 a i b i 2 a i c i 2 (4.2.41) Ahora bien, para estimar el valor del ujo en el elemento t 1, lo aproximaremos por el valor de las funciones del ujo en el centroide del triángulo. Así el valor de ujo óptimo a través del elemento será: τ t 1 x = α x (x c, y c ) τ t 1 y = α y (x c, y c ) 3 i=1 3 i=1 a i b i 2 a i c i 2 (4.2.42) Para calcular el valor del ujo en los nodos, dado que el ujo a través de los elementos es constante, obtendremos el ujo en el nodo como el promedio de los ujos de los elementos adyacentes a él. En el caso del ejemplo mostrado por la gura de la página 32. tendremos que: τ x 5 = τ y 5 = 4 i=1 4 i=1 τ t i x τ t i y 31
42 ~ τ x ~ ~ τ x t 4 τ x t 3 ~ τ x t 2 ~ τ x 5 ~ τ x t 1 3 y t 3 5 t 2 2 t 4 t x Figura 4.2.5: Ejemplo de ujo nodal 32
43 Resolución numérica Para encontrar la solución al problema de la página 21 es necesario ensamblar las matrices de cada elemento indicada en de la página 23. Para ello se evalúan para cada elemento las ecuaciones de la matriz de rigidez , y del lado derecho De acuerdo a las condiciones de borde del problema se evaluaran las expresiones para el ujo de cada elemento (ecuaciones , , , , y ), añadiéndolos al lado derecho de la ecuación. Una vez evaluados todos los elementos, se procede a resolver numéricamente el sistema de ecuaciones. Obtenido el resultado, se procederá a evaluar los indicadores del error para cada elemento, que serán usados para determinar cual o cuales de ellos deberán ser incluidos en el renamiento de la malla, para una nueva aplicación del método Problemas parabólicos Planteamiento Sea Ω un dominio cerrado y acotado en R 2 y Γ = Ω su frontera. Sean Γ E y Γ N subconjuntos de Γ tales que: dos Γ = Γ E ΓN y Γ E ΓN = { 0 o, un conjunto finito de puntos sea T = [t 0, t f ) R +, con t 0 0 y t > 0. El software esta diseñado para resolver problemas parabólicos de la forma : en donde: γ t U t x (α x U x ) (α y y U y ) + βu = f en (Ω, T ) U = g E en (Γ E, T ) U n = g N en (Γ N, T ) U = u 0 en t = t o (4.3.1) U, f, γ t, α x, α y y β son funciones de (Ω, t) (R, T ) u 0 es una función de Ω R g E (condición de borde esencial, o de Dirichlet) es una función de (Γ E, t) (R, T ) g N (condición de borde natural o de Neumann) es una función de (Γ N, t) (R, T ) n es el vector unitario normal a la frontera del dominio Discretización del problema Consideremos el espacio de Hilbert H0(Ω) 1 con el producto interno usual: (u, v) H H < u, v >= (u(x, y)v(x, y) + u(x, y) v(x, y))dxdy Ω 33
44 Si multiplicamos la ecuación diferencial por v H, obtenemos: U (γ t t x (α u x x ) y (α u y ) + βu)v = fv (4.3.2) y Podemos luego a integrar la ecuación sobre el dominio Ω, obtenemos: U (γ t Ω t x (α u x x ) y (α u y ) + βu)v dω = fv dω (4.3.3) y Ω Si introducimos de la página 21en : U (γ t t v x (α u x x v) + α u v x x x y (α u y y v) + α u v y + βuv) dω = y y Ω Aplicando el teorema de Green en la ecuación obtenemos: Ω U (γ t t v + α u v x x x + α u v y + βuv) dω = y y Ω Γ Ω fv dω (4.3.4) u fv dω + (α x x vn u x + α y y vn y)dσ De lo anterior podemos concluir que si u es una solución de 4.3.1, entonces u también es solución de la formulación débil (variacional): en donde: Hallar u H tal que v H a(u, v) = l(v) (4.3.5) U a(u, v) = (γ t Ω t v + α u v x x x + α u v y + βuv) dω y y u l(v) = fv dω + (α x x vn u x + α y y vn y)dσ Ω Γ a(u, v) es una forma bilineal y l(v) es lineal. También se cumple que a y l son continuas, y el problema está bien planteado si α x, α y, β, f L 2 (Ω). La formulación variacional del problema que acabamos de presentar es el punto de partida del método de elementos nitos. El método consiste en encontrar una solución aproximada en un subespacio de dimensión nita. Supongamos que Ω R es discretizado en un conjunto de n t elementos triangulares t i, para i = 1,..., n t, que forma una triangulación T h con n v vértices y n l lados L j, para j = 1,..., n l, que forman la discretización de la frontera Γ. Así: Ω = n t n l t i = t i, Γ = L j t i T h i=1 j=1 Luego, si K es un triángulo t i T h cualquiera, podemos plantear el problema variacional como: Encontrar u h V h Ω (γ t U v + α t x u v + α x x y u y v y + βuv) dω = Ω fv dω + Γ (α x u x vn x + α y u y vn y)dσ v V h 34
45 en donde: V h = { v; v es continua en Ω, v K P 1 (K) y v = g E en Γ E } (4.3.6) con P 1 (K) el espacio de los polinomios en K de grado menor o igual a 1; que es el espacio de las funciones lineales en K. V h es un espacio de dimensiones nitas. Consideremos el espacio: Ṽ h = {v; v es continua en Ω, v K P 1 (K) V h, ) se tiene que dim (Ṽh = n v y que {φ i } nv i=1denidas por: φ i P 1 (K) y φ i (x j ) = δ ij i, j = 1,..., n v (4.3.7) es una base de Ṽh. De ahí podemos escribir v Ṽh n v v(x) = η i φ i (x), en donde η i = v(x i ) para x Ω. i=1 Hay que destacar que i = 1,..., n v, el soporte de la función φ i se reduce al conjunto de triángulos de los cuales x i es un vértice (el soporte de una función es el conjunto de puntos en donde la función es diferente de 0). Se puede notar que: n v j=1 Sea Ũ K una solución sobre un elemento K, φ j (x) = 1 x Ω n v Ũ K = a j (t)φ K j (x, y) (4.3.8) j=1 considerando las funciones de forma denidas en base para el espacio V h, tenemos Ũ K Ũ K φ K (γ t K t φk i i + α x x x + α Ũ K φ K i y y y + βũ K φ K i ) dω = fφ K Ũ K Ũ K i dω + (α x Ω Γ x φk i n x + α y y φk i n y )dσ i = 1,... n v (4.3.9) podemos simplicar la expresión si consideramos las componentes para los ujos en x y en y: τ K x = α x Ũ K x (4.3.10) luego como: τ K y = α y Ũ K y (4.3.11) τ K n = τ K x n K x + τ K y n K y = τ K n 35
46 en que el subscriptor n indica que el ujo es hacia el interior de la frontera del elemento. Resumiendo tenemos: luego: n v j=1 n v j=1 ( ( K K Γ Ũ K Ũ K u (α x x φk i n x + α y y φk i n y )dσ = ) γ t φ K j φ K daj (t) i dω + dt φ K j α x x K φ K i x dω + φ K j α y K y fφ K i dω + Matricialmente, estas ecuaciones para el elemento K serían: G K 11 G K 1n v..... G K n v1 G K n vn v en donde: A K ij = Ω da 1 dt. da nv dt + K K τ K nφ K i dσ φ K ) i y dω + βφ K j φ K i dω a j (t) = K τ K nφ K i dσ i = 1,..., n v (4.3.12) A K 11 A K 1n v a A K n v1 A K n vn v a nv G K ij = K φ K j φ K i α x x x dω + φ K j α y K y L K i = fφ K i dω + K = L 1. L nv (4.3.13) γ t φ K j φ K i dω (4.3.14) K φ K i y dω + βφ K j φ K i dω (4.3.15) K τ K nφ K i dσ (4.3.16) la ecuación matricial representa la discretización del problema variacional para un elemento de la triangulación T h. Para resolver el problema es necesario reunir las ecuaciones de cada elemento, en un gran sistema de ecuaciones, o lo que es lo mismo, una gran ecuación matricial. Esto lo haremos para un tipo de elemento especíco Elementos triangulares lineales Corresponden al más simple de todos los elementos que satisfacen las condiciones de convergencia. A continuación procederemos a: Desarrollar las funciones de forma Preparar las integrales para su evaluación numérica Preparar las expresiones para evaluar el ujo 36
47 Funciones de forma para Triángulos lineales C 0 La función solución de prueba del elemento K para el problema 4.3.6, se obtiene de una manera similar a la realizada en la sección de la página 24, y esta es: Ũ K (x, y, t) = 3 a i (t)φ K i (x, y) (4.3.17) i= Evaluación numérica de las integrales En lo que sigue desarrollaré la evaluación de cada una de los distintos tipos de integrales que componen la discretización Estas son: 1. Integrales de rigidez y de carga interior 2. Integrales de ujo a través de los lados del elemento 1 3. Integrales de ujo en el caso de cargas puntuales 2 4. Integrales de ujo en el caso de lados internos 3 Evaluación numérica de integrales de rigidez y de carga interior Para la evaluación numérica de las integrales presentes en , y necesitamos aproximar las funciones γ t, α x, α y, β y f por una constante, válida para el elemento K. Para calcular dichas constantes aproximaremos los valores de las funciones respectivas por el valor que presenten en el centroide del elemento, así: γ K t = γ t (x c, y c ) α K x = α x (x c, y c ) α K y = α y (x c, y c ) β K = β(x c, y c ) con: f K = f(x c, y c ) (4.3.18) x c = 1 3 (x 1 + x 2 + x 3 ) y c = 1 3 (y 1 + y 2 + y 3 ) y x i, y i, i = 1, 2, 3, son las coordenadas de los nodos del elemento K. Si consideramos primero las integrales de rigidez que contienen a α x y a α y en y la función de forma descrita en , de la página 25, tenemos: 1 Este caso es idéntico al desarrollado en de la página 27 2 Este caso es idéntico al desarrollado en de la página 28 3 Este caso es idéntico al desarrollado en de la página 29 37
48 K α x φ K j x φ K i x dω + φ K j α y K y φ K i y dω = αk x b j b i 2 2 K dω + α K y c j c i 2 2 K dω = αk x 4 b jb i + αk y 4 c jc i (4.3.19) Si consideramos la integral de rigidez que contiene a β, dado que las funciones de forma φ K j (x, y) corresponden a las coordenadas de área ζ j de un triángulo y como: ζ1ζ l 2 m ζ n l!m!n! 3 dω = K (l + m + n + 2)! 2 entonces: K βφ K j φ K i dω = β K K φ K j φ K i dω = β K K ζ j ζ i dω = { β K 6 i = j β K 12 i j (4.3.20) Se puede apreciar que la matriz G es similar a la evaluación a la realizada en el cálculo de las integrales de β. De ahí tenemos que: K γ t φ K j φ K i dω = γ K t K φ K j φ K i dω = γ K t K ζ j ζ i dω = { γ K 6 i = j γ K 12 i j (4.3.21) Si consideramos la integral de carga, utilizando la aproximación de la función f a una constante f K denida en , y aprovechando la fórmula de integración basada en las coordenadas de área, tenemos que: fφ K i dω = f K φ K i dω = f K (4.3.22) 3 K Evaluación numérica de integrales de ujo La evaluación numérica de las condiciones de borde del problema sigue el mismo desarrollo que lo revisado en el caso elíptico, ver páginas 29 y siguientes. Obteniendo como resultado las ecuaciones , , , , y Expresión para el ujo De las ecuaciones y prepararé las expresiones del ujo, así si t 1 es un elemento tal que t 1 Ω, tenemos que: τ t 1 x (x, y, t) = α x (x, y, t) Ũ t 1 (x, y, t) x K τ t 1 y (x, y, t) = α y (x, y, t) Ũ t 1 (x, y, t) y (4.3.23) 38
49 Siguiendo un desarrollo similar al mostrado en la sección de la página 31 obtenemos que para estimar el valor del ujo en el elemento t 1, lo aproximaremos por el valor de las funciones del ujo en el centroide del triángulo. Así el valor de ujo óptimo a través del elemento será: τ t 1 x (t) = α x (x c, y c, t) τ t 1 y (t) = α y (x c, y c, t) 3 i=1 3 i=1 a i (t) b i 2 a i (t) c i 2 (4.3.24) Para calcular el valor del ujo en los nodos, el procedimiento es similar al mostrado en la sección , puesto que dado que el ujo a través de los elementos es constante, obtendremos el ujo en el nodo como el promedio de los ujos de los elementos adyacentes a él Resolución numérica Para encontrar la solución al problema de la página 34 es necesario ensamblar las matrices de cada elemento indicada en de la página 36. Para ello se evalúan para cada elemento las ecuaciones de la matriz de rigidez , , , del lado derecho y de las condiciones de borde (página 38)., y se deben añadir al lado derecho de la ecuación. Una vez evaluados todos los elementos, tenemos que la ecuación corresponde a una EDO (Ecuación Diferencial Ordinaria), en el tiempo. Para resolver el problema para un tiempo t = n + 1, dado que ya lo tenemos resuelto para un tiempo t = n (inicialmente t = t 0, condición inicial del problema) procedemos con una aproximación implícita con peso variable para el tiempo, de modo que en la ecuación inicial: G a t + Aa = L aproximaremos en un paso de tiempo t por un parámetro θ tendremos: { } a G + A {a} t θ = {L} θ θ en que θcontrola la estabilidad de la solución y en donde: /θ [0, 1], de modo que θ = t t n 1 t Así las aproximaciones son: con t = t n t n 1 {a} θ (1 θ){a} n 1 + θ{a} n Derivando obtenemos: {L} θ (1 θ){l} n 1 + θ{l} n 39
50 luego, { } a θ t θ t {a} n 1 + θ t {a} n = {a} n {a} n 1 t así, G {a} n {a} n 1 t + A ((1 θ){a} n 1 + θ{a} n ) = (1 θ){l} n 1 + θ{l} n ( ) ( ) 1 1 [G] + θ[a] {a} n = (1 θ){[l]} n 1 + θ{[l]} n + [G] + (1 θ)[a] {a} n 1 t t es la ecuación que resolveremos numéricamente. En este caso se trata de resolver un sistema de la forma: En donde: y [A] {a} n = {f} [A] = 1 [G] + θ[a] t ( ) 1 {f} = (1 θ){[l]} n 1 + θ{[l]} n + [G] (1 θ)[a] {a} n 1 t Obtenido el resultado, se procederá a evaluar los indicadores del error para cada elemento, que serán usados para determinar cual o cuales de ellos deberán ser incluidos en el renamiento de la malla, para una nueva aplicación del método Integración de elementos triangulares curvos. Los bordes curvos se pueden denir por medio de dos mecanismos, circunferencias y splines cúbicas. Las circunferencias, se denen como: (x, y) s curvo si, (x a) 2 + (y b) 2 = r 2 en que (a, b) corresponde al centro de la circunferencia y r corresponde al radio de la misma. Las splines se denen por la mejor spline que representa el borde curvo que se desea representar, para ello los bordes se representan por polinomios de tercer grado en x e y, en que: (x, y) s curvo si, a xt 3 + b x t 2 + c x t + d x = 0 a y t3 + b y t 2 + c y t + d y = 0, t [0, 1] en que los extremos del arco corresponden a los valores extremos de t. 40
51 Capítulo 5 Estimadores del Error 5.1. Descripción Para revisar la calidad de la solución y poder efectuar un adecuado control adaptivo de la malla a usar, se necesita descubrir cuales de los elementos de la malla requieren ser modicado, esto es ser renados o desrenados. Para identicar los elementos, deniré un indicador del error cometido, que permita discriminar la continuidad del procedimiento de renamiento o desrenamiento aplicado. El indicador del error en los elementos de la malla se construye en base a estimadores del error local cometido en la solución. Estos estimadores a posteriori del error, permiten hallar una solución óptima de acuerdo al criterio usado para denir el estimador. Para ello se modica (rena/desrena), la malla usada para resolver el problema planteado. Este trabajo utiliza estimadores locales del error cometido, los que se miden sobre cada uno de los elementos de la malla y no sobre el dominio en su conjunto. De esta forma se realizan los renamientos o desrenamientos en las zonas en que dicho error estimado sea mayor, para así mejorar la solución obtenida A continuación se presentaran el indicador y los estimadores usados Indicador del error Denición Sea I (t)el indicador del error local cometido en el elemento, tal que: I : {t T } [0, 1] si Ω es el dominio del problema, ε el estimador del error local cometido, T la triangulación denida sobre el dominio y t un elemento de la triangulación. Deno: ε m = máx t T ε (t) entonces el Indicador del error sobre los elementos se dene como: I (t) = ε (t) ε m, t T 41
52 Aplicación El indicador del error local cometido se usará para guiar el proceso de modicación de la malla. Para ello se denen valores límite que indicarán que durante el proceso de renamiento si I(t) > limitesup t es refinado y para los desrenamientos: I(t) < limiteinf t es desrefinado 5.3. Estimadores del Error Los estimadores del error local cometido son funciones tales que: ε (t) : {t T } R + En este trabajo se compararon tres estimadores del error diferentes, tanto en su formulación, como en su espíritu o justicación, y se presentan a continuación Estimador basado en la solución del problema. Este estimador se basa en la solución del problema. Sea u r la solución al problema, u e la solución estimada, denimos al estimador como: u r (v i ) u e (v i ) ε (t) = v i t u r (v i ) Estimador basado en la solución obtenida Este estimador se basa en la solución obtenida. Para ello se estima el error cometido de la forma siguiente : 3 ɛ t = At [α Js 2 + βl t ] i=1 en donde A t es el área del triángulo t, J s es una función relacionada con el salto de la derivada normal de la solución en curso de elementos nitos a lo largo del lado s,α y β son constantes y: L t = (βu f) Estimador basado en la suavidad de la solución obtenida t Este estimador se basa en la suavidad de la solución encontrada, busca disminuir la diferencia en la pendiente entre los triángulos con lados comunes. Este estimador se dene en la forma siguiente: Sea t T un elemento de la triangulación denida sobre el dominioω. 42
53 Sean a, b, c los vértices del elemento. Sean T a (t) = {t a T t a t = xy, x, y {a, b, c} x y} los elementos que son contiguos al elemento t y que comparten un lado con este. Entonces deniremos el estimador como: ε (t) = t a T a (t) pendiente (t, t a) en que: pendiente (t, v) = { 0 si t v Ø 2 si cos(v, t) < cos(angmax), t T, v T a (t) 43
54 Capítulo 6 Problemas de prueba Para ilustrar el uso y evaluar el comportamiento del sistema, se presentan cinco problemas de prueba con solución conocida. Se utilizó para estimar el error una medida de la suavidad de la solución obtenida entre elementos vecinos, y así poder construir la solución numérica adaptiva. La solución nal de cada problema mostrado, se ajusta a la solución real de los problemas planteados Problema Descripción En el primer ejemplo se trata de resolver el problema : u = 18π 2 sin(3πx) sin(3πy), en Ω = [0, 1] [0, 1], u = 0 en Ω cuya solución exacta es la función: u = sin(3πx) sin(3πy) Resultados Las guras a 6.1.4, muestran la secuencia de mallas del renamiento 4, 12 y 23, y grácas de la solución hasta la de última malla, en base a un estimador basado en la suavidad de la solución. La tabla 6.1, muestra información estadística sobre la ejecución de la aplicación durante la resolución de este problema. 44
55 1.000 Y X Figura 6.1.1: Malla inicial, problema 1 Iteración Cantidad Errores Valor Duración de nodos Máximo Medio Cuadrático medio máximo real[seg] Cuadro 6.1: Estadísticas ejecución problema 1 45
56 1.000 Y X (a) Malla (b) Solución Figura 6.1.2: Iteración 4, problema 1 46
57 1.000 Y X (a) Malla (b) Soluciones Figura 6.1.3: Iteración 12, problema 1 47
58 1.000 Y X (a) Malla (b) Solución Figura 6.1.4: Iteración 23, problema 1 48
59 6.2. Problema Descripción En el segundo ejemplo se trata de resolver el problema : u = f, en Ω = [0, 1] [0, 1], u = 0 en Ω donde la solución exacta es la función: u = x(x 1)y(y 1)e 100((x 0,5)2 +(y 0,117) 2 ) La Figura muestra la secuencia de mallas del renamiento, y grácas de la solución encontrada en las mallas 5 y 12, con vistas desde dos perspectivas diferentes para apreciar mejor la forma de la solución Resultados Las guras a 6.2.3, muestran la secuencia de mallas del renamiento, y grácas de la solución encontrada para las mallas, en base a un estimador basado en el error real. La tabla 6.2, muestra información estadística sobre la ejecución de la aplicación durante la resolución de este problema. Iteración Cantidad Errores Valor Duración de nodos Máximo Medio Cuadrático medio máximo real[seg] Cuadro 6.2: Estadísticas ejecución problema 2 49
60 1.000 Y X Figura 6.2.1: Malla inicial, problema Problema Descripción En el tercer ejemplo se trata de resolver el problema : u = f, en Ω = {(x, y)/0,5 x 2 + y 2 1}, u = 0 en Ω e = {(x, y)/x 2 + y 2 = 1}, u = cos( π) 8 en Ω i = {(x, y)/x 2 + y 2 = 0,25}, donde la solución exacta es la función: u = cos( π(x2 +y 2 ) 2 ) La Figura muestra la secuencia de mallas del renamiento, y grácas de la solución encontrada en las mallas 2 y 9, con vistas desde dos perspectivas diferentes para apreciar mejor la forma de la solución Resultados Las guras a 6.2.3, muestran las mallas inicial, 2 y 9 del renamiento, y grácas de la solución encontrada para las mallas, en base a un estimador basado en el error estimado. La tabla 6.3, muestra información estadística sobre la ejecución de la aplicación durante la resolución de este problema. 50
61 1.000 Y X (a) Malla (b) Solución Figura 6.2.2: Iteración 5, problema 2 51
62 1.000 Y X (a) Malla (b) Solución Figura 6.2.3: Iteración 12, problema 2 52
63 Iteración Cantidad Errores Valor Duración de nodos Máximo Medio Cuadrático medio máximo real [seg] , , , , , Cuadro 6.3: Estadísticas ejecución problema Y X Figura 6.3.1: Malla inicial, problema 3 53
64 1.000 Y X (a) Malla (b) Solución Figura 6.3.2: Iteración 2, problema 3 54
65 1.000 Y X (a) Malla inicial (b) Solución Figura 6.3.3: Iteración 9, problema 3 55
66 6.4. Problema Descripción En el cuarto ejemplo se trata de resolver el problema : u = 0, en D u(x, y) = 4 en L6 L7 u(x, y) = x en L2 L4 n u(x, y) = y en L5 n x+y (x, y) = 2 y en L1 L3 u n donde la solución exacta es la función: u(x, y) = (x y) + 4 La gura 6.4.1muestra el dominio D y los bordes L1 a L7 en que se aplican las condiciones de borde. Figura 6.4.1: Dominio D La gura muestra la malla inicial y la solución del problema. 56
67 0.400 Y X (a) Malla inicial (b) Solución Figura 6.4.2: Problema 4 57
68 Capítulo 7 Aplicación 7.1. Introducción La aplicación desarrollada opera sobre equipos Linux, utilizando la interfaz gráca XWindows. Como fue desarrollada sobre GLADE, debiera ser sencillo portarla a otras plataformas. La aplicación permite crear, consultar y/o modicar los problemas y sus mallas asociadas, controlar el renamiento, desrenamiento de las mallas y la aplicación del método de elementos nitos sobre un problema. La interacción del usuario con la aplicación se realiza principalmente a través de: Ventana principal o de control Ventana de consulta/creación/modicación del problema Ventana de consulta/creación/modicación de la malla Ventana de presentación de soluciones (3DPlot) Otras ventanas Cada ventana tiene asociadas menús de control y botones para la interacción del usuario Ventana principal La ventana principal presenta un menú superior que está separado de acuerdo a las funciones principales de la aplicación. Estos son: 1. Problema 2. Malla 3. Rena 4. Desrena 5. M.E.F. 6. Spline 58
69 7. Ayuda Además presenta un panel para el despliegue del dominio y su malla asociada, y un segundo panel destinado a presentar información de la aplicación. Figura 7.2.1: Ventana principal aplicación gxmalla El ítem problema dispone de las siguientes opciones: 1. Nuevo 2. Abrir 3. Guardar 4. Guardar como 5. Salir 7.3. Ventana de consulta/creación/modicación del problema La ventana del problema permite consultar, crear o modicar la especicación de un problema. Allí se aprecian los principales parámetros de control del problema, como la descripción de la ecuación diferencial como también de las regiones de renamiento y parámetros que denirán como operará la aplicación al resolver el problema Ventana de consulta/creación/modicación de la malla La ventana de la malla permite consultar, crear o modicar la especicación de una triangulación sobre el dominio especicado. En la práctica, el dominio está representado por la malla de triángulos. A través de esta interfaz, se puede acceder a la información de: 59
70 Figura 7.3.1: Ventana de visualización de un problema 1. los nodos, 2. las triangulaciones 3. los bordes curvos (splines) 4. las condiciones de borde 5. las coordenadas del dominio Figura 7.4.1: Ventana de visualización de una malla 60
ERNESTO ARANDA. 1. Introducción 1 1.1. Ejemplo... 2
EL MÉTODO DE LOS ELEMENTOS FINITOS DESDE UNA PERSPECTIVA MATEMÁTICA ERNESTO ARANDA Índice. Introducción.. Ejemplo....................................... 2 2. Formulación variacional de problemas de contorno

References: resolución 
 resolución 
 Resolución 
 Resolución 
 Resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 Resolución 
 Resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 Resolución 
 Resolución 
 resolución 
 resolución 
 resolución