Adaptive scheme for three-dimensional convection – diffusion problems

Resumen

Se presenta un esquema adaptativo para problemas tridimensionales de convección-difusión discretizados mediante el Método de Elementos Finitos. El esquema adaptativo incluye una estrategia de remallado basada en la imposición de un volumen máximo de los elementos a partir de una malla de referencia. El remallado permite aumentar o disminuir drásticamente el tamaño de los elementos en un solo paso, de forma automática. Se mantiene una calidad de la malla adecuada; y, como consecuencia, el número de iteracciones necesarias para solucionar el sistema de ecuaciones lineales mediante algoritmos iterativos se mantiene constante. Se presentan 2 ejemplos de características muy distintas que permiten analizar la propuesta para un amplio rango de situaciones. Uno es una extensión tridimensional del problema de Smolarkiewicz y el otro es una versión simplificada del problema de transporte de contaminación emitida por un emisor puntual. Los resultados muestran la flexibilidad de la propuesta. Es posible encontrar un valor óptimo de frecuencia de remallado, desde el punto de vista de coste computacional y precisión de los resultados, para ambas tipologías extremas de problemas.

Abstract

We present an adaptive scheme for three-dimensional convection-diffusion problems discretized by the Finite Element Method. The adaptive scheme is based on a remeshing strategy that applies a maximum volume constraint to the elements of a reference mesh. The remeshing can increase or decrease drastically the size of the elements in a single step automatically. With this strategy, the mesh quality does not deteriorate; as a consequence, the number of iterations required to solve the system of linear equations using iterative algorithms is kept constant. Two examples of very different characteristics are presented in order to analyze the proposal for a wide range of situations. The first is a three-dimensional extension of the Smolarkiewicz problem and the second is a simplified version of a point source pollutant transport problem. The results show the flexibility of the proposal. An optimal remeshing frequency, from a computational cost and accuracy of the results point of view, can be defined for both kinds of problems.

Palabras clave

Ecuaciones de transporte ; Estabilización por mínimos cuadrados ; Coste computacional ; Precisión ; Indicador de error ; Transporte de contaminantes

Keywords

Transport equations ; Least-square stabilization ; Computational cost ; Accuracy ; Error indicator ; Pollutant transport

1. Introducción

Los esquemas adaptativos se emplean para resolver numéricamente ecuaciones en derivadas parciales con discretizaciones diseñadas para minimizar o controlar el error de la solución. Mediante la adaptatividad de la discretización espacial se puede ajustar la precisión de la solución numérica según se requiera, optimizando las necesidades de memoria motivadas por el tamaño del sistema de ecuaciones a resolver. Los esquemas adaptativos requieren un indicador o estimador del error que cuantifique la necesidad de incrementar o disminuir la densidad de la discretización, así como una estrategia de remallado eficiente [1] .

Uno de los esquemas de remallado más habituales es la bisección de elementos [2] , [3]  and [4] , que consiste en insertar nodos nuevos en las aristas, las caras o el interior de los elementos que tengan una medida del error superior a una tolerancia. Este algoritmo se aplica de forma recursiva, dividiendo sucesivamente los elementos en un número prefijado de elementos en su interior. La calidad de los nuevos elementos creados es ligeramente inferior a la de los de la malla anterior, pudiendo aparecer elementos de mala calidad en la frontera de las zonas que se deben refinar y las que no [5]  and [6] , especialmente en geometrías complejas o donde la variación fijada del tamaño de los elementos es grande. La disminución de la calidad de la malla provoca el mal condicionamiento de los sistemas lineales que aparecen al resolver el problema, cosa que a su vez hace que sean necesarias más iteraciones para converger cuando se resuelven estos sistemas mediante esquemas iterativos [7]  and [8] . La pérdida de calidad de la malla puede solventarse recolocando los nodos de la discretización mediante técnicas de suavizado [9]  and [10] . Existen métodos que permiten mejorar la calidad de la malla manteniendo el tamaño de los elementos [11] . Sin embargo, al emplear suavizados se pierde una de las ventajas principales de los métodos de bisección que es la facilidad de proyección de la solución entre mallas sucesivas.

Otra técnica de remallado, típicamente empleada para generar elementos anisotrópicos, consiste en generar una nueva malla basándose en una métrica definida a partir del estimador o indicador del error [12] . Se genera la malla de forma que el error presente una distribución uniforme, tomando como criterio de medida la métrica definida. Esta estrategia ha sido utilizada en problemas de convección-difusión para generar mallas de cuadriláteros en problemas bidimensionales [13]  and [14] o tetraedros en tridimensionales [15] .

En este trabajo se propone utilizar un esquema de remallado de este segundo tipo, de forma que se pueda aumentar y disminuir drásticamente, en un solo remallado, la densidad de elementos donde se considere necesario, manteniendo mallas de cálculo de buena calidad. El esquema adaptativo se utiliza para resolver problemas tridimensionales de convección-difusión lineales, discretizados mediante el Método de los Elementos Finitos. En la siguiente sección se introduce el problema matemático, la estrategia de resolución numérica y el esquema adaptativo. Posteriormente, se describe la estrategia de remallado y su aplicación a 2 ejemplos: una extensión tridimensional del problema de Smolarkiewicz [16] y una versión simplificada del problema de transporte de contaminación emitida por un emisor puntual [4] . Los 2 ejemplos tienen soluciones con características muy diferentes. En el primer caso la solución presenta grandes y crecientes variaciones en una zona acotada del dominio. En el segundo la solución tiene variaciones más suaves, pero valores de interés varios órdenes de magnitud inferiores a los impuestos en las condiciones de contorno y se traslada a lo largo del dominio. Se finaliza destacando las principales conclusiones del estudio.

2. Esquema adaptativo

Se considera el siguiente problema genérico de convección-difusión:

${\begin{array}{ll}{\partial }_{t}u+{\boldsymbol {\mbox{a}}}\cdot \nabla u-\nabla \cdot ({\boldsymbol {\mbox{D}}}\cdot \nabla u)=s&en\quad \Omega \times (0,T)\\u({\boldsymbol {\mbox{x}}},0)=u_{0}({\boldsymbol {\mbox{x}}})&en\quad \Omega \\Mu=0&en\quad \partial \Omega \times (0,T)\end{array}}$

( 1)

donde a es la velocidad advectiva, D es el tensor de difusiones, s   el témino fuente y ${\textstyle M}$ las condiciones de contorno, que deben cumplir las condiciones de regularidad necesarias, véase Verfurth [17] .

El problema se integra en el tiempo mediante el esquema de Crank-Nicolson. La discretización espacial se realiza con Elementos Finitos estabilizados mediante la formulación de Mínimos Cuadrados [18] . El sistema lineal de ecuaciones se resuelve mediante Gradientes Conjugados Precondicionados con una factorización incompleta de Cholesky [19]  and [20] .

En el algoritmo 1 se presenta la propuesta, siendo u_n el vector que contiene la solución en los nodos en el instante t_n  = n Δt con Δt  = mδt . El esquema adaptativo se aplica únicamente a la discretización espacial, manteniendo constante el paso de tiempo δt . El valor de m es el número de pasos de tiempo que se calculan con una misma malla de cálculo. Se remalla cada Δt .

Habitualmente los esquemas recalculan la solución de cada bloque de pasos de tiempo, hasta que el indicador o el estimador de error de la solución o la sucesión de mallas generadas cumplen un determinado criterio de convergencia [21]  and [22] . Los esquemas que no recalculan la solución tienden a acumular más error en la solución [15] , en cambio, presentan un coste computacional más reducido. En el esquema propuesto se propone no recalcularla, a excepción del primer bloque de pasos de tiempo. Esta estrategia es adecuada si es la primera malla de cálculo la que determina los resultados posteriores o si la variación espacial de la solución entre remallados no es muy elevada. En otras situaciones puede ser conveniente incorporar el recálculo y la verificación de la convergencia para cada Δt .

Algoritmo 1.

Esquema adaptativo

• while No convergencia do
• Calcular [0, Δt ] (m pasos de δt )
• Remallar
• end while
• Interpolar
• Guardar u₁
• i  = 1
• whilei Δt  < Tdo
• Calcular [i Δt , (i  + 1)Δt ] (m pasos de δt )
• Remallar
• Interpolar
• Guardar u_{(i +1)}
• i  = i  + 1
• end while

La malla de cálculo inicial se utiliza como malla de referencia a lo largo de todo el problema. Es una discretización muy poco densa, comparada con las de cálculo posterior, establecida fundamentalmente con criterios de representación geométrica del dominio. Puede incluir información sobre la solución y su evolución si se dispone de ella.

Tras la convergencia del primer bloque de pasos de tiempo de cálculo, en los siguientes se calcula la solución, se aproxima el error, se genera una nueva malla y se interpola la solución a la nueva malla.

El error se aproxima mediante un indicador del error. En la literatura se encuentran distintos indicadores de error para problemas de convección-difusión [23] ; por ejemplo, el valor de la propia solución, su gradiente [4] o su curvatura [24] . La elección es más crítica en esquemas r -adaptativos que en esquemas h -adaptativos, como el planteado en esta propuesta [25] . En los primeros, el indicador de error actúa como criterio en un proceso de optimización que busca el mínimo de error para un número de grados de libertad fijado. Sin embargo, en los segundos se incorporan nuevos grados de libertad con el fin de que la precisión de los resultados aumente. Los primeros están más condicionados, por lo que requieren que el indicador sea más preciso.

Los indicadores definidos a partir del gradiente o de la máxima diferencia de la solución en los elementos son habituales; por ejemplo:

${\eta }_{e}=\Vert \nabla u{\Vert }_{{\Omega }^{e}}$

( 2)

siendo ${\textstyle \Vert \cdot {\Vert }_{{\Omega }^{e}}}$ la norma de · en los elementos Ω^e . Este tipo de indicadores tienden a producir mallas refinadas donde el gradiente de la solución aproximada es mayor. Esto es útil en diversos problemas, por ejemplo en los que se quiere localizar una capa límite [26] , pero no para otros. En problemas de calidad de aire con emisores puntuales, la solución, la concentración de las especies consideradas, presenta valores varios órdenes de magnitud inferiores a los de emisión. Tanto ciertos resultados de interés, como los valores de inmisión, como las oscilaciones introducidas por la resolución numérica se encuentran en zonas en las que la solución presenta valores relativos muy bajos. Interesa determinar estas zonas para minimizar la presencia de errores y oscilaciones. Un indicador de error adecuado para estas situaciones es:

${\eta }_{e}={\begin{array}{ll}\Vert \nabla log(u){\Vert }_{{\Omega }^{e}}=\Vert \nabla u/u{\Vert }_{{\Omega }^{e}}&si\quad u>{tol}_{u}\\0&otros\end{array}}$

( 3)

La tolerancia tol_u fija el límite inferior de los valores de la solución sobre los que ya no se refina.

3. Remallado e interpolación

El remallado se implementa a partir de la malla de referencia, donde se impone, en cada elemento, un volumen máximo que deben cumplir los nuevos elementos que se definan en su interior. Este volumen máximo depende del indicador de error y del volumen de los elementos de la malla utilizada previamente. Cada elemento de la malla de referencia se considera una región de remallado.

Sea ${\textstyle V_{r}^{n+1}}$ el volumen que se impone a la región r para la malla de cálculo n  + 1. Si el indicador de error de algún elemento de la región r es superior a una tolerancia, la expresión del volumen máximo que pueden tener los nuevos elementos de la región se define como:

$V_{r}^{n+1}={\underset {e\in S_{r}}{min}}\left({\frac {V_{e}}{1+\alpha {\eta }_{e}}}\right)$

( 4)

donde S_r es el conjundo de elementos de la malla que pertenece a la región r cuyo indicador es superior a una tolerancia dada, y V_e y η_e son, respectivamente, el volumen y el indicador de error del elemento e de la malla de cálculo. El parámetro α  > 0 controla el grado de refinamiento que se quiere obtener. Mediante esta estrategia de refinamiento se consigue que el tamaño de los elementos sea menor en las regiones donde el error es mayor, y que el error esté distribuido uniformemente en los elementos de la malla [27] . Se impone un volumen mínimo de los elementos, para no refinar en exceso y obtener un número de nodos mayor del deseado.

Si los errores de todos los elementos de la región r son inferiores a la tolerancia fijada, se impone un nuevo tamaño de elemento, mayor al impuesto en el remallado anterior:

$V_{r}^{n+1}=min\left({\underset {e\in S_{r}}{min}}\left({\frac {V_{e}}{1+\alpha {\eta }_{e}}}\right),\beta V_{r}^{n}\right)$

( 5)

donde β  > 1 es un parámetro que controla el desrefinamiento.

La estrategia de remallado propuesta permite aumentar y reducir fuertemente el tamaño de los elementos de una malla a la siguiente. Los parámetros α y β permiten modular ambos efectos. Posteriormente, se muestra su influencia tanto en las mallas de cálculo como en la solución del problema.

Al imponer un tamaño de elemento constante para cada región, la capacidad de adaptación de la malla a la solución está limitada por la definición de la malla de referencia. Esto puede limitar la efectividad de la estrategia adaptativa. En esos casos puede ser útil modificar la malla de referencia para ajustarla a las particularidades del problema, ya sea desde el inicio, manteniéndola constante a lo largo del problema, o actualizándola a lo largo de la evolución del mismo.

El remallado se ha implementado mediante el mallador Tetgen[28]  and [29] . Se ha asegurado que las mallas generadas son de calidad adecuada imponiendo no solo las restricciones de volumen definidas, sino también requisitos sobre la calidad de los elementos generados.

Una vez definida la nueva malla es necesario interpolar la solución. Entre 2 mallas consecutivas solo se puede asegurar que los nodos de la malla de referencia son comunes a ambas. El esquema de interpolación implementado en este trabajo es lineal.

4. Ejemplos

El esquema adaptativo propuesto se aplica a 2 ejemplos. El primero es una versión tridimensional del problema de Smolarkiewicz [16] . La formulación original del problema es bidimensional y se ha utilizado para evaluar métodos de resolución de problemas convectivos transitorios. Aquí, la formulación del problema se extiende a un dominio tridimensional. Es un ejemplo con solución analítica, con fuertes y crecientes variaciones de la solución pero que se desarrolla en una zona acotada y constante del dominio.

El segundo ejemplo es un problema simplificado de emisión y transporte de contaminantes en la atmósfera [4] . La difusión vertical y la velocidad dependen de la altura sobre el terreno. Se incluye la discretización del emisor puntual en la malla de referencia. La emisión se activa solo en un intervalo de tiempo. A diferencia del primer ejemplo, en este la solución presenta variaciones más suaves, pero se traslada por el dominio. Los valores de interés de la solución son varios órdenes de magnitud inferiores a los impuestos en la emisión.

En ambos casos se presenta y analiza la solución para un intervalo de remallado m fijo, y se estudia su influencia en los resultados. En el segundo caso se estudia también la influencia de las constantes de remallado α y β .

4.1. Problema de Smolarkiewicz

El problema se define en el dominio Ω = [0, L ] × [0, L ] × [0, L /4]. Se fija D  = 0, s  = 0 y el campo de velocidades:

${\boldsymbol {\mbox{a}}}=(Aksin(kx)sin(ky),Akcos(kx)cos(ky),0)$

( 6)

con ${\textstyle k={\frac {4\pi }{L}}}$ , A  = 8 y L  = 100. La condición inicial es, para todo z , una pirámide circular de radio 15 y altura máxima uno, centrada en el dominio x – y . En la figura 1 se representa un esquema bidimensional de las trayectorias y de la condición inicial. Las condiciones de contorno son Dirichlet nulas en los contornos de entrada y Neumann nulas en las salidas y en ambos extremos z  = 0 y z  = L /4.

Figura 1. Esquema bidimensional del problema de Smolarkievicz: En trazo continuo las líneas de corriente y en trazo discontinuo las isolíneas de la condición inicial.

La solución del problema se puede obtener teniendo en cuenta que el flujo es únicamente convectivo y, como consecuencia, la concentración se conserva a lo largo de las trayectorias. Las trayectorias se encuentran integrando analíticamente el campo de velocidades [30] . La solución no es continua entre las distintas celdas. El período de las trayectorias cerradas es proporcional a la distancia al centro de los vórtices, lo que provoca fuertes y crecientes gradientes en la solución. El tiempo T  = 2.637, 6 coincide con el momento en que la solución presenta 200 máximos relativos en y  = L /2. El tiempo final de cálculo es T   /50 y el paso de tiempo ${\textstyle {\frac {T/50}{100}}}$ .

Se ha utilizado el indicador de error de la ecuación (2) , con un valor umbral para activar el remallado de los elementos igual a 0.05. Se han fijado las constantes de remallado α  = 30 y β  = 2.5. En primer lugar se ha aplicado el esquema adaptativo considerando un único bloque de pasos de tiempo para el remallado, es decir con m  = 100 y Δt  = T /50. Por tanto, repitiendo el cálculo completo, hasta convergencia, con sucesivas mallas, refinadas a partir de la malla inicial de referencia y adaptadas a la solución en el tiempo final de integración. Posteriormente, se ha analizado la influencia de variar la frecuencia de remallado, variando m .

4.1.1. Resultados

Las figuras 2 a y 2 c muestran la solución para ${\textstyle t={\frac {T}{50}}}$ en z  = L /4 y z  = L /8 respectivamente. Los resultados obtenidos son comparables, cuando no mejores, a los presentados en la literatura para el problema en 2 dimensiones, por ejemplo en Hundsdorfer y Spee [31] , donde se soluciona el problema usando splitting de segundo orden y volúmenes finitos de tercer orden, y en en Sarma [32] y Drange y Bleck [33] , donde se usa un esquema upstream . La solución obtenida es también mejor que la presentada en Bott [34] , donde se utiliza una discretización espacial demasiado gruesa. Las diferencias de la solución según el eje z , entre las figuras 2 a y 2 c, son reducidas. Son debidas a las oscilaciones que aparecen en la solución numérica y del mismo orden a las que aparecen en las otras dimensiones espaciales.

Figura 2. Solución y mallas de cálculo y referencia (trazo fino y trazo grueso) en el contorno superior z  = L /4, (a) y (b), y para z  = L /8 en (c) y (d).

Las figuras 2 b y 2 d muestran las mallas de cálculo y referencia (trazos fino y grueso respectivamente) para los 2 valores de z . La primera es una vista del contorno superior del dominio y la segunda una sección de la malla tridimensional. El esquema adaptativo genera mallas con un tamaño de elemento reducido en las zonas donde el gradiente de la solución es elevado. Se puede apreciar como todos los elementos de la malla de cálculo que están en el interior de un mismo elemento de la malla de referencia tienen un volumen parecido. Las excepciones son debidas a las condiciones de calidad de la malla impuestas, que suavizan las transiciones en los tamaños de elemento.

Por último se destaca que la calidad de la malla se mantiene adecuada a lo largo de todo el problema y no afecta la resolución del mismo. El sistema lineal de ecuaciones se resuelve en menos de 6 iteraciones en todos los pasos de tiempo.

4.1.2. Frecuencia de remallado

Una vez resuelto el problema, se ha analizado la influencia de la frecuencia de remallado, variando m . La figura 3 a muestra el coste computacional relativo al de m  = 100. Puede observarse que al reducir m el coste computacional baja hasta m  = 20, a partir del cual el coste aumenta. En principio, si se remalla con más frecuencia es de esperar que el coste computacional sea mayor, debido a que se realizan más veces los procesos de cálculo del indicador, generación de malla, interpolación de la solución y factorización de las matrices del sistema lineal de ecuaciones. Pero, por otro lado, al remallar antes, la malla se adapta antes a la solución y son necesarias menos iteraciones para converger en el primer bloque de pasos de tiempo de cálculo. La combinación de ambos efectos explica el resultado obtenido.

Figura 3. Evolución del coste computacional (a), la norma L2 del error en ${\textstyle t={\frac {T}{50}}}$ (b), el valor máximo de la solución en el dominio (c), y el número de nodos de la última malla de cálculo (d), en función del período de remallado m .

Las figuras 3 b y 3 c muestran, respectivamente, los valores de la norma L² del error de la solución y el valor máximo de la misma, ambos en función del período de remallado. El valor máximo de la solución analítica es uno, por tanto la figura 3 c es, como la figura 3 b, una medida del error de la solución numérica. Se observa valores bajos de m el error es mucho mayor valores de m igual a 20 y mayores. Este hecho es debido al error que se introduce en la interpolación. La interpolación tiende a suavizar la solución y a reducir los valores máximos. La figura 4 muestra el error introducido al interpolar la solución entre 2 mallas consecutivas. Con el suavizado disminuye el indicador de error, ecuación (2) . Además, como los indicadores propuestos no detectan los errores de interpolación, con periodos de remallado pequeños la influencia del error de interpolación es mayor y las mallas generadas son menos densas, como muestra la figura 3 d.

Figura 4. Solución en una malla de cálculo y resultado de la interpolación en una nueva malla.

A la vista de los resultados se puede concluir que el esquema adaptativo propuesto presenta una frecuencia de remallado óptima para problemas en que la solución se concentra en una parte del dominio y presenta fuertes variaciones. Remallar menos veces que el número óptimo es computacionalmente más caro, sin ofrecer mejoras significativas en los resultados. Remallar más veces también es computacionalmente desfavorable, introduciendo, además, errores de interpolación que no pueden ser corregidos mediante el esquema adaptativo.

4.2. Ejemplo meteorológico simplificado

El segundo ejemplo consiste en una emisión puntual en el interior de un dominio Ω = [0, 96000] × [0, 48000] × [0, 3000] m³ con valores no constantes de velocidad y difusión vertical [4] . El emisor puntual está situado en las coordenadas (6.000, 24.000, 250) m. Se ha discretizado como una esfera de radio 5 m. La difusión es diagonal, con valores ${\textstyle {\mbox{D}}_{xx}={\mbox{D}}_{yy}=}$$50{\frac {m^{2}}{s}}$ y D_zz función de la altura. La velocidad es paralela al eje x y de módulo variable también en altura. La figura 5 presenta las variaciones de ambos para condiciones atmosféricas neutras [35] . La condición inicial es nula en todo el dominio. El término fuente, s , también es nulo. Las condiciones de contorno son:

${\begin{array}{cc}u=0&{\mbox{en}}{\mbox{ }}{\Gamma }_{ext}^{in}\times (0,{\mbox{T}})\\{\boldsymbol {n}}\cdot {\boldsymbol {D}}\cdot \nabla u={\mbox{g}}({\mbox{t}})&{\mbox{en}}{\mbox{ }}{\Gamma }_{int}\times (0,{\mbox{T}})\\{\boldsymbol {n}}\cdot {\boldsymbol {D}}\cdot \nabla u=0{\frac {\mbox{g}}{{\mbox{m}}^{2}{\mbox{s}}}}&{\mbox{en}}{\mbox{ }}{\Gamma }^{out}\times (0,{\mbox{T}})\end{array}}$

con Γ^out el contorno con flujo saliente, ${\textstyle {\Gamma }_{ext}^{in}}$ el contorno exterior con flujo entrante y Γ_int . el contorno del emisor. La emisión es dependiente del tiempo según:

$g(t)={\begin{array}{ll}100\quad {\frac {g}{m^{2}s}}&si\quad 0<t<1.800\quad s\\0\quad {\frac {g}{m^{2}s}}&si\quad t\geq 1.800\quad s\end{array}}$

Se fija δt  = 1 s, y se utiliza el indicador de la ecuación (3) con tol_u  = 10⁻⁶${\textstyle {\frac {g}{m^{3}}}}$ .

Figura 5. Perfil del módulo de la velocidad y de la difusión vertical en función de la altura.

Primero se aplica el esquema adaptativo con m  = 180 (10 remallados durante la emisión) y con las constantes de remallado fijas, después se varían las constantes de remallado α y β manteniendo m  = 180, y posteriormente se analiza la influencia de m .

4.2.1. Resultados

La figura 6 muestra la solución obtenida para distintos instantes de tiempo. La escala de las concentraciones es logarítmica. La emisión forma una pluma que se desarrolla en la dirección de la velocidad. Se aprecia el efecto de la variación de la difusión en altura: en las zonas donde es menor, la variación vertical de la solución es también menor. Cuando finaliza la emisión, para t  > 1.800 s, el penacho avanza y tiende a diluirse en el medio. Con el esquema adaptativo propuesto la solución no presenta oscilaciones apreciables para valores mayores a 10⁻⁵ .

Figura 6. Corte de la solución del problema de emisón en la dirección del transporte para distintos instantes de tiempo.

La figura 7 a presenta la malla de referencia y las figuras Figura 7 b,d, 7 f y 7 h las mallas utilizadas en t  = 900, t  = 1.800, t  = 2.700 y t  = 3.600 s respectivamente. La malla de referencia está adaptada a la geometría del problema, que incluye explícitamente el emisor. Las mallas de cálculo tienen entre 6 y 13 veces el número de elementos de la malla de referencia, y se adaptan a la evolución de la solución. La densidad de elementos es alta en las zonas donde la solución toma valores altos, y es suficientemente fina en los valores bajos, evitando que aparezcan oscilaciones. Para los valores menores al umbral de remallado, inferiores a 10⁻⁶ , la malla es muy poco densa y básicamente igual a la de referencia.

Figura 7. Malla de referencia (a), y mallas de cálculo para distintos tiempos y distintos valores de las constantes del esquema adaptativo, de (b) a (i).

4.2.2. Constantes del esquema de remallado

La influencia de la constante que determina el grado de refinamiento, α , se ha analizado en los primeros 1.800 s del problema, en los que domina el refinamiento debido a que se forma la pluma; se han mantenido constantes m y β . La de la constante de desrefinamiento, β , se ha estudiado en los 1.800 s posteriores, en los que cesa la emisión, el penacho se diluye en el medio y se espera un desrefinamiento en la zona entre el penacho y el emisor.

Las diferencias en la malla de cálculo según el valor de α se aprecian comparando las figuras 7 b y 7 d, calculadas con α  = 30, con las figuras 7 c y 7 e, calculadas con α  = 10. Valores mayores de α conducen a mallas con mayor número de elementos, entre un 20 y un 30% más para un factor 3 en α . La solución para ambos casos es muy parecida. No se observan diferencias en los valores altos de la solución, pero en los valores bajos aparecen menos oscilaciones en la malla más fina.

A partir de la solución obtenida con α  = 30 y β  = 2.5, se han calculado los 1.800 s posteriores con 2 valores de β . Esta constante puede considerarse una velocidad de desrefinamiento, ya que el nuevo volumen que se impone es el producto de β con el volumen impuesto en el remallado anterior. Las figuras 7 f y 7 h se han calculado con β  = 2.5 y las figuras 7 g y 7 i con β  = 10. Ambas soluciones son parecidas. El desrefinamiento crece con β , sobre todo en la zona situada entre el penacho y el emisor. No se modifica la malla de referencia, por ello en ambos casos el tamaño de los elementos cercanos al emisor se mantiene reducido. Valores mayores de β generan mallas con menos elementos, un 10% menos para valores de β 4 veces mayores, en este ejemplo. Ambas constantes producen el efecto esperado. A continuación se cuantifica más detalladamente su influencia.

En la figura 8 se muestra la norma L² del error, el número de nodos y el tiempo de cálculo para α de 1 a 10² , con β  = 2.5 y para el cálculo hasta t  = 1.800 s. El valor de α es relativo (multiplica al valor del indicador de error utilizado, véase la ecuación (4) ). Se puede observar que con valores elevados de α disminuye el error, tomando como referencia la solución obtenida con una malla muy densa. El número de nodos de la malla de cálculo aumenta con el valor de α , y como consecuencia también lo hace el coste computacional. En media, un aumento de α de 2 órdenes de magnitud duplica el tamaño del problema y multiplica por 10 el coste computacional, pero reduce a la mitad el error.

Figura 8. Evolución de la norma L2 del error en t  = 1.800 s (a), el número de nodos en t  = 1.800 s (b), y el coste computacional (c), en función del parámetro α .

Por otro lado, la constante β está directamente relacionada con la frecuencia de remallado; véase la ecuación (5) : β multiplica al tamaño de elemento impuesto en el remallado anterior; por tanto, en las zonas en que realizan repetidos desrefinamientos, el efecto de β aplica de forma repetida. En la figura 9 se presentan los resultados para los últimos 1.800 s, variando β entre 1.5 y 10² , con α constante y m  = 180 (10 desrefinamientos como máximo en determinadas zonas del dominio). Se pueden apreciar 2 comportamientos diferenciados: para β  < 10 el error disminuye con β y el número de nodos y el coste computacional aumentan, en cambio, para valores de β  > 10 no se aprecia una diferencia significativa en los resultados. La malla de referencia limita el desrefinamiento para valores elevados de β . Por otro lado, en el límite, β  = 1, no se realiza desrefinamiento, y en este caso se obtiene la mitad de error al doble de coste, en relación a las mallas más desrefinadas. En general, la reducción de error será como máximo la diferencia entre el de la malla de referencia y el de la malla de cálculo que se desea desrefinar.

Figura 9. Evolución de la norma L2 del error en t  = 3.600 s (a), el número de nodos en t  = 3.600 s (b), y el coste computacional de los últimos 1.800 pasos de tiempo (c), en función el parámetro β .

4.2.3. Frecuencia de remallado

Por último, se analiza la influencia de m . Se considera el cálculo hasta el tiempo final t  = 1.800 s, período de tiempo con la emisión activa. La figura 10 muestra la solución y la malla de cálculo utilizada en t  = 1.800 s para m igual a 90, 180, 300, 600, 900 y 1.800 (de 20 a 1 remallados). El número de elementos para todos los casos es parecido. La diferencia es menor del 5% y se concentra en la zona de avance del frente, para valores de m elevados (300, 600 y 900). El tamaño y el número de elementos en las zonas de desarrollo inicial de la pluma y lejanas al avance de la misma coinciden en los distintos casos.

Figura 10. Solución y malla de cálculo en t  = 1.800 s para distintos valores de m .

Los resultados obtenidos con m  = 90, 180 y 1.800 son parecidos: tanto la solución como las mallas de cálculo utilizadas en t  = 1.800 s. La solución con más remallados intermedios presenta un frente algo más suave, con la posición de la isosuperficie 10⁻⁵ ligeramente más avanzada que el caso recalculado completamente. En cambio, para valores de m intermedios, m  = 300 a 900, el resultado es cualitativamente diferente. El frente no se determina con precisión, la solución se difunde en los grandes elementos presentes en la zona sin refinar, y las oscilaciones son visibles en el rango de valores de interés. Esto se debe a que la zona con valores significativos de la solución varía mucho entre entre 2 remallados consecutivos. Por ejemplo, en el caso m  = 900, figuras 10 i y 10 j, el esquema adaptativo itera el remallado y el cálculo hasta t  = 900 δt . A continuación, se calculan los 900 pasos de tiempo siguientes con la malla fija. En estas condiciones el avance del frente en la zona de elementos grandes no puede ser simulada con precisión. Esto puede solventarse introduciendo el recálculo del problema en los bloques de tiempo posteriores al inicial, pero como ya se ha indicado, eso aumenta el tiempo de cálculo y hace que el coste pase a ser parecido al remallado y recálculo completo (m  = 1.800).

A diferencia del problema de Smolarkievicz, el primer ejemplo, donde se obtenían mejores resultados con periodos de remallado grandes en relación al tiempo final de cálculo (pocos, hasta 5 remallados), en este caso se producen con periodos menores (hasta 20 remallados). Esta diferencia es atribuible a que en este caso la zona de interés de la solución se modifica sustancialmente a lo largo del problema. En este caso, además, la solución es mucho más suave, por lo que la influencia del error de interpolación es mucho menor, y efectuar más remallados no influye significativamente en los resultados.

5. Conclusiones

El esquema adaptativo presentado puede aplicarse con éxito a distintas tipologías de problemas tridimensionales de conveción-difusión. Se ha aplicado a un problema con fuertes variaciones de la solución, concentradas en una parte del dominio, y a otro en el que la solución es más suave pero va afectando a diversas partes del dominio a lo largo de su evolución temporal.

El esquema se basa en una estrategia de remallado en la que se impone un volumen máximo de elemento sobre una malla de referencia, que se mantiene constante a lo largo del problema. El remallado permite refinar y desrefinar de forma muy flexible, aumentando y disminuyendo mucho la densidad de elementos en una única iteración. Las mallas obtenidas son de calidad suficiente como para no afectar la eficacia del esquema iterativo utilizado para resolver los sistemas lineales de los problemas discretizados. Dos parámetros permiten modular el proceso. El refinamient se impone de forma relativa a los valores del indicador de error. El desrefinamiento, mediante un término de velocidad que modula la transición de la malla de cálculo hasta la malla de referencia.

El esquema adaptativo propuesto solo recalcula el primer bloque de pasos de tiempo de cálculo. Esto permite reducir el coste computacional en relación con la estrategia de recalcular el problema completamente o recalcular en todos los remallados intermedios. En función del problema, el número de remallados óptimo varía. En el caso con fuertes variaciones en la solución pero que se desarrolla en solo una parte del dominio, el número de remallados debe ser reducido para que los errores de la interpolación no afecten la calidad de la solución. En cambio, en el caso con solución más suave pero que aumenta progresivamente el dominio de interés, el número de remallados debe ser tal que la solución no avance mucho más rápido que la adaptación de la malla; si no, puede ser conveniente imponer recálculo, hasta convergencia, junto con el remallado.

Ante un problema nuevo, se propone utilizar inicialmente un periodo de remallado alto (pocos remallados), y disminuirlo paulatinamente, mientras el coste computacional disminuya. De esta forma se consigue llegar a una solución aceptable con el número de grados de libertad más bajo posible. Mientras el número de remallados sea reducido, el error introducido por la interpolación también lo será. Con esta estrategia puede controlarse el efecto del error de interpolación, comparando sucesivas soluciones para periodos de remallado decrecientes.

El esquema puede ser utilizado con diversos indicadores de error. Se ha verificado la utilidad de un indicador típico para problemas de convección-disfusión, basado en el gradiente de la solución, así como una variante calculada como el gradiente del logaritmo de la solución. Este segundo indicador funciona correctamente para problemas en los que los valores de interés de la solución son varios órdenes de magnitud inferiores a los valores impuestos en las condiciones de contorno.

Agradecimientos

Se agradece el apoyo del Ministerio de Ciencia y Tecnología de España a través de los proyectos con número de referencia CGL2008-06003-C03-02 y CGL2011-29396-C03-00.

References

1. [1] A. Huerta, A. Rodríguez-Ferran, P. Díez, J. Sarrate; Adaptive finite element strategies based on error assessment. Int. J. Numer. Meth. Engng., 46 (10) (1999), pp. 1803–1818
2. [2] R. Biswas, R.C. Strawn; A new procedure for dynamic adaption of three-dimensional unstructured grids; Appl. Numer. Math., 13 (6) (1994), pp. 437–452
3. [3] R. Löhner; An adaptive finite element scheme for transient problems in CFD; Comput. Method. Appl. M., 61 (3) (1987), pp. 323–338
4. [4] S. Ghorai, A.S. Tomlin, M. Berzins; Resolution of pollutant concentrations in the boundary layer using a fully 3D adaptive gridding technique; Atmos. Environ., 34 (18) (2000), pp. 2851–2863
5. [5] D.N. Arnold, A. Mukherjee, L. Pouly; Locally adapted tetrahedral meshes using bisection; SIAM J. Sci. Comput., 22 (4) (2000), pp. 431–448
6. [6] D.N. Arnold, A. Mukherjee, L. Pouly; Adaptive finite elements and colliding black holes; ,in: D. Griffiths, D. Higham, G. Watson (Eds.), Numerical analysis 1997 (Dundee), vol. 380, Longman (1998), pp. 1–15
7. [7] L.A. Freitag, C. Ollivier-Gooch; Tetrahedral mesh improvement using swapping and smoothing; Int. J. Numer. Meth. Engng., 40 (21) (1997), pp. 3979–4002
8. [8] J.R. Shewchuk; What is a good linear finite element? interpolation, conditioning, anisotropy, and quality measures; 11th International Meshing Roundtable, Sandia National Laboratories (2002), pp. 115–126
9. [9] R. Montenegro, G. Montero, J. Escobar, E. Rodríguez, J. González-Yuste; Tetrahedral mesh generation for environmental problems over complex terrains; ,in: P.M.A. Sloot, C.J.K. Tan, J.J. Dongarra, A.G. Hoekstra (Eds.), Computational Science - ICC 2002, vol. 2330, Springer (2002), pp. 335–344
10. [10] A. Oliver, G. Montero, R. Montenegro, E. Rodríguez, J.M. Escobar, A. Pérez-Foguet; Adaptive finite element simulation of stack pollutant emissions over complex terrains; Energy, 49 (2013), pp. 47–60
11. [11] J. Sarrate, A. Huerta; An improved algorithm to smooth graded quadrilateral meshes preserving the prescribed element size; Commun. Numer. Meth. En., 17 (2) (2001), pp. 89–99
12. [12] N. Roquet, P. Saramito; An adaptive finite element method for Bingham fluid flows around a cylinder; Comput. Method. Appl. M., 192 (31-32) (2003), pp. 3317–3341
13. [13] S. Sun, M.F. Wheeler; Anisotropic and dynamic mesh adaptation for discontinuous Galerkin methods applied to reactive transport; Comput. Method. Appl. M., 195 (25-28) (2006), pp. 3382–3405
14. [14] M. Picasso, V. Prachittham; An adaptive algorithm for the Crank-Nicolson scheme applied to a time-dependent convection-diffusion problem; J. Comput. Appl. Math., 233 (4) (2009), pp. 1139–1154
15. [15] F. Alauzet, P.L. George, B. Mohammadi, P. Frey, H. Borouchaki; Transient fixed point-based unstructured mesh adaptation; Int. J. Numer. Meth. Fl., 43 (6-7) (2003), pp. 729–745
16. [16] P.K. Smolarkiewicz; The multi-dimensional Crowley advection scheme; Mon. Weather Rev., 110 (1982), pp. 1968–1983
17. [17] R. Verfurth; Robust a posteriori error estimates for stationary convection-diffusion equations; SIAM J. Numer. Anal., 43 (4) (2006), pp. 1766–1782
18. [18] J. Donea, A. Huerta; Finite Element Methods for Flow Problems; Wiley (2003)
19. [19] A. Rodríguez-Ferran, M.L. Sandoval; Numerical performance of incomplete factorizations for 3D transient convection-diffusion problems; Adv. Eng. Softw., 38 (6) (2007), pp. 439–450
20. [20] C.J. Lin, J.J. Moré; Incomplete Cholesky factorizations with limited memory; SIAM J. Sci. Comput., 21 (1) (2000), pp. 24–45
21. [21] J.M. Cascón, L. Ferragut, M.I. Asensio; Space-time adaptive algorithm for the mixed parabolic problem; Numer. Math., 103 (3) (2006), pp. 367–392
22. [22] A. Lozinski, M. Picasso, V. Prachittham; An anisotropic error estimator for the Crank-Nicolson method: Application to a parabolic problem; SIAM J. Sci. Comput., 31 (4) (2009), pp. 2757–2783
23. [23] V. John; A numerical study of a posteriori error estimators for convection-diffusion equations; Comput. Method. Appl. M., 190 (5) (2000), pp. 757–781
24. [24] E.M. Constantinescu, A. Sandu, G.R. Carmichael; Modeling atmospheric chemistry and transport with dynamic adaptive resolution; Comput. Method. Appl. M., 12 (2) (2008), pp. 133–151
25. [25] F. Garcia-Menendez, M.T. Odman; Adaptive grid use in air quality modeling; Atmosphere, 2 (3) (2011), pp. 484–509
26. [26] R. Montenegro, G. Montero, G. Winter, L. Ferragut; Aplicación de métodos de elementos finitos adaptativos a problemas de convección-difusión; Rev. Int. Métodos Numér. Cálc. Diseño Ing., 5 (4) (1989), pp. 535–560
27. [27] A.B. Díaz-Morcillo, J.V. Balbastre, L. Nuño; New recovery error indicator for adaptive finite-element analysis on waveguiding structures; IEEE T. Microw. Theory, 51 (5) (2003), pp. 1467–1475
28. [28] H. Si; On refinement of constrained Delaunay tetrahedralizations; 15th International Meshing Roundtable (2006), pp. 61–69
29. [29] H. Si; Constrained Delaunay tetrahedral mesh generation and refinement; Finite Elem. Anal. Des., 46 (1) (2010), pp. 33–46
30. [30] A. Staniforth, J. Côté, J. Pudykiewicz; Comments on «Smolarkiewiczs deformational flow»; Mon. Weather Rev., 115 (4) (1987), pp. 894–902
31. [31] W. Hundsdorfer, E.J. Spee; An efficient horizontal advection scheme for the modeling of global transport of constituents; Mon. Weather Rev., 123 (12) (1995), pp. 3554–3564
32. [32] A. Sarma; Application of the multidimensional positive definite advection transport algorithm (MPDATA) to environmental modelling on adaptive unstructured grids; Int. J. Numer. Meth. Fl., 50 (2006), pp. 1247–1268
33. [33] H. Drange, R. Bleck; Multidimensional forward-in-time and upstream-in-space-based differencing for fluids; Mon. Weather Rev., 125 (4) (1997), pp. 616–630
34. [34] A. Bott; A positive definite advection scheme obtained by nonlinear renormalization of the advective fluxes; Mon. Weather Rev., 117 (5) (1989), pp. 1006–1016
35. [35] J.H. Seinfeld; Atmospheric chemistry and physics of air pollution; Wiley (1986)