## 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

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] .

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

 ${\displaystyle {\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 un el vector que contiene la solución en los nodos en el instante tn  = 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 .

## Algoritmo 1.

• while No convergencia do
• Calcular [0, Δt ] (m pasos de δt )
• Remallar
• end while
• Interpolar
• Guardar u1
• 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.

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

 ${\displaystyle {\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:

 ${\displaystyle {\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 tolu fija el límite inferior de los valores de la solución sobre los que ya no se refina.

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:

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

donde Sr es el conjundo de elementos de la malla que pertenece a la región r cuyo indicador es superior a una tolerancia dada, y Ve 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:

 ${\displaystyle 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.

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:

 ${\displaystyle {\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 xy . 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}}}$ .

Las figuras 2 a y 2 c muestran la solución para ${\textstyle t=}$${\displaystyle {\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.

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=}$${\displaystyle {\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 L2 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.

El segundo ejemplo consiste en una emisión puntual en el interior de un dominio Ω = [0, 96000] × [0, 48000] × [0, 3000] m3 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}=}$${\displaystyle {\mbox{D}}_{yy}=50{\frac {m^{2}}{s}}}$ y Dzz 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:

 ${\displaystyle {\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:

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

Se fija δt  = 1 s, y se utiliza el indicador de la ecuación (3) con tolu  = 10−6${\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 .

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−5 .

 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−6 , 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 L2 del error, el número de nodos y el tiempo de cálculo para α de 1 a 102 , 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 α .

 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 β .

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 .

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.

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)

### Document information

Published on 30/03/17

Licence: Other

### Document Score

0

Views 0
Recommendations 0