SciELO - Scientific Electronic Library Online

 
vol.25 número1DETECCIÓN DE LA ISLA URBANA DE CALOR MEDIANTE MODELADO DINÁMICO EN MEXICALI, B.C., MÉXICODISEÑO DE FÁRMACOS ANTICANCEROSOS DERIVADOS DE CIS-PLATINO(II) MEDIANTE LA TÉCNICA QSAR, BASADA EN EL FUNCIONAL DE LA DENSIDAD índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Revista

Articulo

Indicadores

Links relacionados

Compartir


Información tecnológica

versión On-line ISSN 0718-0764

Inf. tecnol. vol.25 no.1 La Serena  2014

http://dx.doi.org/10.4067/S0718-07642014000100016 

ARTÍCULOS 

SOLUCIÓN NUMÉRICA DE UNA ECUACIÓN DEL TIPO ADVECCIÓN-DIFUSIÓN

NUMERICAL SOLUTION OF AN ADVECTION-DIFFUSION EQUATION

Blanca Bermúdez y Luis Juárez

Facultad de Ciencias  de la Computación, BUAP; 14 Sur y San Claudio, Ciudad Universitaria, Puebla, Puebla-México (e-mail: bbj@solarium.cs.buap.mx)


Resumen

Se analizan dos métodos de diferencias finitas para resolver ecuaciones de advección-difusión no estacionarias y se aplica a algunos ejemplos. Se compara un método iterativo de punto fijo y otro resultante del trabajo con las matrices A y B que resultan de la discretización del laplaciano y del término advectivo, respectivamente. Este segundo método resultó ser mucho más eficiente y puede ser usado para resolver las ecuaciones de Navier-Stokes  y Boussinesq en diferentes formulaciones, con números de Reynolds grandes. El tiempo de procesamiento fue considerablemente más corto al trabajar con las matrices A y B  que con el método iterativo de punto fijo para los ejemplos estudiados.

Palabras clave: ecuaciones advección-difusión, número de Reynolds, diferencias finitas, métodos numéricos


Abstract

Two methods for solving non steady advection-diffusion equations are analyzed and applied to some examples. Two methods are compared:  a fixed point iterative method and another one that is obtained from working with the matrixes A and B that result from the discretization of the laplacian and of the advection term, respectively. This second method has proved to be much more efficient and can be used for solving the Navier-Stokes and Boussinesq equations in different formulations, with large Reynolds numbers. The processing time was found to be considerable lower for the method with matrixes A and B than for the fixed point iterative method, for the examples studied.

Keywords: advection-diffusion equations, Reynolds number, finite differences, numerical methods.



INTRODUCCIÓN

Una ecuación de advección-difusión puede ser resuelta sin ninguna dificultad si el término difusivo domina al advectivo; cuando esto no sucede se tienen graves problemas. Sobre este tema se ha realizado mucho trabajo de investigación; por ejemplo, se han propuesto métodos de Petrov-Galerkin para el control de las oscilaciones de Gibbs observadas a lo largo de capas límite internas para el problema de advección-difusión. Véase por ejemplo, Mitsukami et al (1985), Codina (1993), De Sampaio et al (2001), Do Camo et al (2003), Knobloch (2006), Oñate et al, (2006). Un estudio sobre  la situación actual de métodos para el control de oscilaciones espurias fue hecho en John V. et al (2007). Diversos métodos fueron desarrollados sobre los marcos existentes de métodos globalmente estabilizadores para el control de las oscilaciones de Gibbs en casos de reacción dominante. Entre estos  están los trabajos de Franca et al (2000), Hauke et al (2001), Hauke (2002) y Haukeet al (2007).

En este trabajo, se discuten algunos modelos para resolver ecuaciones tipo  advección-difusión y se analizan dos métodos de diferencias finitas para resolver las ecuaciones no estacionarias aplicándolos a dos ejemplos. En particular centraremos nuestra atención en el problema que se describe en lo que sigue.

Dada una región , T > 0, Ω RN  N=2,3, con frontera Γ, se desea aproximar la solución de un problema parabólico-hiperbólicode la forma:

(1)

donde puede describir la temperatura o concentraciónde algún contaminante del fluido de velocidad w, el cual por ser incompresible satisface . Ecuaciones de este tipo representan un modelo para una clase de fenómenos de difusión en los cuales hay presente flujo;  y  corresponden al transporte de u a través del proceso de difusión y el efecto de convección respectivamente.

Para esta ecuación, en general, no es posible obtener la solución exacta, lo cual implica el tener que recurrir a métodos numéricos para obtener una aproximación a la solución. La solución exacta del problema satisface el principio del máximo:

(2)

Si la magnitud de es muy pequeña en comparación de la de ε, las soluciones numéricas también satisfa-cen dicho principio. Sin embargo, si esto no sucede, el esquema numérico a utilizar debe preservar el prin-cipio del máximo pues en caso contrario ocurren oscilaciones no físicas muy grandes. Una forma de medir qué relación deben tener w y ε es a través del número de Péclet, el cual se define como:

es muy pequeña en comparación de la de ε, las soluciones numéricas también satisfacen dicho principio. Sin embargo, si esto no sucede, el esquema numérico a utilizar debe preservar el principio del máximo pues en caso contrario ocurren oscilaciones no físicas muy grandes. Una forma de medir qué relación deben tener  w y ε es a través del número de Péclet,el cual se define como:

(3)

Una ecuación de advección-difusión puede resolverse numéricamente sin ninguna dificultad siempre y cuando dicho número sea pequeño, o sea, cuando la magnitud de w es pequeña en comparación con la de ε; sin embargo, si esto no sucede, se tienen problemas, a menos que se impongan restricciones muy fuertes con respecto al tamaño de la malla (malla muy fina) o que se recurra a esquemas más sofisticados. Para ilustrar lo anterior, tomemos el siguiente problema unidimensional y estacionario:

(4)

Cuando ε>>w, la solución tiene un efecto de capa límite (boundary layer) de ancho , o sea,cuando o sea, cuando el número de Péclet es grande. Nótese que cuando ε se hace muy pequeño, la ecuación degenera en una ecuación de primer grado, para la cual, la prescripción de 2 condiciones de frontera es inadecuada.

EL PROBLEMA DISCRETO

El  problema en el caso de este trabajo se resuelve mediante el método de diferencias finitas, que consiste, como se sabe, en reemplazar cada una de las derivadas que aparecen en una ecuación diferencial, por una aproximación diferencia-cociente apropiada; ver Burden y Faires (2011). En este caso vamos a discretizar tanto en el tiempo, como en el espacio. Primero discretizaremos en el tiempo usando el esquema de segundo orden.

 Sea ∆t>0, el tamaño de paso en el tiempo, y tn,entonces tenemos lo siguiente:

(5)

con un =u(x, y, tn)

(6)

En este caso, u1 es un valor que hay que generar con un método de un paso como Euler hacia atrás, para poder trabajar con el método de dos pasos arriba descrito. Discretizamos ahora en espacio. Poniéndonos en el nivel de tiempo n+1 y tomando h>0 como tamaño de paso de discretización en el espacio, definimos:

Para cada punto de la red usamos la serie de Taylor en la variable x alrededor de xi, para generar una fórmula de diferencia centrada para aproximar el Laplaciano:

(7)

donde

Ahora, para aproximar el gradiente usamos:

(8)

y similarmente para la parcial con respecto a y, de manera que obtenemos:

(9)

es decir, se tiene el problema

(10)

donde la matriz A es la matriz asociada a la discretización del Laplaciano y la matriz B es la matriz asociada a la discretización del término advectivo.

Entonces, poniéndonos en el nivel de tiempo n+1 tenemos que resolver el siguiente problema:

(11)

SOLUCIÓN DEL SISTEMA DE ECUACIONES

Muchos problemas en la práctica requieren de la solución de sistemas de ecuaciones lineales, en los cuales la matriz A es rala. Sistemas de este tipo surgen por ejemplo, al aplicar métodos de diferencias finitas o de elemento finito para encontrar una solución aproximada de una ecuación diferencial en derivadas parciales. Generalmente los métodos usuales de eliminación, no pueden ser aplicados puesto que sin precauciones especiales tienden a llevar a la formación de matrices intermedias más o menos densas, haciendo que el número de operaciones aritméticas necesarias para encontrar la solución crezca considerablemente. Por estas razones, uno se inclina a escoger un método iterativo para resolver dichos sistemas de ecuaciones.

Algo común para estos métodos es que en una iteración, se requiere de una cantidad de trabajo equivalente a la multiplicación de A por un vector; cantidad muy pequeña si la matriz es rala. Las matrices que surgen al aplicar métodos de diferencias finitas a ecuaciones diferenciales parciales tienen por lo general una estructura de bloques. En particular, para resolver nuestro sistema de ecuaciones se han usado los siguientes métodos: (1) Gauss-Seidel; (2) Gradiente Conjugado; (3) Método iterativo de punto fijo (Nicolás (1991)). El Gradiente Conjugado, como se sabe, puede ser usado siempre y cuando resolvamos un problema como el planteado en la sección siguiente donde tenemos una matriz ) que es simétrica y positiva definida.

Método iterativo de punto fijo

Una vez discretizada la ecuación tanto en tiempo como en espacio, su solución implica resolver en cada nivel de tiempo sistemas lineales de dados por:

C + B U = G (12)

donde , siendo A la matriz que se asocia al laplaciano, y B, la asociada a la discretización del término advectivo.

Supongamos

R U = C + B U - G (13)

El problema es equivalente a:

R U = 0 (14)

Resolvemos el problema anterior mediante el siguiente método iterativo (Nicolás, 1991):

Dado U0=Un resolver hasta convergencia (14)

CUm+1=CUm-PR Um , p>0,

y tomamos

Se tiene entonces un sistema de ecuaciones lineales a resolver en cada iteración del método de punto fijo.

Se tiene entonces un sistema de ecuaciones lineales a resolver en cada iteración del método de punto fijo.

Esquemas de upwind

Una de las mayores dificultades que encontramos en la solución numérica de la ecuación del transporte reside en la discretización del término advectivo. Daremos entonces una breve descripción de esquemas deupwind para diferencias finitas.

Upwind en una dimensión:

Para aproximar el término advectivo, si se usan diferencias centradas (Thomasset, 1982)

en el nodo J (16)

puede haber oscilaciones no físicas muy grandes. Consideremos:

u 0 =u 1 =0

(17)

Si, , la solución tiene un efecto de capa límite (boundary layer) de ancho

Las oscilaciones ocurren cuando .

Notamos que cuando ε  se hace muy pequeño, la ecuación degenera en una ecuación de primer grado para la cual, la prescripción de 2 condiciones defrontera es inadecuada. La solución numérica se ve contaminada dentro del dominio por la condición de frontera en x=1. En este caso, ya se han establecido remedios para las oscilaciones en el contexto de diferencias finitas.

Consideremos el siguiente esquema de diferencias finitas:

, si w>0 , (18)
, si w>0 , (19)

Entonces, la solución no sufre oscilaciones. Este esquema es de primer orden.

Para lograr obtener una buena precisión tomaremos:

(20)

donde para w>0 y para w<0

Upwind en dos dimensiones:

Para aproximar  se toma (Ikeda, 1982):

(21)

Experimentosnuméricos

En el caso específico de este trabajo se resolvieron 2  problemas. Para  el problema 1  se tiene f=1, w=(1,0), y Ω la región del canal.  Para este problema no se conoce la solución exacta, pero en el límite cuando t tiende a infinito y ε tiende a cero, la solución tiende a :

Fig. 1 Región del canal

En la Tabla.1 se muestran los resultados obtenidos con el método de Iteración de punto fijo y Gauss-Seidel. Los resultados comparativos entre los dos métodos son mostrados con diferentes valores de ε. Como puede observarse, el tiempo requerido por el esquema de punto fijo es hasta 40 veces mayor que el requerido al trabajar con las matrices A y B. La figura.2 muestra la gráfica obtenida con h=0.0625 y ε=000001.

Tabla 1: Resultados con el método de Iteración de punto fijo y Gauss-Seidel (H=0.0625)

Fig. 2: Solución para H = 0.0625  y e = 0.000001

El ejemplo 2 es un problema introducido en Zalesak (1979) (véase también Nadukandi etal , 2012 ) que simula la advección de un cuerpo sólido sujeto a un campo de velocidad angular constante. El cuerpo sólido es modelado con una función dedensidad escalar que tiene tres formas, un cilindro ranurado, un cono y una joroba sinusoidal. Los datos del problema son: Ω = 0,1 x [0,1] , w= (0.5-y, x-0.5), ε = 10-30, y f=0.

Definimos un radio R=0.15, un vector de posición arbitrario r = x,y ∈ Ω y un vector específico de posición ra = xa, ya ∈ Ω para algún punto elegido a. La solución inicial puede ser expresada como sigue:

(22)

Donde H esla función de Heaviside definida por:

(23)

r1= 0,5,0,75, r2= 0,5,0,25, r3= 0,25,0,5 el vector posición correspondiente al centro del cilindro ranurado, el cono y la joroba sinusoidal respectivamente. La solución inicial vista a (-20º,20º). Se muestraen la figura.3. Se impone la condición de Dirichlet = 0  en la frontera. En la Tabla.2 se muestran los tiempos para  a) t=π/2, b) t=π, c) t=3π/2 y d) t=2π con los dos métodos y en la figura.4 se muestran los resultados obtenidos con  a)t=π/2, b) t=π, c) t=3π/2 y d) t=2π, vistos a (-20º,20º). En este caso Δt=π/1000 y h=1/192.  

Como puede observarse, todavía para este valor de h se observan algunas pequeñas oscilaciones, pero desafortunadamente ya no se pudo reducir el valor de h, porlas limitaciones del equipo de cómputo del que se dispone en este momento. El esquema de upwind de Thomasset (1982) proporciona bastantes oscilaciones, porlo que los resultados presentados son usando el upwind de Ikeda (1982). Puede observarse de la Tabla.2 que el tiempo requerido por el método iterativo es 6  veces mayor.

Tabla2: Resultados con el método de Iteración de punto fijo y Gauss-Seidel (H=π/1000)

Fig. 3. Datos iniciales para el segundo ejemplo.

Fig. 4. Soluciones obtenidas para el ejemplo 2 para a) t=π/2, b) t=π, c) t=3π/2 y  d) t=2π, vista a (-20º,20º).

CONCLUSIONES

Como puede observarse de las tablas mostradas, para el ejemplo 1 el método iterativo resultó hasta 40 veces más lento y para el segundo ejemplo hasta 6 veces más lento que trabajar con las dos matrices, la matriz A asociada a la discretización del laplaciano y lamatriz B asociada a la discretización del gradiente. En el primer ejemplo se obtuvieron muy buenos resultados con ambos métodos. En el segundo ejemplo, comoya se observó, todavía se presentan algunas pequeñas oscilaciones, por lo que se necesitaría reducir todavía más el valor de h, sin embargo, por las limitaciones del equipo de cómputo de que disponemos ya no fue posible hacerlo.

El esquema iterativo de punto fijo es usado en Nicolás (1991).  La idea de usar el método iterativo era trabajar con una matriz que fuese simétrica y positiva definida.  Dicho esquema funcionó, sin embargo, el tiempo de procesamiento fue en general muy grande, especialmente al disminuir el valor de ε. En este trabajo se utilizó la matriz completa, la cual no es simétrica pero sepuede probar que es estrictamente dominante diagonalmente para Δt suficientemente pequeño; para esto se usó el método de Gauss-Seidel.

En conclusión, el tiempo de procesamiento fue considerablemente más corto al trabajar con las matrices A y B  que con el método iterativo de punto fijo para ambos ejemplos. La idea ahora es resolver las ecuaciones de Navier-Stokes con números de Reynolds grandes usando upwind y trabajando con las matrices A y B.También es posible usar esta idea para el caso de Boussinesq.

 

REFERENCIAS

Bermúdez, B., Nicolás A., Sánchez F. J., Buendía E.,  Operator Splitting and upwinding for theNavier-Stokes equations, Computational Mechanics, 20 (5) pp. 474-477 (1997).         [ Links ]

Bermúdez, B.,Nicolás A, Sánchez F. J., Un esquema numérico para convección natural, InformaciónTecnológica  10(3) pp. 121-127 (1999).         [ Links ]

Bermúdez, B. y Nicolás A., An operator splitting numerical scheme for thermal/isothermal incompressible viscous flows, InternationalJournal for Numerical Methods in Fluids 29 (4) (1999).         [ Links ]

Bermúdez, B. y Nicolás A., Flujos isotérmicos/térmicosincompresibles y viscosos con formulación velocidad-vorticidad, Información Tecnológica 21 (3) pp. 39-49 (2010).         [ Links ]

Burden, R. L., Faires, J.D., Numerical analysis, Brooks/Cole Cengage Learning, Ninth Edition (2011)        [ Links ]

Codina, R. A., A discontinuity-capturingcroswind-dissipation for the finite element solution of theconvection-diffusion equation, comput. Methods Appl. MEch. Engrg. 110, pp.325-342 (1993).         [ Links ]

De Sampaio, P.A.B.. Coutinho A.L.G.A., A natural derivation of discontinuity capturing operator for convection-diffusion problems. Comput. Methods Appl. Mech. Engrg. 190 (46-47) pp. 6291-6308 (2001).         [ Links ]

Do Camo, E. G. D., Alvarez G.B., a new stabilizedfinite element formulation for scalar convection-diffusion problems: thestreamline and approximate upwind/Petrov-Galerkin method, Comput. Methods Appl.Mech. Engrg. 192 (31-32) pp. 3379-3396 (2003).         [ Links ]

Franca, L.P., Valentin F., On an improved unusual stabilized finite element method for the advective–reactive–diffusiveequation, Comput. Methods Appl. Mech. Engrg. 190 (13–14) pp. 1785–1800 (2000).         [ Links ]

Hauke, G.,  Garçia-Olivares A., Variational subgrid scale formulations for the advection–diffusion–reaction equation, Comput. Methods Appl. Mech.Engrg.190 (51–52) pp. 6847–6865 (2001).         [ Links ]

Hauke, G., A simple subgridscale stabilized method for the advection–diffusion–reaction equation,   Comput. Methods  Appl. Mech.  Engrg. 191 (27–28) pp. 2925–2947 (2002).         [ Links ]

Hauke, G., Sangalli G., Doweidar M. H.,Combining adjoint stabilized methods for the advection–diffusion–reaction problem, Math. Models  Methods Appl. Sci.17 (2) pp. 305–326 (2007).         [ Links ]

Ikeda, T., Maximum Principle and Finite Element Models for convection-Diffusion phenomena, Lecture Notes in Numerical and applied analysis Vol. 4, North Holland Publ. Co. Amsterdam-New York-Oxford (1982).         [ Links ]

John, V., Knobloch P., On spurious oscillations at layers diminishing (SOLD) methods for convection–diffusion equations: Part I – A review, Comput. Methods Appl. Mech. Engrg. 196 (17–20)pp. 2197–2215 (2007).         [ Links ]

Knobloch, P., Improvements of the Mizukami–Hughes method for convection–diffusion equations, Comput. Methods Appl. Mech. Engrg. 1996(1–3) pp. 579-594 (2006).         [ Links ]

Mitzukami, A., Hughes T.J. R. , A Petrov Galerkinfinite element method convection dominated flow: An accurate upwinding technique for satisfying the maximum principle. Comput. Methods  Appl. Mech. Engrg. 50, pp. 181-193 (1985).         [ Links ]

Nadukandi, P., Oñate E., García J., A high-resolution Petrov-Galerkin method for the convection-diffusion-reaction problem Part II- Amultidimensional extension. Computer Methods in Applied Mechanics and Engineeringpp. 327-352 (2012).         [ Links ]

Nicolás, A.,  A finite element approach to the Kuramoto-Sivashinski equation, Advances inNumerical Equations and Optimization, Siam (1991).         [ Links ]

Nicolás, A. y Bermúdez B., 2D incompressible viscousflows at moderate and high Reynolds numbers, CMES (6), (5) pp. 441-451 (2004).         [ Links ]

Oñate, E. ,  Zarate F., Idelsohn S.R., Finite elementformulation for convective–diffusive problems with sharp gradients using finitecalculus, Comput. Methods  Appl. Mech.Engrg. 195 (13–16) pp.1793-1825 (2006).         [ Links ]

Thomasset, F., Implementation of Finite Element Methods for Navier-Stokes equations}, Springer Verlag New York Inc (1982).         [ Links ]

Zalesak, S. T., Fully multidimensional flux-corrected transport algorithms for fluids, J. Comput. Phys. 31(3) pp.335-362 (1979).         [ Links ]

 

Recibido May. 2, 2013; Aceptado Jun. 24, 2013; Versión final recibida Sep. 14, 2013