]]>
Figura 1. Localización del dominio de la malla no estructurada Figure 1. Domain localization of the unestructured mesh
UNESTRUCTURED MESH GENERATION FOR NUMERIC MODELS IMPLEMENTATION
FRANCISCO GARCÍA
Ing. Ambiental, M.Sc, Docente, Universidad del Magdalena, francisco.garcia@une.net.co
CARLOS PALACIO
Ingeniería, Ph.D, Docente, Universidad de Antioquia, cpalacio@udea.edu.co
URIEL GARCÍA
Biólogo, Investigador, Universidad del Magdalena, uriel.garcia@yahoo.com
]]> Recibido para revisar marzo 5 de 2008, aceptado octubre 10 de 2008, versión final noviembre 12 de 2008
RESUMEN: Se usó una herramienta numérica para la generación automática de mallas no estructuradas en la implementación del método de los elementos finitos. Las mallas generadas son usadas en el desarrollo de modelos numéricos dos dimensionales (2-D) en la simulación de fenómenos oceánicos. Para ello, se desarrolló un paquete en Matlab basado en el refinamiento triangular no estructurado de Delaunay, complementado por un conjunto de utilidades para chequear y mejorar la calidad de la malla.
PALABRAS CLAVE: Generación de mallas, modelación oceánica, Elementos finitos, Refinamiento de mallas, triangulación de Delaunay, Matlab.
ABSTRACT: A numerical tool was used for the automatic generation of unstructured meshes in the finite element method implantation. The generated meshes are used in the development of two dimensional (2-D) numerical models in the simulation of oceanic phenomena. For it, a package developed in Matlab based on Delaunay unstructured triangular grid refinement algorithms was used and complemented by a set of diagnostic utilities to check and improve grid quality.
KEYWORDS: Mesh generation; Ocean modeling; Finite elements; Bathymetry; Grid refinement; Delaunay triangulation; Matlab.
1. INTRODUCCIÓN
Las ecuaciones diferenciales planteadas para simulaciones de fenómenos físicos, requieren para su solución la implementación de métodos numéricos, muchos de los cuales asumen que el dominio de interés puede ser dividido en pequeños elementos cuya unión y continuidad forman una malla [1]. La discretización del
dominio permite que las ecuaciones de gobierno puedan ser resueltas numéricamente con aproximaciones a las ecuaciones de conservación que dan un estimativo de las variables de interés. Estas ecuaciones diferenciales son linealizadas dentro de los elementos de control, donde las mismas son aplicadas [2]. Hay dos tipos de mallas; las estructuradas y las no estructuradas [1]. Una malla estructurada se define como aquella donde cada elemento de control tiene el
]]> mismo número de elementos vecinos [2-3], mientras que en una malla no estructurada estos tienen diferentes números de elementos vecinos [2]. La solución de las mallas estructuradas comúnmente se implementan mediante el método de las diferencias finitas, mientras que el método de los elementos finitos se usa generalmente para las mallas no estructuradas. En las estructuradas, que ofrecen ciertas ventajas sobre las no estructuradas, cada vértice de la malla excepto en las fronteras, tiene un vecino isomorfo. Su implementación es simple y muy conveniente para la solución de diferencias finitas, requieren menos memoria de computo y ofrecen la posibilidad de controlar directamente la forma y el tamaño de cada elemento. La gran desventaja de este tipo de mallas, es su poca flexibilidad para ajustarse a dominios con geometrías complicadas, lo que las hacen poco atractivas para efectuar simulaciones en lugares donde abunden islas, estrechos y en general líneas costeras irregulares. Cuando las mallas estructuradas se basan en coordenadas geográficas y son implementadas para modelos en las cercanías del polo norte aparecen serios problemas de estabilidad [4-5]. A partir de estas dos complicaciones expresadas por las mallas estructuradas, aparecen las ventajas del método de elementos finitos y las mallas no estructuradas, usadas en el mismo. Los elementos finitos tienen una inherente habilidad para tratar con mallas no estructuradas, estas permiten fácilmente refinamientos locales para dar resoluciones altas en regiones de interés, sin perder precisión. En ellas se logra una adecuada representación de líneas costeras en geometrías irregulares, lo cual sin duda es una ventaja frente al uso de mallas estructuradas en el método de las diferencias finitas [5].Crear mallas no estructuradas requiere la selección de los puntos internos de la malla (vértices de nodos) y su triangulación [6]. La resolución debe ser tal que permita a ésta, tener pequeños elementos para resolver los detalles finos de la geometría e igualmente debe poseer elementos grandes en el área del dominio donde sea posible, para reducir el número total de nodos. Para la solución numérica es necesario que los elementos de la malla tengan una alta calidad a fin de asegurar una adecuada discretización y buenas condiciones en el sistema lineal [6].
En los últimos años se han desarrollado diferentes algoritmos para la generación de mallas no estructuradas, detalles de estos se pueden encontrar en diferentes publicaciones [6-11]. Uno de los más populares es el algoritmo de Delaunay, en el cual se parte de una triangulación inicial y este la refina hasta que la calidad de los elementos que se forman es suficientemente alta.
Numerosos autores han demostrado que una buena discretización del domino es un paso esencial para lograr una mejor aproximación al fenómeno físico del flujo de marea y circulación en el dominio [12-15]. Se encontró que la calidad de la solución de los métodos computacionales usados, tiene una fuerte relación de dependencia de la resolución (el número de elementos que representan el área de interés) de la malla [16]. La resolución resultante de la discretizacion del dominio no puede ser separada de la respuesta computacional del modelo [17], de allí que existe la necesidad de crear mallas que representen adecuadamente el dominio y que también puedan contar con una calidad de elementos que favorezca la solución computacional.
El método de elementos finitos ha evolucionado rápidamente en las últimas dos décadas para permitir el control de la calidad de los elementos de la malla en el proceso de refinamiento. Las metodologías de refinamiento de mallas no estructuradas están bien documentadas en la literatura, para ello puede revisarse [5], [18-23]. Legrand (2000) mostró que es posible refinar la malla en función de la propagación de las variables de interés en la aplicación de modelos hidrodinámicos de gran escala [5]. Wang (2007), presentó los resultados de un algoritmo que usa la generación de frentes de avance circular para originar mallas de buena calidad [23].
En este artículo se muestran los resultados de la aplicación de un algoritmo para la construcción de mallas no estructuradas para la simulación hidrodinámica en zonas costeras mediante la implementación del método de los elementos finitos.
2. METODOLOGIA
2.1 Dominio ]]>
Figura 1. Localización del dominio de la malla no estructurada
Figure 1. Domain localization of the unestructured mesh
La longitud de la curva G en el espacio físico y en el espacio dual, respectivamente, están dadas por las expresiones:
De acuerdo a ello, la meta de un generador de mallas triangulares de elementos finitos es lograr que los lados de los triángulos que conforman los elementos sean de la misma longitud. Para lograr esta meta los vértices deben ser creados, removidos y reemplazados hasta satisfacer esta condición.
Los puntos fijados en el dominio de simulación son triangulados con el algoritmo de Delaunay [5]. La conexión entre cada par de nudos tiene una relación de fuerza de desplazamiento 
que depende de sus longitudes actuales 
e iníciales
. En éste sistema las fuerzas externas vienen de las fronteras. En cada nodo de la frontera, hay una fuerza de reacción normal, la magnitud de esta fuerza es lo suficientemente grande para mantener a la frontera moviéndose desde afuera. La posición de los nodos, la cual es la principal variable desconocida, es encontrada resolviendo el equilibrio de fuerzas en estado estático. Se espera que cuando h(x,y)=1 las longitudes de los lados del triángulo en equilibrio tiendan a ser iguales, lo que dará una malla de una buena calidad. Para resolver las fuerzas en equilibrio, se colectan las coordenadas x,y de todos los N puntos de la malla en la matriz p(N,2).
El vector de fuerza 
tiene componentes horizontales y verticales en cada nudo de la malla.
corresponde a las fuerzas internas en los lados de los triángulos y 
las fuerzas externas (en las fronteras). La primera columna contiene los componentes en X de las fuerzas y las segunda corresponde a los componentes de Y. depende de los bordes de los triángulos. En el algoritmo, éstas estructuras son dadas por la triangulación de Delaunay de los puntos de la malla. La triangulación de Delanauy determina la no sobreposición de los triángulos, donde cada borde es compartido por al menos dos triángulos, exceptuando los que conforman las fronteras.
El sistema 
, tiene que ser resuelto por un conjunto de posiciones p; esto es un problema difícil de resolver, en parte debido a la discontinuidad de la función fuerza y a las fuerzas de reacción externas en las fronteras. Una aproximación para resolver el sistema requiere introducir términos artificiales dependientes del tiempo. Así para 
, se considera un sistema de ecuaciones diferenciales parciales.
Si la solución en estado estacionario es encontrada, ésta debe satisfacer el sistema 
. El sistema mostrado por la ecuación 5 es resuelto por el método de Euler (forward). En la discretización del tiempo 
, la solución aproximada 
es actualizada por:
Cuando se evalúan las funciones de fuerza, las posiciones de pn son conocidas así como su verdadera topología. Todos los puntos que quedan fuera de la región durante la actualización de pn a pn+1 son movidos al interior de la fronteras lo que permite que las fuerzas actúen de forma normal a las fronteras. Los puntos pueden moverse al interior y en las frontera, pero no en el exterior de la misma. Hay muchas alternativas para calcular las funciones de fuerza F(p) en cada lado de los triángulos, algunas opciones han sido investigadas por diferentes autores [24-25]. Shimada (1995) calculó las fuerzas de atracción y repulsión F(p) como las fuerzas generadas entre dos burbujas. En una configuración tangencial ideal, la longitud (lados de los triángulos formados) es igual a los centros de dos burbujas adyacentes. Cuando dos burbujas se superponen, se repelen, cuando están lejos la una de la otra, se atraen [24]. Bossen (1996) sugiere que al calcular las funciones de fuerza con la expresión F(p) = (1-d4)*exp(-d4). Se obtienen buenos resultados, siendo d la distancia (m) entre dos puntos. Si en el proceso de actualización de la posición de los nodos, dos de estos se encuentran muy cerca el uno del otro (d<1), se repelen y si por el contario están muy distantes se atraen [25]. El algoritmo propuesto por Persson (1995) y usado para la construcción de las mallas, emplea una respuesta lineal para las fuerzas de repulsión y atracción [5]. Para una revisión detallada del algoritmo y de la forma de calculo de las fuerzas internas y externas, los lectores son referidos a Persson y Strang (2004) [6].
2.3 Discretización del dominio
Donde; q es el criterio de calidad; a, b, y c son los lados de un triángulo; rin es el radio inscrito y rout es el radio circunscrito. Un triángulo
]]> equilátero tiene un q=1; mientras que un triángulo degenerado (área cero) tiene un q=0. Como regla de truncación de error se usó q>0,5 para asegurar buenos resultados.
3. RESULTADOS
3.1 Resolución espacial
Figura 2. Distribución espacial de los ángulos de las mallas
Figure 2. Spatial distribution for angles of meses
Tabla 1. Definición de la resolución de malla para el modelo regional del Caribe
Table 1. Mesh resolution definition for the Caribbean sea model
similar. El análisis de frecuencia del parámetro q para cada elemento de la malla mostró que para las tres resoluciones espaciales la distribución de la calidad es idéntica (figura 3). Estos resultados demuestran que el algoritmo usado propuesto produce mallas con calidades idénticas, no importando la resolución espacial. Más del 60% de los elementos de las mallas tienen un q igual o cercano a 1, lo que significa que los ángulos de los vértices de los elementos son de 60° y que los lados de los triángulos son iguales (criterio de equilateridad). El 90% de los elementos de las mallas tienen un criterio q de calidad entre 0.9 y 1.0.
Figura 3. Histograma de calidad de los elementos
Figure 3. Elements quality histogram
Se encontró que para resoluciones mayores de la malla, el tiempo de procesamiento, es decir, el tiempo de corrida del modelo, es mayor (ver tabla 1). Para éste estudio, el costo computacional en la implementación del modelo depende de la cantidad de elementos en la malla procesada y no de la calidad de los mismos, pues el algoritmo usado produce mallas con calidades similares.
La distribución de los ángulos de los vértices que conforman la malla (figura 4) muestra que la mayor parte de éstos tienen 60 grados (el 40%). Esta distribución es consistente tanto para la malla fina como para la media y gruesa. Se puede notar que entre más fina es la malla (menor tamaño de elementos), hay una mayor tendencia a formar triángulos equiláteros.
Figura 4. Distribución de los angulos en las tres mallas
Figure 4. Distribution of angles in three meshes
3.3 Error de interpolación
El error relativo de interpolación con respecto a 2000 datos de batimetría, se encontró utilizando el estimador del error de la raíz media cuadrática, RMS, (8).
Al interpolar la batimetría en cada una de las mallas se encontró que el error absoluto más bajo correspondió a la malla fina con un 6.60% del promedio de los 2000 datos batimétricos usados. Los errores relativos para las mallas media y gruesa, fueron de 10.23% y 13.13% respectivamente. La malla fina produjo un error de interpolación más pequeño frente a los resultados de las otras dos mallas analizadas. El resultados grafico de la interpolación de la batimetría en la malla fina puede apreciase en la figura 5.
Figura 5. Batimetría en la malla fina
Figure 5. Bathymetry of the fine mesh
4. CONCLUSIONES
]]> La capacidad de los modelos de elementos finitos de trabajar con una malla irregular, sitúan a este método como una de las mejores opciones para el modelado de zonas costeras con geometrías irregulares. Para generar las mallas de este estudio, se utilizó un algoritmo implementado en Matlab y originalmente propuesto por Persson y Strang (2004) el cual permite controlar la calidad de los elementos triangulares que forman la malla no estructurada. Éstos elementos triangulares tienen una tendencia a ser equiláteros, lo cual es una condición deseada para obtener una buena solución en la implementación del método de los elementos finitos.Se probaron tres mallas, diferenciadas en la resolución y cantidad de elementos triangulares. La malla fina, presentó un mayor número de vértices con ángulos de 60 grados, en general, los ángulos de los vértices de los elementos se distribuyeron en un 95% entre 50 y 80 grados para la malla fina, mientras que para la malla gruesa este porcentaje fue del 83%. El algoritmo empleado tiene la propiedad de ajustar a 60 grados la mayoría de los ángulos en los vértices que forman los elementos de las mallas.
El aumento del número de elementos incrementa el tiempo de cálculo para generar una solución estable. El algoritmo implementado tiene la propiedad de producir mallas no estructuradas con calidades similares. Para medir la calidad de la malla se utilizó la relación entre los radios de los círculos inscritos y circunscritos en cada elemento, obteniendo en todos los casos, calidades superiores a 0.9. El costo computacional dependió de la cantidad de elementos de la malla y no de la calidad de los mismos.
Aun que una mayor resolución incrementa el número de elementos, lo que a su vez aumenta el costo computacional, favorece la disminución del error en la interpolación de la batimetría en los nudos de la malla formada, lo cual beneficia la consecución rápida de soluciones estables, en la implementación del método de los elementos finitos.
REFERENCIAS
[1] BERN, M. AND EPPSTEI., D. MESH generation and optimal triangulation. Available: http://www.ics.uci.edu/~eppstein/pubs/BerEpp-CEG-95.pdf [citado 23 Enero de 2008]. [ Links ]
[2] SOUZA, B.A., MATOS, E.M., FURLAN, L.T. AND NUNHEZ, J.R. A simple two-dimensional method for orthogonal and nonorthogonal grid generation, Computers and Chemical Engineering, 31, 800–807, 2007. [ Links ]
[3] MALISKA, C. R. Transferência de calor e mecânica de fluidos computacional. (2nd. ed.). LTC. Rio de Janeiro, Brasil. 2004. [ Links ]
[4] WILLIAMSON, D. Difference approximation for fluid flow on a sphere, En: Numerical methods used in atmospheric models (vol. II), GARP Publication Series no. 17, WMO-ICSU, 1979. [ Links ]
[5] LEGRAND, S., LEGAT, V. AND DELEERSNIJDER, E. Delaunay mesh generation for an unstructured-grid ocean general circulation model, Ocean Modelling, 2, 17-28, 2000. [ Links ]
[6] PERSSON, P. AND STRANG, G. A Simple Mesh Generator in MATLAB, SIAM Review, 46, 329-345, 2004. [ Links ]
[7] JOE, B. Delaunay triangular meshes in convex polygons, SIAM J. Sci. Stat. Comput., 7, 514-539, 1986. [ Links ]
[8] MULLER, J.D., ROE, P.L. AND DECONINCK, H., A frontal approach for internal node generation for Delaunay triangulations, International Journal of Numerical Methods in Fluids, 17, 241–256, 1993. [ Links ]
[9] JOE, B. Tetrahedral mesh generation in polyhedral regions based on convex polyhedron decompositions, Int. J. Numer. Methods Eng., 37, 693-713, 1994. [ Links ]
[10] AHMAD, N., BACON, D., BOYBEYI, Z., DUNN, T., HALL, M., MAYS, M., SARMA, A. AND TURNER, M. A. Solution-Adaptive Grid Generation Scheme for Atmospheric Flow Simulations. Numerical Grid Generation in Computational Field Simulation. Memorias Proceedings of the 6th Intl. Conf. held at Greenwich, England , 327-335. Julio 1998. [ Links ]
[11] BOIVIN, C. AND OLLIVIER-GOOCH, C. Guaranteed-Quality Triangular Mesh Generation For Domains with Curved Boundaries, International Journal for Numerical Methods in Engineering, 55, 1185-1213, 2002. [ Links ]
[12] FOREMAN, M. An Analysis of the Wave Equation Model for Finite Element Tidal Comparisons, J. Comput. Phys,. 52, 290-312, 1983. [ Links ]
[13] WESTERINK, J. LUETTICH, R. AND MUCCINO, J. Modeling Tides in the Western North Atlantic Using Unstructured Graded Grids, Tellus 46A, 178-199, 1994. [ Links ]
[14] LUETTICH, R. AND WESTERINK, J. Continental Shelf Scale Convergence Studies with a Barotropic Tidal Model, Coastal and Estuarine Studies, 47, 349-372, 1995. [ Links ]
[15] LYNCH, D., WERNER, F., CANTOS-FIGUEROLA, A. AND G. PARILLA. Finite Element Modeling of Reduced-Gravity Flow in the Alboran Sea: Sensitivity Studies, Memorias Seminario Sobre Oceanografía Física del Estrecho de Gibraltar, Madrid, Spain, 283-295, Octubre 1988. [ Links ]
[16] GEORGE, L., HECHT, F. AND SALTEL, E. Automatic mesh generator with specified boundary. Comput Methods Appl Mech Eng. , 92, 269–88, 1991. [ Links ]
[17] BLAIN, C.A. AND ROGERS, W.E. Coastal Tide Prediction Using the ADCIRC-2DDI Hydrodynamic Finite Element Model: Model Validation and Sensitivity Analyses in the Southern North Sea/English Channel. Naval Research Laboratory. Available: http://www.adcirc.org/publications/1998/1998_Blain03.pdf [citado febrero20 de 2008]. [ Links ]
[18] BRYAN, K. A numerical method for the study of the circulation of the world, ocean. J. Comput. Phys., 4, 347-376, 1969. [ Links ]
[19] WATSON, D.F. Computing the n-dimensional Delaunay tessellation with application to voronoi polytopes, Comput. J., 24, 167–72, 1981. [ Links ]
[20] RUPPERT, J. A. Delaunay refinement algorithm for quality 2-dimensional mesh generation, Journal of Algorithms, 18, 548–585, 1995. [ Links ]
[21] SHEWCHUK, J.R. Delaunay refinement algorithms for triangular mesh generation. Computational Geometry: Theory and Applications, 22, 21–74, 1996. [ Links ]
[22] FREY, P. AND GEORGE, P. Mesh generation: application to finite elements. Oxford, UK : Hermes Science, 2000. [ Links ]
[23] WANG, W., MING., C. AND LO, S. Generation of triangular mesh with specified size by circle packing, Advances in Engineering Software, 38, 133–142, 2007. [ Links ]
[24] SHIMADA, K. AND GOSSARD, D. C. Bubble mesh: automated triangular meshing of non-manifold geometry by sphere packing, Memorias SMA’95: Proceedings of the Third Symposium on Solid Modeling and Applications, Salt Lake City, U.S.A. , 409–419, Mayo 1995. [ Links ]
[25] BOSSEN, F AND PAUL S. HECKBERT. A pliant method for anisotropic mesh generation. Memorias Proceedings of the 5th International Meshing Roundtable, Pittsburgh, U.S.A., 63–76. Octubre 1996. [ Links ]
[26] FIELD, D. Qualitative measures for initial meshes, International Journal for Numerical Methods in Engineering, 47, 887–906, 2000. [ Links ] ]]>