Análisis de las soluciones numéricas de la ecuación de convección - dispersión - adsorción

Revista hternacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 3,4, 38%409(1987) 1 ANALISIS DE LAS SOLUCIONES NUMERICAS DE LA ECUACION DE CONVECCION - DISPERSION - ADSORCION PATRICIA M. CARPANO* CARLOS A. GRATTONI* SUSANA C. GABBANELLI* Y MIRTHA S. BIDNER** *Departamento de Ingeniería Química Facultad de Ingeniería Universidad Nacional de La Plata 1 esq. 47, 1900 La Plata, Argentina. **Departamento de Ingeniería Quámica Facultad de Ingeniería Universidad de Buenos Aims Pabellón de Industrias, Ciudad Unaversitaria 1427 Buenos Aires, Argentina RESUMEN La ecuación diferencial parcial no lineal que describe el flujo bicomponente miscible a través de medios porosos, con dispersión y adsorción del tipo Langmuir, se resuelve numéricamente. Se aplican cuatro métodos diferentes: explícito, Barakat-Clark, CrankNicolson y ecuaciones diferenciales ordinarias. Se estima el error de truncación de los distintos métodos. Se obtienen las condiciones de estabilidad de los mismos para el caso de adsorción lineal y se infieren las condiciones de estabilidad para la ecuación no lineal. A efectos de comparar los cuatro métodos, se siguen dos caminos. Primero, se define un error global que evalúa las diferencias entre cada una de las soluciones numéricas y la solución analítica que representa el caso de adsorción lineal. Utilizando este error se comparan las soluciones numéricas entre sí. Además se estudia la influencia de los parámetros de la ecuación diferencial y de los incrementos espaciales y temporales de la discretización en las soluciones obtenidas. En especial, se discute el método Barakat-Clark. Segundo, para la ecuación con adsorción no lineal, las cuatro soluciones numéricas se comparan contra resultados experimentales, por medio de un error que evalúa las diferencias entre cada solución y los datos experimentales. SUMMARY The partial differential nonlinear equation which describes the bicomponent miscible flow .through porous media with dispersion and Langmuir equilibrium adsorption is numerically worked out. Four different methods are applied: explicit, Barakat-Clark, CrankNicolson and ordinary differential equations. Truncation error of the methods are estimated. Recibido: Junio 1987 OUniversitat Politecnica de Catalunya (España) ISSN 02 13-1315 SOLUCIONES NUMERICAS DE LA ECUACION DE CONVECCION 391 siendo +, la concentración de soluto; +,, la adosrción del soluto en el medio poroso; Ií, el coeficiente de dispersióin; V, la velocidad intersticial; r, la coordenada temporal; y X, la espacial. La ecuación (1) supone: - medio poroso homogéneo con sección transversal y porosidad constantes; - flujo isotérmico y unidimensional; - velocidad intersticial (obtenida dividiendo la velocidad Darcy por la porosidad), constante; - dispersión del soluto sólo en la dirección longitudinal; - coeficiente de dispersión independiente de la concentración química; - no hay reacción química entre la solución inyectada y la roca o el fluido in-situ. El término % en (1) representa la adsorción de soluto en la roca, considerando que ésta es gobernada por la isoterma de Langmuir, que supone equilibrio instantáneo, as, s,T = - 1 + Ps, El parámetro de adsorción a es adimensional y el parámetro P tiene unidades de inversa de la concentración. Diferenciando 1; ecuación (2) e introduciéndola en la ecuación (l), resulta: as, a2$ as, = IL- - va~ ax2 a~ donde El medio poroso puede suponerse semi-infinito o finito, dando lugar a distintas condiciones de bqrde" En ]nuestra comparación utilizarnos.la siguiente condición de borde: donde es la concentración de inyección. La condición inicial es, donde 392 P.M. CARPANO, C.A. GRATTONI, S.C. GABBANELLI Y M.S. BIDNER l - Se definen las siguientes variables adimensionales: $ X ur C=- x=- t=- $0 I I (6) u y parámetros adimensionales vl NPe = - K b = BIlo (7) siendo Npe el número de Péclet; 1 la longitud del medio poroso y b parámetro de adsorción adimensional. Las ecuaciones (3), (4) y (5) resultan ac 1 a2c ac ~(c)~ = -- - - ax2 ax (8) a = 1 + [ (1 + b~)'] C(0,t) = 1 ; t>O C(x, O) = o SOLUCIONES ANALITICAS Las soluciones analíticas para la ecuación de convección-dispersión (g(C) = 1) y distintas condiciones de borde están dadas en la bibliografía1. Considerando las ecs. (8), (9) y (lo), la solución está representada por C = - erfc - - NP~ x+t l[ ~x~t)+e~pe~erfc($7)] 2 (11) que puede ser extendida al caso de adsorción lineal (g(C) = 1 + a) con sólo cambiar t por t/(l + a). SOLUCIONES NUMERICAS DE LA ECUACION DE CONVECCION 393 SOLUCIONES NUMERICAS La ecuación (8) fue resuelta numéricamente, considerando x = 5 como infinito, en una computadora IBM 4361 por los métodos que se describen a continuación, donde c,u y w representan las aproximaciones numéricas de C. 1. Método explícito5, (E): aproxima la ecuación (8) de una manera simple, en forma centrada alrededor del punto (i, n): 2. Método de Barakat-Clark4, (BC): la ecuación es aproximada por las siguientes expresiones en diferencias finitas. - "Forward" , (F) i u?+, - un - u;+1 + uyfl U~+l - n+l - --[ 'i-i + u:~~~u~] NP~ Az2 1 - [ 2Az (13) - "Backward", (B) La concentración cY+l puede ser evaluada por: un+', la solución de (F), se calcula explícitamente a partir de la ecuación (13) desde x = O en una sucesión de i crecieiites. De la misma. manera, w?+', la solución 394 P.M. CARPANO, C.A. GRATTONI, S.C. GABBANELLI Y M.S. BIDNER de (B), se calcula de la ecuacióii (14) comenzando en x = 5 en una sucesión de .i decrecientes. 3. Método de Crank-Nicolson5, (CN): al aplicarlo a la ecuación (8) resulta 'n+l - n+l - [ <+! 'i-i + '?+í - e?-1 4Ax 4Ax 1 La solución es obtenida resolviendo ecuaciones algebraicas simultáneamente e iterando en cada.paso de tiempo debido a la no linealidad de la ecuación. 4. Método de las Ecuacjones Diferenciales Ordinarias7, (ODE): se discretiza la parte espacial para, reducirla a una ecuación ordinaria. Aproximando las derivadas parciales espaciales y el termino de adsorción como en el método explícito, se obtiene Para realizar la integración numérica se probaron varios métodos: Milne (50. orden, paso variable), Runge-Kutta (40. orden, paso fijo y variable), Stiff (paso variable), Adams (20. orden, paso fijo). Entre ellos se eligi6 el de Runge-Kutta debido a que produce iguales errores que los de 119ilne y Stiff pero utiliza menor cantidad de pasos. El de Adams tiene mayores errores que el de Runge-Kutta de paso fijo para el mismo At. Ante la necesidad de comparar con los otros métodos (E, BC, CN) se utilizó el de Runge-I<iitta de paso fijo, sacrificando la estabilidad incondicional que presenta el de paso variable. Para comparar la exactit,ud de los distintos métodos frente a las diferentes situaciones se usó la siguiente expresión del error, SOLUCIONES NUMERICAS DE LA ECUACION DE CONVECCION 395 debido a que es la más representativa del comportamiento global de la solución, donde 1 es el número de pasos en el espacio y N el número de pasos en el tiempo. Error de truncación Para la ecuación no lineal por expansiones en serie de Taylor se estima el error de truncación (E:)~ de los métodos utilizados; con la excepción del método ODE, para el que se estima el error de truncación en forma empírica. - Explícito: Se observa que la no linealidad de la ecuación (8) no influye sobre el orden de aproximación. - Barakat-Clark: el error de truncación se estima examinando cada una de las ecuaciones (13) y (14), 396 P.M. CARPANO, C.A. GRATTONI, S.C. GABBANELLI Y M.S. BIPNER d3wp 1 d4~? 1 + %:] At2- + [ai" zg(":) - a,zaii 4Np, axat2 4 [a2~? i lat [@w: i - -- -- a3w: ']AxA~+ dxdt Npe Ax ¿3x3 dt 6 Np, dx2 dt 4 d3w: 1 d41u: 1 dx3 6 8x4 12 Np, Debido a que cada una de las ecuaciones (13) y (14) puede ser usada para aproximar la solución de la ecuación (8), cada uno de los términos en las expresiones (20) y (21) tiene aproximadamente la misma magnitud. Por lo tanto al calcular los perfiles de concentración según la ecuación (15), aquellos términos que tienen signos opuestos tenderían a cancelarse y comoconsecuencia el error de truncación sería E: E O(At) + 0(Ax2) P2) Nótese que para el caso de flujo sin adsorción (g(C) = 1) el término de O(At) en las ecuaciones (20) y (21) se anula y por lo tanto el error de truncación sería E; E 0(At2) + 0(Ax2) (23) - Crank-Nicolson n+i n+b n+" [u3ai3 ~(ci '1 + d2ci "+* ac;+* ag(cn+?> :- E i 24 dt2 dt aen+$ 4 1 a3e7+! 1 --- + --] at2+ 8x2 dt2 8 Npe d5 dt2 8 = 0(At2) + 0(Ax2) La expresión anterior muestra como en el caso explícito que la adsorción no modifica el error de truncación. SOLUCIONES NUMERICAS DE LA ECUACION DE CONVECCION 397 - Ecuaciones Diferenciales Ordinarias: se puede decir que a causa de la discretización espacial el error es de 0(Aa2). El método de Runge-Icutta para la integración temporal tiene un error de 0(At5). Por lo tanto: Con respecto al orden de aproximación el método más conveniente es el ODE. Análisis de estabilidad La estabilidad de un esquema en diferencias queda determinada por la no amplificación de los errores de un cierto nivel de tiempo sobre el siguiente, según Peaceman5, Aziz7 y otros8~'. Para comenzar el análisis de la estabilidad se toma el caso sin adsorción con lo cual la ecuación (8) resulta ser una ecuación lineal a coeficientes constantes. En este tipo de ecuaciones es condición necesaria y suficiente para que un esquema bien planteado resulte estable8*' que se verifique el criterio de von Neuman5J0. De la aplicación de dicho criterio se obtuvieron las condiciones siguientes, salvo para el ODE, cuyo análisis se hace empíricamente: - Explícito - Barakat-Clark y "Backward" - Ecuaciones diferenciales ordinarias Las condiciones anteriores pueden ser generalizadas al caso de adsorción lineal (b - 0) con sólo reemplazar At por At/(l+ a). De esta manera la elección de At y Ax es menos restrictiva que en el ca.so de flujo sin adsorción. En base a pruebas computacionales, considerando adsorción tipo Langmuir se observó que las condiciones de estabilidad están comprendidas entre las determinadas para los dos casos: sin adsorción y con adsorción lineal. Entonces, para el caso no lineal, si se satisfacen las ecuaciones (26), (27) y (28) respectivamente, las soluciones son estables. 404 P.M. CARPANO, C.A. GRASTONI, S.C. GABBANELLI Y M.S. BIDNER Figura 5. Variación del error global para las soluciones "forward", (F), "backward", (B), y su promedio, (BC). A: en función del incremento espacial. B: en función del increment,~ temporal. La variación del error local con los parámetros de discretización en el espacio se muestra en la Figura 9-A y en el tiempo en la Figura 9-B, para los cuatro métodos. Cualitativamente el comportamiento es similar al de las Figuras 3 y 4, respectivamente. En la Figura 9-A se visualiza que todos los métodos tienen un comportamiento similar para Ax 2 0.015625. A menores Ax, los métodos E y BC presentan mínimos y aumentos suaves del error. Los otros dos, CN y ODE mantienen el error casi constante. En la Figura 9-B se muestra que para todos los métodos el error disminuye con At tendiendo a un valor común. El mínimo que muestra el método de CN se debe a que éste fue utilizado para la optimización antes mencionada para hallar los parámetros Npe, a, b. SOLUCIONES NUMERICAS DE LA ECUACION DE CONVECCION 405 Figura 6. Perfiles de concentración, pára flujo sin adsorción, de las soluciones analítica, "forward", (F), y "backward", (B). CONCLUSIONES Para resolver la ecuación de convección-dispersión-adsorción se han utilizado cuatro métodos numéricos. Estos fueron comparados con la solución diferencial (para el caso de adsorción lineal) y con datos experimentales (para el caso de adsorción no lineal). De esta comparación surge: . CN es el de menor riesgo, por ser el único incondicionalmente estable, si bien bajo ciertas condiciones puede presentar pequeñas oscilaciones que se amortiguan con el tiempo, al igual que los demás métodos. Muestra un buen comportamiento frente a los parámetros físicos (Npe,a) y de discretización (Ax,At). Lamentablemente requiere mayor tiempo de CPU y en general tiene un error global grande. . ODE tiene como ventajas sobre el CN su mayor orden de aproximación, menor error global y menor tiempo de cálculo, aunque es condicionalmente estable. . BC es un método explícito, muy superior al explícito tradicional. Requiere menor tiempo de cómputo que los demás. Su principal inconveniente es el comportamiento anómalo frente a los parámetros con la aparición de mínimos pronunciados, lo cual puede causar errores mayores que los previsibles, aunque SOLUCIONES NUMERICAS DE LA ECUACION DE CONVECCION 407 éstos sean siempre menores que los de las soluciones "Forward" y "Backwardn que lo componen. E es el más desventajoso en todo sentido. I Figura 9. Variación del error local para todos los métodos. A: en función del incremento espacial. B: en función del incremento temporal. NOMENCLATURA parámetros para el modelo de adsorción de Langmuir, [adimensionales] concentración de solu to, [adimensional] concentración de soluto experimental, [adimensional] aproximaciones numéricas de C, [adimensional] coeficiente de dispersión longitudinal, [L2 / TI error global definido por la ecuación (18) error local definido por la ecuación (31) coeficiente de las ecuaciones (3) y (8) número de bloques en el espacio longitud del medio poroso, [L] número de pasos en el tiempo coeficiente de dispersión, definido por la ecuación (7), [adimensional] 408 P.M. CARPANO, C.A. GRATTONI, S.C. GABBANELLI Y M.S. BIDNER t tiempo, [adimensional] v velocidad intersticial, [L/T] x distancia, [adimensional] Letras Griegas E error de truncación A incremento de la variable /? parámetro para el modelo de adsorción de Langmuir, [L3/ M] S, concentración de soluto, [M/L~] S, cantidad de soluto adsorbido/volumen de fluído [M/L~] T tiempo, [TI X distancia, [L] Sub y Supraindices i índice de la coordenada espacial n índice de la coordenada temporal o de inyección, x = O AGRADECIMIENTOS Los autores agradecen la colaboración brindada por el Señor Juan Ramón Melendi en la realización de las figuras de este trabajo. REFERENCIAS K.H. Coats. y B.D. Smith, "Dead-Eiid Pore Volume and Dispersion in Porous Media", Soc. Pet. Eng. J., Vol. 4, pp. 73-84, (1964). S.P. Gupta y R.A. Greenkorn, "Dispersion During Flow in Porous Media with Bilinear Adsorption", Water Resour. Res., Vol. 9, pp. 1357-1368, (1973). A. Satter, Y.M. Shum, W.T. Adams y L.A. Davis, "Chemical Transport in Porous Media with Dispersion and Rate-Controlled Adsorption", Soc. Pet. Eng. J., Vol. 20, pp. 129-138, (1980). H.Z. Barakat y J.A. Clark, "On the Solution of the Diffusion Equations by Numerical Methods", J. Hent Transfer, Vol. 88, pp. 421-427, (1966). D.W. Peaceman, "Fundamentals of Numerical Reservoir Simulation", Elsevier Scientific Pub. Co., Amsterdam, The Netlierlands, (1978). M.T. Szabo, "Some Aspects of Polymer Retention in Porous Media Using a Cid-Tagged Hydrolized Polyacrylamide", Soc. Pet. Eng. J., Vol. 15, pp. 323-337, (1975). K. Aziz Y A. Settari, "Petroleum Reservoir Simulntion", Elsevier Scientific Pub. Co., Amsterdam, The Netherlands, (1979). A.R. Mitchell, "Computational Methods in Partial Dijferential Equations", John Wiley and Sons, London, Great Britaiii- (1969). SOLUCIONES NUMERICAS DE LA ECUACION DE CONVECCION 409 9. R.D. Richtmeyer y K.W. Morton, '(Diflerence Methods for InitialValue Problems", John Wiley and Sons, New York, U.S.A., (1967). 10. J. von Neuman y R.D. Richtmeyer, "A Method for the Numerical Calculation of Hydrodynamic Shocks" , J. of Appl. Phys., Vol. 21, pp..232-237, (1950). 11. D.U. von Rosenberg, "Methods for the Numerical Solution of Partial Diflerential Equations", Gerald L. Farrar and Associates, Oklahoma, U.S.A., (1977). 12. H.S. Price, R.S. Varga y J.E. Warren, "Application of Oscillation Matrices to Diffusion-Convection Equations" , J. hfath. Phys., Vol. 45, pp. 301-311, (1966). 13. C.A. Grattoni, P.M. Carpano y M.S. Bidrier, "Analysis of the Nonlinear Adsorption Parameters", (En preparación).