Aplicación del método Petrov-galerkin como técnica para la estabilización de la solución en problemas unidimensionales de convección-difusión-reacción

Garzón Alvarado, Diego Alexander; Galeano Urueña, Carlos Humberto; Duque Daza, Carlos Alberto

Services on Demand

Journal

Article

Indicators

Cited by SciELO
Access statistics

Revista Facultad de Ingeniería Universidad de Antioquia

Print version ISSN 0120-6230On-line version ISSN 2422-2844

Rev.fac.ing.univ. Antioquia no.47 Medellín Jan./Mar. 2009

Aplicación del método Petrov-galerkin como técnica para la estabilización de la solución en problemas unidimensionales de convección-difusión-reacción

Application of Petrov-galerkin method in stabilization solution of advection-diffusionreaction unidimensional problems

Diego Alexander Garzón Alvarado^* , Carlos Humberto Galeano Urueña, Carlos Alberto Duque Daza

Universidad Nacional de Colombia, Kr. 30 No. 45-03. Edificio 453. Oficina 401, Bogotá, Colombia

Resumen

El presente artículo estudia el método Streamline Upwind Petrov Galerkin como técnica de estabilización de la solución numérica de las ecuaciones diferenciales de advección-difusión-reacción; se analiza el método a la luz de la naturaleza no auto adjunta del operador diferencial convectivo y de las transformaciones necesarias para la estabilización de la solución por medio de la eliminación del efecto no autoadjunto inducido por el término convectivo. Se desarrollaron seis diversos ejemplos numéricos, los cuales incluyen problemas de coeficientes variables, altamente convectivos, fuertemente reactivos, sistemas de ecuaciones diferenciales y soluciones transitorias. Se encuentra un excelente desempeño de esta técnica de estabilización para todos los casos anteriormente mencionados, exceptuando los problemas con términos reactivos fuertes.

Palabras clave: Petrov-Galerkin, advección, difusión, funciones de perturbación, soluciones inestables.

Abstract

This paper examines the Streamline Upwind Petrov Galerkin method as a stabilizing technique for the numerical solution of differential equations of advection-diffusion-reaction; it analizes the method taking into account the non self-adjoint nature of the convective diferential operator and the necessary transformations for the solution stabilization through the elimination of non self-adjoint effect induced by the convective term. Presents six different numerical examples, which include problems of variable coefficients, high convective problems, highly reactive systems and transitional solutions. This method presents an excellent performance of this stabilization technique for all the cases mentioned above, except for the problems with strong reactives terms.

Keywords: Petrov-Galerkin, advection, diffusion, perturbation functions, unstable solutions.

Introducción

El problema de la advección-difusión-reacción en estado estable, consisten en la determinación de una función de campo escalar φ , la cual debe satisfacer la ecuación diferencial:

así como unas condiciones de contorno definidas por las expresiones (2) y (3):

donde k ≥ 0 es el coeficiente difusivo, es el campo de velocidad asociado al proceso advectivo, s es el coeficiente fuente ( s >0 implica producción y s > 0 significa disipación) es la función de generación, g es la función que define el valor del campo escalar φsobre la frontera y h es la función que define el valor del flujo sobre la frontera . La simplificación de esta expresión permite la obtención de ecuaciones específicas de amplio uso en diferentes campos de la física, tales como:

Ecuación de Helmholtz: empleada en los problemas de acústica, electromagnetismo y sismología, la cual se obtiene considerando el término de producción ( s >0 ) y eliminando el término advectivo [1].

Ecuación de Advección-Difusión ( s = 0): utilizada ampliamente en problemas de dispersión de gases contaminantes (sin reacción, especialmente en dispersión de productos de la combustión altamente estables).

Ecuación de Advección-Reacción ( k = 0 ): usada en el modelado de dispersión de contaminantes acuíferos en aguas con alta velocidad.

Ecuación de Difusión-Reacción : la cual tiene aplicaciones en problemas morfogénesis, donde las moléculas difunden a través de los tejidos y reaccionan con otras sustancias directamente y por señalización a través de las células.

Todas estas expresiones son empleadas de forma recurrente en el modelado de diferentes tipos de problemas de transporte de masa y energía, destacándose especialmente el problema del estudio del flujo de fluidos a través de las ecuaciones linealizadas de Navier-Stokes. Otros usos se encuentran en campos como: el financiero, para la estimación de los valores de derivados financieros a través de la ecuación de Black-Scholes [2], el ambiental para el modelado de los procesos de contaminación ambiental [3] o para el análisis de la propagación de incendios [4], el demográfico para el modelado de la distribución de poblaciones y de la transmisión de la información genética [5], y el biológico para modelado de los mecanismos de morfogénesis [6], entre otros muchos campos.

Además de la enorme utilidad para el estudio de fenómenos físicos, algunas formas de la expresión (1) se caracterizan por la dificultad en la implementación de soluciones numéricas y el requerimiento de técnicas especiales para construir una aproximación.

La solución de la ecuación (1) exhibe dos formas generales. La primera se caracteriza por un régimen exponencial (creciente o decreciente), alcanzado cuando las raíces de la ecuación característica son reales, es decir, cuando el discriminante se hace mayor o igual que cero . La segunda se denomina generalmente régimen de propagación y se caracteriza por un comportamiento sinusoidal modulado por una función exponencial. Para casos de coeficientes variables en la ecuación, se pueden alcanzar soluciones que varían su comportamiento de una zona a otra.

En numerosos casos las soluciones numéricas que se han implementado para el desarrollo de este tipo de ecuaciones presentan oscilaciones falsas que no se ajustan al comportamiento real del fenómeno modelado [7], de modo que el desarrollo de técnicas numéricas exactas y confiables para el desarrollo de este tipo de expresiones, se ha convertido en un campo de permanente investigación en los últimos treinta años.

En el presente artículo se presenta un estudio del método Streamline Upwind Petrov-Galerkin (SUPG) para la solución por elementos finitos de diferentes formas de la ecuación de adveccióndifusión- reacción, tanto en régimen transitorio como estacionario, prestando especial atención a los problemas con coeficientes variables. En la primera parte del texto se hace referencia al método, listando los aportes más importantes en los últimos años. Posteriormente se explica y analiza la formulación Petrov-Galerkin en elementos finitos, examinando el tipo y las condiciones requeridas para la estabilidad de la solución. En la tercera sección del documento se presentan diversos casos de estudio que permiten concluir acerca de la exactitud de la herramienta numérica para diversas formas de la ecuación (1), incluyendo casos con coeficientes variables y análisis en estado transitorio. En la última parte del artículo se presentan las conclusiones del trabajo.

Método streamline upwind petrov-galerkin

La solución del problema de advección-difusión- reacción empleando una formulación convencional de elementos finitos (formulación de Bubnov-Galerkin) lleva frecuentemente a soluciones poco exactas o con oscilaciones erróneas (falsas), especialmente en las zonas aguas abajo donde se presentan fuertes variaciones del campo escalar, tal como en los problemas de capa límite 7], exceptuando los casos en los que el término relacionado con la difusión predomina numéricamente sobre los demás (también denominados problemas de difusión dominante). La presencia de estas oscilaciones esta relacionada con los valores de los coeficientes de la ecuación, es decir con la naturaleza del fenómeno predominante. El método de Galerkin presenta significativa divergencia, con respecto a la solución real, para todos los casos en los que se cumple al menos una de las dos condicionesdefinidas en (4) y (5) [8]:

donde Pe se denomina número de Peclet, σ es el número de Damkohler y h es la longitud característica del dominio. La primera condición define todos los casos en donde el fenómeno advectivo predomina sobre el difusivo, en tanto que la segunda acota los casos en los que el problema posee una componente reactiva superior a la difusiva. De acuerdo con lo anterior, una forma de estabilizar la solución puede ser reduciendo la longitud característica del elemento, disminuyendo por tanto el número de Pe, sin embargo para condiciones de advección fuerte esta técnica puede no ser la más apropiada desde el punto de vista de costo computacional. Numéricamente se observa como para números Pe pequeños, los valores propios discretos de la matriz de rigidez alcanzada con la formulación Bubnov-Galerkin, son reales y la solución del problema se puede interpretar como la minimización de un funcional. Sin embargo, para números Pe elevados estos valores propios son imaginarios (en su mayoría), es decir el problema no puede ser derivado de un principio variacional por tratarse de un problema no auto adjunto.

Algunas de las formulaciones de elementos finitos más comunes para la solución de estos problemas se basan en la modificación de la función de peso de la ecuación de residuos ponderados, adicionando una función de perturbación que logra estabilizar el método (éste es denominado método de Petrov- Galerkin), ponderando con mayor peso la información proveniente de los nodos ubicados aguas arriba (de ahí el término streamline upwind). La primera referencia al uso de estas funciones de peso modificadas aparece en 1975 [9], las cuales más tarde fueron empleadas formalmente por Christie [10] y Zienkiewicz [11]. Posteriormente Hughes [12], retoma el fundamento del método streamline upwind para problemas de convección dominante; este método fue reformulado en 1982 por Brooks & Huhges [13], quienes presentan una formulación general del método SUPG como una forma de aproximar la solución del problema generalizado de Stokes. El método SUPG es adoptado en 1984 por Johnson 14] para la solución de problemas que involucran ecuaciones diferenciales hiperbólicas. En 1986 Huhges [15], asocia el método SUPG con el cumplimiento de las condiciones de Babuska- Brezzi en los problemas de dinámica de fluidos. Posteriormente, Huhges [16] emplea los mínimos cuadrados para la obtención de la función de perturbación adecuada en la estabilización de los elementos finitos (Streamline Petrov-Galerkin con mínimos cuadrados, GLS). Franca [17], emplea el método SUPG a la solución de problemas advectivos-difusivos. En 1993 Baiocchi [18] adiciona el empleo de funciones burbuja para el enriquecimiento de las funciones polinomiales, en tanto que Harari [19] en 1994, incluye técnicas de optimización para mejorar el desempeño de los métodos de estabilización, específicamente en problemas de advección-difusión-producción. En 1996 Russo [20] aplica la estabilización con funciones burbuja en la solución de la ecuación linealizada de Navier-Stokes, mientras que Codina [21] en 1998 presenta una comparación de técnicas de estabilización de elementos finitos para la solución del problema de advección-difusión- reacción. Un nuevo método para la solución de este tipo de problemas, es presentado por Cockburn [22], en donde se hace uso de un método de Galerkin empleando funciones discontinuas, el cual brinda la posibilidad de paralelizar los algoritmos desarrollados, emplear mallas irregulares con nodos no necesariamente coincidentes y emplear funciones polinomiales de diferente grado para diferentes elementos, lo cual permite el desarrollo de métodos adaptativos. Araya [23, 24] incorpora al método de estabilización clásico, un estimador de error a posteriori, que permite la implementación de técnicas de mallado adaptativo. Por último, resulta interesante mencionar el trabajo de Xin [25], el cual emplea un término de viscosidad artificial, que permite la estabilización de la formulación del método de los elementos finitos cambiando virtualmente la naturaleza de la expresión.

Planteamiento del problema de advección-difusión-reacción

Al revisar la expresión (1) en su forma unidimensional, resulta claro que esta no se trata de una ecuación diferencial auto adjunta, debido a la presencia del término advectivo. De manera que si se desea plantear éste caso como un problema variacional (minimización de un funcional), se requiere adicionar una función conveniente p a la expresión de residuos ponderados [26], tal como se muestra en (6).

donde W es la función de peso y L la longitud del intervalo. Así, con la debilitación de esta última expresión, resulta claro que una forma de convertir esta, en una expresión auto adjunta, es eliminar el primer término de la integral en (7), es decir, hacer que se cumpla la ecuación (8):

lo cual permite definir la forma adecuada de la función p [26], tal como se muestra en (9).

De aquí que la expresión de residuos ponderados (6) toma ahora la forma de la ecuación (10):

lo que resulta claramente como un problema auto adjunto.

Asumiendo una aproximación en cada elemento (E) definida por la combinación de dos funciones de forma lineales (N_j), tal como se expresa en (11),

se puede reescribir la ecuación de residuos ponderados como se muestra en (12).

La cual se puede expresar por medio de la ecuación (13):

Así, la i-esima expresión de la formulación de Galerkin estándar esta dada por:

donde φj-1,φ j ,φ j+,1 son tres nodos consecutivos. Es de esperarse que una solución alcanzada empleando la ecuación (14) sea una respuesta estable, es decir no oscilante.

Estabilización empleando el método de Petrov-Galerkin

El uso de la expresión de residuos ponderados (15) para la ecuación de advección-difusión-reacción, suponiendo coeficientes constantes y una función de peso lineal de la forma W_i = N_i, permite obtener, en coordenadas locales (ξ), una ecuación de ensamble con la forma de la expresión (17).

Se observa en esta última expresión una diferencia entre los coeficientes que acompañan a las variables φi-1 y φi+1 , esto se entiende como una asimetría en la matriz de rigidez producida por la presencia del término advectivo en la ecuación diferencial (ecuación no auto adjunta), lo que lleva a la obtención de valores propios mayoritariamente imaginarios y una solución inestable. Adicionalmente se puede ver como la ponderación sobre la información proveniente de los nodos ubicados aguas abajo es mayor que la aplicada sobre la información del punto j -1 (aguas arriba), mientras que para casos con advección nula la información de ambos puntos ( j -1, j +1) tiene el mismo peso, lo que lleva a una matriz de rigidez simétrica (este ejercicio es más sencillo si se obvia el término fuente σ = 0 ). El sentido de la formulación de Petrov-Galerkin consiste en modificar la función de peso, adicionándole una función de perturbación, de manera que se logre equilibrar la ponderación sobre los términos φ j-1 y φ j+1, mejorando la simetría de la matriz de rigidez (logrando valores propios reales) y estabilizando las oscilaciones producidas con Bubnov-Galerkin. La función de peso a emplear tiene la forma mostrada en (18), siendo una de las funciones de perturbación más trabajadas la definida en (19).

donde α es un escalar (coeficiente de perturbación) empleado para fortalecer el peso sobre la información del nodo ubicado aguas arriba. De esta forma, la expresión de la matriz de rigidez para un elemento típico, ahora esta definida por (20), mientras que la ecuación de ensamble, encontrada con una función de peso perturbada por la expresión (19), se muestra en (21).

Del análisis del operador diferencial de la expresión (15), resulta claro que la inestabilidad es producida por el término convectivo, más no por la parte reactiva. De este modo, y con el fin de analizar las condiciones que convierten la solución alcanzada mediante la ecuación (21) en una aproximación estable, se pueden comparar las expresiones (14) y (21) para σ = 0 . Así, se puede encontrar fácilmente el valor de α requerido para hacer coincidir las dos expresiones, o en otras palabras, para que la estabilización planteada al introducir una función de perturbación lleve al problema a una forma auto adjunta, con una matriz de rigidez definida con valores propios reales. Así, el valor del coeficiente de perturbación óptimo (α_op) que permite alcanzar soluciones nodales exactas, esta dado por:

De otro lado, la máxima perturbación que puede inducirse en la solución, sin que la matriz de rigidez deje de ser diagonal dominante, se obtiene eliminando el término que multiplica la información nodal φi+1 en la ecuación (14), es decir con un valor α dado por la expresión (23),

Consideraciones acerca de la aplicación del método SUPG

Considerando el caso de la expresión de residuos ponderados para la ecuación de advección-difusión- reacción de coeficientes variables (24), se puede observar como la debilitación de esta expresión, empleando funciones de forma lineales, lleva a la ecuación (25), la cual es análoga a la del caso de coeficientes constantes presentada en (16).

Donde K’ (ξ), u’ (ξ), s’ (ξ) son las funciones que definen los coeficientes de la ecuación diferencial, pero expresados en términos de las coordenadas locales. Empleando la función de peso modificada (19), la matriz de rigidez para este caso se puede escribir así:

Y la ecuación de ensamble resultante estará definida como se muestra en (27).

en donde los coeficientes se encuentran definidos dentro de la expresión (28).

El cálculo de estos términos se puede efectuar empleando integración numérica por cuadratura gaussiana, lo cual permite considerar el efecto de la variación de los términos dentro de cada uno de los elementos, así como el empleo de valores óptimos para el coeficiente de perturbación α en cada uno de los puntos de integración empleados para el elemento, procedimiento que garantiza la estabilidad del método, independientemente del comportamiento del número de Peclet a lo largo del dominio.

SUPG empleando funciones de perturbación de orden superior

Así entonces, la solución estabilizada SUPG se consigue perturbando (o modificando) la función de peso, de manera que se pueda mejorar la ponderación sobre los nodos aguas arriba, lo cual permite un mejor comportamiento en la solución del sistema de ecuaciones lineales resultante. En la Figura 1 se muestran las funciones de peso originales y las funciones de peso perturbadas de acuerdo a (19) para α =1, se observa como en este último caso el área bajo la curva de la función W₂ se ha aumentado respecto a la función original sin perturbar, en tanto que para la función W₁ perturbada dicha área se muestra reducida. Estas áreas son una medida de la ponderación generada por cada una de estas funciones sobre la información nodal (W₂ sobre el nodo j -1 y W₁ sobre el nodo j +1). Se pueden construir diversos tipos de funciones de peso perturbadas que cumplan este comportamiento upwind, sin embargo es necesario que las mismas satisfagan el principio de complementariedad, de modo que se garantice que el método sea global y localmente conservativo. Para una función de peso, como la descrita en (18), construida a partir de funciones de forma lineales y una función de perturbación cuadrática, como se muestra en la figura 2 y se expresa en la ecuación (30), la tendencia a la estabilización es más fuerte que en el caso de la figura 1, tal como se puede deducir de la comparación de las áreas bajo la curva en cada caso. Introduciendo las funciones de perturbación definidas en (29) y (30), dentro de la expresión (26), se obtienen los coeficientes de la ecuación de ensamble, los cuales se muestran en (31).

Figura 1 Función de peso lineal convencional (izq.) y función de peso perturbada con α =1 (der.)

Análisis de casos

Caso 1: Problema de advección-difusión con coeficientes constantes

En este caso se plantea un problema de advección- difusión con coeficientes constantes y una función de generación cúbica. El problema se describe matemáticamente en la ecuación (32)

En la figura 3 se muestra el resultado de la implementación del método de estabilización de Petrov-Galerkin, comparándolo con el comportamiento mostrado por la solución encontrada mediante el método convencional (Bubnov-Galerkin) y la solución analítica exacta. Se observa como la solución Petrov-Galerkin implementada con 10 elementos (Pe = 20, α_opt α_crit = 0,95 en todo el dominio) y un coeficiente de perturbación α = α_opt, es una solución nodalmente exacta y libre de oscilaciones, en tanto que la solución convencional presenta fuertes oscilaciones alrededor de la solución exacta (debido a que P_e > 1). Sin embargo, al aumentar el número de elementos a 20 (P_e = 10, α_opt α_crit = 0,9 ) y 100 (P_e = 2, α_opt ≈ 0,54, α_crit = 0,5 ), se observa un reducción en la amplitud de las oscilaciones y una convergencia hacia la solución exacta debido a la reducción en el número de Peclet.

Figura 2 Función de perturbación cuadrática (arriba-izq.), función de peso perturbada con α =1 (arriba-der.), α con ½ (abajo-izq.), α = 3&/2 (abajo-der.)

Se puede ver como para la malla de 100 elementos la aproximación por Galerkin logra una respuesta adecuada, excepto en las zonas aguas abajo con cambios muy fuertes en el comportamiento de φ (cerca de x = 1), en donde la respuesta numérica aún exhibe una tendencia oscilatoria localizada. La figura 4 muestra el comportamiento del error absoluto acumulado porcentual a lo largo del dominio para cada una de las tres mallas empleadas con el método SUPG, casos en los que se empleó un valor óptimo en el coeficiente de perturbación.

Caso 2: Problema de advección-difusiónreacción con coeficientes constantes

En este problema, planteado en la expresión (33), se analiza la estabilidad del método Petrov-Galerkin incluyendo términos reactivos en la ecuación diferencial y considerando coeficientes constantes en todo el dominio analizado.

Como se puede apreciar en la Figura 5, para una malla formada por 40 elementos el comportamiento de la solución Bubnov-Galerkin es inestable dado que el problema esta caracterizado por valores altos en el número de Peclet (P_e = 7,5), lo cual hace que el problema pueda clasificarse como predominantemente advectivo. Las pruebas numéricas realizadas con diferentes valores en el coeficiente de perturbación (α) muestran que se requiere un valor cercano a a = 0,9 para lograr estabilizar la solución empleando el planteamiento de Petrov-Galerkin con funciones de forma lineales y una función de perturbación como la descrita en la expresión (19). Lo anterior coincide con procedimiento utilizado para la obtención de la expresión (22), el cual lleva a un valor . Adicionalmente, se observa como para valores en el coeficiente de perturbación inferiores al valor crítico, la solución no alcanza la estabilización; en tanto que para valores superiores a α_crit se obtienen aproximaciones estabilizadas pero poco exactas, la cuales se pueden denominar sobre perturbadas. Al igual que en el caso anterior, resulta interesante observar como la amplitud de las respuestas no estabilizadas se hacen más fuertes en las zonas aguas arriba en las que la variación de φ se hace más alta. Así, a medida que se incrementa el valor del coeficiente de perturbación, también se puede ver como se reducen progresivamente las oscilaciones sobre la respuesta, siendo la región cercana a x = 3 la última en alcanzar la estabilidad. En la Figura 6 se observa el comportamiento del error acumulado porcentual para tres valores del coeficiente de perturbación, resulta claro de esta gráfica como para valores diferentes de α =α_opt la aproximación puede perder estabilidad (solución sub-perturbada) o exactitud (solución sobre perturbada).

Figura 3 Soluciones obtenidas para el caso de estudio 1 empleando el método de estabilización Petrov-Galerkin con 10 elementos (arriba-izq), 20 elementos (arriba-der.) y 100 elementos (abajo)

Figura 4 Comportamiento del error absoluto acumulado porcentual alcanzado con el método SUPG para tres tipos de malla diferentes

Caso 3: Problema de advección-difusiónreacción con coeficientes variables

En este caso se presenta la solución implementada para un problema de advección-difusión-reacción con coeficientes variables, el cual es planteado en la ecuación (34).

La figura 7 muestra el comportamiento de los números de Peclet y Damkohler para tres tipos de discretización diferentes (10, 100 y 1000 elementos). Como primer punto se debe resaltar como los valores de 2 σ P_e son muy inferiores que los valores del número de Peclet, por lo que el fenómeno físico modelado mediante la ecuación diferencial puede ser catalogado, por encima de las componentes difusivas y reactivas, como predominantemente advectivo. Se observa claramente, que solo la discretización con 1000 elementos permite una solución estable y exacta empleando un formulación Bubnov-Galerkin (P_e < 1 y 2σ P_e < 1). Para las demás mallas, en ninguna zona se presentan valores de Peclet inferiores a la unidad y por tanto se prevén aproximaciones oscilantes, tal como se puede verificar en la figura 8, en donde se muestran oscilaciones generalizadas sobre todo el dominio para 10 elementos y una oscilación localizada (cercana a x = 3 ) para 100 elementos.

Figura 5 Soluciones obtenidas para el caso de estudio 2 empleando Petrov-Galerkin con 40 elementos y diferentes valores de coeficiente de perturbación α

Figura 6 Comportamiento del error absoluto acumulado porcentual alcanzado con el método SUPG para el caso de estudio 2 para diferentes valores de coeficiente de perturbación α

Adicionalmente, en la Figura 8, se observan los resultados alcanzados para este problema empleando una formulación modificada (Petrov- Galerkin), la cual empleó funciones de forma lineales y perturbaciones obtenidas a partir de la expresión (19), así como un coeficiente de perturbación α = α_opt calculado para cada punto de integración. Se observa como las soluciones estabilizadas alcanzan valores nodales exactos y comportamientos no oscilatorios con pocos elementos, aún en las zonas con altos valores en el gradiente de φ, gracias a que se garantiza en cada punto del dominio un valor óptimo para el coeficiente de perturbación. Esta estrategia de adaptación implementada muestra excelentes resultados frente a la solución convencional (Bubnov-Galerkin) para este tipo de problemas.

Figura 7 Comportamiento de los valores de Peclet, Damkohler, α_opt y α_crit para el caso de estudio 3, empleando 10 elementos (arriba-izq.), 100 elementos (arriba-der.) y 1000 elementos (abajo)

Caso 4: Problema de advección-difusiónreacción con altos valores de reacción

Con este ejemplo se busca observar el comportamiento de problemas con un alto componente reactivo, lo cual aunque no aporta inestabilidad en la solución (debido a que su presencia no hace el operador diferencial no auto adjunto), si incorpora problemas de exactitud en los resultados numéricos alcanzados con el método Bubnov-Galerkin. La ecuación diferencial que modela el problema se encuentra enunciada en la expresión (35).

En este caso, a diferencia del caso de estudio 2, el número 2σ P_e es muy superior al número de P_e, por lo que el fenómeno descrito por la ecuación diferencial es predominantemente reactivo, tal como se puede apreciar en la Figura 9, en donde se ilustra el comportamiento de estos números a lo largo del dominio del problema para cada una de las tres mallas empleadas (10, 100 y 1000 elementos).

Figura 8 ción obtenida para el caso de estudio 3 empleando Petrov-Galerkin con 10 elementos (arriba-izq), 100 elementos (arriba-der.) y 1000 elementos (abajo)

Al igual que en el caso de estudio 3, para la solución Petrov-Galerkin se emplearon valores óptimos de α calculados para cada punto de integración. Sin embargo, como se puede apreciar en la figura 10, la aproximación construida con este método, aunque estable no logra valores exactos, alcanzando desviaciones tan grandes como las obtenidas con la formulación convencional de elementos finitos. El análisis de la naturaleza del método de estabilización SUPG permite concluir que la presencia de un término reactivo fuerte en la ecuación diferencial no modifica las condiciones para la estabilización del método, es decir, los valores de α_opt y α_crit , lo cual se evidencia en el comportamiento de las soluciones alcanzadas en este caso. No obstante, se observa como la corrección upwind no logra mejorar el método en cuanto a la exactitud alcanzada para fenómenos con reacción dominante, problema que se produce por diferencia en las magnitudes de diversos términos en la matriz de rigidez.

Figura 9 Comportamiento de los valores de Peclet, Damkohler, α_opt y α_crit para el caso de estudio 4, empleando 10 elementos (arriba-izq), 100 elementos (arriba-der.) y 1000 elementos (abajo)

Solución obtenida para el caso de estudio 4 empleando Petrov-Galerkin con 10 elementos (arriba-izq), 100 elementos (arriba-der.) y 1000 elementos (abajo)

Caso 5: Solución de un modelo de advección-difusión-reacción

Dado que muchos fenómenos físicos deben ser modelados por sistemas de ecuaciones diferenciales, resulta interesante evaluar el desempeño del método SUPG para problemas de adveccióndifusión- reacción formulados con sistemas de ecuaciones. En este ejemplo se plantea un sistema de dos ecuaciones diferenciales acopladas, con coeficientes constantes y términos de generación, las cuales se escriben en la expresión (36).

Se observa como el término advectivo predomina ampliamente en ambas ecuaciones, por lo que se prevé que las aproximaciones obtenidas con la formulación convencional puedan mostrar inestabilidades. Las soluciones alcanzadas con una malla con 60 elementos se muestra en la figura 11.

Figura 11 Solución obtenida con Petrov-Galerkin para un sistema de dos ecuaciones diferenciales

En estas gráficas se observa un comportamiento estable de las aproximaciones para φ₁ y φ₂ construidas empleando el método Petrov-Galerkin con funciones de forma lineales, una perturbación constante y un coeficiente α = α_opt . Por otro lado, las oscilaciones esperadas en la solución Bubnov-Galerkin, se hacen más notorias aguas abajo sobre la función φ₁, que es la región en donde se presenta una acumulación del error favorecida por la presencia de un alto gradiente.

Caso 6: Solución del modelo de advección-difusión-reacción en estado transitorio

En este ejemplo se estudia el comportamiento de la solución SUPG para problemas transitorios. La ecuación (37) define el problema advectivo difusivo dinámico planteado para en este caso.

La solución planteada se desarrollo con 50 elementos e involucró, al igual que en todos los casos anteriores, un planteamiento Bubnov-Galerkin y una aproximación upwind Petrov-Galerkin para el problema espacial, así como una formulación de integración en el tiempo backward-Euler, para la dimensión restante. Se observa en las soluciones alcanzadas, las cuales se ilustran en la Figura 12, como las inestabilidades exhibidas por la solución convencional se transmiten a lo largo del tiempo, obteniéndose finalmente una aproximación estacionaria inestable, que muestra oscilaciones que no guardan relación con la naturaleza de la ecuación diferencial, como tampoco con las condiciones iniciales del problema. Entre tanto la solución SUPG muestra una estabilización de la solución para cada paso de tiempo, llevando a una respuesta estacionaria muy exacta y libre de falsas oscilaciones.

Figura 12 Solución transitoria alcanzada con el método Bubnov-Galerkin (izq) y SUPG (der.) para el caso de estudio 6

Conclusiones

(i) Se analizó el método Streamline Upwind Petrov Galerkin como técnica para la estabilización de la solución de las ecuaciones diferenciales de advección-difusión-reacción, se observó como esta técnica es equivalente a la transformación del operador diferencial de la ecuación de residuos ponderados, el cual pasa a ser un operador auto adjunto al eliminar el término convectivo presente inicialmente y causante de la inestabilidad en las soluciones por elementos finitos. Por otro lado, se analizó la naturaleza del método a la luz de la modificación de las funciones de peso, por medio de la adición de expresiones de perturbación, lo cual lleva a una sobre-ponderación de la información proveniente de los nodos aguas arriba, una matriz de rigidez con valores propios reales y la obtención de un problema de origen variacional, lo que en últimas es equivalente a la modificación del operador diferencial mencionada anteriormente. (ii) El análisis de los ejemplos numéricos mostró que el uso de coeficientes de perturbación iguales o superiores a α_opt permite la estabilización de la función, aunque valores superiores a α_crit reducen drásticamente la exactitud del método. Para el caso de ecuaciones diferenciales con coeficientes variables, el uso de coeficientes de perturbación óptimos calculados para cada uno de los puntos de integración del elemento, permitió alcanzar soluciones altamente estables y exactas de forma independiente del tipo de malla empleada. Este último método adaptativo presentó excelentes resultados tanto para ecuaciones diferenciales independientes, como para sistemas acoplados y problemas transitorios. Sin embargo, para problemas con términos convectivos fuertes (valores de 2σ P_e muy superiores al número P_e), no se lograron obtener soluciones exactas aunque si estabilizadas. El origen de la inexactitud se atribuye a la diferencia en las magnitudes de los términos que conforman la matriz K , efecto que no esta asociado al problema convectivo, el cual se observa superado con el empleo de las funciones de perturbación. Por último, las pruebas realizadas con problemas transitorios mostraron que mientras las inestabilidades se transmiten en el tiempo hasta la solución estacionaria en el caso de aproximaciones convencionales Bubnov- Galerkin, la solución Petrov-Galerkin estabiliza la respuesta en cada uno de los pasos de tiempo, obteniéndose comportamientos estables y exactos, tanto para cada paso de tiempo como para el sistema en estado estable.

Referencias

1. I. Babuska, F. Ihlenburg, E. Paik, S. Sauter. “A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution”. Computer Methods in Applied Mechanics and Enginnering. Vol. 128. 1995. pp. 325-359.

2. R. Smith. “Optimal and near-optimal advection diffusion finite-difference schemes iii. Black-Scholes equation”. Proceedings: Mathematical, Physical and Engineering Sciences. Vol. 456. 2000. pp. 1019-1028.

3. M. Sanín, G. Montero. “A finite difference model for air pollution simulation”. Advances in Engineering Software. Vol. 38. 2007. pp. 358–365.

4. L. Ferragut, M. Asensio, S. Monedero. “A numerical method for solving convection–reaction–diffusion multivalued equations in fire spread modelling”. Advances in Engineering Software. Vol. 38. 2007. pp. 366–371.

5. O. Richter. “Modelling dispersal of populations and genetic information by finite element methods”. Environmental Modelling & Software. Vol. 23. 2008. pp. 206–214.

6. D. Dan, C. Mueller, K. Chen, J. Glazier. “Solving the advection-diffusion equations in biological contexts using the cellular Potts model”. Physical Review. Vol. 72. 2005. pp. 041909.

7. O. Zienkiewicz, R. Taylor. Finite Element Method. Ed. Butterworth-Heinemann College. Vol. 3. 2000. pp. 5- 150. [ Links ]

8. M. Stynes. “Steady-state convection-diffusion problems”. Acta Numerica. Vol. 14. 2005. pp. 445- 508.

9. O. Zienkiewicz, R. Gallagher, P. Hood. Newtonian and non-Newtonian viscous impompressible flow. Temperature inducedflows and finite elements solutions. The Mathematics of Finite Elements and Applications. Ed. Academic Press. New York. 1975. pp. 13-63. [ Links ]

10. I. Christie, D. Griffiths, O. Zienkiewicz. “Finite element methods for second order differential equations with significant first derivatives”. International Journal for Numerical Methods in Engineering. Vol. 10. 1976. pp. 1389-1396.

11. O. Zienkiewicz, J. Heinrich, P. Huyakorn, A. Mitchel. “An upwind finite element scheme for two dimensional convective transport equations”. International Journal for Numerical Methods in Engineering. Vol. 11. 1977. pp. 131-144.

12. T. Hughes, A. Brooks. “A multidimensional upwind scheme with no crosswind diffusion”. Finite Element Method for Convection Dominated Flows (ASME). Vol. 34. 1979. pp. 19-35.

13. A. Brooks, T. Hughes. “Stream upwind Petrov- Galerkin formulations for convection dominated flow with particular emphasis on the incompressible Navier Stokes equations”. Computer Methods in Applied Mechanics and Engineering. Vol. 32. 1982. pp. 199- 259.

14. C. Johnson, U. Navert, J. Pitkaranta. “Finite element methods for linear hyperbolic problems”. Computer Methods in Applied Mechanics and Engineering. Vol. 45. 1984. pp. 285-312.

15. T. Hughes, L. Franca, M. Balestra. “A new finite element formulation for computational and fluid dynamics: V. Circumventing the Babuska-Brezzi condition: a stable Petrov-Galerkin formulation of the Stokes problema accommodating equal-order interpolations”. Computer Methods in Applied Mechanics and Engineering. Vol. 59. 1986. pp. 85-99.

16. T. Hughes, L. Franca, G. Hulbert. “A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective diffusive equations”. Computer Methods in Applied Mechanics and Engineering. Vol. 73. 1989. pp. 173- 189.

17. L. Franca, S. Frey, T. Hughes. “Stabilized finite element methods: I. Application to advective difusive model”. Computer Methods in Applied Mechanics and Engineering. Vol. 95. 1992. pp. 253-276.

18. C. Baiocchi, F. Brezzi, L. Franca. “Virtual bubbles and Galerkin/least-square type methods”. Computer Methods in Applied Mechanics and Engineering. Vol. 105. 1993. pp. 125-141.

19. I. Harari, T. Hughes. “Stabilized finite element methods for steady advection-diffusion with production”. Computer Methods in Applied Mechanics and Engineering. Vol. 115. 1994. pp. 165-191.

20. A. Russo. “Bubble stabilization of the finite element method for the linearized incompressible Navier-Stokes equation”. Computer Methods in Applied Mechanics and Engineering. Vol. 132. 1996. pp. 335-343.

21. R. Codina. “Comparison of some finite element methods for solving the diffusion-convection-reaction equation”. Computer Methods in Applied Mechanics and Engineering. Vol. 156. 1998. pp. 185-210.

22. B. Cockburn, G. Karniadakis, C. Shu. Discontinuos Galerkin Methods, Theory, Computational and Application. Ed. Springer. 2000. pp. 20-80. [ Links ]

23. R. Araya, E. Behrens, R. Rodríguez. “An adaptive stabilized finite element scheme for the advection– reaction–diffusion equation”. Applied Numerical Mathematics. Vol. 54. 2005. pp. 491–503.

24. R. Araya, E. Behrens, R. Rodríguez. “Error estimators for advection-reaction-diffusion equations based on the solution of local problems”. Journal of Computational and Applied Mathematics. Vol. 206. 2007. pp. 440- 453.

25. J. Xin, J. Flaherty. “Viscous stabilization of discontinuous Galerkin solutions of hyperbolic conservation laws”. Applied Numerical Mathematics. Vol. 56. 2006. pp. 444-458.

26. G. Guymon, V. Scott, L. Herrmann. “A general numerical solution of the two dimensional diffusionconvection equation by the finite element method”. Water Resources Research. Vol. 6. 1970. pp. 1611- 1617.

(Recibido el 2 de abril de 2008. Aceptado el 6 de noviembre de 2008)

^*Autor de correspondencia: teléfono: + 57 + 1 + 316 50 00 ext. 14062, fax: + 57 + 1 + 316 53 33, correo electrónico: dagarzona@unal.edu.co. (D. Garzón)