Modelización y simulación numérica de campos de viento mediante elementos finitos adaptativos en 3-D

D. LUIS MAZORRA MANRIQUE DE LARA, SECRETARIO DEL DEPARTAMENTO DE INFORMÁTICA Y SISTEMAS DE LA UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA, CERTIFICA, Que el Consejo de Departamento en su sesión de fecha 11 de junio de 2004 tomó el acuerdo de dar el consentimiento para su tramitación, a la tesis doctoral titulada "Modelización y simulación numérica de campos de viento mediante elementos finitos en 3-D" presentada por el doctorando D. Eduardo Rodríguez Barrera y dirigida por los Doctores D. Gustavo Montero García y D. Rafael Montenegro Armas. Y para que así conste, y a efectos de lo previsto en el Art° 73.2 del Reglamento de Estudios de Doctorado de esta Universidad, firmo la presente en Las Palmas de Gran Canaria, a 11 de junio de 2004. ra Manrique de Lara UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA Departamento: Informática y Sistemas Programa de Doctorado: Sistemas Inteligentes y Aplicaciones Numéricas en Ingeniería Modelización y simulación numérica de campos de viento mediante elementos finitos en 3-D Tesis Doctoral presentada por D. Eduardo Rodríguez Barrera Dirigida por los Doctores D. Gustavo Montero García y D. Rafael Montenegro Armas El Director El Director El Doctorando /^./LJ: A. / Fdo:^ustavo Montero García Fdo: Rafael Montenegro Armas Fdo: Eduardo Rodríguez Barrera Las Palmas de Gran Canaria, a 26 de mayo de 2004 Universidad de Las Palmas de Gran Canaria Departamento de Informática y Sistemas 27 Tesis Doctoral >u /f AB!A LAS PALMAS DE G CANARIA N" Copia Modelización y simulación numérica de campos de viento mediante elementos finitos adaptativos en 3-D Autor Director Director Eduardo Rodríguez Barrera Gustavo Montero García Rafael Montenegro Armas Las Palmas de Gran Canaria, Mayo de 2004 A Nancy Agradecimientos Sirvan las siguientes líneas como reconocimiento a las personas que de una forma u otra han contribuido al desarrollo de este trabajo. En primer lugar quiero agradecer a mis directores de tesis Gustavo Montero García y Rafael Montenegro Armas su dedicación y apoyo durante todo el proceso, desde sus clases en los cursos de doctorado, hasta las largas horas pasadas delante de los ordenadores buscando "bichos" en los códigos. Tanto los conocimientos que ambos acumulan como sus talantes personales han sido imprescindibles para la realización de este trabajo. Quisiera hacer una mención especial a José María Escobar Sánchez y a José María González Yuste, miembros del grupo de investigación formado por las divisiones GANA y DDA del lUSIANI. El primero me ayudó a comprender algunos de los entresijos de la triangulación de Delaunay en 3D y el suavizado y desenredo de mallas, y compartió con todos los miembros del grupo el código desarrollado en su tesis doctoral. El segundo aportó el código de refinamiento adaptativo, en el que sigue trabajando para completar su tesis doctoral. A todos los miembros del Departamento de Matemáticas de la Universidad de Las Palmas de Gran Canaria, donde trabajo como informático, por darme ánimos. En este apartado debo hacer una especial mención a Angelo Santana del Pino, director del Departamento, porque siempre que acudí a él encontré soluciones. Quiero agradecer de una forma muy especial a mi esposa, Nancy, su infinita paciencia; han sido muchas las horas que le he hurtado para dedicárselas a esta tesis. Su apoyo ha sido constante, sobre todo en los momentos en que la tarea parecía ser más fuerte que yo. También quiero expresar mi gratitud a Maly, mi madre, y Juan José, mi hermano, por estar siempre ahí. VI Esta tesis ha sido desarrollada en el marco del proyecto subvencionado por el Ministerio de Ciencia y Tecnología y FEDER, REN2001-0925-C03-02/CLI, titulado Modelización numérica de transporte de contaminantes en la atmósfera. índice general Introducción 1 1.1. Estado del arte 3 1.2. Justificación 7 1.3. Objetivos 8 1.4. Metodología 8 Generación de mallas pstra dominios definidos sobre orografías irregulcires 13 2.1. Introducción 13 2.2. Refinamiento y desrefinamiento en mallas bidimensionales .... 16 2.3. Algoritmo de Delaunay en 3-D 18 2.3.1. Aspectos generales 19 2.4. Punción de espaciado vertical 23 2.5. Definición de la nube de puntos 25 2.6. Construcción de la malla 31 2.7. Optimización de la malla. Suavizado y desenredo simultáneo ... 33 2.7.1. Funciones objetivo 35 2.7.2. Funciones objetivo modificadas 37 2.7.3. Optimización de las funciones objetivo modificadas .... 41 2.7.4. Problema test 44 2.8. Aplicaciones del generador de mallas 45 2.8.1. Sur de la Isla de La Palma. Comportamiento de las cuatro estrategias de generación de puntos 45 2.8.2. Emplazamiento de una chimenea en el dominio 53 2.8.3. Detalle de la Isla de La Palma. Análisis del proceso de optimización de la malla 56 2.9. Enlace entre códigos 59 viii índice general 3. Modelización de campos de viento 63 3.1. Modelo de masa consistente 63 3.2. Construcción del campo inicial 66 3.2.1. Análisis de las medidas de las estaciones 66 3.2.2. Interpolación horizontal 67 3.2.3. Extrapolación vertical 68 3.2.3.1. Estratificación atmosférica 68 3.2.3.2. Estabilidad atmosférica 70 3.2.3.3. Perfil vertical de velocidades de viento 70 3.2.4. Corrección de la componente vertical en la trayectoria de la pluma 75 3.3. Discretización mediante elementos finitos 83 3.4. Resolución del sistema de ecuaciones 85 3.4.1. Precondicionamiento 86 3.4.1.1. Precondicionador de Jacobi 87 3.4.1.2. Precondicionador SSOR 88 3.4.1.3. Precondicionador ILL^(°) 89 3.4.1.4. Precondicionador diagonal óptimo 89 3.4.2. Algoritmo de gradiente conjugado precondicionado .... 89 3.5. Refinamiento adaptativo 91 3.6. Experimentos numéricos 99 3.6.1. Problema test 99 3.6.2. Efecto de una chimenea en el campo de velocidades .... 103 4. Estimación de pEirámetros 107 4.1. Definición del problema 107 4.2. Algoritmos genéticos 109 4.3. Experimentos numéricos 113 4.3.1. Influencia del viento en la estimación de parámetros .... 114 4.3.2. Influencia de la topografía en la estimación de parámetros 116 4.3.3. Influencia de la malla en la estimación de parámetros . . . 124 5. Simulación numérica en la Isla de Gran Canaria 133 5.1. Estimación de los parámetros 133 5.2. Validación de los resultados en otra secuencia de mallas adaptativas 143 índice general ix 6. Conclusiones y líneas futuras 151 í- <' i A ^ y'> -^'^•' ^ÍMS: .•• ' 'ñ '-'''' í - "yii ^^V: - • .^J' • y' " ,/y . - ^^J- 2 Introducción molinos de viento para la elevación de agua y la molienda de grano, y se mantienen hasta bien entrado el siglo XIX, donde su desarrollo se interrumpe con la revolución industrial y la utilización masiva de vapor, electricidad y combustibles fósiles como energía. Sin embargo, en la segunda mitad del siglo XIX aparece el popular molino multipala tipo americano, utilizado para bombeo de agua prácticamente en todo el mundo, y cuyas características habrían de sentar las bases para el diseño de los modernos generadores cólicos. El siglo XX ha visto el desarrollo de los modernos aerogeneradores y cómo la superficie del planeta se va cubriendo paulatinamente de más y más parques cólicos que se han convertido en una alternativa viable a las centrales térmicas. Por otra parte, la revolución industrial ha traído como consecuencia el vertido masivo a la atmósfera de sustancias contaminantes. Existen cada vez más evidencias de que la contaminación atmosférica tiene graves repercusiones que provocan la alteración de las condiciones medioambientales del planeta. Las consecuencias de la contaminación van desde la lluvia acida, hasta el aumento de trastornos respiratorios y alérgicos de la población, pasando por el preocupante efecto invernadero, de graves consecuencias a largo plazo. En los últimos años se ha producido un crecimiento notable de la producción de energía eléctrica de origen cólico. La Conferencia de Madrid (marzo de 1994) considera viable que las energías renovables contribuyan con un 15 % a la demanda total de energía primaria en la CE antes del año 2010. España ocupa un lugar destacable en el panorama cólico comunitario, con el quinto puesto por potencia cólica instalada, detrás de Dinamarca, Alemania, Reino Unido y Holanda. Las empresas del sector demandan herramientas cada vez más sofisticadas que les permitan hacer frente a las demandas de un mercado cada vez más competitivo y exigente. Así, la sociedad ha ido adquiriendo conciencia de los problemas medioambientales y valora cada vez más el uso de las energías renovables. Esta creciente inquietud social tiene una gran importancia desde el punto de vista político (ya no hay formación política que se sustraiga a los problemas ecológicos y no los incluya en su programa electoral) y económico (las empresas emplean más dinero en estudios de impacto ambiental y en publicidad para alardear de sus valores ecológicos, sean reales o no). Los modelos de viento son herramientas que permiten el estudio de diversos problemas relacionados con la atmósfera, tales como, el estudio de cómo afecta el Estado del arte viento a una determinada estructura, la dispersión de contaminantes, el estudio del emplazamiento de parques eólicos o la propagación de incendios. Los modelos de viento son, pues, herramientas cada vez más importantes para afrontar con solvencia una amplia gama de problemas de interés social, político y económico, y cada vez se exige más de ellos. 1.1. Estado del arte Tradicionalmente los modelos meteorológicos se dividen en modelos físicos y modelos matemáticos. Los primeros emplean túneles de viento sobre representaciones a pequeña escala del terreno en estudio. En los modelos matemáticos se emplean técnicas algebraicas y de cálculo para resolver ecuaciones meteorológicas. Los modelos matemáticos se dividen a su vez en analíticos y numéricos; son estos últimos los que ofrecen mejores perspectivas, ya que la gran complejidad de las ecuaciones que describen la atmósfera hacen muy difícil la resolución exacta en dominios irregulares por métodos analíticos. Los modelos de viento pueden clasificarse según la extensión del dominio de estudio, así hablamos de modelos de macroescala cuando el área de estudio abarca desde un continente hasta el globo terráqueo completo. La mesoescala se refiere a extensiones que van desde unos pocos kilómetros hasta alrededor de cien, y microescala a regiones locales que tienen como máximo alrededor de un kilómetro. Los modelos de pronóstico se basan en la solución de ecuaciones hidrodinámicas y termodinámicas que dependen del tiempo (llamadas también ecuaciones primitivas porque derivan directamente de los principios de conservación) modificadas para su aplicación a la atmósfera. Estos modelos también se denominan dinámicos [Lalas et al., 1988] para indicar la inclusión explícita de las ecuaciones dinámicas. Sin embargo, la solución del conjunto completo de ecuaciones sigue siendo un tarea costosa. Además, cuanto más elaborado es el modelo, más fiables han de ser los datos de entrada para usar las ventajas ofrecidas. Desafortunadamente, con frecuencia estos datos no están disponibles. Autores como Lalas et al. [1988] incluyen en los modelos dinámicos algunos códigos que introducen aproximaciones a las ecuaciones primitivas, a la vez que desprecian su dependencia del tiempo. Estos códigos, denominados modelos JH, se basan en una propuesta realizada por Jackson y Hunt [1975] y son ampliamente utilizados. Como ejemplo podemos citar el modelo empleado en el Atlas Euro- 4 Introducción peo de Viento [Troen y Petersen, 1989]. Otros autores, sin embargo, consideran estos modelos dentro de los de diagnóstico precisamente porque eliminan de sus ecuaciones la dependencia del tiempo. Los modelos de diagnóstico deben su nombre a que, como observó Pielke [1984], no se utilizan para realizar previsiones a través de la integración de las relaciones conservativas. Por esta misma razón se llaman también cinemáticos [Lalas et al., 1988]. Estos modelos generan un campo de viento que satisface algunas restricciones físicas. Si la única ecuación que se impone es la de continuidad, que impone la conservación de la masa, el modelo se define como de masa consistente. La relativa simplicidad de los modelos de diagnóstico los hacen atractivos desde el punto de vista práctico, en el sentido de que no requieren muchos datos de entrada y son fáciles de usar. Estos modelos utilizan los datos disponibles de forma eficiente para generar un campo de viento que satisface algunas restricciones físicas. Autores como Pennel [1983] han comprobado que en algunos casos los modelos de masa consistente mejorados, tales como NOABL y COMPLEX, mejoraron los resultados de modelos dinámicos más complicados y costosos. Sin embargo, hay que tener en cuenta que los modelos de diagnóstico (tanto los de masa consistente como los de tipo JH) no consideran los efectos térmicos ni los debidos a cambios de gradientes de presión. Como consecuencia, ñujos tales como brisas marinas, vientos en pendiente y otros tales como los de separación a favor del viento no pueden simularse con estos modelos, a no ser que éstos se incorporen en los datos de viento inicial a partir de observaciones realizadas en lugares apropiados a tal efecto [Kitada et al., 1983; Moussiopoulos et al., 1988]. Así, los modelos de diagnóstico están diseñados específicamente para predecir los efectos de la orografía sobre el flujo medio de viento considerado de manera estacionaria, esto es, flujos promediados en intervalos de tiempo de entre 10 minutos y 1 hora. Algunos autores [Troen, 1996] afirman que deben emplearse modelos JH en lugar de los de masa consistente, porque incorporan más restricciones físicas y no únicamente la ecuación de continuidad. Sin embargo, los modelos de masa consistente tienen la ventaja de que incorporan de manera natural los valores de viento observados en diferentes puntos del terreno, mientras que los JH generalmente, aunque no siempre, describen las modificaciones producidas por la orografía sobre el flujo no perturbado considerado. Es más, las pendientes muy pronunciadas Estado del arte afectan más a los modelos JH que a los de masa consistente. NOABL [Phillips y Traci, 1978] es un modelo meteorológico que proporciona una representación precisa del terreno gracias a una transformación de la coordenada vertical en la que la coordenada más baja es conforme a la superficie del terreno. Posteriormente, diversos autores [Lalas, 1983, 1985; Lalas et al., 1988; Tombrou y Lalas, 1990] realizaron algunas modificaciones en la inicialización de NOABL con el fin de que el modelo considerara el efecto de la rugosidad del terreno sobre el perfil de viento, de forma que el modelo dispusiera de perfiles más realistas que los del código original y tuviera en cuenta el cambio debido a la fuerza de Coriolis en la dirección del viento en la capa límite atmosférica. Los códigos resultantes fueron bautizados como NOABL* [Lalas et al., 1988] y EOLOS [Tombrou y Lalas, 1990], aunque actualmente es más conocido por AIOLOS, por la palabra griega aioXocr. Más tarde se introdujeron modificaciones que describen con más precisión los perfiles de viento en diferentes condiciones de estabilidad. Este nuevo código [Ratto et al., 1990] se llamó WINDS {Wind-field Interpolation by Non Divergent Squemes). Ambos modelos, AIOLOS y WINDS, usan datos de estaciones situadas en tierra y, opcionalmente, datos observados tanto de perfiles verticales como de viento geostrófico. También utilizan coordenadas conformes al terreno y se basan en la minimización de la diferencia de cuadrados entre un viento inicial obtenido por interpolación de los datos observados y el viento a ajustar, sujeto a la restricción de que el campo de viento ajustado ha de tener divergencia nula [Sasaki, 1970]. Además de los ya citados, existe toda una gama de modelos de diagnóstico que la comunidad científica ha venido utilizando en problemas relacionados con la meteorología y/o con la contaminación atmosférica. Dentro del grupo de los más conocidos podemos citar el ATMOS-1 [Davis et al., 1984; King y Bunker, 1984], que calcula campos de vientos basándose en una minimización del error de la conservación de la masa y apoyándose en medidas de viento. Utihza coordenadas conformes al terreno y un sistema de malla expandida verticalmente que asegura la solución cerca de la superficie. Otro modelo que se basa en la conservación de la masa es el DWM [Douglas y Kessler, 1988], capaz de generar campos de viento 3D sobre terreno complejo a partir de un número limitado de observaciones; incorpora observaciones tanto a nivel del terreno como de capas altas del aire, si estuvieran disponibles, y proporciona información de flujos de aire generados en el terreno en regiones donde no hay observaciones locales. Utiliza un sistema de coordenadas Introducción conformes al terreno. Al igual que los anteriores, MASCÓN [Dickerson, 1978] y MATHEW están basados en la ecuación de continuidad. El primero utiliza técnicas de cálculo variacional para ajustar los flujos horizontales observados y es muy similar al segundo. MATHEW por su parte [Sherman, 1978; Rodríguez et al., 1985] proporciona un campo de viento medio 3D ajustado por mínimos cuadrados, de acuerdo con la función integral E{u,v,'w,X) = Jv al(u - uo)^ + al{v - VQ)'^ + al{w - 100)^+ fdu dv dw^ \dx dy dz j dxdydz (1.1) donde «(x, y, 2), v{x, y, z) y w{x, y, z) son las componentes de viento calculadas por el modelo mediante ajuste; UQ{X, y, z), VQ{X^ y, z) y WQ{X, y, z) son las componentes del campo interpolado a partir de las observaciones; X{x, y, z) es el multiplicador de Lagrange y ai, 02 y 03 son los módulos de precisión de Gauss. La solución {u, V, w) se determina minimizando E en la ecuación (1.1). La minimización da una fórmula para {u, v, w) y una ecuación diferencial para A, que puede resolverse si se proporcionan las condiciones de contorno adecuadas. Por último, cabe mencionar MINERVE [Geai, 1987], modelo de masa consistente que reduce la divergencia mediante un procedimiento iterativo que puede tener en cuenta la estabilidad atmosférica. Entre los modelos de pronóstico cabe destacar el NMM [Pielke et al., 1983] que proporciona simulaciones meteorológicas de mesoescala en regiones de orografía compleja (zonas costeras y terrenos complejos). El modelo simula circulaciones 3D con intervalos de malla horizontal de 1 km a 10 km y considera una atmósfera incompresible, hidrostática y sin condensación. Ha sido utilizado ampliamente por los investigadores y se ha comportado razonablemente bien en simulaciones básicas de circulación de brisas marinas y otros regímenes de flujo de mesoescala generados topográficamente. El modelo ARAMS es el resultado de incorporar al modelo NMM el modelo de nubes no hidrostático desarrollado por Cotton. ARAMS {Advanced Regional Atmospheric Modeling System) es un modelo numérico predictivo que proporciona diversa información que ha ido evolucionando a lo largo de los años y que representa la fusión de tres modelos distintos, dos hidrostáticos y uno de nubes no hidrostático. Utiliza un esquema de mallas encajadas que le permite tratar los Justificación 7 fenómenos de gran escala con mallas más "groseras" mientras que utiliza mallas más "finas" para otros que así lo requieren. 1.2. Justificación En terrenos con orografía compleja, como los de las Islas Canarias, se hace imprescindible disponer de mallas de alta calidad para discretizar los dominios de estudio. Los modelos analizados anteriormente utilizan mallas regulares encajadas para discretizar el dominio. Con este tipo de mallas, en un dominio con orografía compleja el tamaño de elemento ha de ser excesivamente pequeño para captar la información digitalizada del terreno, lo que provoca que en zonas donde no se requiera tanto detalle exista también una innecesaria concentración de nodos, dando lugar por tanto a sistemas de ecuaciones de alto orden cuya resolución puede suponer un coste computacional impracticable. Frente a estas limitaciones nosotros proponemos el uso de mallas no estructuradas adaptativas, que permiten un tamaño pequeño de elemento allí donde es necesario, manteniendo tamaños de elemento mayores en otras zonas del dominio donde no se requiera gran precisión. Además proponemos el refinamiento local de la malla en función de los requerimientos de la solución numérica. Asimismo, el mallador debe concentrar los elementos en las proximidades del terreno, que es donde se necesita mayor precisión. Debemos disponer además de un método de suavizado que asegure una buena calidad de la malla y su eventual desenredo, si fuera necesario. Ya hemos mencionado que los modelos de masa consistente son los más utilizados por su mayor economía. Sin embargo, se les critica porque sus resultados dependen fuertemente de los parámetros que los definen. En los modelos anteriores, los valores de los parámetros se determinan de forma empírica. Nosotros proponemos una herramienta para la estimación de estos parámetros, dando así mayor robustez al modelo. Aunque en la actualidad existen modelos de viento que emplean elementos finitos con mallas regulares, la mayoría de los programas comerciales usan diferencias finitas. El método de elementos finitos, y especialmente cuando se utilizan técnicas adaptativas en 3D, mejoran sustancialmente los resultados que se obtienen con los modelos mencionados anteriormente. Introducción 1.3. Objetivos El objetivo que se plantea en esta tesis es desarrollar un modelo de viento tridimensional de tipo masa consistente y su correspondiente implementación. Para cubrir este objetivo se deben contemplar los siguientes aspectos: • Un generador automático de mallas de tetraedros que se adapte a una orografía irregular con una precisión prefijada, concentrando mayor densidad de puntos cerca del terreno. • Un procedimiento de optimización de las mallas que incluya la posibilidad de suavizado y desenredo simultáneos. • Un método de refinamiento local de mallas de tetraedros encajadas para mejorar la solución numérica en aquellas zonas del dominio donde exista mayor error. • Un modelo de viento de masa consistente que contemple aspectos de la física atmosférica y medidas experimentales. • Un método adaptativo de elementos finitos para la aproximación de la solución del modelo de viento. • Una técnica para la modificación del campo de velocidades de viento que incorpore el efecto de la emisión de gases a través de chimeneas. • Un método basado en algoritmos genéticos que permita la estimación automática de los parámetros principales que intervienen en el modelo de viento. Se incluye también como objetivos la aplicación del modelo a problemas con datos reales y el análisis gráfico de los resultados obtenidos, haciendo uso de herramientas de visualización de datos. 1.4. Metodología Para alcanzar el objetivo asociado a la generación de la malla tridimensional se parte de un generador de mallas bidimensional [Ferragut et al., 1994] que permite refinar y desrefinar triángulos utilizando el algoritmo 4T [Rivara, 1987]. Por 3 Metodología 9 otra parte, se dispone de un código capaz de generar mallas de tetraedros mediante el algoritmo de triangulación de Delaunay [Escobar y Montenegro, 1996]. La idea inicial para construir la malla 3D adaptada a la topografía del terreno intenta combinar estos dos códigos programados en FORTRAN. En primer lugar, se utiliza el código bidimensional para generar una malla adaptada de triángulos que aproxime la orografía del terreno con una precisión preestablecida. Esta precisión vendría determinada por el parámetro de desrefinamiento introducido en [Ferragut et al., 1994]. Una vez definida la distribución de puntos sobre el terreno, es necesario introducir una función de espaciado vertical para generar la nube de puntos en el resto del dominio. Es fundamental que esta función concentre una mayor densidad de puntos cerca del terreno, debido a que es en esta zona donde se producen mayores variaciones en el campo de velocidades de viento. Esta nube de puntos sirve de base para la construcción de una malla conforme de manera que los puntos de la nube se conviertan en los nodos de la malla. Si bien existen diferentes posibilidades para abordar este objetivo, puesto que se dispone del código de triangulación de Delaunay [Escobar, 1995] se plantea su utilización. Sin embargo esta elección conlleva problemas relativos a la conformidad de la triangulación con el terreno. Por ello, se piensa en diferentes formas de atacar este problema y se opta por la transformación de la nube de puntos a un paralelepípedo auxiliar y construir en él la triangulación. Deshaciendo la transformación anterior y manteniendo la topología de la malla construida en el paralelepípedo, se obtiene la malla resultante en el dominio real conforme con la superficie del terreno [Montenegro et al., 2002a,b; Montero et al, 2003b]. De esta forma queda resuelto el problema de la conformidad de la malla con el terreno, pero eventualmente surge un importante problema de enredo de tetraedros. Para resolverlo se desarrolla un algoritmo capaz de suavizar y desenredar la malla simultáneamente [Escobar et al., 2003; Montenegro et al., 2003]. Desde el punto de vista de la implementación, es necesario modificar el programa existente de refinamiento/derefinamiento en 2D para la introducción de la información topográfica digitalizada y las diferentes estrategias de generación de puntos sobre la vertical de cada uno de los nodos de la malla bidimensional adaptada que aproxima el terreno. Asimismo el programa de triangulación de Delaunay debe modificarse atendiendo a las características particulares del dominio. Estos dos códigos se integran en un único programa escrito en lenguaje PERL, que permite generar la malla con una mínima intervención del usuario. Por otra parte. 10 Introducción se desarrolla un nuevo código en lenguaje C para el suavizado y desenredo de la malla resultante, si fuera necesario. Todos los detalles relativos a la generación y optimización de las mallas tridimensionales se presentan en el capítulo 2 de esta tesis, así como aplicaciones en diversos problemas test y regiones definidas en el Sur de la Isla de La Palma. El origen del modelo de viento propuesto en esta tesis data de las primeras experiencias de nuestro grupo de investigación en problemas bidimensionales de ajuste de campos de viento. Dichos modelos primitivos no consideraban la orografía del terreno y simplemente construían un campo inicial mediante la interpolación de medidas de las estaciones que dependía únicamente de la distancia a los nodos. Este campo se ajusta después mediante un modelo de masa consistente en 2D que utiliza elementos finitos adaptativos [Winter et al., 1995]. Con el desarrollo del generador de mallas de tetraedros, se pasa a un modelo más sofisticado que tiene en cuenta la orografía del terreno [Montero et al., 1998; Rodríguez et al., 2002]. Así, se construye una fórmula de interpolación que considera la distancia horizontal y la diferencia de cotas entre puntos. El perfil vertical de viento tiene en cuenta la rugosidad del terreno. Por otro lado, el carácter 3D del modelo permite integrar consideraciones físicas de la atmósfera y su estratificación. Una vez construido el campo tridimensional inicial, mediante una interpolación horizontal de las medidas y la correspondiente extrapolación vertical, el problema elíptico asociado al modelo de ajuste se resuelve mediante elementos finitos, utilizando técnicas de refinamiento local basadas en la subdivisión en 8-subtetraedros [Lohner y Baum, 1992; Liu y Joe, 1996; González-Yuste et al., 2003] para mejorar la solución numérica del problema. Este método de refinamiento local acota la degeneración de los elementos, introduce una mínima propagación de la zona refinada por conformidad y es razonablemente rápido. Para su desarrollo se opta por el paradigma de la programación orientada a objetos y se implementa en el lenguaje de programación C-|—t-. Para determinar los elementos de la malla que deben ser refinados se considera el uso de un simple indicador de error basado en el gradiente de la solución numérica y el diámetro de cada elemento. Los sistemas de ecuaciones asociados, de matriz simétrica y definida positiva, se resuelven con el método de gradiente conjugado precondicionado. El código dispone de varios precondicionadores, como Jacobi, SSOR, ILL"^ y diagonal óptimo [Montero et al, 2003a, 2004a]. Asimismo, el modelo se diseña con la intención de que pueda proporcionar un Metodología 11 campo de velocidades a un modelo de dispersión de contaminantes en la atmósfera. La capacidad del mallador para discretizar chimeneas situadas sobre el terreno nos motivó para transformar el campo de velocidades de viento con el objeto de incluir el movimiento de los gases emitidos por las chimeneas debido a la velocidad de salida y a las diferencias de temperatura con el aire atmosférico [Montero et al., 2004b]. Esta modificación se llevó a cabo incorporando un modelo de pluma gaussiana al modelo de viento [Boubel et al., 1994]. Así, se comprueba que si dicha aportación al campo de velocidades se lleva a cabo en el campo interpolado, el modelo de masa consistente nos proporciona un campo ajustado que incluye el efecto de las chimeneas y sigue siendo incompresible. Un campo de estas características tiene mucho interés a la hora de resolver las ecuaciones de convección-difusiónreacción de un modelo de contaminación atmosférica [Seinfeld, 1998; Winter et al., 2004]. La implementación del modelo de viento, así como de los módulos de resolución de sistemas de ecuaciones y de precondicionadores, se realiza utilizando el lenguaje C. Es importante, dado el tamaño considerable que pueden alcanzar algunos problemas, que se haga una gestión dinámica de la memoria RAM. También conviene usar estructuras de datos adecuadas a cada aspecto del código; así, se emplea almacenamiento compacto tipo morse para las matrices y listas encadenadas para facilitar el ensamblaje del sistema de ecuaciones. En definitiva, la programación debe tener como meta un uso eficiente de la memoria sin penalizar el tiempo de ejecución. Además, la entrada de datos al modelo debe ser lo más cómoda posible para el usuario, evitando el uso de ficheros crípticos donde únicamente figuran cifras. En este sentido, se propone el uso de ficheros de texto con una sintaxis sencilla que permita al usuario introducir la información de manera natural y donde se puedan insertar comentarios descriptivos. En el capítulo 3 se recogen los detalles relacionados tanto con el modelo de viento, como con el algoritmo de refinamiento local de tetraedros, y se aplican a un problema test y a un dominio definido en la Isla de La Palma que incluye una chimenea. La mayor crítica que se hace a los modelos de masa consistente es su fuerte dependencia de los parámetros que lo gobiernan, de ahí la importancia de contar con un método que permita su estimación. Se desarrolla un código en lenguaje C que hace uso de Algoritmos Genéticos [HoUand, 1992] para la estimación automática de cuatro de los parámetros principales del modelo [Rodríguez et al., 2002; 18 Generación de mallas para dominios definidos sobre orografías irregulares cota del nodo estudiado y el valor interpolado de las cotas correspondientes a los dos nodos extremos de su lado entorno, es decir, el lado en que ese nodo fue introducido en su punto medio durante el proceso de refinamiento. Si esa diferencia es menor que el parámetro de desrefinamiento e, entonces el nodo podría ser eliminado, aunque en algunos casos deberá permanecer por razones de conformidad. Destacamos que la malla bidimensional obtenida puede ser modificada al construir la triangulación de Delaunay en el dominio tridimensional, puesto que lo único que necesitamos y conservamos es la posición de sus nodos. También nos interesa tener presente el nivel en que cada nodo es propio, para proceder a la generación de nodos en el interior del dominio. Este último aspecto se utiliza en las estrategias propuestas. 2.3. Algoritmo de Delaunay en 3-D La triangulación de Delaunay es un método ampliamente utilizado en la construcción de mallas tanto en 2-D como en 3-D, para la aplicación del método de los elementos finitos por la buena calidad de las mallas resultantes. De hecho, en el caso bidimensional, de todas las triangulaciones de un conjunto de puntos del plano, la de Delaunay es la que hace máximo el mínimo ángulo de cualquier triángulo. Se sabe que el condicionamiento de las matrices asociadas al MEF depende, entre otros aspectos, del mínimo ángulo de todos los triángulos de la triangulación, lo que justifica el uso de la triangulación de Delaunay en el MEF. Sin embargo, no existe una generalización de esta propiedad para dimensión mayor que dos, al menos en términos de alguna medida angular de los símplices. A pesar de ello, los resultados experimentales revelan que las mallas tridimensionales, construidas a partir de la triangulación de Delaunay, dan buenas medidas de calidad cuando la distribución de puntos es suficientemente regular. Para crear la triangulación de Delaunay existen algoritmos de tipo incremental, como el de Bowyer-Watson [Bowyer, 1981; Watson, 1981], que construyen la triangulación de un conjunto X por adición de puntos uno a uno. Una de las principales dificultades que plantea la utilización de la triangulación de Delaunay para generar mallas es conseguir que la triangulación respete la frontera del dominio de definición, esto es, que no haya ningún tetraedro que corte sus superficies. Cuando esto ocurre decimos que la triangulación es conforme con la frontera del dominio. Este problema se resuelve, bien imponiendo que la triangu- Algoritmo de Delaunay en 3-D 19^ lación contenga todas las aristas que definen el contorno del dominio (Constrained Delaunay Triangulation), bien colocando puntos sobre las aristas del dominio en posiciones adecuadas de manera que éstas sean la unión de aristas de la triangulación {Conforming Delaunay Triangulation). Hay algoritmos (véase por ejemplo [George et al., 1991; Weatherhill y Hassan, 1992]) basados en el intercambio de aristas y caras que consiguen la conformidad con la frontera del dominio. Una vez que se ha definido una malla conforme con la frontera del dominio, una buena estrategia para adaptar la malla atendiendo a la solución numérica, asegurando que en todo momento se mantiene dicha conformidad, es utilizar métodos de refinamiento/desrefinamiento de mallas encajadas. Esta estrategia se ha desarrollado para triangulaciones bidimensionales [Ferragut et al., 1994; Plaza et al., 1992], pero en la actualidad es una línea de investigación abierta en tres dimensiones. La tesis doctoral [González-Yuste, 2004] está siendo desarrollada en esta línea y ya se tienen resultados en el refinamiento [González-Yuste et al., 2003]. 2.3.1. Aspectos generales Sea X un conjunto de puntos distintos, no todos coplanarios, del espacio euclídeo tridimensional E^. Se define el poliedro de Voronoi, V{xi)., asociado al punto Xi ^ X como el conjunto de puntos del espacio E^ que están más cercanos a Xi (a) Diagrama de Voronoi (b) Triangulación de Delaunay Figura 2.1: Diagrama de Voronoi y triangulación de Delaunay de 16 puntos aleatorios.. í --^ ) "••;' -. \- u\'v-:- i-' •', ,-'• 20 Generación de mallas para dominios definidos sobre orografías irregulares que a cualquier otro punto de X, V{xi) = {xeE^ / d{x, Xi) < d{x, Xj), Vxj e X, j j^ i} (2.1) donde d{x, y) es la distancia euclídea entre los puntos x ey.El poliedro de Voronoi V{xi) sería un conjunto de puntos no acotado si y solo si Xi se encuentra sobre la frontera de la envolvente convexa de X. Conforme a lo anterior, definimos el diagrama de Voronoi, V{X), como el conjunto de poliedros de Voronoi asociados a los puntos de X. El diagrama de Voronoi así definido cubre todo el espacio E^, en el sentido de que cualquier punto de E^ pertenece, al menos, a un poliedro de Voronoi. Siempre es posible encontrar subconjuntos R con cuatro o más puntos de X, no todos coplanarios, tal que existe un punto UR G E^, equidistante a todos los puntos de R que satisface d{i'R, Xk) < d{vR, Xi), donde Xk E Ry Vxi E X — R. La envolvente convexa de R recibe el nombre de politopo de Delaunay (poliedro convexo acotado), D{R). A partir de la definición anterior, se deduce que todos los puntos de R están sobre una misma esfera y que no hay ningún punto de X en su interior. Normalmente, todos los politopos de Delaunay contienen solamente cuatro puntos de X, en este caso, los politopos serán tetraedros y se dirá que los puntos de X están en posición general [Fortune, 1992]. En el caso que el politopo D{R) contenga a más de cuatro puntos de X, se dice que está degenerado. La triangulación de Delaunay del conjunto de puntos X, D{X), sería el conjunto formado por todos los politopos de Delaunay. Sólo cuando este conjunto esté compuesto por tetraedros, la triangulación de Delaunay será una triangulación propia tridimensional. Sin embargo, si D{X) es una triangulación impropia, siempre será posible descomponer los politopos degenerados en tetraedros. Se conoce como algoritmos de triangulación increméntales a aquellos que se basan en la adición de puntos uno a uno sobre triangulaciones previas. Se describe a continuación un algoritmo incremental que construye una triangulación propia de Delaunay, basado en el algoritmo de Bowyer-Watson [Bowyer, 1981; Watson, 1981] que ha sido modificado por varios autores. Los principales pasos del algoritmo son: 1. Triangulación inicial: El algoritmo se simplifica si construimos un tetraedro inicial, TQ, o un prisma descompuesto en tetraedros, que contenga a todos los puntos de X. 2. Adición de Xi^i en la triangulación: Algoritmo de Delaunay en 3-D 21 Supongamos que ya ha sido construida la triangulación D{Xi), de los primeros i puntos de X. Entonces buscamos el conjunto D\ de tetraedros cuyas esferas circunscritas contienen al nuevo punto Xj+i: D{ = {T,- G D{X,) I x,+i G 5(T,.)} (2.2) Sea C]-£,j el conjunto formado por las caras de la frontera de D\, esto es, caras que sólo pertenecen a un tetraedro de D\. Podemos definir la triangulación £>(Xi+i) que incluye el punto Xj+i como: D(X,+i) = {D{X,) - D\) U D\ (2.3) donde D\ = \} {Cj , Xj+i} es el conjunto de tetraedros formados por el punto j Xj+i y la cara genérica Cj de C}jj^. Si las circunsferas de los tetraedros de D\ están evaluadas exactamente, entonces los tetraedros de D2 estarán bien construidos, en el sentido de que las caras de C)j^^ serán visibles con respecto al punto Xi-^-l, esto es lo mismo que decir que D\ es en forma de estrella con respecto a Xj+i. Como ejemplo en la figura 2.2 se muestra la triangulación del conjunto Xg = {xi,... ,Xs} (figura 2.2(a)) y la triangulación que resulta al añadir el punto Xg (figura 2.2(d)). El conjunto de caras de la frontera de Df está resaltado con líneas gruesas en la figura 2.2(c). Cuando los puntos del conjunto X no están claramente en posición general, los problemas debidos a errores de redondeo cometidos al usar representación de números en coma flotante pueden hacer que no funcione correctamente el algoritmo presentado anteriormente, pudiendo generarse tetraedros planos o cruces entre ellos. Este problema ha sido estudiado y resuelto satisfactoriamente en [Escobar, 1995]. El objeto o dominio Í2 que deseamos mallar deberá ser definido por puntos distribuidos en su interior y en su superficie. Los puntos deben estar situados de forma que la triangulación del objeto, T(f2), sea un subconjunto de la triangulación global, D{X), en la que en el conjunto X se incluyen los puntos que definen el objeto y los puntos que pertenecen al prisma inicial descompuesto en tetraedros. Si se satisface la condición anterior, entonces diremos que la triangulación T{Q) es conforme con la frontera de íl. En este caso, T(Q) resulta directamente al eliminar 22 Generación de mallas para dominios definidos sobre orografías irregulares (a) Inclusión del punto Xg en la triangulación D{X^) ^5 Xi x-rY (b; ícer XjV ) Conjunto Z)í X5 \ (líneas gri XA Xg (XQ ~so» (c) Frontera de D\ (líneas gruesas) (d) Triangulación final D{XQ) Figura 2.2: Triangulación del conjunto Xg = {xi,... ,^8} y proceso de adición de un punto xg.. de la triangulación global D{X) los tetraedros que no están incluidos en el objeto. Como se mencionó anteriormente, la triangulación de Delaunay presenta problemas en lo relativo a la conformidad con la frontera, puesto que crea una malla de la envolvente convexa de los puntos que definen el dominio. En nuestro caso concreto, la conformidad con la frontera es un problema particularmente difícil de tratar debido a la complejidad de las orografías que delimitan el dominio que se pretende mallar. En la sección 2.6, se detalla la técnica que hemos utilizado y que permite el uso de la triangulación de Delaunay sin recurrir a técnicas tan complejas como la triangulación de Delaunay restringida {Constrained Delaunay Función de espaciado vertical 23 Triangulation) o la triangulación de Delaunay conforme (Conforming Delaunay Triangulation), y al mismo tiempo evita los problemas de conformidad con la frontera del dominio. 2.4. Función de espaciado vertical Como ya se ha indicado, interesa generar una nube de puntos con mayor densidad en la zona cercana al terreno. Para ello, cada nodo va estar situado atendiendo a una función del tipo, Zi = ai'' + b (2.4) tal que a medida que aumenta el exponente a > 1, proporciona una mayor concentración de puntos cerca de la superficie del terreno; Zi es la cota correspondiente al z-ésimo punto insertado, del tal manera que para z = O se obtiene la cota del terreno y para Í = n, la cota del último punto introducido que debe coincidir con la altura h del plano superior que delimita el dominio a discretizar. En estas condiciones el número de puntos definidos en la vertical sería n -I-1 y la función de espaciado vertical (ver figura 2.3) se puede expresar como Zi = -i'^ + zo ; i = 0,1,2,...,n (2.5) 71°' En ocasiones conviene expresar la altitud de un punto en función de la del punto anterior, evitando así tener que conservar en memoria el valor de ZQ, Zi = ^i-i + ^.^_~(/l"\)a r -{i1)1 ; i = l,2,...,n (2.6) A partir de las ecuaciones (2.5) o (2.6), los puntos quedan perfectamente definidos una vez fijados los valores de o; y n. No obstante, también puede ser interesante fijar la distancia del primer punto insertado (z = 1) a la superficie del terreno con el fin de mantener unos parámetros mínimos de calidad en la malla tridimensional que se pretende generar. Esto reduciría el número de grados de libertad a uno, bien sea a o bien n. Consideremos fijado y conocido el valor de esa distancia d tal que d = zi — ZQ; véase figura 2.3. Sustituyendo en la ecuación (2.5), d = zi-zo = —^^ (2.7) n" 24 Generación de mallas para dominios definidos sobre orografías irregulares Figura 2.3: Distribución de n + 1 puntos sobre el eje de ordenadas mediante la función de espaciado vertical. Si fijamos a y dejamos libre el valor de n, de (2.7) se obtiene, n d l/a (2.8) No obstante, en la práctica se aproximará el valor de n al número natural más cercano. En cambio, si fijamos el valor de n y dejamos libre a, resulta, h-zn a — log n (2.9) En ambos casos, dado uno de los dos parámetros, se calcula el otro mediante las expresiones (2.8) o (2.9), respectivamente. De esta forma, la distribución de puntos en la vertical respeta la distancia d entre Zi y ZQ. Si además fijamos la distancia entre los dos últimos puntos introducidos, esto es D = Zn—Zn-i (ver figura 2.3), entonces los parámetros ayn quedan perfectamente determinados. Supongamos que a es definido por la ecuación (2.9). Para i = n — 1, la ecuación (2.5) resulta h — ZQ ^\a . Zn-\ = —r— \Ji-l) + zo (2.10) Definición de la nube de puntos 25 y por tanto, usando la ecuación (2.9), log (n - 1) _ log ^^^^ log n log h2¡í (2.11) A partir de las características con que se ha definido la malla buscada a priori se puede afirmar que h — Zo>D>d>0. Por tanto, el valor de n estará acotado, 2 < n < ^^^, y el de o; no podrá ser inferior a 1. Asimismo, para llegar a introducir al menos un punto intermedio entre la superficie del terreno y la frontera superior del dominio, se tiene que verificar que d + D < h — ZQ. 1 h—zo~D Llamando k = °^ h-z^ , se puede comprobar fácilmente que O < A; < 1. De esta forma, la ecuación (2.11) se transforma en n = l + n'' (2.12) Si denotamos g{x) = 1 + x'^, se puede comprobar que g{x) es contractiva en [2, ^^^] con constante de Lipschitz C — jir^ y además está acotada, 2<,(.)<l+(^)'<^ (2.13) En virtud del teorema del punto fijo, podemos asegurar que la ecuación (2.12) tiene solución única, y puede obtenerse numéricamente, por ejemplo, mediante el método del punto fijo, puesto que éste converge para cualquier aproximación inicial escogida en el intervalo [2, ^^^] • No obstante, en general, la solución no tomará valores enteros. Consecuentemente, si aproximamos su valor al número natural más próximo, la condición impuesta con la distancia D no se cumplirá exactamente, sino de forma aproximada. 2.5. Definición de la nube de puntos Cualquiera que sea la estrategia a seguir, la generación de puntos se realizará en tres etapas. En la primera se define una malla regular bidimensional con la densidad de puntos que se desea obtener sobre la frontera superior del dominio. En segundo lugar, y sobre esta malla ri, se lleva a cabo un proceso de refinamiento global y desrefinamiento en función de la topografía del terreno para obtener la malla r^, que define la distribución de puntos en la superficie del terreno. Para 26 Generación de mallas para dominios definidos sobre orografías irregulares ilustrar estas dos primeras etapas, en la figura 2.4 se muestra la distribución de puntos sobre ambas superficies para un problema test. Una vez definida la distribución de puntos en la superficie del terreno y en la frontera superior del dominio, se procede a la generación de la nube de puntos distribuida entre estas dos capas. Para ello, sobre la vertical de cada nodo P de la malla del terreno r^/ situaremos puntos atendiendo a la función de espaciado vertical y al nivel j en el que P es propio, siendo 1 < j < m'. La función de espaciado vertical quedará determinada mediante la estrategia utilizada para definir los siguientes parámetros: la cota topográfica ZQ de P; la altitud h de la frontera superior del dominio; el máximo número posible de puntos n+1 en la vertical de P, incluyendo el propio P y el de la frontera superior del dominio, caso de existir; el grado de la función de espaciado a; la distancia entre los dos primeros puntos generados d = zi — ZQ; Y Isi distancia entre los dos últimos puntos generados D = Zn — Zn-iDe esta forma, la cota del z-ésimo punto generado sobre la vertical de P viene dada por, -2i = ^^-^í" + ^o ; i = l,2,...,n-l (2.14) ib Independientemente de la función de espaciado vertical definida, utilizaremos el nivel j en el que P es propio para determinar el número definitivo de puntos que se generan sobre la vertical de P excluyendo el terreno y la frontera superior. Distinguiremos: 1. Si j = 1, es decir, si el nodo P es propio de la malla base ri, se generan a partir de la ecuación (2.14) para i = 1, 2,..., n — 1. 2. Si 2 < j < m' — 1, generamos nodos para z = 1,2,..., min{m' — j,n — 1). 3. Si j = m', esto es, el nodo P es propio del nivel más fino r^,, entonces no se genera ningún nodo nuevo. Este proceso tiene su justificación, ya que la malla r^, corresponde al nivel más fino de la secuencia de mallas encajadas T' = {ri < TJ < ... < r^/}, obtenida mediante el algoritmo de refinamiento y desrefinamiento, y por tanto el número de puntos introducidos decrece suavemente con la altura y además resultan distribuidos eficientemente con el fin de construir la malla tridimensional en el dominio. Definición de la nube de puntos 27 Figura 2.4: Representación tridimensional de una distribución de puntos sobre el terreno (ABCD) y la frontera superior del dominio (A'B'C'D'). Estrategia 1: número de capas y grado de espaciado vertical fijos En este caso, se considera impuesto el mismo valor de a y n para todo punto P de T^,, con el fin de generar los puntos a partir de la ecuación (2.14). Obsérvese que en la práctica, el parámetro n permite fijar el número de capas (?^+l) que se desea generar en el dominio, sobre las cuales se van a distribuir los puntos. Por otro lado, el valor de a determina el grado de concentración de capas hacia el terreno. En concreto, para a = 1 la distancia entre dos capas consecutivas es constante sobre la vertical de cada punto de la superficie del terreno. En cambio, si escogemos valores de a superiores, la concentración de capas es mayor cerca del terreno. Con esta elección, los valores de o: y n son introducidos como datos y, por tanto, la nube de puntos queda definida completamente. Consecuentemente, se pierde la libertad de fijar las distancias d y D. En general, esto conduce a que aquellos elementos con algún vértice sobre la superficie del terreno o sobre la frontera superior del dominio puedan tener una baja calidad. Con este procedimiento, la regularidad de las funciones que definen las capas va aumentando de una capa a otra superior, siendo la más irregular la correspondiente a la superficie del terreno y, la más regular, la capa horizontal correspondiente a la frontera superior del dominio. En base a esta regularidad, esta estrategia se ha diseñado de tal forma que, además de 34 Generación de mallas para dominios definidos sobre orografías irregulares cambia el signo de su determinante jacobiano. Para resolver este problema se puede proceder tal y como se propone en [Preitag y Knupp, 1999, 2002; Preitag et al., 2000]. En una primera etapa se desenredan los elementos invertidos usando un algoritmo que maximice sus determinantes jacobianos negativos [Freitag y Knupp, 2002]; en una segunda etapa se suaviza la malla resultante usando otra función objetivo basada en alguna medida de calidad de los tetraedros de N (v) [FVeitag et al., 2000]. Más adelante se presentan dos de estas funciones objetivo. Tras el proceso de desenredo la malla tiene una calidad muy pobre porque la técnica no busca la creación de elementos de buena calidad. Como queda dicho en [Freitag y Knupp, 1999], para optimizar la función objetivo asociada al proceso de suavizado no se puede aplicar un algoritmo basado en el gradiente porque ésta no está definida en todo M.^. Es necesario entonces utilizar otras técnicas para evitar este problema. Algunos autores [Tinoco-Ruiz y Barrera-Sánchez, 1998, 1999] han propuesto el desplazamiento de las barreras (singularidades) mediante la introducción de parámetros globales calculados en función del área, con el propósito de obtener mallas convexas estructuradas en dominios bidimensionales. Nosotros proponemos una alternativa a estas técnicas, de forma que el desenredo y el suavizado de la malla se realicen de manera simultánea. De esta forma se obtienen mallas con mejor calidad tras el desenredo y además se necesitan menos iteraciones de suavizado para alcanzar una cierta calidad. Para ello se usa una modificación de la función objetivo de forma que esté definida sobre K^, pero preservando las posiciones en que se encuentran los mínimos, de forma que cuando existe una región factible (subconjunto de E^ donde puede situarse v siendo N (v) una submalla válida) los mínimos de la función original y la modificada están muy cercanos. Cuando ésta región no existe, el mínimo de la función objetivo modificada es tal que tiende a desenredar A'' (v). Esto último ocurre, por ejemplo, cuando la frontera fija de N (v) está enredada. De esta forma se pueden usar métodos de optimización estándar para localizar el mínimo de la función objetivo modificada (véanse por ejemplo [Dennis y Schnabel, 1996], [Gilí et al., 1981] y [Bazaraa et al., 1993]). En nuestros trabajos se han aplicado las modificaciones propuestas para dos funciones objetivo distintas derivadas de las medidas de calidad algebraicas estudiadas en [Knupp, 2001], aunque sería posible aplicarlas también a otras funciones objetivo que presentaran singularidades, tales como las estudiadas en [Knupp, 2000]. Optimización de la malla. Suavizado y desenredo simultáneo 35 2.7.1. Funciones objetivo Las funciones objetivo se pueden construir partiendo de alguna medida de calidad de los tetraedros [Dompierre et al., 1998]. Sin embargo, aquellas que se obtienen mediante operaciones algebraicas son especialmente adecuadas para nuestro propósito ya que el coste computacional requerido para su evaluación puede ser muy bajo. Sea T un tetraedro en el espacio físico cuyos vértices están dados por x^ = {xk^VkiZkY ^ ^^) ^ = 0,1,2,3 y sea TR el tetraedro de referencia con vértices Uo = (0,0,0)^, ui = (1,0,0)^, U2 = (0,1,0)^ y U3 = (0,0,1)^. Si se toma XQ como vector de traslación, la aplicación afín que transforma TR en T es X =AU -\- XQ, donde A es la matriz jacobiana de la aplicación afín referida al nodo XQ, y cuya expresión es '^ Xi - xo X2Xa XzXQ^ A= Vi-yo Vi-yo 2/3 - yo (2-15) \ ^1 — ^O Z2 — ZQ Z3 — ZQ J Sea ahora T/ un tetraedro equilátero con todas sus aristas de longitud uno y vértices en VQ = (0,0,0f, Vi = (l,0,Of, V2 = (1/2,^3/2,0)^ y V3 = (1/2, -s/S/G, v^/x/S) . Sea v —Ww. la aplicación afín que convierte TR en T/ y siendo su matriz jacobiana W = ( 1 1/2 1/2 \ O \/3/2 V3/6 V O O \/2/V3 / (2.16) La aplicación afín que transforma T/ en T viene dada por x =ÁW~^\^ -V XQ, y su matriz jacobiana es 5 = AW~^. La matriz 5' es independiente del nodo que se toma como referencia y se dice que es invariante respecto a los nodos [Knupp, 2001]. Para construir medidas de calidad algebraicas de T puede usarse alguna norma matricial o la traza de S. Por ejemplo, la norma de Frobenius de 6*, definida por \S\ = -y^tr (6'^S'), es especialmente adecuada porque es fácilmente 2 computable. Así, se demuestra en [Knupp, 2001] que tanto q»^ (5) = y|¡lcomo q^iS) = is\\s-i\ son medidas de calidad algebraicas de T, donde a = det (S). El valor máximo de estas medidas es la unidad y corresponde al tetraedro equilátero, mientras que cualquier tetraedro plano o degenerado tiene medida nula. A partir de estas medidas de calidad se pueden obtener funciones objetivo. 36 Generación de mallas para dominios definidos sobre orografías irregulares Así, sea x = {x, y, z) la posición del nodo libre v, y Sm la matriz jacobiana correspondiente al m-ésimo tetraedro de A^" {v). Se define la función objetivo de x, asociada al m-ésimo tetraedro como \Sr, 3(7^ (2.17) Pm| I'-', -1| m |5„ _ \'~>m\ l^m I _ I'-'mi |^m| /r, -, Q\ l^m — 7, — o \^-^°) O ÓUm donde E^ = (^mS^ es la matriz de los adjuntos de SmEntonces, las funciones objetivo correspondientes de N [v) pueden construirse usando la p-norma de (^I,?72,---,?7M) O (Ati,/í2,---,/íM) como l^.l,(x) = M T.'Ci^) m=l (2.19) l^«L(x)-^ M E<w m=l (2.20) donde M es el número de tetraedros de N (v). La función objetivo ¡K^^l^ es deducida y utilizada en [Bank y Smith, 1997] para suavizar y adaptar mallas 2-D. La misma función se utiliza en [Djidjev, 2000] para el suavizado de mallas 2-D y 3-D como resultado de un método forcé-directed. Las funciones li^T^I están propuestas y ampliamente analizadas en [Preitag y Knupp, 1999; Preitag et al., 2000] y finalmente ambos tipos de funciones, entre otras, están estudiados y comparados en [Knupp, 2000]. Todas las funciones citadas únicamente pueden usarse para suavizar mallas válidas, es decir, aquellas en las que se cumple que 0"^ > O, Vm = 1,..., M. Estas funciones objetivo son suaves en aquellos puntos donde N (v) es una submalla válida, pero se vuelven discontinuas cuando el volumen de cualquier tetraedro de N (v) se hace cero. Esto se debe a que \Kn\ y ¡K^L tienden a infinito cuando cr^ tiende a cero, ya que sus numeradores están acotados. Se puede demostrar que IS'ml y l^ml alcanzan sus mínimos, con valores estrictamente positivos, cuando v está situado en el centro geométrico de la cara fija del 7?i-ésimo tetraedro. Las posiciones de v tales que N (v) sea válida, es decir, la región factible, se Optimización de la malla. Suavizado y desenredo simultáneo 37 corresponden con el interior del poliedro P definido como M P= f]H„ (2.21) m=l donde H^ son los semiespacios definidos por cr^ (x) ^ O (la región sombreada de la figura 2.7a). Este conjunto puede estar vacío, por ejemplo, cuando la frontera fija de N (v) está enredada (figura 2.7c). En esa situación, las funciones \Kjj\^ y ¡K^lp dejan de útiles como funciones objetivo. Por otra parte, si existe región factible, esto es, int P 7^ 0, las funciones objetivo tienden a infinito a medida que v tiende a la frontera de P. Debido a esas singularidades, se forma una barrera que impide que se alcance el mínimo usando algoritmos basados en el gradiente cuando se parte de un nodo libre que se encuentra fuera de la región factible (véase figura 2.7b). En resumen, no podemos optimizar una malla enredada A'' {v) trabajando con estas funciones objetivo, ni con algoritmos de tipo gradiente, aunque exista región factible. (a) (b) (c) Figura 2.7: (a) Malla válida con el nodo libre en el interior de la región factible, (b) Malla enredada con el nodo libre fuera de la región factible, (c) Malla enredada con la frontera fija enredada. 2.7.2. Funciones objetivo modificadas A continuación proponemos una modificación de las funciones objetivo (2.19) y (2.20), de manera que se elimine la barrera asociada a las singularidades y se consiga que la nueva función sea suave en M^. Un requisito esencial es que los 38 Generación de mallas para dominios definidos sobre orografías irregulares h{a) • / / / • • / • / / / / / / s a Figura 2.8: Representación de la función h{a). mínimos de la funciones objetivo originales y las modificadas sean casi idénticos cuando int P ^$. Nuestra modificación consiste en sustituir a en (2.19) y (2.20) por la función creciente y positiva 1 h(a) =-{a + Va^+W) donde b = h{0). En la figura 2.8 está representada la función h{a). Así, las nuevas funciones objetivo propuestas aquí están dadas por (2.22) 1^:1 (x) = M E (^^)' (^) m=l (2.23) donde \K\A^) = M E (<y w m=l \Sr, y 19 I IV I * |^7n| \^-'m\ '^^~ 3h{am) son las funciones objetivo modificadas para el tetraedro m-ésimo. (2.24) (2.25) (2.26) El comportamiento de h{a) en función del parámetro 5 es tal que, lím h{a) = a, Optimización de la malla. Suavizado y desenredo simultáneo 39 Va > O y lím h{a) = O, Ver < 0. Así, si int P ^ ^, entonces Vx G int P tenemos Gm (x) > O, para m = l,2,...,My, a medida que vayamos eligiendo valores de 5 más pequeños, h{(7m) se va pareciendo más a a^, de manera que las funciones originales y sus correspondientes versiones modificadas son muy próximas en la región factible. Así, en dicha región, \K^\ y \K'^\p convergen puntualmente a \KJ y \KA cuando ¿ ^ 0. Además, considerando que Va > O, lím h'{(j) = 1 y lím h'^'^^a) = O, para n > 2, es fácil comprobar que las derivadas de las funciones objetivo verifican la misma propiedad de convergencia. Como resultado de estas consideraciones, podemos concluir que las posiciones de v que minimizan las funciones objetivo originales y las modificadas son casi idénticas cuando el valor de 6 es pequeño. En realidad, S se selecciona en función del punto v bajo consideración, haciéndolo tan pequeño como sea posible pero de manera que la evaluación del mínimo de la función objetivo modificada no presente ningún problema computacional. En particular, sea 7 el épsilon de la máquina (O < 7 << 1), entonces para evitar divisiones por cero al calcular las funciones objetivo modificadas (2.25) y (2.26) hay que imponer que h{a) > 7 en todos los tetraedros de N{v). Como h{a) es una función creciente el peor caso corresponde a a = a^¿„, donde amm = min {(Tm)- Se puede demostrar que para los siguientes valores de 5, m=l,...,M 5>5..„ = í ^^^^""-"^ sia..„<7 (2.27) [ O SI amin > 7 se verifica la condición h{a) > 7. En la práctica se ha multiplicado 7 por un factor de seguridad. Hay que destacar que cuando el nodo libre está claramente dentro de la región factible {amin ^ 7)) nuestras funciones objetivo modificadas coinciden con las originales. Supongamos ahora que, int P = 0, entonces las funciones objetivo originales, \K.n\ y \KJ, no son adecuadas para nuestro propósito ya que no están correctamente definidas. Sin embargo, las funciones objetivo modificadas están correctamente definidas y tienden a resolver el enredo. Podemos razonar esto desde un punto de vista cualitativo considerando que los términos dominantes en l/í*! o I 'I \p \K*\p son aquellos que están asociados a tetraedros con valores de a más negativos y, por ello, la minimización de estos términos implica el incremento de estos valores. Hay que destacar que h{a) es una función creciente y que tanto \K*\ 40 Generación de mallas para dominios definidos sobre orografías irregulares como \K^\p tienden a oo cuando el volumen de cualquier tetraedro de A'' (v) tiende a —oo, ya que lím h (a) = 0. cr—>—oo En resumen, mediante las funciones objetivo modificadas, podemos desenredar la malla a la vez que mejoramos su calidad. Obviamente, la modificación propuesta aquí puede aplicarse fácilmente a otras funciones objetivo del mismo tipo. Para comprender mejor el comportamiento de las funciones objetivo y su modificación, se propone el siguiente problema ejemplo. Considérese una malla en 2-D formada por tres triángulos, vBC, vCA y vAB, donde se han fijado A{0, —1), -B(-\/3, 0), C(0,1), y v{x, y) es el nodo libre (véase la figura 2.9). En este caso, la región factible es el interior del triángulo equilátero ABC. En la figura 2.10(a) se muestra \Krj\^ (línea continua) y \K*\ (línea discontinua) como una función de x para un valor fijo y = O (la coordenada y de la solución óptima) y «5 = 0.1, 0.2, 0.3. La figura 2.11 muestra las gráficas equivalentes para \KK\2 (línea continua) y \K*\2 (líneas discontinuas). Puede observarse que las funciones originales presentan múltiples mínimos locales y discontinuidades, al contrario de lo que les ocurre a las funciones modificadas. Además, las funciones objetivo originales alcanzan sus mínimos absolutos fuera de la región factible. Las asíntotas verticales en las funciones objetivo originales corresponden a posiciones del nodo libre para las que cr = O para alguno de los triángulos de la malla local. Como cabría esperar en este ejemplo, la solución óptima para ambas funciones modificadas es Í;(\/3/3, 0). La funciones modificadas y las originales son casi idénticas en la proximidad de este punto, véase figura 2.10(a) y figura 2.11 (a). Consideremos ahora la malla enredada obtenida al cambiar la posición del punto5(x/3,0)por5'(-\/3,0); véase la figura 2.9(b). Aquí, la malla está constituida por los triángulos vB'C, vCA y vAB', donde vB'C y vAB' están invertidos. En esta nueva situación no existe región factible. Las gráficas de las funciones iKj^l^ y |ií*| están representadas en la figura 2.10(b). En la figura 2.11(b) están las gráficas correspondientes a \K^\2 y |-fírKl2Si bien la malla no puede ser desenredada, se obtiene Í;(—\/3/3,0) como posición óptima del nodo libre utilizando las funciones objetivo modificadas. Para esta posición los tres triángulos están "igualmente invertidos" (tienen iguales valores de a). En este ejemplo se habrían obtenido los mismos resultados si se hubiera maximizado el mínimo valor de a en la malla, tal como se propone en [Freitag y Knupp, 1999] y [Freitag y Plassmann, 2000]. Optimización de la malla. Suavizado y desenredo simultáneo 41 (a) (b) Figura 2.9: Ejemplo test: (a) malla válida y (b) malla enredada. 2.7.3. Optimización de las funciones objetivo modificadas Los algoritmos convencionales de optimización, como el descenso mínimo o el algoritmo de Newton, necesitan evaluar el gradiente y, en algunos casos, el hessiano de la función objetivo. Por eso en esta sección se van a obtener las expresiones de las derivadas de rj* y K* con respecto a parámetros arbitrarios a y /? que representan cualquiera de las coordenadas x, y y z del nodo libre. Considérese el producto interno de dos matrices, Ry S, como {R, S) = triR^S) (2.28) de forma que, la norma de Frobenius de S es l^l = •\/{S, S). Si llamamos da al operador derivada parcial respecto de a, tal que d^R — [darij] para una matriz n X n R= [rij], 1 <i,j <n, se puede demostrar que la derivada de r/* es daV* = '¿V* {dgS, S) daO Z^/^^TW (2.29) Para K* se tiene (Xy/C /í \daS,S) (<9aI],S) daO I o \sx 7^2+4^2 (2.30) 42 Generación de mallas para dominios definidos sobre orografías irregulares y 250 \ \200 \\ \ «O \ \ \ 10b 7/ V/ -1 (a) 1 \ \ \ \ w V 500 400 \ 30pf 200 / / A/ / j^ -3 -2 -1 5=0 6=0.1 5=0.2 5=0.3 — 5=0 --- 5=0.1 5=0.2 - 5=0.3 (b) Figura 2.10: (a) Corte transversal de IK^I^ (línea continua, 5 = 0) y \K*\^ (líneas discontinuas, S = 0.1,0.2,0.3j para el problema ejemplo representado en la figura 2.9(a); (b) las mismas funciones objetivo para la malla enredada de la figura 2.9(b). Para obtener las segundas derivadas hay que considerar que 5, E y cr son funciones lineales de x, y, z, y así dapS, da/s^ y dapo' son cero. Con estas consideraciones tenemos dc^pT]* = daV*d/3V* rj* + 2r}* {daS, dpS) 2 {d^S, S) {d^S_^ ^ ad^udpa_ \S\ \s\' 3 (a2 + 4¿2) (2.31) Optimización de la malla. Suavizado y desenredo simultáneo 43 \ 1 VTi M ; / > ' / í ! ! II I %J I 5=0 5=0.1 5=0.2 5=0.3 (a) \ \ \ \ 4000 3000 2000 f \ 1000 I / / / / y -3 -1 6=0 5=0.1 5=0.2 5=0.3 (b) Figura 2.11: (a) Corte transversal de \Kfi\2 (linea continua, S — 0) y \K*\2 (líneas discontinuas, 5 = 0.1,0.2,0.3^ para el problema test representado en la figura 2.9(a); (b) las mismas funciones objetivo para la malla enredada de la figura 2.9(b). y para K* se tiene OafiH = 1K* I ,0 \SY I r, (2.32) \S\ \A ((72 + 452)§_ Las ecuaciones (2.31) y (2.32) pueden simplificarse dado que {daS,daS) = |, 50 Generación de mallas para dominios definidos sobre orografías irregulares (a) Malla generada (b) Malla optimizada Figura 2.16: (a) malla generada mediante la estrategia 3 y (b) malla resultante después de cinco pasos del proceso de optimización. Aplicaciones del generador de mallas 51 (a) Malla generada (b) Malla optimizada Figura 2.17: (a) malla generada mediante la estrategia 4 y (b) malla resultante después de cinco pasos del proceso de optimización. 52 Generación de mallas para dominios definidos sobre orografías irregulares 30000 40000 50000 60000 (a) Estrategia 1 (b) Estrategia 2 10000 20000 30000 40000 50000 60000 (c) Estrategia 3 (d) Estrategia 4 Figura 2.18: Curvas de calidad de las mallas generadas y optimizadas para las distintas estrategias. estrategia. Además, con esta última estrategia se obtiene una distribución de puntos cuasi-óptima en el dominio real atendiendo a los tamaños de los elementos que se crean sobre la superficie del terreno y la parte superior del dominio. Para evitar la aparición de tetraedros invertidos se ha aplicado de forma eficiente la técnica propuesta en la sección 2.7. De hecho, la peor medida de cahdad de los tetraedros de las mallas optimizadas para las distintas estrategias está entorno a 0.2. En general, el número de parámetros necesarios para definir la malla resultante es muy reducido, así como el coste computacional. Aplicaciones del generador de mallas 53 2.8.2. Emplazamiento de una chimenea en el dominio Para la modelización del transporte de contaminantes en la atmósfera en el entorno de fuentes contaminantes, tales como centrales térmicas, debemos añadir la geometría de la chimenea a los datos topográficos y luego aplicar nuestro generador de mallas tridimensionales. Como la malla debe ser capaz de detectar los detalles de la chimenea, si decidimos un tamaño de elemento del orden de unos pocos metros en la chimenea, partiendo de una malla 2-D uniforme ri del área rectangular de la base del dominio, con un tamaño de elemento típico del orden de kilómetros, tendríamos que realizar un número considerable de pasos de refinamiento global usando el algoritmo 4-T [Rivara, 1987], lo que podría provocar problemas de falta de memoria en computadores convencionales. Para evitarlo, realizamos unos pocos pasos de refinamiento global sobre TJ y, posteriormente, un número suficiente de refinamientos locales de los elementos situados dentro de la zona que define la chimenea. Seguidamente, aplicamos el algoritmo de desrefinamiento desarrollado en [Ferragut et al., 1994] y [Plaza et al., 1996] con un parámetro de desrefinamiento e adecuado, tal que los nodos situados dentro de la chimenea no deben ser eliminados. La chimenea queda definida con el radio de la base, Tg, y el radio del extremo superior, por donde salen los gases, rj, (rj < Tg) (ver figura 2.19). A los nodos que están en el interior del círculo definido por TJ se les asigna z — Zc-\- hchim, siendo Zc la cota del terreno correspondiente al centro del círculo y hchim la altura de la chimenea. A los nodos interiores a la corona circular definida por los radios TJ y Te se les asigna un valor de z dado por la función (2.33), donde de es la distancia del punto considerado al centro de la chimenea. z = Zc + hchim e^í^c (2.33) Figura 2.19: Radios de definición de la chimenea. 54 Generación de mallas para dominios definidos sobre orografías irregulares Figura 2.20: Zona de la Isla de La Palma que incluye la discretización de una chimenea. Como ejemplo de lo expuesto, se considera una central térmica test situada en una zona rectangular de 22.8 x 15,6 km de la Isla de La Palma, donde las cotas varían desde O a 2279 m, con una chimenea de 200 m de altura sobre el terreno y de diámetro 20 m en la superficie de salida de los gases y 40 m en su base. La topografía viene dada por una digitalización donde los datos de altura están equiespaciados sobre una malla con un tamaño de paso de 200 m en las direcciones X e y. Para que la malla capte los detalles de la chimenea, si decidimos un tamaño de elemento aproximadamente igual a 2 x 2 m en la chimenea, partiendo de una malla 2-D uniforme TI del área rectangular, con un tamaño de elemento de unos 2 X 2 km, tendríamos que realizar diez pasos de refinamiento global usando el algoritmo 4-T de Rivara [1987]. Sin embargo, sólo es necesario realizar cinco pasos de refinamiento global sobre TI y, posteriormente, cinco refinamientos locales de los elementos situados dentro de la zona que define la chimenea. Seguidamente, aplicamos el algoritmo de desrefinamiento con un parámetro de desrefinamiento e = 40 m, manteniendo todos los nodos situados dentro de la chimenea. La distribución de nodos de TI se toma para la frontera superior del dominio. La figura 2.20 Aplicaciones del generador de mallas 55 Figura 2.21: Detalle de la malla tridimensional de la figura 2.20 con una chimenea cerca de la esquina inferior derecha. 56 Generación de mallas para dominios definidos sobre orografías irregulares representa el dominio utilizado en este experimento. Se observa la discretización adaptada a la topografía y la distribución de puntos en la vertical siguiendo un función de espaciado que permite mayor concentración a medida que nos aproximamos al terreno. En la figura 2.21 se muestra un detalle de la malla donde se puede distinguir la ubicación de la chimenea considerada en este estudio. 2.8.3. Detalle de la Isla de La Palma. Análisis del proceso de optimización de la malla Como aplicación práctica del mallador y del procedimiento de optimización se considera el mismo dominio del apartado anterior, siendo en este caso su altura h = 6 km. Se parte de una malla uniforme 2-D ri del área de estudio con un tamaño de elemento de aproximadamente 2x2 km para realizar a continuación seis etapas de refinamiento global usando el algoritmo 4-T de Rivara [1987]. Una vez interpolados los datos topográficos sobre esta malla refinada, se hace uso del algoritmo de desrefinamiento desarrollado en [Ferragut et al., 1994] y [Plaza et al., 1996] con un parámetro de desrefinamiento e = 25 m. De esta forma se asegura que la malla se adapta al terreno con un error menor que ese valor. La distribución de nodos de Ti es la considerada en la frontera superior del dominio. La malla generada tiene Ntet = 81068 tetraedros y Nnod = 16504 nodos, con una valencia máxima de 36 (véase la figura 2.22(a)). La malla inicial tiene 574 tetraedros invertidos con una calidad media de ^^ = 0.626; véase la figura 2.23. La distribución de los nodos se modifica significativamente tras diez iteraciones del proceso de optimización usando l-íí'^L • Tras la primera iteración deja de haber tetraedros invertidos y la calidad media se incrementa a ^^ = 0.706. Esta medida tiende a estabilizarse rápidamente; en la quinta iteración se tiene q^ = 0.732 y en la décima ^^ = 0.734. Al finalizar el proceso de optimización la peor de las cahdades de los elementos de la malla optimizada es de g™" = 0.112. Hay que destacar que el coste computacional es bastante bajo; la complejidad del algoritmo de refinamiento/desrefinamiento 2-D es lineal [Plaza et al., 1996], y además, los resultados experimentales obtenidos revelan que nuestro algoritmo de triangulación de Delaunay en tres dimensiones tiene una complejidad lineal en función del número de puntos [Escobar y Montenegro, 1996]. Esta aplicación se realizó en un ordenador equipado con un procesador Xeon de Intel y se emplearon unos pocos Aplicaciones del generador de mallas 57 (a) (b) Figura 2.22: Área rectangular de la Isla de La Palma: (a) malla inicial con 574 tetraedros invertidos y (b) malla válida resultante tras diez iteraciones del proceso de optimización. 58 Generación de mallas para dominios definidos sobre orografías irregulares 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 O Malla original Malla optimizada Jl I I L j_ _L 10000 20000 30000 40000 50000 60000 70000 80000 Figura 2.23: Curvas de calidad de la malla generada y la malla obtenida tras diez iteraciones del método de optimización. La función ^«(e) es una medida de calidad del tetraedro e. segundos en construir la malla. La complejidad de cada paso de optimización es también lineal; en la práctica también fueron pocos segundos los empleados en hacer diez iteraciones del proceso de optimización, usando el método BFGS [Bazaraa et al., 1993] para minimizar la función objetivo. Finalmente, la tabla 2.1 muestra los resultados de la aplicación del proceso de desenredo y suavizado en tres nuevas mallas enredadas de este mismo dominio. A partir de la malla de la figura 2.22(b), se mueven aleatoriamente un cierto número de nodos interiores una distancia acotada h/3 para obtener Ninv = 431, 3441 y 10051 tetraedros invertidos, respectivamente. Se observa que solamente se necesitan unas pocas iteraciones, lunt, para desenredarlas. Además, la calidad de las mallas tras 5 iteraciones de suavizado son prácticamente idénticas. En general, la calidad final de una malla tras un número suficiente de iteraciones del proceso de optimización depende más de su topología que del estado inicial de enredo. De hecho, siempre existe un límite de calidad definido por las conectividades de los nodos de la malla. Enlace entre códigos 59 Malla inicial ^^nod 16504 Ntet 81068 N431 3441 10051 QK 0.723 0.667 0.547 SUS ^unt 2 3 4 „min 0.112 0.112 0.118 QK 0.735 0.735 0.734 Tabla 2.1: Resultados comparativos del uso de la técnica de desenredo y suavizado simultáneo para el detalle de La Palma, usando 5 iteraciones de suavizado tras lunt iteraciones de desenredo. 2.9. Enlace entre códigos La implementación de las técnicas de mallado descritas anteriormente se ha realizado aprovechando códigos que ya estaban disponibles. En concreto, utilizamos un generador de mallas 2-D, el programa de refinamiento/desrefinamiento en 2-D (NEPTUNO) y un generador de malla 3-D basado en la triangulación de Delaunay. Los anteriores programas fueron creados con el lenguaje de programación FORTRAN 77 y no fueron pensados para interactuar entre ellos. A pesar de que han sido modificados a lo largo del tiempo, bien para incorporar ciertas mejoras o bien para adaptarlos a las sucesivas plataformas informáticas donde se han ido ejecutando, siguen teniendo algunas características que los hacen muy poco amigables con el usuario: • Los ficheros de entradas son complejos y hay que atenerse a un formato tan estricto como incómodo para el usuario. Una de las mayores dificultades era la adaptación manual de la salida de un programa para que sirviera como entrada para el siguiente, en la secuencia de generación de la malla 3-D. • Las salidas que proporcionan carecen de flexibilidad; como ejemplo, baste decir que no se puede elegir el nombre ni el directorio del fichero de salida de ninguno de los programas. • Debido a que FORTRAN 77 no tiene capacidad de memoria dinámica, en ocasiones es necesario recompilarlos para adaptarlos a la capacidad de la máquina donde se ejecutan y a los requerimientos del problema concreto. Todo esto justifica la creación de un programa que, enlazando todos los anteriores, presente al usuario una interfaz coherente y cómoda. Los detalles engorrosos. 66 Modelización de campos de viento 3.2. Construcción del campo inicial Para la construcción del campo inicial partimos de los valores de la velocidad del viento y de su dirección obtenidos en las estaciones de medida. El campo inicial vo se construye en tres etapas. En primer lugar, se calcula mediante interpolación horizontal el valor de VQ en los puntos del dominio situados a la misma altura Zs (sobre el terreno) que las estaciones de medida. Con esta información se realiza una extrapolación vertical para definir el campo de velocidades en todo el dominio. Finalmente, la componente vertical del campo de velocidades es corregida en el entorno de posibles fuentes de emisión de contaminantes (chimeneas) con el fin de simular el movimiento de salida de los gases. 3.2.1. Análisis de las medidas de las estaciones Los datos de viento se toman de estaciones de medida ubicadas en el dominio de estudio. Cada estación de medida proporciona la velocidad (en m/s) y dirección del viento a una altura Zg sobre el nivel del terreno (típicamente 10 metros). La dirección del viento viene dada en grados sexagesimales medidos en sentido horario y tomando como referencia la dirección norte. Así el norte se corresponde a O grados, el sur a 180 grados, el este a 90 grados y el oeste a 270 grados. A efectos de cálculo en el modelo, es necesario obtener el ángulo medido en sentido antihorario, tomando como referencia el semieje positivo horizontal. Por otro lado, como las estaciones miden el viento en intervalos discretos de tiempo, en general es necesario interpolar las medidas para calcular el viento en un instante concreto. El cambio se realiza de la siguiente manera (ver figura 3.1), a = ^ ~ <p9 a' = 270-a o de una manera más simplificada a' = 180 + e' + (l)g donde, d es el ángulo tomado por las estaciones de medida en grados norte en sentido horario y 6' es su complementario. Asumimos que el viento sobre la superficie gira un ángulo (f)g con respecto a la dirección del viento geostrófico Vg,a es Construcción del campo inicial 67 Norte Vest Figura 3.1: Ángulo de desviación del viento sobre la superficie respecto al viento geostrófico. el ángulo respecto a la vertical y a' es el ángulo que se desea obtener en sentido antihorario sobre la horizontal. 3.2.2. Interpolación horizontal La técnica más común de interpolación se formula en términos de la inversa de la distancia al cuadrado entre el punto y la estación de medida [Winter et al., 1995]. Sin embargo, otros autores usan simplemente la altitud de los puntos de medida [Palomino y Martín, 1995]. Aquí se propone una fórmula que tiene en cuenta ambas consideraciones, 1=1 "n N E Vo{ze) = e ^^ + (1 - e)— — y — y —— (3.14) El valor de Vn corresponde a la velocidad observada en la estación n, A'' es el número de estaciones utilizadas en la interpolación, d„ es la distancia horizontal desde la estación n al punto donde estamos calculando la velocidad del viento. 68 Modelización de campos de viento \Ahn\ es la diferencia de altura entre la estación n y el punto en estudio, y e es un parámetro de peso que toma valores entre O y 1. Cuando e —> 1 aumenta la importancia de la distancia horizontal desde cada punto a las estaciones de medida. Esta aproximación se emplea en problemas con una orografía regular o en análisis bidimensionales. De manera análoga, si £ -^ O es entonces la diferencia de altura entre cada punto y las estaciones de medida la que resulta determinante, en detrimento de la distancia horizontal. Esta segunda aproximación es la que se usa cuando la orografía del terreno es irregular. En la práctica, las regiones geográficas estudiadas suelen combinar zonas de orografía irregular con otras de orografía mas regular, por lo que tomar un valor intermedio para e suele ser lo más apropiado. 3.2.3. Extrapolación vertical El viento se desarrolla, en primer lugar, como consecuencia de diferencias espaciales en la presión atmosférica. Estas diferencias de presión normalmente son causadas por una diferente absorción de la radiación solar. En un plano horizontal, el viento fluye de las zonas de alta presión a zonas de baja presión y verticalmente de zonas de baja presión a zonas de alta presión. La velocidad del viento es proporcional al cambio de presión por unidad de distancia o gradiente de presión. Las zonas con presiones similares se representan en los mapas meteorológicos unidas mediante líneas imaginarias denominadas isóbaras. Cuanto más juntas están unas isóbaras, mayor será la fuerza del viento. Un segundo factor que afecta el movimiento del aire es la fuerza de Coriolis, debida a la rotación terrestre. El parámetro / = 29 sen ¡pi se denomina parámetro de Coriolis, siendo © = 7.292 x 10~^ s~^ la velocidad de rotación de la Tierra y (pi la latitud. Se considera positiva en el hemisferio norte, nula en el ecuador y negativa en el hemisferio sur. En tercer lugar puede aparecer una aceleración centrípeta, cuando el viento gira en torno a un centro. Por último, aparece la fricción debida al desplazamiento del aire. Los vientos influenciados por el gradiente de presión y la fuerza de Coriolis se denominan vientos geostróficos. 3.2.3.1. Estratificación atmosférica En este modelo consideraremos una división de la capa más baja de la atmósfera en distintas subcapas, en las que la extrapolación vertical de las velocidades Construcción del campo inicial 69 de viento se realiza de forma diferente, como puede observarse en la figura 3.2. Así, la capa límite planetaria está situada a una altitud Zpu sobre el nivel del terreno, y es la capa de la atmósfera, situada por debajo de la atmósfera libre, que está afectada directamente por la fricción de la superficie de la tierra (conocida también como capa límite atmosférica). La altitud de la capa límite planetaria Zpbi sobre el terreno se ha tomado tal que la dirección e intensidad del viento es constante a partir de esa altura [Sempreviva, 1996], zpu = ^ (3.15) siendo 7 una constante comprendida entre 0.15 y 0.45 que depende de la estabilidad de la atmósfera y t/* la velocidad de fricción que será definida más adelante a partir de los valores obtenidos en la interpolación horizontal. La capa de mezcla, también llamada capa límite convectiva, es la capa límite atmosférica sujeta a fenómenos convectivos causados por el calor superficial. El aire está bien mezclado, es decir el viento y el potencial de temperatura son prácticamente constantes con la altura. La altitud de la capa de mezcla hm se considerará igual a Zpu para condiciones neutras e inestables. En condiciones estables se aproxima por hm - y ;/^ (3.16) / donde usualmente se toma 7' — 0.4 [Zannetti, 1990] y L es la longitud de MoninObukov, que se calcula a través de la fórmula de Liu [Ratto, 1996], \ = a^o (3.17) con a y b, definidas por la clase de estabihdad de Pasquill (ver tabla 3.2). La capa superficial, localizada a una altura Zgi sobre la superficie, es la capa baja, dentro de la capa límite planetaria, inmediatamente adyacente a la capa de la superficie de la tierra, en la que la fuerza de arrastre de fricción es dominante. Conocido el valor de la altura de la capa de mezcla hm, la altitud de la capa superficial se suele fijar en [Zannetti, 1990], ^si = ^ (3.18) 10 Modelización de campos de viento 3.2.3.2. Estabilidad atmosférica El concepto de estabilidad atmosférica está relacionado tanto con la turbulencia atmosférica como con el gradiente vertical de temperatura y las situaciones de inversión térmica. La estabilidad atmosférica nos proporciona una medida cualitativa de las variaciones de la densidad del aire, debidas a los cambios de presión y temperatura y que influyen en determinados movimientos atmosféricos. Las condiciones atmosféricas pueden clasificarse como: • Estable: Si una masa de aire sube se encontrará rodeada de aire más caliente y, por tanto, menos denso que ella, lo que la hará bajar; y si baja, se encontrará rodeada de aire más frío (más denso), y tenderá a subir. Esta tendencia que tiene el aire de permanecer en la misma capa es lo que se denomina estabilidad de la estratificación atmosférica. • Inestable: En condiciones inestables la temperatura potencial disminuye con la altura, incrementándose los movimientos verticales, es decir si el aire sube se encontrará rodeado de aire más frío y denso que él, y tenderá a seguir subiendo; y si baja se encontrará con aire más caliente y ligero, y tenderá a seguir bajando. • Neutra: Si un volumen de aire (después de un desplazamiento vertical en una capa atmosférica sin mezclar con el aire circundante) experimenta una fuerza neta vertical nula, los movimientos ascensionales no se verán perturbados por el gradiente térmico, entonces la capa atmosférica se asume neutralmente estratificada. Bajo tales condiciones, dicho volumen ni tiende a volver a su posición original (estratificación estable) ni acelera alejándose de ella (estratificación inestable). La estabilidad atmosférica puede ser caracterizada mediante la tabla (3.1) definida por Pasquill. 3.2.3.3. Perfil vertical de velocidades de viento Como se muestra en la figura 3.2, se considera un perfil logarítmico-lineal [Lalas y Ratto, 1996] en la capa límite planetaria, que tiene en cuenta la interpolación horizontal [Montero et al., 1998], el efecto de la rugosidad en la intensidad y dirección del viento, y la estabilidad del aire (neutra, estable o inestable) según la clasificación de Pasquill. En la capa superficial se construye un perfil logarítmico Construcción del campo inicial 71 Clase de estabilidad de Pasquill Insolación Noche Velocidad del viento en la superficie (ni/s) Fuerte Moderada Ligera Cubierto ó > 4/8 < 3/8 nubes nubes A A-B B C C A-B B B-C C-D D B C C D D - E D D D - F E D D <2 2-3 3-5 5-6 >6 Para A-B, tomar la media de los valores de A y B, etc. Tabla 3.1: Clases de estabilidad de Pasquill según la velocidad del viento en la superficie y la insolación. Insolación fuerte corresponde al mediodía soleado de mitad de verano en Inglaterra; insolación ligera a condiciones similares en mitad del invierno. La noche se refiere al periodo que va desde una hora antes de ponerse el sol hasta una hora después de salir. La clase neutra D debería ser usada también, a pesar de la velocidad del viento, para cielos cubiertos durante el día o la noche, y para cualquier condición del cielo durante las horas precedente y siguiente de la noche definida anteriormente. Clase de estabilidad de Pasquill A (Extremadamente Inestable) B (Moderadamente Inestable) C (Ligeramente Inestable) D (Neutra) E (Ligeramente Estable) F (Moderadamente Estable) a -0.08750 -0.03849 -0.00807 0.00000 0.00807 0.03849 b -0.1029 -0.1714 -0.3049 0.0000 -0.3049 -0.1714 Tabla 3.2: Coeficientes a y b para el cálculo de la longitud de Monin Obukov según la clase de estabilidad de Pasquill. de velocidades de viento definido por, Mz) = i-( log —-$r, Zo < Z < Zsl (3.19) donde VQ es la velocidad del viento. A; ~ 0.4 es la constante de von Karman y z es la altura sobre el terreno del punto estudiado. El término v* representa la velocidad de fricción. En el flujo turbulento atmosférico las fuerzas que se oponen al movimiento están caracterizadas por la acción que ejercen las rugosidades o 12 Modelización de campos de viento Viento Geostrófico Capa de Mezcla ^^o(^) = PÍZ) vo{zsi) + [1 - p{z)]vg «oW = f(^-*".) Figura 3.2: Perfil vertical de viento definido sobre cada capa de la estratificación atmosférica. asperezas propias de la orografía del terreno. La velocidad de fricción se obtiene en cada punto a partir de las medidas interpoladas a la altura de las estaciones {interpolación horizontal), k Vo{Ze) lnf^-$^(ze) (3.20) Asimismo, ZQ corresponde a la longitud de rugosidad de la zona. El concepto de longitud de rugosidad viene a definir una altura por encima del terreno diferente de z = O, donde, en teoría de la capa superficial, la velocidad del viento es cero. El valor de ZQ depende de las características del terreno. Una forma de estimarla es mediante valores estándar para diferentes tipos de terreno [McRae et al., 1982]; ver figura 3.3. Otros autores la definen como 2:0 = ^) donde e es la altura media de los obstáculos existentes en la zona de estudio. Por último, ^rn es una función que depende de la estabilidad del aire [Zannetti, Construcción del campo inicial 73 fe? ^m) ^ Centros de ciudades con edificios muy altos ^ Centros de grandes poblaciones, ciudades Centros de pequeñas poblaciones Alrededores de poblaciones ^ Muchos arboles, setos, pocos edificios Muchos setos Áreas escarpadas o montañosas - Bosques Reglón de nivel medio de bosque - Pocos arboles, verano Arboles aislados Hierba sin cortar Pocos arboles, invierno Hiertia cortaela ( 3 cm) > Tienes de cultivo Hierba alta( 60cms) campos cultivados Aeropuertos Llanuras de hierbas medianas _ Superficie natural nevada (tierras de cultivo) Viento mar adentro en zonas costeras \ Mar abierto en calma Desierto (llanura) > Grandes extensiones de agua Planicie cubierta de nieve Terreno ondulado Hielo, planicie (llanura) enlodada. Figura 3.3: Longitud de rugosidad: valores aproximados de ZQ para distintos tipos de terreno definidos por McRae (1982) 7^ Modelización de campos de viento 1990], $m = O (neutra) z $TO = -5— (estable) 1J *m = log donde ^m + l^ /'^m + 1 7¡- — 2 arctan 9m-\- — (inestable) Bm = il-16^)^/^ (3.21) El viento geostrófico es una buena aproximación al viento real con flujo uniforme en la alta atmósfera (atmósfera libre), donde la fricción y aceleraciones no son importantes. La forma general de la expresión usada para calcular el viento geostrófico es la bien conocida ley de resistencia geostrófica (geostrophic drag low) Ratto [1996]. Los valores de los coeficientes A y B tienden a ser ~ 1.8 y ~ 1.5 respectivamente, que son los valores aceptables para condiciones neutras de estabilidad atmosférica. El viento en la superficie se supone que gira un ángulo (¡)g con respecto a Vg, dado por la relación En nuestro modelo, desde Zgi hasta Zpu se realiza una interpolación lineal en p{z) con el viento geostrófico Vg Vo{z) = p{z) Vo{zsi) + [1 - p{z)]vg con Zsi < z < Zpbi (3.24) donde p(z) es p(^) = 1 _ (J-ZIILVÍS _ 2 ^~^^' ) (3.25) \Zpbl - ZslJ V ^pbl ~ ^sl J Finalmente, este modelo considera VQ{Z) = Vg si z> Zpbi (3.26) vo{z) = 0 si z < ZQ (3.27) Construcción del campo inicial 75^ 3.2.4. Corrección de la componente vertical en la trayectoria de la pluma La idea es incorporar al perfil inicial de velocidades de vientos, donde usualmente sólo se calculan las componentes horizontales de la velocidades debido a la ausencia de medidas de la componente vertical en las estaciones, una componente vertical no nula en la trayectoria de una posible pluma de contaminantes originada por una fuente emisora (chimenea). De esta forma se pretende simular el campo de velocidades del fluido compuesto por dos aportaciones: la del viento y la de expulsión del gas contaminante de la chimenea. Los actuales modelos de pluma gaussiana permiten aproximar los valores de la altitud efectiva ZH de la pluma y la distancia horizontal df desde el centro de la superficie de salida de la chimenea hasta donde se alcanza ZH, en función de las características de la emisión, del viento y de la estabilidad atmosférica [Boubel et al., 1994]. Los gases que salen de una chimenea alcanzan una altura superior a la de la chimenea cuando estos son de menor densidad que el aire del entorno (elevación por flotación) o bien son expulsados a una velocidad suficiente que les proporciona una energía cinética (elevación por momento). La elevación por flotación se denomina en ocasiones elevación térmica ya que la causa más común de disminución de la densidad es el aumento de temperatura. Para estimar la altura efectiva de la pluma se utilizan las ecuaciones de Briggs [Briggs, 1969, 1971, 1972, 1973, 1975]. Consideraremos el viento VQ resultado de la extrapolación anterior conocido en todo el dominio. La altitud física Zc de la chimenea se sustituye en la práctica por la altitud z'^, algo menor que la primera cuando la velocidad de salida de los gases de la chimenea Wc es menor que 1.5 veces la velocidad del viento {Stack Downwash), z'c = Zc si Wc > 1.5 \vo {xc, Ve, Zc)\ (3.28) z'^^Zc + 2Dc [{wj \VQ {XC, Ve, Zc)\) - 1.5] si Wc < 1.5 \vo {xc, Ve, Zc)\ (3.29) siendo (xcVcZc) y De, las coordenadas del centro y el diámetro de la superficie de salida en la chimenea, respectivamente. Consideremos en primer lugar que predomina el fenómeno de elevación por w flotación, es decir, -zr-7 rr < 4. Para ello es necesario definir el parámetro \VQ[XC, Ve, Zc)\ 82 Modelización de campos de viento .-..w^t^=2(Zjj-z;) •^ t Figura 3.6: Curvas de la evolución de componente vertical de la velocidad con t para los casos relativos a la figura 3.5. y la aceleración, ÜQ = -Wc tf (3.66) Por tanto, la velocidad vertical en un punto de altitud z en función del tiempo vendrá dada por, (3.67) wo{t) =IÍ;C í 1- — siendo ya que t = tf[l-Jl 2 {z - O Wctf 1 t Z = Zc + Wct 11- - — (3.68) (3.69) En este caso, la componente vertical se modifica en los puntos de un cilindro recto de base la superficie de salida de gases en la chimenea y altura ZH — ZcPor tanto, sólo se consideran aquellos puntos {xo,yo,zo), con Zc < ZQ < ZH, que cumplan la Discretización mediante elementos finitos 83 " ^ ^ y^ -..w^t^=2(z^-z;) Figura 3.7: Curvas de la aceleración en los casos relativos a la figura 3.5. condición, ^J{xc - xof + {Vo - yof < -^ (3.70) a los cuales se le asigna el valor ifo(ío), siendo ÍQ el valor de í en ^ = ^o3.3. Discretización mediante elementos finitos Para la discretización mediante elementos finitos de la formulación clásica del problema dada en (3.9), (3.10) y (3.11) se ha utilizado una malla de tetraedros, generada mediante las técnicas ya descritas, e interpolación lineal. Nótese que en la formulación variacional del problema, las integrales de contorno en la parte de la frontera con condición Neumann se cancelan utilizando la ecuación (3.11) y las de tipo Dirichlet se eliminan anulando la correspondiente función test. Esto conduce a un conjunto de matrices elementales de dimensión 4x4 asociadas al elemento 0,^, siendo -^j la función de forma correspondiente a su z-ésimo nodo, i = 1,2,3,4, definidos en el elemento de referencia íle y I J| el jacobiano de 84 Modelización de campos de viento H c j ifcM Tita ni^^Sgtg^^Sj^^BBBs i /// // ^" . — jjff IfaV^ Figura 3.8: Zona de corrección de la componente vertical de la velocidad del fluido con elevación por flotación. la transformación de ííe a. Cíe, ,di¡)idC d'^idr] d-^id^p d^jdi dipj drj dipj d<f Jíle + + •X- + + )+ + + d^ dx drj dx díp dx d$, dx dr] dx díp dx' (3.71) )} • |J| d^drjdif a^ae dAdv di^dip Mdi di^dv d¿_d^ ^ d^ dy dr] dy d^ dy^^ d^ dy dr] dy dip dy' + :)(- + + Th ^ di dz dr] dz d^p dz'^ d^ dz dr] dz dcp dz y de vectores elementales de 4 x 1, Resolución del sistema de ecuaciones 85 Figura 3.9: Zona de corrección de la componente vertical de la velocidad del fluido con elevación por momento. {b-}i / J_r rdipjd^ dipidí] dj^idip 'Th ^"°^ di dx dridx d<f dx^ a^a£ d^dji d^id<p d^ dy dr) dy dip dy (3.72) 3.4. Resolución del sistema de ecuaciones La aplicación del método de elementos finitos en este tipo de problemas implica la resolución de grandes sistemas de ecuaciones, que se caracterizan porque la matriz de coeficientes es escasa (sparse), ya que muchos de los términos de la matriz de los coeficientes son cero. Ambas características, tamaño y elevado nú- 86 Modelización de campos de viento mero de términos nulos, justifican la elección de métodos iterativos de resolución de sistemas de ecuaciones. Se ha optado por representar las matrices en un formato compacto denominado morse, lo que permite un gran ahorro de memoria en comparación con la representación estándar. En este formato se utiliza un vector para almacenar los valores no nulos de la matriz, un vector de enteros con tantos elementos como valores no nulos tiene la matriz donde se almacenan las columnas a las que pertenece cada valor, un vector de enteros con tantos elementos como la dimensión de la matriz más uno y un entero que contiene la dimensión original de la matriz. Así, el tamaño de memoria necesario para la representación de una matriz escasa es del orden del número de elementos no nulos que contiene, mientras que la representación estándar sería de orden m?, siendo n la dimensión de la matriz. En nuestro problema tenemos dos tipos de condiciones de contorno. Las de tipo Neumann se introducen en la formulación variacional. En cambio, aquellos nodos con condición de tipo Dirichlet pueden ser eliminados del sistema de ecuaciones. Para ello eliminamos de la matriz del sistema la fila y la columna correspondiente al nodo ya que el valor de la condición Dirichlet es, en este caso, cero. De esta manera el sistema se convierte en simétrico y, por tanto, puede emplearse el método del gradiente conjugado para resolverlo, con el consiguiente ahorro de tiempo de ejecución, ya que la aplicación del gradiente conjugado resulta más eficiente que otros algoritmos basados en los subespacios de Krylov para sistemas no simétricos. Esto es debido a que utiliza un sólo producto matriz-vector frente a los dos de los otros algoritmos. Además se han implementado una serie de precondicionadores que permiten mejorar la convergencia de aquellos sistemas que estuvieran mal condicionados. Una vez calculada la solución del sistema "reducido", se completa con las aportaciones de los nodos con condición Dirichlet, obteniendo de esta manera la solución completa del sistema de ecuaciones. Una vez obtenido (f) como solución del sistema de ecuaciones, se calcula el campo el campo de viento usando la ecuación (3.7). 3.4.1. Precondicionamiento La convergencia de los métodos basados en los subespacios de Krylov, y en particular la del gradiente conjugado, mejora con el uso de las técnicas de precondicionamiento. Éstas consisten generalmente en cambiar el sistema original Resolución del sistema de ecuaciones 87 Ax = b por otro de idéntica solución, de forma que el número de condicionamiento de la matriz del nuevo sistema sea menor que el de A, o bien tenga una mejor distribución de autovalores. Generalmente, se considera una matriz de precondicionamiento M~^, siendo M una aproximación de A, esto es, M-^Ax = M-^b tal que, K (M~^A) < K (A). El menor valor corresponde aM. = A, K (A~^A) == 1, que es el caso ideal y el sistema convergería en una sola iteración, pero el coste computacional del cálculo de A~^ equivaldría a resolver el sistema por un método directo. Se sugiere que M que sea una matriz lo más próxima a A sin que su determinación suponga un coste elevado. Por otro lado, la matriz M debe ser fácilmente invertible para poder efectuar los productos M"^ por vector que aparecen en los algoritmos precondicionados sin excesivo coste adicional. Dependiendo de la forma de plantear el producto de M~^ por la matriz del sistema obtendremos distintas formas de precondicionamiento. Estas son, M~^Ax = M~^b (Precondicionamiento por la izquierda) AM~^Mx — b (Precondicionamiento por la derecha) (3.73) M]'^AM^^M2X — M^^b (Precondicionamiento por ambos lados) si M puede ser factorizada como M = M1M2. El campo de posibles precondicionadores es muy amplio. Algunos de los más usados son el de Jacobi, SSOR, ILL^(O) y el Diagonal Óptimo. 3.4.1.1. Precondicionador de Jacobi Surge comparando la fórmula de recurrencia para la solución que resulta de aplicar el método de Richardson, cuya relación de recurrencia viene dada por Xj+i = x¿ + /x(b — Axj), con /^ > O, al sistema precondicionado con la fórmula correspondiente que se obtiene aplicando el método de Jacobi al sistema sin precondicionar. De la aplicación del método de Richardson al sistema precondicionado M~^Ax = M~^b, se obtiene para el cálculo de los sucesivos valores de la solución, Xj+i ^Xi + iJ, (M~^b - M~^Ax¿) 88 Modelización de campos de viento Multiplicando por la matriz de precondicionamiento M, queda, Mxi+i = Mxi + ^ (b - Ax¿) (3.74) Por otro lado, descomponiendo la matriz del sistema enA = D — E — F, (D matriz diagonal formada por los elementos de la diagonal de A y E y F matrices triangulares inferior y superior respectivamente), y utilizando el método de Jacobi para la resolución del sistema, se obtiene, x¿+i = D-^ (E + F) x¿ + D-^b que multiplicando por D y operando resulta. Dxj+i = Dxí + (b - Axj) (3.75) Comparando las expresiones de recurrencia finales de ambos métodos, se observa que el método de Jacobi aplicado al sistema sin precondicionar, equivale al de Richardson, con a = 1, menos robusto y más simple, cuando este se aplica al sistema precondicionado con la matriz diagonal D = diag{A). 3.4.1.2. Precondicionador SSOR Si aplicamos el método SSOR al sistema sin precondicionar, considerando la descomposición de la matriz AenA = D — E — F, como en el caso anterior, y siendo u el parámetro de relajación, se obtiene, U}{2-Lü) 1 (D-CJE) D-' (D-a;F) x^+i (D-uE) D' (B-LvF) Xi + (b - Ax^) uj{2-u) con lo que resulta como matriz de precondicionamiento, que en el caso de sistemas simétricos, podemos expresarla como. M = (D-a;E)D-^/^ y/uji2-üj) (3.76) (3.77) (3.78) Resolución del sistema de ecuaciones _55 3.4.1.3. Precondicionador ILL^^*^) Resulta de la aproximación de A por una factorización incompleta LL^, conservando las mismas entradas nulas en las matrices triangulares L y L^, A =« ILL^(O) = M (3.79) donde kj = Iji son las entradas de las matrices triangulares incompletas tal que, hj = 0 si Uij = O (3.80) {A-LL'^},. = 0 si a¿j7¿0 (3.81) Es decir, que los elementos nulos de la matriz del sistema siguen siendo nulos en las posiciones respectivas de las matrices triangulares para no incrementar el coste computacional. 3.4.1.4. Precondicionador diagonal óptimo Se trata de un caso particular de inversa aproximada de estructura diagonal. Si se ha precondicionado por la izquierda, el mejor precondicionador diagonal del sistema resulta, M = diag I -^, -^,..., -^^ I (3.82) IP^AII IIP^AII lle^AII 1*^1-^112 11*^2-^112 Il*^n-^ll2 ||MA -1||' =n-y —^ (3.83) 1=1 \\*^i -^112 De forma similar se puede definir el precondicionador por la derecha. 3.4.2. Algoritmo de gradiente conjugado precondicionado El método del Gradiente Conjugado está basado en una técnica de proyección ortogonal sobre el subespacio de Krylov JCk (A; r^) donde ro es el residuo inicial, y fundamentado en el algoritmo D-Lanczos de la siguiente forma. Una aproximación Xj+i puede ser expresada, Xj+i = Xj- + ajPj (3.84) 90 Modelización de campos de viento Algoritmo 3.1 Gradiente Conjugado precondicionado (PCG). Aproximación inicial XQ. TQ = b — AXQ; Resolver MZQ = TQ, po = ZQ; Mientras || TJ || / || ro ||> £ {j = 0,1,2,3,...) Hacer ' JAp^.,p,->' Xj+i - Xj + Q:ÍPJ; rj+i=rj-ajAp^.; Resolver Mzj+i= r^+i; ^ (rj+i,Zj+i). Pj+i = Zi+i + ^jPj; Fin Mientras entonces, los vectores residuos deben satisfacer que, rj+i = r^- - ajAp^. (3.85) Estos vectores residuos son ortogonales, (rj-QíjAp^.,rj) =0 por tanto, «.• = 7^^ (3-86) Como las siguientes direcciones Pj+i son una combinación lineal de r^+i y Pj, después de reescalar los vectores p apropiadamente, se tiene que Pj+i = Tj+i + /3jPj (3.87) con lo cual, (Apj,rj> = (Ap^.,p,- - /?j-iPj-i) = (Ap^-,p,) ya que Ap^ es ortogonal a Pj-i. Entonces, teniendo en cuenta (3.86) y (3.87) obtenemos que, (rj+i,Ap^.) ^' (P.Ap,) Refinamiento adaptativo 91^ y como de (3.85), podemos escribir, Ap,.---(r,+i-r,) (3.88) Oij entonces, Pi = 1 (rj+i, (fj+i - rj)) _ (rj+i,rj+i) "j (Apj>Pi) (rj,rj) 3.5. Refinamiento adaptativo En la actualidad, la mayor parte de los programas que utilizan el método de elementos finitos se apoyan en técnicas adaptables basadas en una estimación del error cometido con nuestra solución numérica, o al menos en indicadores de error fiables que nos señalen los elementos que deben ser refinados o desrefinados en la malla. En la generación de mallas adaptables podemos considerar dos aspectos diferentes: la discretización del dominio atendiendo a su geometría o a la solución numérica, mediante estimadores de error o mediante indicadores de error. Existen muchas formas de abordar estos aspectos. La primera cuestión es: ¿mallas estructuradas o no estructuradas?. En este sentido, está claro que el uso de mallas no estructuradas nos proporciona más flexibilidad a la hora de mallar geometrías complejas utilizando un número óptimo de nodos. En este caso, los métodos más clásicos para la obtención de triangulaciones tridimensionales se basan fundamentalmente en algoritmos de avance frontal [Lohner y Baum, 1992] o en algoritmos basados en la triangulación de Delaunay [George et al., 1991; Escobar y Montenegro, 1996]. Una vez que se ha discretizado la geometría del dominio, la malla debe adaptarse atendiendo a las singularidades de la solución numérica. Este proceso implica la introducción (refinamiento) o eliminación (desrefinamiento) de nodos de la malla actual. Por esta razón, se debe definir una nueva malla. Los cambios pueden afectar a la malla actual de forma local o global, dependiendo del método de triangulación elegido. Diferentes estrategias de refinamiento han sido desarrolladas para triangulaciones en 2-D, y han sido generalizadas a 3-D. Si se ha optado por un refinamiento que afecte localmente a la malla actual, cabe plantease otra cuestión: ¿mallas encajadas o no encajadas?. La respuesta en este caso no es tan clara. El uso de mallas encajadas tiene varias ventajas importantes. Podemos con- 98 Modelización de campos de viento (a) Tipo I (b) Tipo II (c) Tipo Illa (d) Tipo Ill.b (e) Tipo IV (f) TipoV Figura 3.10: Clasificación de las subdivisiones de un tetraedro en función de los nuevos nodos (círculos huecos). Experimentos numéricos 99 3.6. Experimentos numéricos 3.6.1. Problema test El caso estudiado corresponde a un problema test de ajuste de campos de viento correspondiente a una zona cuadrada de 10000 x 10000 m?, donde la topografía del terreno viene dada por la función, Z Z"n {"i^y (y-vc \ «y )1 (3.90) siendo en nuestro ejemplo z^nax — 1500, Xc = yc = 5000, Sx = 1000 j Sy = 800. La situación de las estaciones y sus respectivas medidas de las velocidades de viento se indica en la tabla 3.3. Se ha considerado que la altitud de las estaciones es de 10 m sobre el terreno y que la componente vertical de la velocidad medida es nula. Estación m 1 2 3 4 5 6 7 8 9 Xm en m 0.0 5000.0 10000.0 0.0 5000.0 10000.0 0.0 5000.0 10000.0 Vm enm 0.0 0.0 0.0 5000.0 5000.0 5000.0 10000.0 10000.0 10000.0 Um en m/s 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Vm en m/s 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 Tabla 3.3: Situación de las estaciones y medidas consideradas. Otros valores de los parámetros del modelo considerados en esta aplicación son E = 0.5, 7 = 0.3, Y = 0.4, a = 0.1, k = 0.4, 0 = 28.6°, Ug = 0.0 m/s, Vg = 10.0 m/s, así como una atmósfera ligeramente estable. La figura 3.11 representa la malla inicial obtenida con el código desarrollado en el capítulo 2, que se ha utilizado para resolver este problema de viento. A partir de este resultado, utilizando el indicador de error (3.89) con p = 1, se realizaron dos etapas de refinamiento adaptativo de la malla con ^ = 0.8. Las sucesivas mallas refinadas obtenidas en este proceso se ilustran en las figuras 3.12 y 3.13. Se observa una mayor concentración de nodos en la zona de la montaña donde la velocidad del viento experimenta mayores cambios, tanto en dirección como en módulo. 100 Modelización de campos de viento Figura 3.11: Malla original. Figura 3.12: Malla correspondiente a la primera etapa de refinamiento. Experimentos numéricos 101 Figura 3.13: Malla correspondiente a la segunda etapa de refinamiento. Figura 3.14: Velocidades de viento a 10 m de altitud. 102 Modelización de campos de viento Figura 3.15: Velocidades de viento a 500 m de altitud. Figura 3.16: Velocidades de viento a 1500 m de altitud. Experimentos numéricos 103 Las figuras 3.14, 3.15 y 3.16 muestran las velocidades de viento obtenidas después del último refinamiento sobre diferentes planos horizontales. Concretamente a 10 m (cota de las 8 estaciones de medida a menor altitud), a 500 m y a 1500 m (altitud de la cima de la montaña). Como era previsible, para el valor de a utilizado, el viento tiende a bordear la montaña horizontalmente. 3.6.2. Efecto de una chimenea en el campo de velocidades Como aplicación práctica se estudia el efecto que produce la emisión de gases en el campo de velocidades de viento calculado en un dominio. Se parte de la malla construida en la sección 2.8.2 sobre una zona de la Isla de La Palma, donde se incorpora una chimenea a la discretización del domino. Se considera una velocidad de salida de los gases Wc — 30 m/s a una temperatura Te = 500° K y una temperatura ambiente de T = 293° K. La chimenea tiene altura h — 200 m y el diámetro de la superficie de salida de los gases es De = 20 m. Con estos valores, el modelo gaussiano utilizado en esta tesis genera una pluma en la que la elevación de los gases se produce por flotación. Tras aplicar seis pasos de refinamiento en la trayectoria de la pluma sobre la malla de tetraedros, la malla alcanza 34626 nodos y 184017 tetraedros (véase la tabla 3.4). En la figura 2.21 puede verse la situación de la chimenea sobre la zona de estudio. Las figuras 3.17 y 3.18 muestran el efecto de la emisión de la chimenea a distintas escalas. Etapa 0 1 2 3 4 5 6 Nodos 28387 28652 29996 33277 33322 34551 34626 Tetraedros 153085 154595 160960 177473 177685 183659 184017 Tabla 3.4: Refinamiento del dominio para la adaptación a la pluma. 104 Modelización de campos de viento Figura 3.17: Modificación del campo de velocidades por la salida de gases contaminantes. Experimentos numéricos 105 Figura 3.18: Detalle del campo de velocidades en el entorno de la chimenea. Capítulo 4 Estimación de parámetros La eficiencia de los modelos de masa consistente para ajuste de campos de viento depende en gran medida de ciertos parámetros que aparecen en las distintas etapas del proceso, especialmente de algunos de los que intervienen en la construcción del campo de viento inicial y de los módulos de precisión de Gauss. En general, los valores de estos parámetros se toman usando una serie de reglas empíricas. Nosotros hemos planteado su estimación de manera automática, tal que las velocidades observadas en las estaciones de medida sean regeneradas de la forma más exacta posible por el modelo, dando lugar por tanto a un problema inverso. Existen diversos métodos de resolución de problemas inversos relacionados con la estimación de parámetros. De entre ellos, se han elegido los algoritmos genéticos, por ser una herramienta robusta y flexible, que puede ser competitiva ya que los cálculos pueden paralelizarse. 4.1. Definición del problema Para el cálculo automático de ciertos parámetros del modelo de viento se plantea el siguiente problema inverso: de las A'^ estaciones de medida disponibles se toman A^^ como referencia; el resto, se utilizan para el cálculo del viento. El viento así obtenido se compara con el medido en las N^. estaciones de referencia. Para la estimación de los parámetros del modelo se procede a minimizar la diferencia entre los resultados obtenidos y las medidas observadas en las estaciones de referencia. Esta técnica una mejora sustancial sobre la propuesta en [Barnard et al., 1987] para la estimación de uno solo de los parámetros por un procedimiento de ensayo y error. Nosotros proponemos extenderlo a un total de cuatro de los parámetros 108 Estimación de parámetros que intervienen en el modelo y automatizar su cálculo. En primer lugar se considera el parámetro de estabilidad « = — = ^/;7^ (4.1) que se deriva del funcional (3.3) y cuyo mínimo no varía si se divide por a^. Hay que señalar que para a >> 1 predomina el ajuste de viento en la dirección vertical, mientras que para a « 1 el ajuste tiene lugar predominantemente sobre el plano horizontal. Por lo tanto la elección de a determina que el viento tienda a rodear los obstáculos o a sobrepasarlos. Diversos experimentos numéricos han hecho patente que el comportamiento de los modelos de masa consistente depende sensiblemente de la elección de los valores de a, por lo que se presta particular atención a este problema. Diversos autores han estudiado cómo parametrizar la estabilidad debido a que la dificultad en la determinación de los valores de a han limitado el uso de modelos de masa consistente en terrenos de orografía compleja. En [Sherman, 1978], [Kitada et al., 1983] y [Businger y Arya, 1974], los autores proponen tomar a = 10"'^, o sea, proporcional a la magnitud de w/u. Otros, como Ross et al. [1988] y Moussiopoulos et al. [1988], relacionan a con el número de Fronde, mientras que Geai [1985], Lalas et al. [1988] y Tombrou y Lalas [1990], proponen que el parámetro a varíe en la dirección vertical. Finalmente, Barnard et al. [1987] proponen un procedimiento para obtener a en cada simulación del campo de viento. La idea es usar N velocidades de viento observadas para obtener el campo de viento y usar las restantes A',, como referencia. Entonces se realizarían diversas simulaciones con distintos valores de a. El valor que más acerque el viento estimado al observado en las estaciones de referencia es el que se toma como valor del parámetro de estabilidad. Este método proporciona valores de a que sólo son válidos para cada caso particular y, por tanto, no proporciona valores válidos a priori para otras simulaciones. Aquí se estudia una versión del método propuesto en [Barnard et al., 1987], utilizando algoritmos genéticos como herramienta de optimización que permite una selección automática de a. El segundo parámetro que va a estimarse es el coeficiente de peso e (O < e < 1) de la ecuación (3.14), correspondiente a la interpolación horizontal de las medidas de viento observadas. Cuando £ —> 1 adquiere más importancia la distancia horizontal de cada punto a las estaciones de medida, mientras que para £ —>• O se da más peso a la distancia vertical entre cada punto y las estaciones Experimentos numéricos 115 (3 entre 0° y 45° en incrementos de 5° (véase la figura 4.3). En cada variación de (5 se realiza la estimación de los parámetros a, e, 7 y 7'. Como los AG tienen una naturaleza estocástica, repetimos cada uno de los experimentos de AG 10 veces, utilizando semillas aleatorias distintas en caxla ejecución con el fin de asegurarnos de que los resultados obtenidos no se deben al azar, sino que son consistentes. La repetición de cada experimento permite obtener la media (x), desviación típica (a) y coeficiente de variación (cr/x) de la estimación de cada parámetro. Que las medidas de dispersión sean pequeñas indica que los valores obtenidos están razonablemente alrededor de la media y, por tanto, que la estimación es consistente y el método de AG es adecuado para realizar el estudio. Observando las figuras 4.4 y 4.5 se puede concluir lo siguiente: • Los valores de la mejor evaluación dependen de la dirección del viento medida en las estaciones, mejorando los resultados de la estimación a medida que ¡3 se acerca a 45°. La mejora se produce de forma gradual, sin saltos bruscos. • Los valores obtenidos para a tras el ajuste dependen de la dirección del viento. • En la gráfica correspondiente a e con /5 = 0°, la desviación típica y el coeficiente de variación son muy grandes (este último está en torno al 50%). Esto se explica porque en este caso, al ser la velocidad del viento Í/„ medida en todas las estaciones idéntica en módulo y dirección, e se cancela en la ecuación (3.14) y no tiene influencia en el modelo. Por tanto, cualquier valor de e sería factible y la estimación de un valor concreto no tiene sentido. • Con esta excepción, tanto 7 como e permanecen prácticamente constantes independientemente de la dirección del viento. • La variación del parámetro 7' con la dirección del viento sigue un cierto paralelismo con la de a. A medida que el ajuste de viento es predominantemente horizontal (valores pequeños de a), la altura de la capa de mezcla disminuye (menores valores de 7'). Asimismo, cuando el ajuste comienza a dar más peso a la componente vertical del viento (aumento de a), la altura de la capa de mezcla aumenta. Estos resultados están de acuerdo con la definición dada para la capa de mezcla en el apartado 3.2.3.1. 116 Estimación de parámetros Gaussiana de 1500 m de altura Variación de Best Eva!, frente al ángulo /? del viento 0.14 0.12 0.1 0.08 •ü 0.06 co 0) m 0.04 0.02 - I r valor medio • desv. tiplea • coef. variac. • 00 05 10 15 20 25 30 Ángulo j3 del viento Figura 4.4: Mejores evaluaciones del primer problema test. 4.3.2. Influencia de la topografía en la estimación de parámetros Para averiguar si la topografía influye en la estimación de parámetros repetimos el experimento anterior con las mallas correspondientes a las topografías de 50 y 500 m. Obsérvese que el número de nodos y tetraedros de las tres mallas consideradas es del mismo orden (ver tabla 4.1). Los resultados obtenidos están representados gráficamente en las figuras 4.6, 4.7, 4.8, 4.9 y 4.10. Se deduce de estos resultados que los valores de los parámetros dependen de la topografía y observando las gráficas se puede deducir lo siguiente: • Como era de esperar, la estimación de parámetros resulta más complicada en topografías irregulares que en superficies más suaves. Por otro lado, esta dificultad se incrementa para las configuraciones menos naturales de viento en las estaciones, correspondientes en este caso a valores pequeños de /?. Este efecto también se produce en los experimentos del apartado anterior. • Según se observa en la figura 4.7 los rangos de variación del parámetro a en cada experimento son distintos para cada orografía. Obsérvese que, no obstante, los rangos de variación obtenidos para este parámetro son pequeños si se comparan con los propuestos por los diversos autores. • Para /3 = O, cualquier valor de e es factible ya que no afecta al cómputo de viento interpolado. En el caso de la gaussiana de altura 1500 m, e es Experimentos numéricos 117 Gaussiana de 1500m de altura Varíadón de a frente al ángulo 6 def viento 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 1 1 1 - 1 1 1 valwmec&o ^^^— desv. típica --^— ~ coef. variac. 1 - - - - - - a 1 « o s i 1 lA) 1 0.9 0.8 0.7 0.6 0.5' 0.4 0.3< 0.2 0.1 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 i- •s s ° •s ^ H k rí 0.5 0.45 0.4 0.35 n.1 0 2.5 0.2 0.15 01 0.05 0 15 20 25 Ángulo fJ del viento Gaussiana de 1500m de altura Varíadón de e frente at ángulo ¡3 del viento - - - . ^ \ V! d co lor medie sv. típica - - - 15 20 25 30 Ángulo ¡3 del viento Gaussiana de 1500m de altura Varíadón de 7 frente al ángulo 3 del viento w or medio di sv. típica co f, variac. 15 20 25 30 Anguk) ¡3 del vierrto Gaussiana de 1500m de altura Variación de 7/ frente al ángulo /? del viento w or medio • di sv. típica - co f. variac- - 15 20 25 30 Ángulo fj del viento Figura 4.5: Resultados del primer problema test. 118 Estimación de parámetros prácticamente igual a 1, lo que implica una interpolación horizontal en la que predominamente van a influir las distancias horizontales de cada punto a las estaciones de medida. Esta situación empieza a cambiar a medida que disminuimos la altura y aumentamos /5. Así, por ejemplo para 50 m, en la interpolación horizontal se da más peso a la diferencia de cota entre los puntos y las estaciones. Esta tendencia se ve acentuada para valores altos de p. Esto parece indicar que el valor de e depende no sólo de la orografía del terreno sino además de la dirección en que esa orografía es atacada por el viento. Por tanto, queda en parte justificada la idea de ponderar tanto el efecto de la distancia horizontal como la diferencia de cotas en la interpolación horizontal. • 7 alcanza el valor máximo permitido y permanece prácticamente invariable enlos tres experimentos, incluso para distintos ángulos de ataque. • Para sacar conclusiones sobre 7' habría que realizar más experimentos. No obstante, el paralelismo con a apuntado en el experimento anterior se sigue manteniendo para las otras dos alturas del obstáculo. Experimentos numéricos 119 0.14 í 0.12 0.1 c o 5 0.08 (O ^ 0.06 (O <D CQ 0.04 0.02 O' Gaussiana de 50m de altura Variación de Best Eval. frente al ángulo /3 del viento T" "~i r valor medio • desv. tiplea - coef. variac. (dividido por 10) - 00 05 10 15 20 25 30 Ángulo /? del viento 35 40 45 0.045 0.04 0.035 0.03 0.025 0.02 m 0.015 0.01 0.005 •'"••v,. . o UJ 00 05 Gaussiana de 500m de altura Variación de Best Eval. frente al ángulo /3 del viento - - " 1 1 1. 1 Mwraryss; 1 ==*=: lili valor medio ^^— desv. tiplea ••»•• - coef. variac. (dividido por 10) I / t 11 ' i - i i f 1 1 1 1 i t i ~ -i ^—T t—í 10 15 20 25 30 Ángulo p del viento 35 40 45 íu <0 -I > <D ffi 0.14 0.12 0.1 Ü.Ü8 II OH 0.04 0.02 Gaussiana de 1500 m de altura Variación de Best Eval. frente al ángulo /3 del viento 1 1 1 valor medio ^^^ desv. tiplea ••»•• coef. variac. 1 00 05 10 15 20 25 30 Ángulo P del viento 35 40 45 Figura 4.6: Comparativa de las mejores evaluaciones para las tres orografías. 120 Estimación de 'parámetros i 0.5 0.45 0.4 0.35 0.3 0.25 Gaussiana de 50m de altura Variación de a frente al ángulo /3 del viento 1 r valor medio • desv. tipica • coef. variac. • 15 20 25 30 Ángulo /? del viento 35 40 45 E i 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 - 00 Gaussiana de 500m de altura Variación de a frente al ángulo ¡3 del viento M^S 05 10 15 20 25 30 Ángulo /3 del viento 1 r valor medio ' desv. tipica • coef. variac. • / / - 40 45 Gaussiana de 1500m de altura Variación de a frente al ángulo /? del viento 15 20 25 30 Ángulo j3 del viento Figura 4.7: Comparativa del parámetro a para las tres orografías. Experimentos numéricos 121 I I % i Gaussiana de 50m de altura Variación de e frente al ángulo ¡3 del viento 0.4 0.35 0.3 0.25- < 0.2 0.15 0.1 0.05 0 - K \ \ - \ 1 1 1 1 1 1 1 1 valor medio —^— desv. tipica -•»-• - coef. variac. (dividido por 10) | • •• ->" « -• .: III 00 05 10 15 20 25 30 Ángulo /? del viento 35 40 45 Gaussiana de 500m de altura Variación de e frente al ángulo (i del viento .1 0.8 -g 0.7 1 0.6 o E 0.5 <s 1 0.4 o í 0.3 •2 0.2; ^ 0.1' - - \ lili valor medio —— desv. tipica --^— ~ coef variac. (dividido por 10) | - ^ 00 05 10 15 20 25 30 Ángulo /? del viento 35 40 45 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 Gaussiana de 1500m de altura Variación de e frente al ángulo ¡3 del viento w d sv. tipies ;o if variac lor medie 00 05 10 15 20 25 30 Ángulo (5 del viento 35 40 45 Figura 4.8: Comparativa del parámetro e para las tres orografías. 122 Estimación de parámetros o c: <D > (D T3 O <D O fc <S fo «5 1- "í" (0 0.5 0.45 (14 0.35 0.3 0.25 0.2 0.15 0.1 0.05 Gaussiana de 50m de altura Variación de 7 frente al ángulo 0 del viento v« ormedio ái sv. tipies co f. variac 00 05 10 15 20 25 30 Ángulo ¡3 del viento 35 40 45 E 2 •55 E •s 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 O Gaussiana de 500m de altura Variación de 7 frente al ángulo /3 del viento 00 05 10 15 20 25 30 Ángulo ¡3 del viento Vi or medio sv. tipica co f. variac. / - E i Gaussiana de 1500m de altura Variación de 7 frente al ángulo (3 del viento u.o 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 - - - - - Vé d co or medio sv. tipies f. variac ^ II1 ~ —1_~ - - - - - - - 00 05 10 15 20 25 30 Ángulo ¡3 del viento 35 40 45 Figura 4.9: Comparativa del parámetro 7 para las tres orografías. Experimentos numéricos 123 E E fi Gaussiana de 50nn de altura Variación de 7/ frente al ángulo /3 del viento 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0^ - ; / / t \ \ 1 vaÉor medio d< sv. típica co f. variac. 1 2*. ^^^^^ - - - - - 00 05 10 15 20 25 30 Ángulo P del viento 35 40 45 p I 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 - 0.1 0.05 Gaussiana de 500m de altura Variación de 7/ frente al ángulo P del viento 00 05 10 15 20 25 30 Ángulo /3 del viento —I r valor medio • desv. lipica - coef. variac. - •> <D I Gaussiana de 1500m de altura Variación de 7/ frente al ángulo /? del viento 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 1 1 I Ve di co 1 or medio ^^— sv. tiplea • ~ f. variac. 1 - - - - - - 00 05 10 15 20 25 30 Ángulo /3 del viento 35 40 45 Figura 4.10: Comparativa del parámetro 7' para las tres orografías. 124 Estimación de parámetros 4.3.3. Influencia de la malla en la estimación de pcirámetros El siguiente experimento tiene como objeto estudiar la influencia del grado de refinamiento de la malla en la estimación de los parámetros del modelo. Para ello se generan dos nuevas mallas a partir de la original, ri, correspondiente a la gaussiana de 1500 m, realizando sobre ella dos refinamientos globales sucesivos. Tras el primer refinamiento global obtenemos una malla de 11787 nodos y 61160 tetraedros a la que llamaremos T2. Refinando globalmente T2 obtenemos T3, con 87865 nodos y 489280 tetraedros. La figura 4.11 muestra una vista de las tres mallas TI, T2 y T3. En sendos experimentos numéricos se estiman los cuatro parámetros a, e, 7 y 7' para un viento con ¡3 = 45°, que es el que se corresponde con la mejor estimación de parámetros para ri, utilizando las nuevas mallas T2 j T3. Los resultados obtenidos se pueden observar en la tabla 4.2 y de ellos se pueden sacar las siguientes conclusiones: • Los parámetros varían con el refinamiento y por lo tanto habría que aj listarlos en cada malla. • Las estimaciones realizadas en este problema sobre las mallas más finas mejoran ligeramente los valores de la función objetivo. Este resultado es coherente con el hecho de que la aproximación de la solución del problema depende de la discretización del dominio. Sin embargo, al aumentar el número de elementos de la malla también crece el tiempo de cómputo necesario para estimar los parámetros. En general, parece conveniente llegar a un compromiso entre grado de refinamiento y coste computacional en relación con la posible mejora de la función objetivo. • En nuestros experimentos hemos comprobado que tras estimar adecuadamente los parámetros, si se altera uno de ellos de manera arbitraria dándole a e 7 7' F. Objetivo Ti 0.3400 0.9980 0.4980 0.2023 0.0047 T2 1.1235 0.7422 0.4944 0.1505 0.0042 T3 1.7952 0.8079 0.3112 0.1549 0.0023 Tabla 4,2: Estimación de parámetros sobre las mallas TI, T2, T3. Experimentos numéricos 131 (a) Malla ra con P(ri) (b) Malla rg con P(T2) Figura 4.15: Distribución espacial de las zonas de mayor sensibilidad en la malla Capítulo 5 Simulación numérica en la Isla de Gran Canaria En este capítulo se desarrolla una aplicación numérica completa utilizando las técnicas descritas en los capítulos 2, 3 y 4. Se pretende calcular un campo de viento sobre la zona noroeste de la isla de Gran Canaria. Los datos de que disponemos son, una topografía digitalizada del terreno y un mapa meteorológico del viento previsto por el Instituto Nacional de Meteorología para la zona de Canarias. La aplicación numérica empieza por discretizar el dominio construyendo una malla de tetraedros utilizando las técnicas descritas en el capítulo 2; a continuación se mejora la calidad de esta malla inicial y se estiman los parámetros del modelo. Con estos parámetros se realiza un refinamiento local que adapta la malla en función de la solución numérica, y se realiza una nueva estimación de parámetros, esta vez sobre la malla refinada. Esta última fase refinamiento-estimación se repite hasta que la malla alcanza el número máximo de elementos con los que se desea trabajar. Finalmente se representa el campo de viento correspondiente a la mejor evaluación, utilizando otra estrategia adaptativa, y se comparan los resultados. 5.1. Estimación de los parámetros Planteamos el problema sobre un dominio de 16.5 x 9.5 x 7 km situado en la zona noroeste de la isla de Gran Canaria. Disponemos de una digitalización del terreno con una resolución de 25 x 25 m y un error máximo de 5 m en la cota. Con esta información se genera una malla tridimensional utilizando el programa MALLA (ver sección 2.9). Los parámetros empleados en la generación de la 134 Simulación numérica en la Isla de Gran Canaria malla son los siguientes: • El tamaño de los catetos de los triángulos de la malla grosera 2D sea aproxima a 3000 m. • Se realizan 8 pasos de refinamiento global sobre esta malla grosera 2D. • La nube de puntos se genera con la estrategia 1, fijando el grado de espaciado a = 2, y 8 capas incluyendo el terreno y el plano superior del dominio (n = 7). • El parámetro de desrefinamiento es £ = lOm. La malla resultante, que llamaremos r¿ y que se muestra en la figura 5.1, tiene 44886 nodos y 216043 tetraedros, con una calidad media g^ = 0.395. Tras diez iteraciones de suavizado la calidad media pasa a ser g„ = 0.749. En la tabla 5.1 pueden verse las calidades mínimas, medias y máximas obtenidas en cada paso del proceso de suavizado. En la figura 5.3 se muestra cómo mejoran las curvas de las calidades de todos los tetraedros con el suavizado; llamaremos TQ a la malla obtenida tras el último paso de suavizado realizado sobre la malla inicial r¿ (ver figura 5.2). En cuanto al coste computacional, la generación de la malla TQ se realizó en 47 segundos y el suavizado en 10.5 minutos, ejecutándose sobre un ordenador con dos procesadores Intel Xeon a 2.1 GHz con 4 Gb de memoria RAM. Como no disponemos de datos de estaciones meteorológicas de la zona de estudio optamos por obtenerlos de los mapas de previsión que el Instituto Nacional de Meteorología publica en sus páginas web (http://www.inm.es/puertos/mapas. html). En la figura 5.4 puede verse el que se ha utilizado en esta aplicación, correspondiente a la previsión realizada el día 27/04/04 a las 00 horas UTC, con horizonte de previsión a 48 horas, esto es, para el 29/04/04 a la misma hora. Destacamos que esta configuración de vientos del NO es atípica en esta región, siendo la más frecuente la de los alisios, con vientos del NE. La escala de colores indica la intensidad del viento y la dirección viene dada por las flechas. Con el programa de manipulación de imágenes GIMP (http://www.gimp.org) se aproximaron ambos datos. Decidimos utilizar tres medidas de viento repartidas al norte de la isla y una al oeste, todas sobre el mar (véase la figura 5.5). Los valores de viento asignados a las estaciones ficticias son los correspondientes a sus localizaciones en el mapa meteorológico y pueden verse en la tabla 5.2. La intensidad y dirección del Estimación de los parámetros 135 viento geostrófico se aproximó aplicando la ley de resistencia geostrófica, según la expresiones (3.22) y (3.23) respectivamente, resultando un vector de componentes (0.47, -3.19, 0) en m/s. Una vez fijados los valores de viento, se realizó una primera estimación de los parámetros del modelo con algoritmos genéticos, resultando los siguientes valores: a = 5.3727, e — 1, 7 = 0.1501, 7' = 0.15. La mejor evaluación de la función objetivo fue 0.0254. A continuación se reaUzaron tres iteraciones de cálculo del viento aplicando refinamiento local de la malla, con el indicador de error de la expresión (3.89) y un parámetro de refinamiento 9 = 0.45, tras las que se obtiene la malla ri con 66248 nodos y 335865 tetraedros. Iteraciones de suavizado Malla inicial 1 2 3 4 5 6 7 8 9 10 Mínima 0.079 0.099 0.155 0.185 0.195 0.200 0.202 0.203 0.204 0.204 0.204 Calidades Media 0.395 0.533 0.616 0.665 0.696 0.716 0.729 0.738 0.743 0.747 0.749 Máxima 0.987 0.994 0.998 0.999 0.998 0.999 0.998 0.999 0.998 0.998 0.998 Tabla 5.1: Evolución del suavizado de la malla inicial r¿ definida para la aplicación numérica sobre Gran Canaria. Con esta nueva malla ri como soporte se realizó una nueva estimación de los parámetros del modelo, resultando a = 4.9344, e ^ 0.9999, 7 = 0.15 y 7' = 0.1502; el valor de la función objetivo fue 0.0315. Estos nuevos parámetros se utilizaron para calcular nuevamente el campo de velocidades de viento y realizar un paso de refinamiento local sobre Ti. El parámetro de refinamiento en este caso fue 9 = 0.5, resultando la malla T2 con 103744 nodos y 545420 tetraedros. Finalmente se realizó otra estimación de parámetros sobre r2, en la que resultó Oí = 5.059, e = 1, 7 = 0.1501 y 7' = 0.1501, con un valor de la función objetivo de 0.032. Se calcularon después las sensibiUdades (ver ecuación (4.3)) correspondientes 136 Simulación numérica en la Isla de Gran Canaria Figura 5.1: Vista de la malla TQ antes del proceso de suavizado. a las combinaciones de parámetros de la tabla 5.3. Puede comprobarse que los valores de las sensibilidades son pequeños, ya que los parámetros estimados son bastante parecidos. En las gráficas 5.6, 5.7, 5.8, 5.9, 5.10y5.11 pueden observarse algunos detalles de la distribución espacial de las zonas de mayor sensibilidad en cada uno de los casos. Como se ve, las zonas donde el modelo es más sensible a la variación de los parámetros coinciden en todos los casos. Esto viene a ratificar que el estudio de la sensibilidad puede ser útil para el diseño de una red de medida que permita la evaluación del potencial eólico de una zona, en el sentido de que colocando más estaciones en las zonas de mayor sensibilidad se espera que se reduzca la sensibilidad del modelo. Estimación de los parámetros 137 Figura 5.2: Vista de la malla TQ obtenida tras el suavizado de TL 138 Simulación numérica en la Isla de Gran Canaria en Calidad Inicial Calidad tras un paso de suavizado Calidad tras seis pasos de suavizado Calidad tras diez pasos de suavizado 50000 100000 e 150000 200000 Figura 5.3: Gráficas de calidad de los elementos de la malla durante el proceso de suavizado. Estación El E2 E3 E4 v\ 3.75 3.75 3.75 3.00 direc. 155 155 155 178 X 433926 439460 443333 430053 y 3116868 3116868 3116868 3113068 Tabla 5.2: Valores de viento asignados a las estaciones y localización de las mismas en coordenadas UTM. Malla Ti Ti Ti T2 T2 T2 Parámetros Pin) P{r2) Piro) Pin) Pin) Piro) Pin) Pir2) F. Objetivo 0.0255 0.0256 0.0255 0.0316 0.0315 0.0315 0.0320 0.0320 0.0320 Max. Sensib. 2.4304 1.6919 0.7188 0.1786 0.9055 0.4151 Tetraedro 55475 55475 166981 154743 349208 349208 Tabla 5.3: Valores de las sensibilidades. P(TJ) es el conjunto de los cuatro parámetros estimados sobre la malla r,-. Estimación de los parámetros 139 SAUDA DELMOOELO OELDIA Z7;D4A)4 A LA5 DO UTC de INM / PE Canpa da Viamo (m^) afaa 00 QMT<ialdia29,'04<t>4.) Horizomi Pradiocian = 46 ham MMU» devieni» eo supesfic»? (codo de ei*j» ' C'iiícción d? pioivjíocort ^ Figura 5.4: Mapa de viento del I.N.M. usado para asignar valores de viento a las estaciones. E4 ¡f. Jtri-^: ', I'J Figura 5.5: Localización de las cuatro estaciones sobre el dominio. 140 Simulación numérica en la Isla de Gran Canaria Figura 5.6: Sensibilidad de TQ con -P(ri). Zona con valores de sensibilidad por encima de 0.024. Figura 5.7: Sensibilidad de TQ con P{T2). Zona con valores de sensibilidad por encima de 0.016. Estimación de los parámetros 141 Figura 5.8: Sensibilidad de TI con P{TQ). Zona con valores de sensibilidad por encima de 0.016. Figura 5.9: Sensibilidad de TX con P{T'¿). Zona con valores de sensibilidad por encima de 0.005. w..: .«^'