Apuntes_sny_completos.pdf

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS

CONTENIDO.

1. Introducción. 2. Formulación de problemas de flujo de fluidos en medios porosos. 3. Aproximación de ecuaciones diferenciales en diferencias finitas. 4. Simulación numérica de flujo monofásico en yacimientos. 5. Modelos de pozos en la simulación numérica de yacimientos. 6. Solución de sistemas lineales de ecuaciones. 7. Introducción a la simulación numérica de flujo multifásico. 8. Aspectos prácticos de la simulación numérica de yacimientos.

BIBLIOGRAFIA. 1. Aziz, K. and Settari, A.: Petroleum Reservoir Simulation. Elsevier Applied Science Publishers, New York, 1979. 2. Peaceman, D.W: Fundamentals of Numerical Reservoir Simulation, Elsevier Scientific Publishing, Co. Amsterdam, 1977. 3. Mattax, C.C. and Dalton, R.L.: Reservoir Simulation, SPE Monograph Series, SPE Richardson, Tex., 1990. 4. Artículos Técnicos Selectos.

1

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO

DIVISIÓN DE ESTUDIOS DE POSGRADO DE LA FACULTAD DE INGENIERÍA.

SECCIÓN INGENIERÍA PETROLERA

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS

“APUNTES”

Dr. Femando Rodríguez de la Garza M. en I. Agustín P. Galindo Nava

ENERO 2000

2

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS CONTENIDO DEL CURSO 1. Introducción a la Simulación Numérica de Yacimientos, (SNY). 1.1 Enfoque de la SNY. 1.2 Tipos de simuladores. 2. Formulación de Problemas de Flujo de Fluidos en Medios Porosos. 2.1 Ecuaciones de flujo multifásico composicional. 2.1.1 Consideraciones y Antecedentes. 2.1.2 Ecuaciones de Flujo. 2.2 Flujo Multifásico Pseudocomposicional: Fluidos Tipo Beta y Beta-Modificado. 2.2.1 Flujo de Fluidos Tipo Beta-Modificado. 2.2.2 Flujo de Fluidos Tipo Beta. 2.3 Flujo Monofásico. 2.3.1 Flujo de Aceite. 2-3.2 Flujo de un Aceite de Compresibilidad Pequeña y Constante. 2.3.3 Flujo de Gas. 2.4 Condiciones iniciales y de frontera. 2.4.1 Condiciones Iniciales. 2.4.2 Condiciones de Frontera. 3. Aproximación de Ecuaciones Diferenciales en Diferencias Finitas. 3.1 Diferencias Finitas. 3.2 Aproximaciones a la Primera Derivada. 3.2.1 Diferencias progresivas. 3.2.2 Diferencias regresivas, y ; 3.2.3 Diferencias centrales. ∂ ⎛ ∂p ⎞ 3.3 Aproximación de los Términos de la Forma: ⎜ λ ⎟ ∂x ⎝ ∂x ⎠ ∂p 3.4 Aproximación de los Términos de la Forma: ∂t 3.5 Notación de las Ecuaciones Aproximadas en Operadores en Diferencias. 3.6 Los Esquemas Explícitos, Implícito y Crank-Nicholson de discretización en Tiempo. 3.6.1 Esquema Explícito. 3.6.2 Esquema Implícito. 3.6.3 Esquema Crank-Nicholson. 3.7 Conceptos de Consistencia, Convergencia y Estabilidad de una Aproximación Numérica. 3.8 Construcción de Mallas 3.8.1 Malla Uniforme de Nodos Distribuidos. 3.8.2 Malla uniforme de Bloques Centrados. 3.8.3 Definición de la malla radial en coordenadas cilíndricas. 3.8.3.1 Malla Radial de nodos distribuidos.

3

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

3.8.3.2 Malla Radial de bloques centrados. 4. Simulación Numérica del Flujo Monofásico Unidimensional. 4.1 Flujo Monofásico Unidimensional. 4.2 Aproximación de las ecuaciones diferenciales en deferencias finitas. 4.3 Acoplamiento de las Condiciones de Frontera. 4.4 Cálculo de las Transmisibilidades en las fronteras: Ti + 1 , T 1 . 2

i−

4.4.1 Evaluación de la permeabilidad absoluta en las fronteras, k 4.4.2 Evaluación de la movilidad en las fronteras, λ 4.5

i±

i±

2

1 2

1 2

Métodos de Linealización de las Ecuaciones de Flujo en Diferencias Finitas. 4.5.1 Método de Newton-Raphson: El Método General. 4.5.2 Estructura Matricial del Sistema Lineal de Ecuaciones. 4.5.3 Casos Particulares del Método de Newton. 4.5.3.1 Método de Linealización Directa. 4.5.3.2 Método Semi-Implícito Linealizado.

5. Simulación Numérica del Flujo Monofásico Multidimensional. 5.1 Flujo Monofásico Bidimensional. 5.1.1 Aproximación de las ecuaciones diferenciales en diferencias finitas. 5.1.2 Acopiamiento de las Condiciones de Frontera. 5.1.3 Método de Linealización de las Ecuaciones de Flujo en Diferencias Finitas. Método do Newton-Raphson. El Método General. 5.1.4 Estructura Matricial del sistema lineal de ecuaciones. 5.2 Ordenamiento Normal. 5.3 Ordenamiento D4. 5.4 Flujo Monofásico Tridimensional. 5.4.1 Aproximación de las ecuaciones diferenciales en diferencias finitas. 5.4.2 Acoplamiento de las Condiciones de Frontera. 5.4.3 Método de Linealización de las Ecuaciones de Flujo en Diferencias Finitas. Método de Newton-Raphson. El Método General. 5.4.4 Estructura Matricial del sistema lineal de ecuaciones. 5.5 Anotaciones para concluir. 6. Modelos de Pozos. 6.1 Modelo de pozo de Peaceman para un yacimiento homogéneo e isotrópico. 6.2 Modelo de pozo de Peaceman para un yacimiento heterogéneo y anisotrópico. 6.3 Modelo de Hoo-Jeen Su para pozos fuera de centro. 7. Solución del Sistema de Ecuaciones Lineales. 7.1 Métodos Directos. 7.1.1 Solución de Sistemas Tridiagonales: Algoritmo de Thomas. 7.1.2 Algoritmos de matrices dispersas: Descomposición LDU. 7.2 Métodos Iterativos. 7.2.1 Método de Jacobi. 4

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

7.2.2 Método de Gauss-Seidel. 7.2.3 Método de Sobrerelajamiento en Punto o PSOR. 7.2.4 Método de SOR en Línea o LSOR. 7.2.5 Algoritmo para el cálculo del parámetro de sobrerelajación. 7.2.6 Método de SOR en Bloque o BSOR. 7.2.7 Procedimiento iterativo Implícito de Dirección Alternada o ADIP. 7.2.8 Procedimiento Fuertemente Implícito o SIP. 7.3 Solución de Sistemas de Ecuaciones estructurados en Bloques. 7.4 Métodos Directos Versus Métodos Iterativos. 8. Introducción a la Simulación Numérica del Flujo Multifásico. 8.1 Flujo Multifásico Unidimensional de Fluidos Tipo Beta. 8.2 Forma de las Ecuaciones de Flujo en Diferencias Finitas. 8.3 Cálculo de las Transmisibilidades en las Fronteras. 8.4 Acoplamiento de las Condiciones de Fronteras a las Ecuaciones Aproximada. 8.5 Análisis de no-linealidades y tratamiento. 8.5.1 Método IMPES: Stone y Garder, 1961. 8.5.1.1 Estructura matricial del sistema de ecuaciones. 3.5.2 Totalmente Implícito. 8.5.2.1 Estructura matricial del sistema de ecuaciones. 8.6 Extensiones al Flujo Multifásico Bidimensional. 9. Aspectos Prácticos de la Simulación Numérica de Yacimientos. 9.1 Información necesaria para un estudio. 9.2 Construcción de la malla de cálculo. 9.3 Ajuste de historia del yacimiento. 9.4 Predicción del comportamiento futuro del yacimiento. BIBLIOGRAFIA. 1. Aziz, K. and Settari, A.: Petroleum Reservoir Simulation. Elsevier Applied Science Publishers, New York, 1979. 2. Peaceman, D.W: Fundamentals of Numerical Reservoir Simulation, Elsevier Scientific Publishing, Co. Amsterdam, 1977. 3. Mattax, C.C. and Dalton, R.L.: Reservoir Simulation, SPE Monograph Series, SPE Richardson, Tex., 1990. 4. Artículos Técnicos Selectos.

5

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

1. Introducción a la Simulación Numérica de Yacimientos, (SNY). El objetivo de la Ingeniería de Yacimientos es obtener la máxima recuperación, económicamente posible, de hidrocarburos de un yacimiento petrolero. Para lograr este objetivo, el ingeniero de yacimientos emplea principios científicos para desarrollar modelos del yacimiento. Estos se usan para simular el comportamiento del yacimiento ante diversas opciones de producción y recuperación de hidrocarburos a lo largo de su vida productiva.

1.1 Enfoque de la Simulación Numérica dé Yacimientos. Los modelos son básicamente de dos tipos: 1. Modelos físicos. 2. Modelos matemáticos. En estas notas, los modelos matemáticos son los que ocuparán nuestra atención. El modelamiento matemático en la ingeniería de yacimientos, se refiere a la representación de los procesos de transferencia de masa, y en algunas instancias de energía, que ocurren en el medio poroso, el yacimiento, a través de un conjunto de ecuaciones diferenciales y a su solución matemática. Las ecuaciones diferenciales, constitutivas de un modelo, se obtienen básicamente de la aplicación de los principios de conservación de masa y de energía en un volumen elemental, representativo del medio poroso. Las ecuaciones diferenciales se complementan con: 1. Ecuaciones de estado, que describen el comportamiento volumétrico y de fases de los fluidos. 2. Ecuaciones de movimiento de las fases en el medio poroso, como son la ecuación de Darcy ó la ecuación de Forcheimmer en el caso más general para flujo no Darciano. 3. Ecuaciones adicionales apropiadas. En el caso general, las ecuaciones diferenciales de flujo de fluidos en medios porosos son ecuaciones no lineales. En situaciones particulares, como es el caso de flujo monofásico de un fluido ligeramente compresible y de compresibilidad constante, las ecuaciones diferenciales adquieren formas lineales de manera que pueden resolverse mediante métodos analíticos clásicos. Generalmente las ecuaciones diferenciales deben resolverse numéricamente. Esto ha dado origen a la disciplina de Simulación Numérica de Yacimientos. La simulación numérica de yacimientos consiste en el desarrollo de técnicas y métodos para resolver numéricamente las ecuaciones diferenciales de flujo de fluidos en medios porosos, y en la aplicación de esos modelos numéricos en el estudio del comportamiento de yacimientos. El proceso de solución numérica de las ecuaciones diferenciales consiste básicamente en obtener una representación aproximada de las ecuaciones en puntos específicos del espacio y del tiempo: esto se logra mediante el empleo de métodos finitos, tales como diferencias finitas o elementos finitos. Así, las ecuaciones diferenciales de flujo, cuyo dominio de aplicación en espacio y tiempo es continuo, son reemplazadas por un sistema algebraico de ecuaciones. Este sistema se

6

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

genera al aplicar en un cierto nivel de tiempo, las ecuaciones de flujo aproximadas en puntos predeterminados, nodos, de la malla de cálculo que discretiza al yacimiento, de tal manera que cuando se usa un simulador es necesario subdividir al yacimiento en una serie de celdas o bloques. Idealmente estas celdas deben ser lo suficientemente pequeñas para evitar que los errores de truncamiento sean grandes, y en donde en cada una de ellas se consideran constantes las propiedades del yacimiento y de los fluidos, ver Fig. 1.1.

Fig. 1.1 Discretización del yacimiento.

Generalmente, el sistema de ecuaciones aproximado es no lineal. Existen diversos métodos para resolver ese sistema, cada uno con características propias de estabilidad numérica. Desde el punto de vista practico; la estabilidad se mide por la capacidad de un simulador de resolver problemas difíciles de flujo, como es el caso de variaciones fuertes de presión y de saturación en alguna región dada del yacimiento. En problemas de flujo multifásico, el método IMPES ocupa el extremo inferior de implicitud, mientras que el Método Totalmente implícito ocupa el extremo superior. La estabilidad tiene sin embargo un precio, que es un mayor requerimiento de memoria y tiempo de cómputo para resolver una etapa de tiempo dada. En contraparte, el tamaño del paso de tiempo, Δt , con que puede avanzar la simulación es considerablemente mayor. Existen incluso problemas donde el método IMPES no puede aplicarse, como es el caso de problemas de conificación, o de fuerte segregación gravitacional. Independientemente del método empleado en la linealización del sistema de ecuaciones aproximadas, el resultado final es un sistema algebraico de ecuaciones lineales. La solución del sistema lineal puede básicamente obtenerse mediante métodos directos y métodos iterativos. Los métodos directos realizan un numero fijo de operaciones para resolver un sistema dado de ecuaciones; de particular interés para la simulación numérica de yacimientos son los métodos de matrices bandadas y de matrices dispersas. Los métodos iterativos se basan en la aplicación de algoritmos cíclicos, en espera en que a medida que las iteraciones progresen, el método converja a la solución.

7

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Los simuladores numéricos de yacimientos han sido ampliamente usados debido a que ellos pueden resolver problemas que de otra manera sería imposible hacerlo. No existe otro método que pueda describir el flujo multifásico (gas-aceite-agua) - multidimensional (ID a 3D) en un yacimiento real. Simulación numérica es la única herramienta que puede describir cuantitativamente el flujo multifásico en un yacimiento heterogéneo a partir de un programa de producción, el cual no solo es determinado de las propiedades del yacimiento sino que también por demandas del mercado, de las estrategias de inversión y por las regulaciones gubernamentales. Debe quedar claro que la respuesta de un simulador numérico en un estudio de simulación siempre es cualitativa por todos los factores que intervienen. Simulación numérica frecuentemente es el mejor método de solución, debido a que puede dar una respuesta inmediata, barata y/o más creíble que otros métodos. Lo anterior no quiere decir que siempre simulación sea el mejor método de análisis por si solo para un problema de ingeniería de yacimientos en particular. Algunos métodos de análisis que complementan con simulación numérica de yacimientos, son las pruebas de variación de presión, las observaciones directas en el campo, las pruebas de laboratorio, las pruebas piloto en el campo, los análisis matemáticos sencillos (soluciones analíticas) y la extrapolación del comportamiento de otros yacimientos. Los problemas deben resolverse con el método más sencillo y barato posible, pero que permita simular el proceso de desplazamiento con suficiente realismo que conduzca a decisiones adecuadas. El ingeniero de yacimientos debe primero determinar el apropiado nivel de simplificación para que en base a ese nivel pueda seleccionar el método de análisis técnicamente factible. Dentro de los factores que se deben considerar para esta adecuada selección están los siguientes: 1. Tipo y complejidad del problema (Geometría del sistema, heterogeneidad del yacimiento, tipo de fluidos presentes y tipo de proceso de depresionamiento a considerarse). 2. Calidad de la respuesta necesaria para las decisiones administrativas del yacimiento. 3. Tiempo disponible para terminar el estudio de simulación. 4. Factores económicos. 5. Disponibilidad y calidad de los datos requeridos para el estudio. 6. Capacidades del simulador numérico y características del sistema de cómputo disponibles. Simulación numérica de yacimientos puede ofrecer al personal involucrado dentro de un estudio de simulación, lo siguiente: 1. Al Ingeniero un rápido método matemático compresivo del análisis del comportamiento del yacimiento. 2. Al Encargado del proyecto un estudio completo del yacimiento que puede incluir comparaciones adecuadas de alternativos planes operativos. 3. Al ingeniero inexperto una efectiva herramienta de entrenamiento. 4. A terceras personas (gobierno) un sistemático análisis del comportamiento del yacimiento, creíble que incluye restricciones de operación que se pueden especificar con detalle.

8

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

1.2 Tipos de simuladores. 1) Modelo tanque (dimensión cero): es muy útil cuando se necesitan respuestas rápidas y el comportamiento de la presión promedio del yacimiento se considera como el único factor importante al realizar las decisiones operativas o de inversión. Los gradientes de presión deben ser pequeños.

2π

h h w L re

rw Modelo de 0D (horizontal) Sistema coordenado Cartesiano

Modelo de 0D (horizontal) Sistema coordenado Cilíndrico

2) Modelo de 1D: estos modelos no pueden emplearse en estudios de yacimientos con un gran ancho y/o espesor debido a que no pueden modelar la eficiencia de barrido areal o vertical. Estos modelos pueden usarse para investigar la sensibilidad del comportamiento del yacimiento en la variación de ciertos parámetros del mismo. Estos modelos son muy útiles en la evaluación de la influencia de las heterogeneidades en la dirección preferencial de flujo. 2π

Δxi

h

i-1

i

i+1 L

Modelo de 1D (horizontal) Sistema coordenado Cartesiano

h

w

radial rw

re

Modelo de 1D (horizontal) Sistema coordenado Cilíndrico

9

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

3) Modelo de 2D (Areal): estos modelos son usados cuando el patrón de flujo areal domina el comportamiento del yacimiento. Se deben emplear estos modelos para estudiar yacimientos completos, pero con un espesor no muy importante o bien la mayoría de los estudios con estos modelos areales en 2D incluyen el uso de pseudofunciones para tomar en cuenta el flujo de fluidos en la dirección vertical. También se emplean en la optimización de los ritmos de producción, de los requerimientos de los pozos, de las facilidades en la superficie, en la localización óptima de los pozos productores así como de los pozos inyectores de gas y/o agua y del tiempo óptimo de las instalaciones de un equipo de producción artificial o de la modificación de las facilidades en la superficie. Estos modelos son empleados para estimar la recuperación final (calculo de reservas), para probar diferentes plataformas de producción y en las estrategias operativas en la recuperación final.

Δxi Pozo inyector de gas

0

.

.

.

.

. . . Ly

.

.

. . .

.

.

. i-1, j

. .

.

.

.

. .

.

. .

. .

.

. .

. .

. i+1, j

. .

. .

.

. .

i, j-1

. i, j

i, j+1

.

. .

h

. .

.

Δyi

. Pozo productor de aceite

0

Lx Modelo de 2D (areal). Sistema Coordenado Cartesiano

Δθj

.

i+1, j i,. j . i-1, j

. .

i, j-1.

.

. . i, j+1

. .

. . .

.

.

cima . h

. . . .

base

radial rw

re

Δri Modelo de 2D (areal). Sistema Coordenado Cilíndrico.

10

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

4) Modelo de 2D (Sección Transversal): estos modelos son empleados para desarrollar las funciones de pozos o pseudofunciones para usarlas en los modelos areales de 2D y de esta forma reducir los requerimientos computacionales y el tiempo de proceso. También son empleados en la simulación de la inyección periférica, inyección de gas en el casquete u otros procesos en los cuales las velocidades frontales hacia los pozos productores son altamente uniformes. Estos modelos a su vez pueden emplearse para evaluar la interacción de las fuerzas gravitacionales, capilares y viscosas y el efecto resultante que tienen sobre las eficiencias de desplazamiento y barrido verticales. En los modelos de sección en 2D, pero considerando una geometría cilíndrica (r-z) se pueden usar para representar el flujo convergente o divergente en una región radialmente simétrica del yacimiento, así como evaluar el comportamiento del pozo cuando los efectos verticales dominan el comportamiento: como la conificación de gas y/o agua. w

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

i, k-1

i+1, k

i, k

i-1, k

i, k+1

. . . . . . .

. . . . . . .

Δz k

Δx i

Modelo de 2D (Sección Transversal). Sistema Coordenado Cartesiano.

.

.

. . . .

.

. . . . . . . . . . . . .

2π

. .

.

.

.

. . . .

. .

. i, k-1 .

. .

i-1, k

i, k

. .

.

.

. .

.

.

Cima

i+1, k

. i, k+1

Vertical

Angular

Δz k

. Base

radial rw

Δri

re

Modelo de 2D (sección). Sistema Coordenado Cilíndrico.

11

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Una importante aplicación de los modelos en 2D, es el evaluar el efecto de los ritmos de producción en el comportamiento del yacimiento y en la recuperación final. 5) Modelo en 3D: estos modelos son empleados cuando la geometría del yacimiento puede ser demasiado compleja para reducirla a modelos en 2D areales o de sección transversal. Por ejemplo, yacimientos con grandes extensiones de barreras de flujo pero con zonas permeables donde ocurre flujo cruzado y que son difíciles de modelar con modelos en 2D. O bien cuando los mecanismos de los fluidos del yacimiento pueden ser complejos para analizarlos con modelos en 2D, por ejemplo, yacimientos con un gran estado de depresionamiento caen centro de esta categoría. Estos modelos también son usados cuando el desplazamiento a ser estudiado está dominado por el flujo vertical, por ejemplo, cercas del pozo en donde, tanto del casquete de gas como del acuífero, puede ocurrir una conificación. Ocasionalmente, los modelos en 2D pueden ser mas problemáticos y caros que los simuladores en 3D debido a que el modelamiento de algunos yacimientos que son arealmente complejos y altamente estratificados pueden requerir docenas o cientos de conjunto de pseudofunciones.

Ly . .

0

. .

. . .

0 Pozo inyector de gas . . . . . o. . . . . . . . . . . . i, j-1,k. . . . . . i-1, .j,k i, j,k. i+1, .j,k . . . . . . i, j+1,k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i, j,k-1

. .

. .

. .

.

.

.

. . i-1, j,k i, j,k . .

i, j,k+1

.

.

. . i+1, j,k . . .

.

. .

. .

.

.

.

Δxi

. . . . . .

o. . . o. . . aceite Pozo productor de . . . .

. . . .

. .

. . .

.

. .

. . . . .

. .

.

.

.

. .

. .

. . . .

i, j,k

.

.

i, j,k

.

.

i,j+1,k

. . .

i, j,k

.

. .

.

.

. . .

i, j-1,k.

. .

. .

.

h

Δz k

.

.

Δy i

Lx

Modelo de 3D. Sistema Coordenado Cartesiano

12

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Angular

Modelo en 3D. Sistema Coordenado Cilíndrico.

Habiendo presentado una breve introducción a la simulación numérica de yacimientos, se pretende ahora, en los diversos capítulos que componen estas notas, presentar con mayor detalle los aspectos teóricos que sirven como base a esa disciplina. En el capitulo 2 se hace una revisión de la formulación de problemas de flujo multifásico en medios porosos; se parte del problema general de flujo multifásico composicional para después obtener, como casos particulares, las ecuaciones de flujo multifásico de fluidos tipo Beta modificado, tipo Beta y de flujo monofásico. En el capitulo 3 se presentan los fundamentos de diferencias finitas. Este material es básico para la aproximación de las ecuaciones diferenciales de flujo. En este capitulo se introducen los esquemas de aproximación explícito, implícito y Crank-Nicholson. Se revisan los conceptos de consistencia, convergencia y estabilidad de un algoritmo numérico. En el capitulo 4 se aborda el problema de simulación numérica de problemas de flujo monofásico unidimensional y multidimensional. Se presenta el método de Newton-Raphson para linealizar el sistema de ecuaciones de flujo aproximadas y se obtienen otros métodos de linealización como casos particulares de éste. Se analiza la estructura de los sistemas de ecuaciones generados por problemas uni-, bi- y tridimensionales. En el capitulo 5 se revisa la simulación numérica de flujo multifásico multidimensional. Se presentan los métodos IMPES y de Newton-Raphson, mejor conocido en la literatura petrolera como el método Totalmente Implícito. Se analiza la estructura de los sistemas de ecuaciones y las similitudes con el problema de flujo monofásico. Finalmente, en el capitulo 6 se discuten las técnicas de solución de sistemas lineales de ecuaciones.

13

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

2. Formulación de Problemas de Flujo de Fluidos en Medios Porosos. En este capítulo se presenta primeramente la formulación de las ecuaciones que describen el flujo multifásico isotérmico de fluidos composicionales en un yacimiento petrolero. Posteriormente se consideran casos particulares de esas ecuaciones, tales como el flujo multifásico de fluidos tipo beta y beta modificado, así como el flujo monofásico de líquidos y gases y el tratamiento del flujo de alta velocidad. Se revisa el procedimiento para definir las condiciones iniciales, con base en los conceptos de equilibrio gravitacional y capilar así como las diversas condiciones de frontera de interés. Se revisan también algunos conceptos relativos al flujo multifásico, como son: el cálculo de permeabilidades relativas en sistemas bifásicos y trifásicos.

2.1 Ecuaciones de Flujo Multifásico Composicional. En esta sección se verá la formulación de problemas de flujo multifásico composicional.

2.1.1 Consideraciones y Antecedentes. Considérese el flujo multifásico isotérmico de aceite, gas y agua en un yacimiento. Las fases aceite y gas forman una mezcla multicomponente, preponderantemente de hidrocarburos, véase la tabla mostrada a continuación, en equilibrio termodinámico en el yacimiento, en todo punto del dominio en espacio y en tiempo. Considérese también que la fase agua no intercambia su masa con las otras fases. Tabla 2.1 Distribución de componentes en las Fases.

COMPONENTE C1 C2 C3 C4

FASE Aceite

Gas

Cnc Agua

Agua

Definición de Conceptos. La composición total de la mezcla de hidrocarburos, a una presión y temperatura dadas, es determinada por las fracciones molares definidas por: z1, z2, ..., znc, de sus componentes, es decir: Moles mM zm = m = 1, 2, ..., nc …(2.1) Moles M Donde MolesmM son los moles del componente m en la mezcla, y Moles M = ∑ J =1 Moles jM . nc

14

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Puesto que los moles del componente m en la mezcla son la suma de los moles en las fases aceite y gas: MolesmM = Molesmo + Molesmg , se puede entonces escribir la ecuación anterior como:

zm =

Moles mo + Moles mg Moles M

...(2.2)

Empleando las siguientes definiciones de las fracciones molares de los componentes en las fases aceite y gas, x m y y m , y las fracciones molares de las fases aceite y el gas en la mezcla, L y V, es decir: Moles mo xm = ...(2.3) Moleso Moles mg ym = ...(2.4) Moles g

L=

V =

Moleso Moles M Moles g

Moles M

...(2.5) ...(2.6)

donde: nc

Moles o = ∑ Moles jo

...(2.7)

j =1

nc

Moles g = ∑ Moles jg

...(2.8)

j =1

Las Ecs. 2.7 y 2.8 son los moles en las fases aceite y gas, respectivamente. Se puede entonces escribir la Ec. 2.2 como:

z m = x m L + y mV

...(2.9)

Nótese que son aplicables las siguientes ecuaciones de restricción. nc

∑ xm = 1

...(2.10)

∑ ym = 1

...(2.11)

∑ zm = 1

...(2.12)

L +V =1

...(2.13)

m =1 nc

m =1 nc

m =1

Y

15

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

2.1.2 Ecuaciones de Flujo. Las ecuaciones que describen el flujo multifásico composicional en un yacimiento, en condiciones isotérmicas, conforme a las consideraciones previamente establecidas, se obtienen a partir de:

a)

Ecuaciones de conservación de materia: una para cada uno de los componentes que constituyen la mezcla de hidrocarburos y una más para la fase agua.

b)

Ecuaciones de movimiento de las fases que fluyen en el medio poroso.

c)

Ecuaciones de equilibrio termodinámico entre las fases aceite y gas.

d)

Relaciones de capilaridad, que establecen la relación entre las presiones de las fases.

e)

Ecuaciones de restricción.

f)

Es necesario además tener expresiones, o correlaciones, para el cálculo de propiedades físicas de las fases: densidad, ρ p , viscosidad, μ p , y tensiones interfaciales, σ p , p = o, g , w , en función de su composición, de su presión, p, y de la temperatura del yacimiento, Ty.

g)

Así como también de ecuaciones, o datos experimentales, para las permeabilidades relativas Kr's.

Con el fin de obtener las ecuaciones de flujo multifásico composicional, considérese un volumen elemental representativo del medio poroso, donde existe flujo unidimensional, ver Fig. 2.1. Con base en el principio de conservación de materia, se puede establecer el siguiente balance molar para el componente m, m=1,2,…,nc, en el elemento:

A φ, k

(xm ρˆ oν ox )x

(y

m

(xm ρˆ oν ox )x + Δx

ρˆ gν gx )x

(y

(ρˆ wν wx )x

m

ρˆ gν gx )x + Δx

(ρˆ wν wx )x + Δx

x

x + Δx

Fig. 2.1 Volumen de control representativo del medio poroso

16

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Ritmo de Entrada de Moles de m al Elemento, a través de x Ritmo de Salida de Moles de m al Elemento, a través de x + Δx Ritmo de Producción / Inyección de Moles de m en el Elemento Ritmo de Acumulación de Moles de m en el Elemento.

+ = ...(2.14)

Ahora bien, asumiendo que cada uno de los componentes de la mezcla de hidrocarburos están presentes simultáneamente en las fases aceite y gas, se tiene que los ritmos de entrada y salida de moles del componente m en el elemento son: Ritmo de Entrada de Moles de m = Ritmo de Entrada de Moles de m Contenidos en la Fase Aceite + Ritmo de Entrada de Moles de m Contenidos en la Fase Gas =

( Axm ρˆ oν ox )x + (Aym ρˆ gν gx )x

...(2.15)

y Ritmo de Salida de Moles de m del Elemento =

( Axm ρˆ oν ox )x+Δx + (Ay m ρˆ gν gx )x+Δx

...(2.16)

donde: A es el área transversal del elemento, expuesta al flujo. ρˆ o y ρˆ g son las densidades molares de las fases aceite y gas,

ν o y ν g son las velocidades microscópicas de las fases aceite, o, y gas, g, en la dirección x. El ritmo de acumulación de moles del componente m en el elemento, está dado por el ritmo de cambio de los moles de m, contenidos en las fases aceite y gas en el espacio poroso, esto es: Ritmo de Acumulación de Moles de m en el Elemento = ∂ ⎡∂ ⎤ AΔx ⎢ (φS o ρˆ o xm ) + (φS g ρˆ g y m )⎥ ∂t ⎣ ∂t ⎦

...(2.17)

El ritmo de producción/inyección de moles de m del elemento, definido como: q~ p* = gasto volumétrico de la fase p = o, g a condiciones de yacimiento por unidad de volumen de roca, es, Ritmo de Prod. / Iny. de Moles de m = AΔx(xm ρˆ o q~o* + ym ρˆ g q~g* )

...(2.18)

Substituyendo las Ecs. 2.15 a 2.18 en el postulado de conservación de masa, Ec. 2.14 y dividiendo la expresión resultante entre el volumen del elemento de control, AΔx, rearreglando y tomando el límite cuando Δx → 0, se llega a:

−

∂ (x m ρˆ oν ox ) − ∂ ( y m ρˆ gν gx ) + (x m ρˆ o q~o* + y m ρˆ g q~g* ) = ∂x ∂x ∂ [φ (xm ρˆ o So + ym ρˆ g S g )] ∂t

...(2.19)

17

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

La ecuación anterior puede extenderse al caso de flujo multidimensional, por lo que la Ec. 2.19 se escribe en su forma más general, en términos del operador diferencial nabla, ∇, como:

− ∇ ⋅ [xm ρˆ oν o + ym ρˆ gν g ] + (xm ρˆ o q~o* + ym ρˆ g q~g* ) = ∂ [φ (xm ρˆ o So + ym ρˆ g S g )] ∂t

...(2.20)

Un balance molar para el agua conduce similarmente a la siguiente ecuación de conservación: ∂ ...(2.21) − ∇ ⋅ [ρˆ oν w ] + ρˆ w q~w* = (φρˆ w S w ) ∂t Nótese que el operador diferencial ∇ se define en el sistema de coordenadas cartesiana como: ∂ ∂ ∂ ∇= i+ j+ k ...(2.22) ∂x ∂y ∂z También, el vector velocidad de la fase p = o, g , w,ν p , se define como:

ν p = ν px i + ν py j + ν pz k

...(2.23)

Suponiendo flujo Darciano en el yacimiento, se puede expresar la velocidad de las fases como:

νp =−

kkr p

μp

(∇p p − γ p ∇D )

...(2.24)

donde, γ p es el peso específico de la fase p, definido como: γ p = ρ p g gc , D es la profundidad con respecto a un plano de referencia, y k es, en el sentido más general, un tensor de permeabilidades absoluta. Substituyendo la Ec. 2.24 en las Ecs. 2.20 y 2.21 se obtienen:

⎡ ⎤ kkrg kkr (∇pg − γ g ∇D )⎥ + ∇ ⋅ ⎢ xm ρˆ o o (∇po − γ o∇D ) + ym ρˆ g μo μg ⎣ ⎦ (xm ρˆ o q~o* + ym ρˆ g q~g* ) = ∂ [φ (xm ρˆ o So + ym ρˆ g S g )] ∂t

...(2.25)

m = 1, 2,…, nc y

⎡ kkr ⎤ ∂ ∇ ⋅ ⎢ ρˆ w w (∇pw − γ w∇D )⎥ + ρˆ w q~w* = (φρˆ w S w ) μw ∂t ⎣ ⎦

...(2.26)

Ahora bien, la condición de equilibrio termodinámico entre las fases gas y aceite se estipula a través de la igualdad de las fugacidades de los componentes en las fases, es decir:

18

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

f mo = f mg

m= 1, 2,…, nc

...(2.27)

Siendo la fugacidad del componente m en la fase p = o, g , función de la presión, de la temperatura y de la composición de la fase; esto es: f mo = f mo po , T y , x1 , x 2 ,..., x nc ...(2.28)

(

)

(

f mg = f mg p g , T y , y1 , y 2 ,..., y nc

)

...(2.29)

Se tienen además las siguientes relaciones y ecuaciones de restricción: Pc go S g = p g − po

…(2.30)

Pc wo (S w ) = po − p w

…(2.31)

So + S g + S w = 1

…(2.32)

( )

nc

∑ xm = 1

…(2.33)

∑ ym = 1

…(2.34)

m =1 nc m =1

Se puede verificar que las Ecs. 2.25, 2.26, 2.27 y 2.30 a 2.34 definen un conjunto de 2nc + 6 ecuaciones con 2nc + 6 incógnitas ó variables primarias, que son: po , p g , p w , S o , S g , S w , x1 , x 2 ,..., x nc , y1 , y 2 ,..., y nc Como puede apreciarse, el tratamiento composicional del flujo multifásico involucra un número grande de ecuaciones y de incógnitas, que crecen con el número de componentes a considerar. Consecuentemente, el estudio del comportamiento de yacimientos composicionales requiere de mayores recursos computacionales y por ende de mayores recursos económicos. El tratamiento composicional es necesario cuando la composición de las fases varía importantemente durante la vida del yacimiento, como es el caso de yacimiento de aceite, volátil y de gas y condensado, y en procesos de recuperación mejorada. Coats ha demostrado que el agotamiento natural de yacimientos composicionales puede describirse aproximadamente a través de formulaciones pseudo-composicionales más simples, como se verá a continuación.

Ec. de Conservación de la materia Fugacidad Capilaridad Restricción

nc + 1 nc 2 3 __________ Total = 2nc + 6

19

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

2.2 Flujo Multifásico Pseudocomposicional: Fluidos tipo Beta y Beta-Modificado. En yacimientos de aceites ordinarios, ó del tipo beta, la composición de las fases aceite y gas no experimentan cambios importantes durante la vida del yacimiento. El comportamiento termodinámico de las fases y sus propiedades físicas pueden consecuentemente ser descritas únicamente en términos de la presión, a la temperatura correspondiente del yacimiento. Esto es, a través de Bo, Bg, y Rs. Desde el punto de vista composicional, el gas y el aceite obtenidos en la superficie, después de los separadores, se consideran ser los pseudocomponentes que constituyen a las fases aceite y gas del yacimiento. EI pseudocomponente aceite esta presente solo en la fase aceite, mientras que el pseudocomponente gas esta presente tanto en la fase aceite como en la fase gas. El tratamiento del comportamiento de fases de fluidos tipo beta se ha extendido al tratamiento aproximado de fluidos composicionales. Se ha mostrado en la literatura que el comportamiento de fase de un aceite volátil y de un gas y condensado, en estudios de agotamiento natural, puede ser descrito aproximadamente extendiendo las ideas aplicadas a fluidos tipo beta. Esto es, a través de dos pseudocomponentes hidrocarburos y del uso de factores de volumen y de relaciones de solubilidad, funciones que dependen solo de la presión de las fases. Los pseudocomponentes aceite y gas, obtenidos en la superficie después de una ó más etapas de separación, se consideran presentes tanto en la fase gas como en la fase aceite. Este tratamiento es conocido en la literatura como tipo beta-modificado. Coats presenta una metodología para obtener los pseudocomponentes y para el cálculo de las propiedades termodinámicas de las fases.

2.2.1 Flujo de Fluidos Tipo Beta-Modificado. La composición de fluidos tipo beta-modificado, conforme a lo expresado anteriormente, se muestra en la Tabla 2.2. Tabla 2.2 Distribución de componentes en las Fases.

COMPONENTE

FASE

Aceite del separador

Aceite

Gas del separador

Gas

Agua

Agua

Las ecuaciones de flujo multifásico de fluidos tipo beta-modificado pueden obtenerse como un caso particular de las ecuaciones generales del flujo multifásico composicional. La Ec. 2.25, para los componentes aceite y gas, es por lo tanto: Para el componente aceite, m = 0, ⎡ ⎤ kkr kkro (∇po − γ o ∇D ) + y o ρˆ g g ∇p g − γ g ∇D ⎥ + ∇ ⋅ ⎢ xo ρˆ o μo μg ⎢⎣ ⎥⎦ ...(2.35) ∂ xo ρˆ o q~o* + y o ρˆ g q~g* = φ xo ρˆ o S o + y o ρˆ g S g ∂t

(

(

)

[(

)

)]

20

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

y para el componente gas, m = g, ⎡ kkr kkro (∇po − γ o ∇D ) + y g ρˆ g g ∇p g − γ g ∇D ∇ ⋅ ⎢ x g ρˆ o μg μo ⎢⎣ ∂ x g ρˆ o q~o* + y g ρˆ g q~g* = φ x g ρˆ o S o + y g ρˆ g S g ∂t

(

(

)

[(

⎤

)⎥ + ⎥⎦

)]

...(2.36)

La ecuación de flujo de agua se mantiene igual, esta es:

⎡ kkrw ⎤ (∇p w − γ w ∇D )⎥ + ρˆ w q~w* = ∂ (φρˆ w S w ) ∇ ⋅ ⎢ ρˆ w μw ∂t ⎣ ⎦

...(2.37)

La condición de equilibrio termodinámico entre las fases puede ser alternamente expresada a través de las constantes de equilibrio de los componentes en las fases, es decir:

(

)

(

)

y k o po , p g , T y , xo , y o = o xo

...(2.38)

y,

k g po , p g , T y , x g , y g =

yg xg

Ahora, las ecuaciones de restricción son: xo + x g = 1

yo + y g = 1

...(2.39)

...(2.40) ...(2.41)

y

So + S g + S w = 1 Las relaciones de presión capilar están dadas por las Ecs. 2.30 y 2.31: Pc go S g = p g − p o

( )

Pc wo (S w ) = p o − p w

...(2.42)

...(2.43) ...(2.44)

Se puede verificar que las Ecs. 2.35 a 2.44 constituyen un conjunto de 10 ecuaciones con 10 incógnitas: po , p g , p w , S o , S g , S w , xo , x g , y o , y g . Se puede sin embargo, expresar las fracciones molares de los componentes; xo , x g , y o y y g , en términos de propiedades volumétricas, funciones de po y p g , y consecuentemente reducir el número de ecuaciones y de incógnitas, como se muestra a continuación:

21

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Empleando la definición de xo , ver Ec. 2.3, y utilizando la relación entre moles y masa, se tiene que:

xm =

Molesoo M o moo = Moleso M o mo

...(2.45)

Donde: Moles oo = Moles de aceite en la fase aceite,

M o = Peso molecular del componente aceite, M o = Peso molecular de la fase aceite, moo = Masa del componente aceite en la fase aceite y

mo = moo + m go = Masa total de la fase aceite.

A partir de la definición de densidad másica ρ = m /V , y de su relación con la densidad molar, ρˆ = ρ / M , así como de la definición de Bo , se puede re-escribir la Ec. 2.45 de la siguiente forma: ρˆ o, cs M o ρ o, cs Voo, cs ...(2.46) xo = = M o ρ o Vo, cy ρˆ o Bo Similarmente se tiene que: Moles go M o ρ g , cs V go, cs Voo, cs ρˆ g , cs Rs = = xg = Moles o M g ρ o Voo, cs Vo, cy ρˆ o Bo

yg =

yo =

Moles gg Moles g Moles og Moles g

=

=

ρˆ g , cs M g ρ g , cs V gg , cs = M g ρ g V g , cy ρˆ g B g

M g ρ o,cs Vog ,cs V gg ,cs ρˆ o,cs rs = M o ρ g V gg ,cs V g ,cy ρˆ g B g

...(2.47)

...(2.48)

...(2.49)

Donde: rs = Relación de solubilidad del aceite en el gas, concepto definido como la Relación entre los volúmenes de los componentes aceite y gas, medidos a condiciones de superficie, obtenidos al llevar a la superficie un cierto volumen de la fase gas del yacimiento, rs =

V og , cs

...(2.50)

V gg , cs

Este concepto adquiere sentido cuando se trata de representar el fenómeno de condensación retrograda en yacimientos de gas y condensado ó en el caso de representar el comportamiento de fase de yacimientos de aceite volátil. Substituyendo las Ecs. 2.42 a 2.45 en las Ecs. 2.31 a 2.40 se obtiene:

22

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Para el componente aceite, ⎡ kkro ⎤ kkrg (∇p o − γ o ∇D ) + ∇⋅⎢ rs ∇p g − γ g ∇D ⎥ + μ g Bg ⎢⎣ μ o Bo ⎦⎥

(

q o* + rs q *g =

)

...(2.51)

∂ ⎡ ⎛⎜ S o S g ⎞⎟⎤ ⎢φ rs ⎥ + ∂t ⎢ ⎜⎝ Bo B g ⎟⎠⎥ ⎣ ⎦

Para el componente gas, ⎤ ⎡ kkro kkrg R s (∇p o − γ o ∇D ) + ∇⋅⎢ ∇p g − γ g ∇D ⎥ + μ g Bg ⎢⎣ μ o Bo ⎥⎦ S g ⎞⎤ ∂ ⎡ ⎛S ⎟⎥ R s q o* + q *g = ⎢φ ⎜ o R s + B g ⎟⎠⎥ ∂t ⎢ ⎜⎝ Bo ⎣ ⎦

...(2.52)

La ecuación del agua se expresa en términos volumétricos como sigue: ⎡ kkrw ⎤ ⎛ S ⎞ (∇p w − γ w ∇D )⎥ + q *w = ∂ ⎜⎜φ w ⎟⎟ ∇⋅⎢ ∂t ⎝ B w ⎠ ⎣ μ w Bw ⎦

...(2.53)

(

)

Donde q o* , q *g y q *w son los gastos a condiciones de superficie (c.s.) de aceite, gas y agua por unidad de volumen de roca. Las Ecs. 2.51 a 2.53 juntamente con las Ecs. 2.42 a 2.44 definen un conjunto de 6 ecuaciones con 6 incógnitas: po , p g , p w , S o , S g y S w . Se puede ahora, establecer la relación entre las densidades de las fases, ρ p , y sus propiedades volumétricas. Para esto se emplean las restricciones estipuladas por las Ecs. 2.40 a 2.41, es decir: ρˆ o, c.s. + ρˆ g , c.s. Rs ...(2.54) xo + x g = 1 = ρˆ o Bo Lo que conduce a: ρˆ o, c.s. + ρˆ g , c.s. Rs ...(2.55) ρˆ o = Bo

Similarmente, yo + y g = 1 =

ρˆ o,c.s. rs + ρˆ g ,c.s. ρˆ g B g

y para la densidad de la fase gas se tiene: ρˆ o, c.s.rs + ρˆ g , c.s. ρˆ g = Bg

...(2.56)

...(2.57)

23

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Nótese que las constantes de equilibrio se expresan en términos de las relaciones de solubilidad Rs y rs de la siguiente manera. Substituyendo las Ecs. 2.46 y 2.49 en la Ec. 2.38, se tiene: ρˆ B ...(2.58) K o = o o rs ρˆ g B g Ahora bien, substituyendo las Ecs. 2.55 y 2.57 en la Ec. 2.58 y rearreglando, se llega a:

Ko =

1 + ρˆ g , c.s. Rs ρˆ o, c.s. 1 + ρˆ g , c.s. ρˆ o, c.s. rs

...(2.59)

Similarmente para la fase gas, se tiene que: 1 + ρˆ o, c.s. ρˆ g , c.s. Rs 1 Kg = 1 + ρˆ o, c.s. rs ρˆ g , c.s. Rs

...(2.60)

2.2.2 Flujo de Fluidos Tipo Beta. La composición de fluidos tipo beta, conforme a lo expresado anteriormente, se muestra en la tabla 2.3. Tabla 2.3 Distribución de componentes en las Fases.

COMPONENTE

FASE

Aceite del separador

Aceite

Gas del separador

Gas

Agua

Agua

Nótese que al no existir el componente aceite en la fase gas, la relación de solubilidad del aceite en el gas es cero, rs = 0 y las ecuaciones de flujo multifásico de fluidos tipo beta se obtienen como un caso particular de las ecuaciones anteriores. Estas son: Para el componente aceite, ⎡ ⎛ ⎡ kkro ⎤ ⎞⎤ (∇po − γ o∇D )⎥ + qo* = ∂ ⎢φ ⎜⎜ So ⎟⎟⎥ ∇⋅⎢ ∂t ⎣ ⎝ Bo ⎠⎦ ⎣ μo Bo ⎦

...(2.61)

Para el componente gas, ⎤ ⎡ kkro kkrg Rs (∇po − γ o ∇D ) + ∇⋅⎢ ∇p g − γ g ∇D ⎥ + μ g Bg ⎥⎦ ⎣⎢ μ o Bo

(

S g ⎞⎤ ∂ ⎡ ⎛S ⎟⎥ Rs qo* + q*g = ⎢φ ⎜ o Rs + B g ⎟⎠⎥ ∂t ⎢ ⎜⎝ Bo ⎣ ⎦

)

...(2.62)

24

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

La ecuación del agua:

⎡ kkrw ⎤ ⎛ ⎞ (∇pw − γ w∇D )⎥ + qw* = ∂ ⎜⎜φ S w ⎟⎟ ∇⋅⎢ ∂t ⎝ Bw ⎠ ⎣ μ w Bw ⎦

...(2.63)

Con las siguientes relaciones adicionales:

So + S g + S w = 1

( )

…(2.64)

Pc go S g = p g − po

…(2.65)

Pc wo (S w ) = po − p w

…(2.66)

Como en el caso de aceites tipo beta-modificado, se puede apreciar que las Ecs. 2.61 a 2.66 forman un conjunto de 6 ecuaciones con 6 incógnitas: po , p g , pw , S o , S g y S w .

2.3 Flujo Monofásico. 2.3.1 Flujo de Aceite. Se considera ahora un yacimiento bajosaturado donde existe flujo monofásico de aceite, esto es, en presencia de una saturación inmóvil de la fase agua. La fase agua se considera compresible y su saturación, por lo tanto, esta sujeta a cambios debido a las variaciones de presión. Las ecuaciones que describen el flujo de aceite, en estas condiciones, se obtienen a partir de las ecuaciones presentadas en la sección anterior, estas son: Para el componente aceite, ⎡ k ⎤ ∂ ⎡ ⎛ S ⎞⎤ ∇ ⋅ ⎢ oc (∇po − γ o∇D )⎥ + qo* = ⎢φ ⎜⎜ o ⎟⎟⎥ ∂t ⎣ ⎝ Bo ⎠⎦ ⎣ μ o Bo ⎦

...(2.67)

Nótese que k oc es la permeabilidad efectiva de la fase aceite, medida a la saturación

irreducible de agua, koc = kkro (S wc ). Para la fase agua se tiene:

∂ ⎡ ⎛ S wc ⎞⎤ ⎟⎥ = 0 ⎢φ ⎜ ∂t ⎣ ⎜⎝ Bw ⎟⎠⎦ Se tienen las siguientes relaciones adicionales: S o + S wc = 1

p w = po − Pcwo (S w )

…(2.68)

...(2.69) …(2.70)

25

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

2.3.2 Flujo de un Aceite de Compresibilidad Pequeña y Constante Considérese el flujo monofásico de un aceite de compresibilidad pequeña y constante en un medio poroso homogéneo e isotrópico. Supóngase que la viscosidad de: aceite es aproximadamente constante y que al agua y la roca son incompresibles. Los efectos gravitacionales son también ignorados. Tomando en cuenta las consideraciones anteriores en la Ec. 2.61, haciendo p = p o y empleando la porosidad de hidrocarburos, φ hc = φSo , se tiene que: Ecuación de Flujo:

⎡ 1 ⎤ μ φ μ ∂⎛ 1 ⎞ ⎟ ∇ ⋅ ⎢ ∇p ⎥ + o qo∗ = hc o ⎜⎜ k oc ∂t ⎝ Bo ⎟⎠ ⎣ Bo ⎦ k oc Ahora bien, la compresibilidad de un aceite bajosaturado se define como: 1 dBo co = − po dp

…(2.71)

…(2.72)

Y para un aceite de compresibilidad pequeña y constante se tiene, al integrar la Ec. 2.72, la siguiente expresión para el factor de volumen: Bo = Bob exp[co ( pb − p )] …(2.73) Introduciendo la Ec. 2.73 en la Ec. 2.71 y operando, se obtiene: μ φ μ c ∂p ∇ 2 p + co (∇p )2 + o qo∗ Bo = hc o o k oc k oc ∂t

…(2.74)

Considerando que los gradientes de presión son pequeños, de manera que co (∇p )2 ≈ 0, se tiene que: μ φ μ c ∂p ∇ 2 p + o qo∗ Bo = hc o o …(2.75) k oc k oc ∂t

2.3.3 Flujo de Gas. Considérese el flujo monofásico de un gas en un medio poroso homogéneo e isotrópico. Supóngase nuevamente que el agua y la roca son incompresibles y que los efectos gravitacionales en el flujo son ignorados.

26

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Las ecuaciones de flujo multifásico se reducen en este caso a: ⎡ 1 ⎤ 1 ∗ φhc ∂ ⎛ 1 ⎜ qg = ∇⋅⎢ ∇p ⎥ + ⎜ Bg μ B k k t ∂ gc gc ⎝ ⎣⎢ g g ⎦⎥

⎞ ⎟ ⎟ ⎠

…(2.76)

Donde p = p g y φ hc = φ S g = φ (1 − S wc ) Para gases reales se tiene la siguiente ecuación de estado:

pV = znRT

…(2.77)

Esta ecuación puede usarse en la definición del factor de volumen del gas, para dar la siguiente relación: Vg ,cy T y p s z = …(2.78) Bg = Vg ,cs Ts p Substituyendo la Ec. 2.78 en la Ec. 2.76 y rearreglando:

⎡ p ⎤ pT φ ∂ ⎛ p⎞ ∇⋅⎢ ∇p ⎥ + s y q ∗g = hc ⎜ ⎟ k gc ∂t ⎝ ' z ⎠ ⎢⎣ μ g z ⎥⎦ Ts k gc

…(2.79)

La Ec. 2.79 se reduce a una forma casi lineal empleando la transformación de pseudopresión, o el potencial de gas real, de Al-Hussainy y Colaboradores. Para esto considérese: 1 dBg z d ⎛ p ⎞ …(2.80) cg = − = ⎜ ⎟ Bg dp p dp ⎝ z ⎠ Usando la regla de la cadena en la derivada del tiempo de la Ec. 2.79 y considerando la Ec. 2.80, se puede reescribir la Ec. 2.79 como: ⎡ p ⎤ pT φ μ c p ∂p ∇⋅⎢ ∇p ⎥ + s y q ∗g = hc g g …(2.81) k gc μ g z ∂t ⎢⎣ μ g z ⎥⎦ Ts k gc Al-Hussainy y colaboradores, definieron el potencial de gas real, m( p ), como:

m( p ) = 2

p

∫

pref

p' dp' μg z

Introduciendo la Ec. 2.82 en la Ec. 2.81 se obtiene: 2p T φ μ c ∂m( p ) ∇ 2 m( p ) + s y q g∗ = hc g g Ts k gc k gc ∂t

…(2.82)

…(2.83)

La Ec. 2.83 es aún no lineal ya que μ g y c g dependen de la presión.

27

FI-UNAM

2.4

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Condiciones Iniciales y de Frontera.

La formulación de un problema de flujo de fluidos no esta completa si no se definen las condiciones iniciales y de frontera. Las condiciones iniciales definen en el caso más general, de flujo multifásico composicional, la distribución de presiones, de saturaciones y la composición de las fases al tiempo cero de nuestro problema. Por otra parte, las condiciones de frontera introducen la información correspondiente a la forma en que el yacimiento interactúa con sus alrededores durante su vida productiva.

2.4.1 Condiciones Iniciales. Considérese que al tiempo cero, t = 0, existe equilibrio gravitacional y capilar en el yacimiento. Supóngase también que se conoce la presión po, ref , de la fase aceite, a una cierta profundidad del yacimiento, z ref y que así como la profundidad del contacto agua-aceite, zcwo . En el caso de un yacimiento inicialmente saturado, donde se conoce la posición del contacto gas-aceite, zcgo , se puede establecer zref = zcgo , sabiendo que la presión en el contacto gas-aceite corresponde a la presión de burbuja, po, ref = pb . La suposición de equilibrio gravitacional y capilar significa que no existe flujo de las fases al tiempo cero; esto es, los gradientes de potencial de las fases en cualquier punto y en cualquier dirección son cero. Esto es, conforme a la ley de Darcy, Ec. 2.24, se tiene que: kkrp ∇ p p − γ p ∇D = 0 ...(2.84) μ pBp

(

)

A lo largo de un plano horizontal ubicado en una posición cualquiera en z, se tiene que: ∂p p =0 ...(2.85) ∂x

∂p p ∂y

=0

...(2.86)

Las Ecs. 2.85 y 2.86 indican que en condiciones de equilibrio, las presiones, p p en el plano horizontal (x, y) ubicado en z es constante. En estas condiciones, la Ec. 2.84 también indica que: ∂p p −γ p = 0 ...(2.87) ∂z Lo que significa que en condiciones de equilibrio, la distribución vertical de presiones, p p ( z ), esta dado por el peso de la columna de los fluidos p, o sea:

28

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

p p ( x, y, z , t = 0) = p p 0 ( z )

...(2.88)

Conocida la presión, po, ref , medida a z ref , es posible calcular la presión po en cualquier profundidad z dentro del yacimiento. Integrando la Ec. 2.87 se obtiene lo siguiente:

gc g

dp 'p

pp

∫ p

p , ref

ρp

=

z

dz ∫ z

...(2.89)

ref

Lo que da como resultado,

p p = p p, ref +

g ρ ( z − zref ) gc p

...(2.90)

Donde:

ρ p = ρp(p p) Y p p puede aproximarse como el promedio aritmético de las presiones

...(2.91)

p p y pref . Se

puede ver que la Ec. 2.90 contiene implícitamente a la presión p p , ya que ρ p depende de esa variable. La solución de p p de la Ec. 2.90 puede obtenerse a través del método de NewtonRaphson, por ejemplo: Conocida po,ref , a zref , se calcula po en cualquier z del yacimiento a partir de la Ec. 2.90, escrita para la fase aceite, esto es: g po = po, ref + ρ ( z − zref ) ...(2.92) gc o Conocida la presión en el contacto agua-aceite, po, zcwo se puede calcular la presión de la fase agua en el contacto, pw, zcwo . Dado que:

Pcwo ( S w = 1.0) = Pcwo, e = po, zcwo − pw, zcwo

...(2.93)

Donde Pcwo, e es la presión de entrada de los poros. De la Ec. 2.92 se puede obtener la presión de la fase agua en el contacto agua-aceite, pw, zcwo = po, zcwo − Pcwo, e ...(2.94) Similarmente, cuando existe capa de gas en el yacimiento, se puede obtener la presión de la fase gas en el contacto gas-aceite mediante: p g , zcgo = Pc go, e − po, zcgo ...(2.95) La distribución de presiones de las fases gas, p g (z ), y agua, pw (z ), pueden entonces calcularse a partir de la Ec. 2.90. Es decir:

29

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

g ρ ( z − z zcgo ) ...(2.96) gc g g pw = pw, zcwo + ρ ( z − z zcwo ) ...(2.97) gc w Conocida la distribución vertical de presiones de las fases en el yacimiento, po (z ), p g (z ), p g = p g , zcgo +

y pw (z ), es posible calcular la distribución de presiones capilares. Esto es; Y

Pc go ( S g ( z )) = p g ( z ) − po ( z )

...(2.98)

Pcwo ( S w ( z )) = po ( z ) − pw ( z )

...(2.98)

Mediante un proceso de interpolación inversa, se puede obtener la distribución inicial de saturaciones de las fases gas y agua, S g (z ), y S w (z ). La saturación de aceite se obtiene a partir de la Ec. 2.64, S o ( z ) = 1 − S g ( z ) − S w ( z ). Nótese que dentro de la zona de aceite existe una cierta posición z arriba de la cual la saturación de agua, S w , se toma irreducible o inmóvil. A partir de este punto y en toda la región donde esto se cumpla, se tiene que el gradiente de presión de la fase agua será el correspondiente a la fase aceite. Un criterio similar se aplica a la región dentro de la capa de gas donde la fase gas se toma inmóvil: el gradiente de presión de la fase gas en esa región será el correspondiente a la fase aceite. La Fig. 2.1 muestra esquemáticamente esta inicialización en equilibrio gravitacional y capilar del dominio de interés.

Z ref

1 − S or P r o f u n d i d a d

Sg

Pg

Cima Gas

Pg = Po − Pcgoe

S g (0% ) = 0Pcgo = Pcgoe

Z cgo

Aceite

S wc Po Pcwo (S w ) = Po − Pw

Sw Pw

S w (100% ) = 1

Z

Z cwo

Pcwo = PcwoeAgua ∴ Pw = Po − Pcwoe Base

Pref Fig. 2.1 Condiciones Iniciales en Equilibrio Gravitacional y Capilar

30

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Para el caso de fluidos composicionales; se requiere también especificar la composición inicial de las fases aceite y gas. Esta información deberá conocerse del análisis en laboratorio de muestras representativas de los fluidos, obtenidas antes de iniciar la explotación del yacimiento.

2.4.2 Condiciones de Frontera. Básicamente existen dos tipos de condiciones de frontera; de interés en problemas de flujo de fluidos en medios porosos: i. Gasto especificado en la frontera y, ii. Presión especificada. 1) Gasto Especificado.

De la ecuación de Darcy, escrita en la frontera para la fase p , se tiene:

(q p ) ⊥ Front = −

kkr p A ⎛ ∂p p ∂D ⎞ ⎟ ⎜ γ − p μ p B p ⎜⎝ ∂η ∂η ⎟⎠ ⊥ Front

...(2.100)

Donde η es dirección de flujo, perpendicular a la frontera. Entonces:

⎛ μ p Bpq p ⎛ ∂p p ∂D ⎞ ⎜ ⎜ ⎟ − γ = − p ⎜ ∂η ⎟ ⎜ kkr p A ∂ η ⎝ ⎠ ⊥ Front ⎝

⎞ ⎟ ⎟ ⎠ ⊥ Front

...(2.101)

2) Presión Especificada.

En este caso se especifica la presión en la frontera, por ejemplo de la fase aceite, a lo largo del tiempo. Es decir: po (Frontera, t ) = po (t )⊥ Front ...(2.102)

3. Aproximación de Ecuaciones Diferenciales en Diferencias Finitas. 3.1 Diferencias Finitas. Suponer que para cada x en el intervalo (a, b) existen la función f(x); y hasta la k esima derivada de la función, o sea: f ( x ), f ' (x ), f ' ' ( x ),..., f ( k ) ( x ). Entonces, la expansión de la función f (x ) al rededor del punto xi contenido en el intervalo será:

( x − xi ) df f ( x ) = f (xi ) + 1! dx

xi

( x − xi ) 2 d 2 f + 2! dx 2

Donde, el error de truncamiento es: ξ = xi + θ ( x − xi )

y

xi

( x − xξ ) k d k f + ... + k! dx k

...(3.1) ξ

0 <θ <1

31

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Nótese que la expansión dada por la Ec. 3.1 contiene k términos solamente y aun es exacta. Reemplazamos la acostumbrada serie infinita de Taylor por una serie finita. Esto es posible ya que como se puede apreciar, el k esimo término esta evaluado no en el punto xi sino en el punto ξi contenido en el intervalo ( x, xi ), que se desconoce, pero sin embargo hace posible la igualdad. La Ec. 3.1 sirve de base en la aproximación de las derivadas que constituyen las ecuaciones de flujo de fluidos en medios porosos que nos ocupan, como se verá a continuación. 3.2 Aproximaciones a la primera derivada. 3.2.1 Diferencias Progresivas.

Considerar en la Ec. 3.1: k = 2 y x = xi + Δx. Esto es:

f ( xi + Δx ) = f ( xi ) + Δx

df dx

+ xi

Δx 2 d 2 f 2! dx 2

...(3.2) ξp

De esta ecuación, se puede obtener la siguiente expresión para la primera derivada: f ( xi + Δx ) − f ( xi ) Δx d 2 f df = − Δx dx xi 2 dx 2 ξ

...(3.3)

p

Nótese que no existe manera de evaluar el último término de la Ec. 3.3, ya que no se tiene información sobre la segunda derivada, f ' ' ( x ), y del punto ξ p donde debe evaluarse; este término acaba despreciándose, y constituye lo que se denomina el error local de truncamiento de la aproximación. Su análisis es sin embargo importante, ya que nos informa sobre el orden de la aproximación. El orden es definido por la potencia de Δx, que multiplica al término despreciado. En este caso la aproximación de la primera derivada de f (x) mediante diferencias progresivas es de primer orden, o sea O(Δx) . Es común reescribir la Ec. 3.3 como:

df dx

= xi

f ( xi + Δx ) − f (xi ) + O p (Δx) Δx

...(3.4)

Siendo O p (Δx) el Error Local de Truncamiento, definido como:

O p ( Δx ) = −

Δx d 2 f 2 dx 2

, para K = 2

...(3.5)

ξp

3.2.2 Diferencias Regresivas. Considerar en la Ec. 3.1: k = 2 y haciendo ahora x = xi − Δx. Con esto, la Ec. 3.2 se escribe como: df Δx 2 d 2 f ...(3.6) f ( xi − Δx ) = f ( xi ) − Δx + dx xi 2! dx 2 ξ r

32

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Y se obtiene la siguiente expresión para la primera derivada: f ( xi ) − f ( xi − Δx ) df = + Or (Δx) dx xi Δx

...(3.7)

Donde, la aproximación de la primera derivada de f(x) mediante diferencias regresivas es también de primer orden. EL error local de truncamiento de la aproximación es: Δx d 2 f , para K = 2 ... (3.8) Or (Δx) = 2 dx 2 ξ r

3.2.3 Diferencias Centrales. Considerar en la Ec. 3.1: k = 3 y escribir la función f(x) en x = xi + Δx y x = xi − Δx , como se muestra a continuación: Δx 2 d 2 f Δx 3 d 3 f df ...(3.9) f ( xi + Δx ) = f ( xi ) + Δx + + dx xi 2! dx 2 x 3! dx 3 ξ p*

i

Y

f ( xi − Δx ) = f (xi ) − Δx

df dx

+ xi

Δx 2 d 2 f 2! dx 2

− xi

Δx 3 d 3 f 3! dx3

...(3.10) ξ

r*

Restando las Ecs. 3.9 y 3.10, se tiene:

df f ( xi + Δx ) − f (xi − Δx ) = 2Δx dx

xi

⎡ Δx 3 ⎢ d 3 f + 3! ⎢ dx 3 ⎣

ξ

p*

⎤ ⎥ ⎥ ξ* r ⎦

d3 f + 3 dx

...(3.11)

La Ec. 3.11 lleva a la siguiente aproximación de la primera derivada mediante diferencias centrales: df f ( xi + Δx ) − f ( xi − Δx ) = + Oc Δx 2 ...(3.12) dx xi 2Δx

( )

( )

Con un error local de truncamiento de segundo orden Oc Δx 2 , es decir:

⎡ Δx 3 ⎢ d 3 f Oc Δx = 3! ⎢ dx 3 ⎣

( ) 2

ξ

p*

d3 f + 3 dx

⎤ ⎥ ⎥ ξ* r ⎦

...(3.13)

( )

Comparando el error local de truncamiento de la aproximación anterior, Oc Δx 2 , con los obtenidos previamente para diferencias progresivas y regresivas, O p (Δx) y Or (Δx) se tiene que:

( )

lim Oc Δx 2 < lim O p (Δx )

Δx → 0

Δx → 0

...(3.14)

Lo que indica que el error de truncamiento de la aproximación en diferencias centrales, para la primera derivada, es menor que el correspondiente a diferencias progresivas o regresivas.

33

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Para abreviar la escritura de las aproximaciones es conveniente usar la notación que a continuación se muestra, en la que se considera que Δx es uniforme. Sea:

xi +1 = xi + Δx, xi −1 = xi − Δx, f i = f ( xi ), f i +1 = f ( xi +1 ), Y

f i −1 = f ( xi −1 ) Empleando esta notación se puede ahora escribir las aproximaciones anteriores de la siguiente manera: Diferencias progresivas,

df dx

= xi

f ( xi + Δx ) − f ( xi ) f −f + O p (Δx) = i +1 i + O p (Δx) Δx Δx

...(3.15)

Diferencias regresivas,

df dx

= xi

f ( xi ) − f ( xi − Δx ) f − f i −1 + Or (Δx) = i + Or (Δx) Δx Δx

...(3.16)

Y diferencias centrales,

df dx

= xi

f (xi + Δx ) − f (xi − Δx ) f −f + Oc Δx 2 = i +1 i −1 + Oc Δx 2 2Δx 2Δx

( )

( )

...(3.17)

34

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

∂ ⎛ ∂p ⎞ ⎜λ ⎟ ∂x ⎝ ∂x ⎠ Suponer que p = p( x, y ) y el coeficiente λ = λ ( p ) . Se desea aproximar el siguiente término diferencial: ∂ ⎛ ∂p ⎞ ...(3.18) ⎜λ ⎟ ∂x ⎝ ∂x ⎠ 3.3 Aproximación de los Términos de la Forma:

Considérese ahora el núcleo típico de celdas mostrado en la Fig. 3.1 para un espacio bidimensional discretizado. Fig. 3.1 Núcleo típico de celdas en Problemas Bidimensionales

j

i

i = 1, 2, …, Nx j = 1, 2, …, Ny k= 1, 2, …, Nz

Frontera = i+1/2, j La cual puede no estar exactamente a la mitad, únicamente sirve para denotar que es frontera.

= i,j-1

Δy j −1 / 2 Δy j

= i-1,j

= i,j

= i+1,j

Δy j +1 / 2 = i,j+1

Δxi Δxi −1 / 2

Δxi +1 / 2

Si se define,

u=λ

∂p ∂x

Se puede entonces escribir la Ec. 3.18 como: ∂ ⎛ ∂p ⎞ ∂u ⎜λ ⎟ = ∂x ⎝ ∂x ⎠ ∂x

...(3.19)

...(3.20)

Empleando diferencias centrales, Ec. 3.17, en la aproximación de la Ec. 3.20 y apoyándose en los puntos i + 1 e i − 1 , que corresponden a las fronteras de la celda i, j en la dirección x, se 2 2 obtiene: u 1 −u 1 i+ , j i− , j ∂u 2 2 ≈ ...(3.21) ∂x i , j Δxi

35

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Sustituyendo u por su definición, Ec. 3.19: ⎡ ⎤ ∂p ∂p ∂u 1 ⎢ ⎥ ≈ λ 1 ...(3.22) −λ 1 ∂x i, j Δxi ⎢ i + 2 , j ∂x i + 1 , j i − 2 , j ∂x i − 1 , j ⎥ ⎢⎣ 2 2 ⎥ ⎦ ∂p en las fronteras de las celdas. Usando de nuevo Se necesita ahora aproximar ∂x diferencias centrales y apoyándose ahora en los nodos i, j , i + 1, j e i − 1, j , se obtienen las siguientes aproximaciones: p i +1, j − p i, j ∂p ...(3.23) ≈ ∂x i + 1 , j Δx 1 i+

2

Y

∂p ∂x

1 i− , j 2

≈

2

p i, j − p i −1, j Δx 1 i−

...(3.24)

2

Substituyendo las Ecs. 3.23 y 3.24 en la Ec. 3.22, se obtiene: ⎡ 1 ⎢⎛ λ ⎞ ∂ ⎛ ∂p ⎞ ⎛ λ ⎞ ≈ p i +1, j − p i, j − ⎜ ⎟ p i, j − pi −1, j ⎜λ ⎟ ⎜ ⎟ ∂x ⎝ ∂x ⎠ i, j Δx i ⎢⎝ Δx ⎠ i + 1 , j ⎝ Δx ⎠ i − 1 , j 2 2 ⎣

(

)

3.4 Aproximación de los Términos de la Forma:

(

⎤ ⎥ ⎥ ⎦

)

...(3.25)

∂p ∂t

Considérese la función p(x,y,t) para la cual se desea obtener aproximaciones en diferencias finitas de su derivada parcial con respecto al tiempo. Se puede para esto aplicar las aproximaciones de diferencias progresivas, regresivas y centrales en tiempo. Antes de mostrar la forma de las aproximaciones; es conveniente resaltar que para fines de representar puntos en el dominio del espacio discretizado, es común el uso de los subíndices i, j, k para las direcciones x, y, z, por ejemplo. Por otra parte; para representar puntos en el tiempo se emplean los subíndices n y n+1, que indican los niveles donde se conoce y se desconoce la solución del problema de interés. Aplicando diferencias regresivas en tiempo, se obtiene:

∂p ∂t

n +1

=

p n +1 − p n i, j

i, j

Δt

i, j

+ O(Δt )

...(3.26)

Empleando diferencias progresivas en tiempo conduce a:

∂p ∂t

n i, j

=

p n +1 − p n i, j

i, j

Δt

+ O(Δt )

...(3.27)

36

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Ahora bien, aplicando diferencias centrales en la aproximación de la derivada resulta en:

∂p ∂t

n+

1 2

=

p n +1 − p n i, j

i, j

Δt

i, j

( )

+ O Δt 2

...(3.28)

Como se vio previamente, la aproximación de la primera derivada mediante diferencias centrales es más exacta que las aproximaciones mediante diferencias regresivas y progresivas. Aunque esto es cierto, el empleo de diferencias centrales en tiempo no es común en la solución de problemas de flujo de fluidos en yacimientos petroleros; esto probablemente se debe a que su aplicación requiere de mayor esfuerzo computacional. Es común, sin embargo, el empleo de diferencias regresivas en tiempo. El empleo de diferencias progresivas se descarta por razones de estabilidad numérica de la aproximación resultante. Como se podrá evidenciar posteriormente, el empleo de cualquiera de estas tres aproximaciones en la discretización de las ecuaciones de flujo, resulta en esquemas numéricos con características totalmente diferentes.

3.5 Notación de las Ecuaciones Aproximadas en Operadores en Diferencias. Con el objeto de abreviar la escritura de las ecuaciones aproximadas, se definen los siguientes operadores de diferencias. Sea el operador de diferencias centrales, Δ, definido conforme a las siguientes propiedades:

Δui, j , k = Δ x ui, j , k + Δ y ui, j , k + Δ z ui, j , k

...(3.29)

donde,

Δ x ui , j , k = u

−u

1 i + , j,k 2

Δ y ui , j , k = u

1 i, j + , k 2

1 i − , j,k 2

−u

1 i, j − , k 2

...(3.30) ...(3.31)

y

Δ z u i , j ,k = u

i , j ,k +

1 2

−u

i , j ,k −

1 2

...(3.32)

Se define también el operador de diferencias regresivas en tiempo, Δ t , como sigue:

Δ t u = u n +1 − u n

...(3.33)

Empleando la notación de estos operadores es posible escribir las aproximaciones mostradas en las Ecs. 3.25 y 3.26 como sigue:

∂ ⎛ ∂p ⎞ ⎟ ⎜ λx ∂x ⎝ ∂x ⎠

≈ i, j

1 ⎞ ⎛λ Δx⎜ x Δx p⎟ Δxi ⎠ i, j ⎝ Δx

...(3.34)

n +1

y

1 ∂p ≈ Δ i pi , j ∂t i, j Δt

...(3.35)

37

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

3.6 Los Esquemas Explícitos, Implícito y Crank-Nicholson de discretización en Tiempo. Se mostrara ahora el uso de los esquemas explícito, implícito y Crank-Nicholson en la aproximación de las ecuaciones diferenciales. Se analizaran las características de las ecuaciones resultantes. Se considera para esto la ecuación diferencial de flujo unidimensional de un fluido de compresibilidad pequeña y constante, Ec. 2.69, esto es:

∂2 p

μB

1 ∂p k η ∂t ∂x 2 Con las siguientes condiciones iniciales y de frontera:

+

q=

...(3.36)

p ( x,0 ) = P 0

...(3.37)

p(0, t ) = P 0

...(3.38)

p ( L, t ) = P 0

...(3.38)

Y

3.6.1 Esquema Explícito. E1 esquema explícito consiste en aproximar la derivada parcial del tiempo de la Ec. 3.36 mediante diferencias progresivas, es decir al nivel de tiempo n. Esto automáticamente define el nivel de tiempo del resto de los términos de la ecuación. Se desea entonces aproximar la Ec. 3.36 en cualquiera de las celdas i, que discretizan el dominio O < x < L, en el nivel de tiempo n, esto es: n

n ⎧⎪ ∂ 2 p ⎫⎪ μB n 1 ⎧ ∂p ⎫ + q = ⎨ 2⎬ ⎨ ⎬ i k η ⎩ ∂t ⎭i ⎪⎩ ∂x ⎪⎭i

...(3.40)

i = 1,2,..., I n = 0,1,2,...

Donde: k φμCT μ = cte.

η=

CT = cte.

φ = uniforme y cte. k = Isotropica. La aproximación del término de flujo mediante diferencias centrales, de la Ec. 3.25 con λ = 1 y Δx = constante, es: n

⎧∂2 p ⎫ pin+1 − 2 pin + pin−1 ...(3.41) ⎨ 2⎬ ≈ 2 Δ x x ∂ ⎩ ⎭i Y la aproximación de la derivada del tiempo mediante diferencias progresivas produce: n 1 ⎧ ∂p ⎫ 1 pin +1 − pin ≈ ⎨ ⎬ η ⎩ ∂t ⎭ i η Δt 2

...(3.42)

38

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Sustituyendo estas aproximaciones, 3.41 y 3.42, en la Ec. 3.40, se tiene: pin+1 − 2 pin + pin−1 μB n 1 pin+1 − pin + qi = Δx 2 k η Δt

...(3.43)

Despejando pin+1 de esta ecuación, se tiene que:

ημB

pin+1 = αpin+1 + (1 − 2α ) pin + αpin−1 +

k

donde:

pin+1 = αpin+1 + (1 − 2α ) pin + αpin−1 +

ημB k

Δtqin

...(3.44)

Δtqin

...(3.45)

Se puede ver que la solución de la Ec. 3.40 mediante un esquema explícito produce una solución explícita para pin+1 .

3.6.2 Esquema Implícito. En un esquema implícito, la derivada parcial del tiempo de la Ec. 3.40 se aproxima en la celda i de la malla de cálculo, mediante diferencias regresivas, en el nivel de tiempo n + 1. Se tiene que: n +1

⎧∂2 p ⎫ ⎨ 2⎬ ⎩ ∂x ⎭i

+

μB k

n +1

q

n +1 i

1 ⎧ ∂p ⎫ = ⎨ ⎬ η ⎩ ∂t ⎭i

...(3.46)

i = 1,2,..., I n = 0,1,2,...

Aproximando el término de flujo mediante diferencias centrales, como en el caso anterior, y el término del tiempo mediante diferencias regresivas, la forma de la Ec. 3.46 en diferencias finitas es: pin++11 − 2 pin+1 + pin−+11 μB n+1 1 pin+1 − pin ...(3.47) + qi = η Δt Δx 2 k Y rearreglando,

αpin−+11 − (2α + 1) pin+1 + αpin++11 = − pin −

ημB

i = 1,2,..., I

k

Δtqin+1

...(3.48)

n = 0,1,2,...

O bien

cpin−+11 − apin+1 + bpin++11 = − d in i = 1,2,..., I n = 0,1,2,...

La aproximación de la Ec. 3.46 mediante un esquema implícito conduce a un sistema algebraico de ecuaciones lineales, cuya solución proporciona la distribución de presiones en el nivel de tiempo n + 1, pin+1 , i = 1,2,..., I .

39

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

3.6.3 Esquema Crank-Nicholson. En el esquema Crank-Nicholson, la derivada parcial del tiempo de la Ec. 3.40 se aproxima mediante diferencias finitas centrales, lo que implica que: n+

⎧⎪ ∂ 2 p ⎫⎪ ⎨ 2⎬ ⎪⎩ ∂x ⎪⎭ i

1 2

+

μB k

n+

qi

1 2

n+

1 ⎧ ∂p ⎫ = ⎨ ⎬ η ⎩ ∂t ⎭ i

1 2

...(3.49)

i = 1,2,..., I n = 0,1,2,...

Donde, n+

⎧⎪ ∂ 2 p ⎫⎪ ⎨ 2⎬ ⎪⎩ ∂x ⎪⎭ i

1 2

1 ⎧⎪ ∂ 2 p ⎫⎪ = ⎨ ⎬ 2 ⎪⎩ ∂x 2 ⎪⎭

n +1

i

1 ⎧⎪ ∂ 2 p ⎫⎪ + ⎨ ⎬ 2 ⎪⎩ ∂x 2 ⎪⎭

n

...(3.50) i

Y

μB k

n+

qi

1 2

=

(

1 μB n +1 n qi + qi 2 k

)

...(3.51)

Substituyendo las Ecs. 3.50 y 3.51 en la Ec. 3.49 y aproximando las derivadas parciales mediante diferencias centrales, se obtiene la siguiente representación de la Ec. 3.40 en diferencias finitas: 1 pin +1 − pin 1 pin++11 − 2 pin +1 + pin−+11 1 pin+1 − 2 pin + pin−1 μB n +1 n ...(3.52) + + q q = + i i 2 2k η Δt 2 Δx 2 Δx 2

(

Rearreglando,

αpin−+11 − 2(α + 1) pin +1 + αpin++11 = −αpin−1 + 2(α − 1) pin − αpin+1 − i = 1,2,..., I

)

ημB k

(

Δt qin +1 + qin

)

...(3.53)

n = 0,1,2,...

O bien,

cpin−+11 + apin +1 + bpin++11 = d in

...(3.54)

Se observa que la aproximación de la Ec. 3.49 mediante el esquema Crank-Nicholson produce nuevamente un sistema algebraico de ecuaciones lineales en las incógnitas pin +1 , i = 1,2,..., I . En esta aplicación simple, se puede apreciar el trabajo adicional que requiere el esquema Crack-Nicholson, comparado con el esquema implícito.

40

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Después de revisar los esquemas de aproximación anteriores, es natural preguntarse: 1. ¿Cómo se comparan las soluciones obtenidas mediante cada uno de los métodos con la solución exacta?, 2. ¿Cuál de los métodos produce la mejor solución?, 3. ¿Cuál es el efecto de Δx y Δt sobre el desempeño numérico de los métodos y sobre la solución? Se pueden obtener respuestas cualitativas empleando los conceptos de: consistencia, convergencia y estabilidad de una aproximación numérica, que se revisan brevemente a continuación.

3.7 Conceptos de Consistencia, Convergencia y Estabilidad de una Aproximación Numérica. Para ejemplificar estos conceptos, se considera el problema de flujo unidimensional de un fluido de compresibilidad pequeña y constante, ecuación de difusividad. Se define el siguiente operador diferencial A . ∂2 1 ∂ A= 2 − η ∂t ∂x De esta manera, la Ec. 3.36 se escribe como: μB ...(3.55) Ap + q=0 k Se considera la aproximación de la Ec. 3.55 producida por el esquema explícito, es decir: p n − 2 pin + pin−1 μB n 1 pin +1 − pin ...(3.56) {Ap}in + μB qin = i +1 + qi − + Oc Δx 2 − O p (Δt ) 2 η k k Δ t Δx

( )

i = 1,2,..., I n = 0,1,2,...

Ahora bien, se define el operador de diferencias Lp como sigue: p n − 2 pin + pin−1 1 pin+1 − pin n L p p i = i +1 − η Δt Δx 2

{ }

La Ec. 3.56 puede entonces escribirse como: {Ap}in + μB qin = {L p p}in + μB qin + Oc Δx 2 − O p (Δt ) = 0 k k

( )

De donde se obtiene:

{Ap}in − {L p p}in

( )

= Oc Δx 2 − O p (Δt ) = Rin

...(3.57)

...(3.58)

...(3.59)

Esta expresión indica que la diferencia entre la representación exacta y aproximada del problema de flujo en cuestión, esta dada por el error local de truncamiento o de discretización, Rin .

41

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Con estos antecedentes, se dice en lo general que un operador de diferencias L es consistente con el operador diferencial A, al cual aproxima, si el orden local de la aproximación es mayor o igual que uno. Que L sea consistente con A, tomando como ejemplo la Ec. 3.59, significa que en el límite, cuando Δx → 0 y Δt → 0 , la norma del vector de errores de discretización también tiende a cero, R → 0. Es claro que una aproximación numérica deberá ser consistente para que tenga algún valor práctico. Además de ser consistente, el operador L debe ser convergente a A. Para analizar la propiedad de convergencia se toma como punto de partida la Ec. 3.58, la cual puede reescribirse como sigue: ...(3.60) {L p p}in + μkB qin = − Rin Como se mencionó previamente, el término Rin se desconoce. Bajo la suposición que es pequeño, se desprecia. Entonces se tiene: n μB n ...(3.61) Lp p∗ i + qi = 0 k

{

}

Nótese que en la Ec. 3.61 se usó p ∗ , para resaltar el carácter aproximado de la solución, comparada con la solución exacta de la Ec. 3.60. Se define ahora el Error Global de la Solución como la diferencia entre las soluciones exacta y aproximada:

ε in = pin − pi∗n

...(3.62)

Entonces se dice que el operador Lp es convergente al operador diferencial A si la norma del vector de errores globales tiende a cero, ε → 0, cuando Δx → 0 y Δt → 0. Existe una relación entre los errores locales de discretización y los errores globales de la solución, ε y Ri . Esta se establece restando las Ecs. 3.60 y 3.61:

{L p p}in − {L p p ∗ }i

n

{

⎛ ⎞ = L p ⎜⎜ p − p ∗ ⎟⎟ ⎝

} = {L s} = −R n

⎠ i

p

n i

n i

...(3.63)

o sea:

ε in++11 − 2ε in+1 + ε in−+11 Δx 2

1 ε in +1 − ε in − = − Rin η Δt

...(3.64)

Se puede ver en la Ec. 3.64 que la ecuación que describe el comportamiento de los errores globales de la solución, ε , es la misma ecuación que gobierna el comportamiento de p , excepto que el término fuente, qi se substituye, en este caso, por el error local de discretización, Rin . Esto se cumple solamente en problemas lineales. La Ec. 3.64 implica que:

42

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

ε →0

cuando

R →0

La convergencia de un operador de diferencias puede probarse indirectamente probando la estabilidad del operador. Esto es posible a partir del Teorema de equivalencia de Lax. Teorema de Equivalencia de Lax: La estabilidad es una condición necesaria y suficiente para que exista la convergencia, cuando la aproximación es consistente.

Estabilidad es un concepto que aplica a problemas dependientes del tiempo. Un algoritmo numérico es estable si cualquier error, introducido en alguna etapa de los cálculos, no se amplifica en cálculos subsecuentes. Existen varios métodos para probar la estabilidad de una aproximación numérica, entre los que se puede mencionar, el método de series de Fourier, el método matricial y el método de la energía. Se ha demostrado que el esquema explícito aplicado a la Ec. 3.35 es condicionalmente estable. La condición de estabilidad es:

α=

η Δt Δx 2

≤

1 2

...(3.65)

Se puede ver, en la Ec. 3.65, que para mantener estable la solución numérica de la Ec. 3.35 mediante un esquema explícito, existe un compromiso entre el tamaño del incremento de tiempo, Δt y el tamaño de las celdas Δx . Análisis de estabilidad de los esquemas implícitos y Crank-Nicholson, aplicados a la solución numérica de la Ec. 3.35, muestran que estos incondicionalmente estables. 3.8

Construcción de Mallas.

Como se evidenció en las secciones anteriores, la solución numérica de las ecuaciones de flujo de fluidos en medios porosos consiste en obtener una representación aproximada de las ecuaciones diferenciales en puntos específicos del espacio y del tiempo, i = 1,2,...,I y n = 0, 1,2,... para el problema unidimensional. El dominio del problema, en espacio y en tiempo, se segmenta o discretiza; se genera así una malla de cálculo, constituida de celdas y nodos, donde se calcula la solución en etapas sucesivas de tiempo. Existen básicamente dos tipos de mallas en la simulación numérica: (a) Mallas de nodos distribuidos. (Uniformes y No uniformes). (b) Mallas de bloques centrados. (Uniformes y No uniformes). Los nodos y los bloques o celdas a su vez pueden ser distribuidos de manera uniforme o no uniforme. Las mallas no uniformes son necesarias en la simulación de problemas con regiones que experimentan cambios fuertes en la presión y en las saturaciones, a lo largo del tiempo. A continuación se indican los métodos para construirlas.

43

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

3.8.1

Malla Uniforme de Nodos Distribuidos:

En un yacimiento lineal de longitud L, y área transversal al flujo A, una malla de nodos distribuidos uniformemente se construye, según lo esquematizado en la Fig. 3.2. Primeramente se colocan nodos en las fronteras del yacimiento y entre ellos se distribuyen, con espaciamiento uniforme, el resto de los nodos. Una vez definida la posición de los nodos, se procede a definir la posición de las fronteras de las celdas en el punto medio entre los nodos.

Fig. 3.2 Malla Uniforme de Nodos Distribuidos

Δxi +1 2

1

2

Fronteras de la celda

I

re

Nodos

Si se considera un total de I nodos, su espaciamiento, Δx, es: L Δxi = ...(3.66) I −1 Nótese que el volumen de las celdas situadas en las fronteras es la mitad del volumen de las celdas internas, es decir: i = 1, I ⎧ AΔxi / 2 ...(3.67) Vri = ⎨ i = 2,3,..., I - 1 ⎩ AΔxi La posición de los nodos es: xi = (i − 1)Δxi i = 1,2,...,I ...(3.68) Y la posición de la frontera de las celdas es: 1 x 1 = (i − )Δxi i = 1,2,...,I ...(3.69) i+ 2 2

3.8.2 Malla Uniforme de Bloques Centrados: Este tipo de mallas se construyen situando celdas de tamaño uniforme en el yacimiento y ubicando posteriormente los nodos en el centro de cada una de ellas, véase la Fig. 3.3: Fig. 3.3 Malla Uniforme de Nodos Centrados

44

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Δxi +1 2

i=

1

2

Fronteras de la celda

I

Δxi

Nodos Si se considera un total de I nodos, el tamaño de las celdas se calcula como: Δxi =

L I

Todas las celdas son de igual volumen, esto es: Vri = AΔxi Las coordenadas de los nodos son: ⎛ 1⎞ xi = ⎜ i − ⎟Δxi ⎝ 2⎠ y las coordenadas de las fronteras de las celdas: x 1 = iΔxi i+

...(3.70)

i = 1,2,3,..., I

...(3.71)

i = 1,2,3,..., I

...(3.72)

i = 1,2,3,..., I

...(3.73)

2

3.8.3 Definición de la malla radial en coordenadas cilíndricas. La discretización de la región de drene de un pozo requiere de una malla no uniforme en la dirección radial, para representar adecuadamente las fuertes variaciones de presiones y de saturaciones en las vecindades del pozo, intervalo disparado. La mejor representación de tales condiciones de flujo se obtiene definiendo el tamaño de las celdas de manera que la caída de presión entré ellas sea de la misma magnitud. Esto es posible mediante el empleo de una malla logarítmica convencional de nodos distribuidos o de bloques centrados. La forma de definir estas mallas se describe a continuación: Considérese flujo radial en régimen permanente con la viscosidad y el factor de volumen del fluido constantes, de tal manera que el gasto está representado por la Ec. 3.74: ∂p 2πkhr ∂p ...(3.74) q=− = −2πhλr = constante μB ∂r ∂r Donde k λ= μB Integrando la Ec. 3.74,

45

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

q = −2πhλ

p − p ref ln(r rref )

...(3.75)

La generación de esta malla radial se establece a partir de este gasto, de la siguiente manera: La representación exacta del gasto en la frontera entre celdas, i + 1 , ver Fig. 3.4, está dada por: 2

Δri +1 2

i-1 rw

i i−1

2

Δri

re

i+ 1

2

Fig. 3.4 Malla radial logarítmica.

q E = q = −2πhλ

p i +1 − ln(ri +1

pi ri )

Donde: q E = Gasto Exacto y la representación numérica del gasto en la frontera entre celdas, está dada como: ⎛ r ⎞ q N = qi + 1 = 2πh⎜ ⎟ λ ( pi +1 − pi ) 2 ⎝ Δr ⎠ i + 1

...(3.76)

...(3.77)

2

Por lo que si se considera que:

q E = − qi + 1

...(3.78)

2

entonces igualando ambas ecuaciones, Ecs. 3.76 y 3.77, se tiene: p − pi ⎛ r ⎞ − 2πhλ i +1 = −2πh⎜ ⎟ λ ( p i +1 − p i ) ln(ri +1 ri ) ⎝ Δr ⎠ i + 1

...(3.79)

2

Simplificando y despejando a ri + 1 de la Ec. 3.79, se establece que los radios de las 2 fronteras de las celdas se calculan como el promedio logarítmico del radio de los nodos, es decir: r −r ...(3.80) ri + 1 = i +1 i n(ri +1 ri ) 2 Por otra parte, el gasto exacto que pasa a través de las fronteras, está dado como: Para la frontera en i + 1 , 2

q Ei + 1 = −2πhλ 2

p i +1 − ln(ri +1

pi ri )

...(3.81)

Para la frontera en i − 1 , 2

46

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

q Ei − 1 = −2πhλ 2

pi − pi −1 ln(ri ri −1 )

...(3.82)

y si este gasto es una constante, entonces satisface la siguiente relación: q Ei + 1 ln(ri ri −1 )( pi +1 − pi ) 2 = =1 q Ei − 1 ln(ri +1 ri )( pi − pi −1 ) 2

...(3.83) Para satisfacer la Ec. 3.83, se tiene que imponer ( pi +1 − pi ) = ( pi − pi −1 ), por lo que simplificándola, se tiene que:

la

siguiente

condición:

ln(ri +1 ri ) = ln(ri ri −1 ) = constante

...(3.84)

Esto implica que: ri +1 r = i = α ó bien ri +1 = αri donde i= 1, 2, 3… I. Valido para cualquier tipo de malla. ri ri −1 Sumando esta constante a lo largo del dominio discretizado, se establece: I −1

∑ ln(ri +1

i =1

ri ) = ( I − 1) ln(ri +1 ri ) = ln(rI r1 )

...(3.85)

3.8.3.1 Malla Radial de nodos distribuidos.

Ahora bien, la distribución de los radios de los nodos, para una malla de nodos distribuidos, ver Fig. 3.5, está dada como sigue: Δri +1 2

i-1 rw

i i−1

2

i+1 Δri

re

i+ 1

2

Fig. 3.5 Malla radial logarítmica de nodos distribuidos.

Si se considera que r1 = rw y rI = re, la Ec. 3.85 se escribe como: I −1

∑ ln(ri +1

i =1

ri ) = ( I − 1) ln(ri +1 ri ) = ln(re rw )

...(3.86)

Por lo que,

(ri +1 ri ) = (re rw )1 /(i −1) = α

...(3.87)

Donde, α es el factor geométrico,

47

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

re es el radio de drene del pozo, rw es el radio del pozo e I es el número de celdas en la dirección radial. Entonces, el radio de los nodos de las celdas esta dado como: ri +1 = α ri

3.8.3.2

...(3.88)

Malla Radial de bloques centrados.

La distribución de los radios de los nodos, para una malla de bloques centrados, ver Fig. 3.4, está dada también por la Ec. 3.88. La diferencia con respecto a la malla de nodos distribuidos está en como calcular ahora el radio del pozo, rw , el radio de drene del pozo, re y el factor geométrico, α , los cuales se calculan para esta malla como sigue: La Ec. 3.85 aplica directamente a este tipo de malla, por lo que simplificándola, se establece que:

(ri +1 ri ) = (rI r1 )1 /( I −1) = α

...(3.89)

Pero de acuerdo con las Ecs. 3.80 y 3.88, se tiene que el radio del pozo, rw , se calcula como:

r r1 − 1

r −r α = (α − 1) r1 r1 = rw = 1 0 = ln α ln α α ln α 2

...(3.90)

Despejando a r1 de la Ec. 3.90, se tiene que:

r1 =

α ln αr (α − 1) w

...(3.91)

En forma similar a rw , se tiene que el radio de drene del pozo, re , se calcula como:

αr − rI (α − 1) r −r rI + 1 = re = I +1 I = I rI = ln α ln α ln α 2

...(3.92)

Despejando a rI de la ec. 3.92, se tiene que:

rI =

ln α r (α − 1) e

...(3.93)

Substituyendo las Ecs. 3.91 y 3.93 en la Ec. 3.89, se obtiene α a partir de los radios del pozo y de drene, así como el número de celdas en que se discretiza la dirección radial, esto es:

⎛ ln α ⎞ re ⎟ ⎜ α = (rI r1 )1 (I −1) = ⎜ α − 1 ⎟ ⎜ ln α αrw ⎟⎟ ⎜ ⎠ ⎝ α −1

1 (I −1)

= (re αrw )1 (I −1)

...(3.94)

48

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Simplificando la Ec. 3.94, se tiene que el factor geométrico se calcula como:

α = (re rw )1 / I

...(3.95)

4. Simulación Numérica del Flujo Monofásico Unidimensional. En este capítulo se revisara la simulación numérica del flujo monofásico unidimensional en medios porosos. En primer término se trata el problema de flujo unidimensional; posteriormente en el capitulo siguiente se extenderán los conceptos al caso de flujo multidimensional.

4.1 Flujo Monofásico unidimensional. Considérese el flujo lineal, en coordenadas cartesianas, de un fluido compresible en un yacimiento. La distribución de presiones al tiempo cero, p ( x, t = 0 ), esta condicionada al equilibrio gravitacional del fluido, ver Ec. 2.62. El yacimiento produce a un gasto q w en la frontera situada en x = 0 . La frontera ubicada en x = L esta cerrada al flujo. Matemáticamente esto se expresa como: ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ∂ ⎛φ ⎞ λ⎜ − γ ⎟⎥ + q* = ⎜ ⎟ ⎢ ∂t ⎝ B ⎠ ∂x ⎣ ⎝ ∂x ∂x ⎠⎦ __________

__

__

(a)

(b)

(c)

...(4.1)

0 < x < L,

t>0

(a) Término de flujo y gravedad. (b) Término de fuente o sumidero en el dominio, no en las fronteras, @c.s. Gasto por unidad de volumen. (c) Término de acumulación. Las condiciones iniciales y de frontera se expresan como: 0< x
q μB ∂D ⎞ ⎛ ∂p =− w ; ⎜ −γ ⎟ kA ∂x ⎠ x = 0 ⎝ ∂x

t>0

...(4.2) ...(4.3)

49

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Y,

∂D ⎞ ⎛ ∂p =0 ; ⎜ −γ ⎟ ∂x ⎠ x = L ⎝ ∂x

t>0

...(4.4)

4.2 Aproximación de las ecuaciones diferenciales en diferencias finitas. Se desea resolver numéricamente el problema definido por las Ec. 4.1 a 4.4. La aproximación en diferencias finitas de la Ec. 4.1 en el nodo i, de la malla de bloques centrados espaciada irregularmente, considerando un esquema implícito, se inicia con: n +1

⎧ ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ⎫ ⎟ ⎬ ⎨ ⎢λ ⎜ − γ ∂x ⎠⎥⎦ ⎭i ⎩ ∂x ⎣ ⎝ ∂x

+ {q *}in +1

n +1

⎧ ∂ ⎛ φ ⎞⎫ = ⎨ ⎜ ⎟⎬ ⎩ ∂t ⎝ B ⎠⎭ i

...(4.5)

i = 1, 2,..., I n = 1, 2,...

La aproximación del término de flujo mediante diferencias finitas centrales, véase la Ec. 3.25, es: n +1 n +1 ⎧ n +1 ⎧ ∂ ⎡ ⎛ ∂p 1 ⎪⎛ λ ⎞ ∂D ⎞⎤ ⎫ ⎡ p − p − (γΔD ) 1 ⎤ ≈ − λ γ − ⎜ ⎟ ⎜ ⎟ ⎨ ⎢ ⎬ ⎨ i +1 i i + 2 ⎥⎦ Δxi ⎪⎝ Δx ⎠ i + 1 ⎢⎣ ∂x ⎠⎥⎦ ⎭ i ⎩ ∂x ⎣ ⎝ ∂x 2 ⎩ ...(4.6) n +1 ⎫ n +1 ⎪ ⎛ λ ⎞ ⎡ p − p − (γΔD ) 1 ⎤ ⎜ ⎟ ⎬ i i − 1 i − 2 ⎥⎦ ⎝ Δx ⎠ i − 1 ⎢⎣ ⎪⎭ 2 La aproximación del término de acumulación de la Ec. 4.5 mediante diferencias regresivas, véase la Ec. 3.26, es: n +1 n +1 n ⎧ ∂ ⎛ φ ⎞⎫ 1 ⎡⎛ φ ⎞ ⎛φ ⎞ ⎤ ⎢ ⎥ ...(4.7) − ≈ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎨ ⎬ B ⎠i ⎥ Δt ⎢⎝ B ⎠ i ⎝ ⎩ ∂t ⎝ B ⎠⎭ i ⎣ ⎦ Sustituyendo las Ecs. 4.6 y 4.7 en la Ec. 4.5 y multiplicando por el volumen de roca de la celda i, Vri = AΔxi , se obtiene: Gasto que pasa a través de I + 1

Gasto que pasa a través de I − 1

644444 474444424 8 T n +11 ⎡ pi +1 − pi − (γΔD )i + 1 ⎤ i + 2 ⎢⎣ 2⎥ ⎦ + qin +1

n +1

644444 474444424 8 − T n +11 ⎡ pi − pi −1 − (γΔD )i − 1 ⎤ i − 2 ⎢⎣ 2⎥ ⎦

Vr = i Δt

n +1

...(4.8)

⎡⎛ φ ⎞ n +1 ⎛ φ ⎞ n ⎤ ⎢⎜ ⎟ −⎜ ⎟ ⎥ ⎝ B ⎠ i ⎥⎦ ⎢⎣⎝ B ⎠ i i = 1,2,..., I

n = 1,2,... Donde T es la transmisibilidad del fluido, dada por:

50

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

⎛ A⎞ Ti + 1 = T ( pi +1 , pi ) = ⎜ ⎟ λ 1 2 ⎝ Δx ⎠ i + 1 i + 2

...(4.9)

⎛ A⎞ Ti − 1 = T ( pi , pi´−1 ) = ⎜ ⎟ λ 1 2 ⎝ Δx ⎠ i − 1 i − 2

...(4.10)

2

2

La Ec. 4.8 puede expresarse en términos de operadores de diferencias, Ecs. 3.29 y 3.32, como sigue: Vr ⎛φ ⎞ T n +11 [Δp − γΔD ]n +11 − T n +11 [Δp − γΔD ]n +11 + qin +1 = i Δ t ⎜ ⎟ ...(4.11) i+ 2 i+ 2 i+ 2 i− 2 Δt ⎝ B ⎠ i O bien, Vr ⎛φ ⎞ ...(4.12) Δ[T (Δp − γΔD )]n +11 + qin +1 = i Δ t ⎜ ⎟ i+ 2 Δt ⎝ B ⎠ i i = 1,2,...I n = 0,1,2,...

4.3 Acoplamiento de las Condiciones de frontera. Se necesita ahora obtener la forma particular que adquiere la ecuación aproximada en las celdas de las fronteras, i = 1 e i = I , debido al acoplamiento de las condiciones de frontera. Para esto, considérese la porción de la malla de cálculo correspondiente a la frontera ubicada en x = 0, mostrada en la figura 4.1, donde se muestra el nodo ficticio i = 0 : Fig. 4.1 Malla de cálculo en la porción de la frontera en x = 0, mostrando el nodo ficticio i = 0 X=0

i=

0

1/2

1

11/2

21/2

2

La aproximación de la condición de frontera, Ec. 4.3, en x = 0, o i = 1 2 , empleando diferencias centrales es: n +1

∂D ⎞ ⎛ ∂p ⎜ −γ ⎟ ∂x ⎠ 1 ⎝ ∂x

2

≈

1 Δx 1

⎡ p − p − (γΔD ) 1 ⎤ 0 ⎢⎣ 1 2⎥ ⎦

n +1

2

n +1

⎛ μB ⎞ =⎜ ⎟ ⎝ kA ⎠ 1

qw

...(4.13)

2

O bien:

q w = −T 1n +1 ⎡ p1 − p0 − (γΔD ) 1 ⎤ ⎢⎣ 2⎥ ⎦ 2

n +1

...(4.14)

51

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

La condición de frontera en x = L, Ec. 4.4 se aproxima similarmente como, véase la Fig. 4.2: Fig. 4.2 Malla de cálculo en la porción de la frontera en x = L, mostrando el nodo ficticio i = I + 1

X=L

i=

I-2

I-11/2

I-1

I-1/2

I

I+1/2

I+1

En x = L, o i = I + 1 2 , se tiene que:

n +1

1 ∂D ⎞ ⎛ ∂p ≈ ⎜ −γ ⎟ ∂x ⎠ I + 1 Δx I + 1 ⎝ ∂x 2 2

⎡p − p I − (γΔD )I + 1 ⎤ ⎢⎣ I +1 2⎥ ⎦

n +1

=0

...(4.15)

=0

...(4.16)

Entonces,

T n +11 ⎡ p I +1 − p I − (γΔD )I + 1 ⎤ I + 2 ⎢⎣ 2⎥ ⎦

n +1

Considérese ahora la Ec. 4.8 escrita en la celda I , i = I, con q In +1 = 0 :

T1n.5+1 [ p 2 − p1 − (γΔD )1.5 ]n +1 − T 1n +1 ⎡ p1 − p0 − (γΔD ) 1 ⎤ ⎢⎣ 2⎥ ⎦

n +1

2

=

n +1 n Vr1 ⎡⎛ φ ⎞ ⎛φ ⎞ ⎤ ⎢⎜ ⎟ −⎜ ⎟ ⎥ Δt ⎢⎝ B ⎠1 ⎝ B ⎠1 ⎥⎦ ⎣

...(4.17)

La condición de frontera, Ec. 4.14, puede ahora acoplarse en esta ecuación. Así la ecuación aproximada para el nodo 1 adquiere la siguiente forma: n n +1 Vr ⎡⎛ φ ⎞ ⎛φ ⎞ ⎤ ...(4.18) T1n.5+1 [ p 2 − p1 − (γΔD )1.5 ]n +1 + q w = 1 ⎢⎜ ⎟ −⎜ ⎟ ⎥ B ⎠1 ⎥ Δt ⎢⎝ B ⎠1 ⎝ ⎣ ⎦ Véase ahora la forma particular que adquiere la Ec. 4.8 en la celda I: T n +11 ⎡ p I +1 − p I − (γΔD )I + 1 ⎤ I + 2 ⎢⎣ 2⎥ ⎦

n +1

− T n +11 ⎡ p I − p I −1 − (γΔD )I − 1 ⎤ I − 2 ⎢⎣ 2⎥ ⎦

n +1

=

VrI Δt

⎡⎛ φ ⎞ n +1 ⎛ φ ⎞ n ⎤ ⎢⎜ ⎟ −⎜ ⎟ ⎥ ⎝ B ⎠ I ⎥⎦ ⎢⎣⎝ B ⎠ I

Acoplando la condición de frontera, Ec. 4.16, esta ecuación se reduce a: n n +1 Vr ⎡ φ n +1 ⎛ ⎞ ⎛φ ⎞ ⎤ n +1 ⎡ I ⎤ ⎢⎜ ⎟ − T 1 p I − p I −1 − (γΔD )I − 1 = −⎜ ⎟ ⎥ I − 2 ⎢⎣ 2⎥ ⎦ Δt ⎢⎝ B ⎠ I ⎝ B ⎠ I ⎥⎦ ⎣

...(4.19)

...(4.20)

52

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Las Ecs. 4.18 para el nodo 1, i = 1, conjuntamente con la Ec. 4.8 para los nodos i = 2,3,..., I − 1 y la Ec. 4.20 para el nodo I , i = I , constituyen un sistema de I ecuaciones algebraicas no lineales con las incógnitas pin +1 , i = 1,2,..., I . Existen varios métodos para resolver este sistema de ecuaciones.

4.4 Cálculo de las Transmisibilidades en las frontera: Ti + 1 , Ti − 1 . 2 2 Antes de revisar los métodos de linealización, es conveniente definir la manera de calcular las transmisibilidades en las fronteras de las celdas, Ti + 1 , Ti − 1 , las cuales están definidas por las 2

2

Ecs. 4.9 y 4.10, respectivamente. Dado que solo se conoce las presiones en los nodos de las celdas en que se discretiza el dominio de interés, i = 1,2 ,...,I, las transmisibilidades deben evaluarse en tal forma que se obtenga la mejor representación posible del flujo entre nodos vecinos.

4.4.1 Evaluación de la permeabilidad absoluta en las fronteras, k i ± 1 : 2 Supóngase que existe una discontinuidad en las propiedades de los nodos i e i + 1 , y que la discontinuidad se localiza en el punto i + 1 2 , la frontera de la celda como se muestra en la Fig. 4.3: Fig. 4.3 Discontinuidad en las propiedades de las celdas

Δxi + 1

2

Δxi +1

Δxi i + 12

El flujo entre los nodos i e i + 1, suponiendo flujo lineal Marciano, es:

⎛ k 1 ⎞ ⎛ ⎟⎟ ⎜ pi + 1 − pi ⎞⎟ = − qi + 1 = A⎜⎜ 2 2 ⎠ ⎝ Δx μB ⎠ i ⎝

⎛ k 1 ⎞ ⎛ ⎟⎟ ⎜ pi +1 − pi + 1 ⎞⎟ A⎜⎜ 2 ⎠ ⎝ Δx μB ⎠ i +1 ⎝

...(4.22)

Se desea entonces tratar el flujo entre los nodos i e i + 1 mediante una λi + 1 , 2 correspondiente a un medio homogéneo equivalente, tal que produzca el mismo qi + 1 , o sea: 2

53

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

⎛ k 1 ⎞ ⎟⎟ ( pi +1 − pi ) − qi + 1 = A⎜⎜ 2 ⎝ Δx μB ⎠ i +1

...(4.23)

Esta ecuación puede escribirse como:

⎛ k 1 ⎞ ⎛ ⎟⎟ ⎜ pi +1 − pi + 1 + pi + 1 − pi ⎞⎟ − qi + 1 = A⎜⎜ 2 2 2 ⎠ ⎝ Δx μB ⎠ i +1 ⎝

...(4.24)

De la Ec. 4.22, se puede despejar a k i + 1 , conduce a: 2

ki + 1 = 2

Δxi + Δxi +1 ⎛ Δx ⎞ ⎛ Δx ⎞ ⎜ ⎟ ⎜ ⎟ ⎝ k ⎠ i ⎝ k ⎠ i +1

...(4.25)

Lo que indica que k i + 1 , debe calcularse como la media armónica de k i e k i +1. 2 Se puede hacer un análisis similar para el caso de flujo radial. Esto conducirá a la siguiente expresión para k i + 1 : 2

ki+ 1= 2

⎛ ⎛ ⎜ ln⎜ ri + 1 ⎜ ⎝ 2 ⎜ k ⎜ ⎝

ln (ri +1 ri ) ⎞ ⎛ ⎞ ri ⎞⎟ ⎟ ⎜ ln⎛⎜ ri +1 ri + 1 ⎞⎟ ⎟ 2 ⎠⎟ ⎠⎟ +⎜ ⎝ ⎟ ⎜ ⎟ k ⎟ ⎜ ⎟ ⎠i ⎝ ⎠ i +1

...(4.26)

4.4.2 Evaluación de la movilidad en las fronteras, k i + 1 : 2 4.5 Métodos de Linealización de las Ecuaciones de Flujo en Diferencias Finitas: 6. 6. 6. 6.

Linealización Directa. Método semi-implícito linealizado. Método de iteraciones sucesivas. Método de Newton Raphson, o de iteraciones Newtoniana.

Estos métodos han sido clásicamente empleados en la solución de sistemas de ecuaciones no lineales. A continuación se presenta un enfoque general de linealización de las ecuaciones aproximadas del flujo de fluidos en medios porosos. El enfoque se basa en el método de Newton, que puede considerarse como el método general: se verán lo métodos de linealización directa, de aproximaciones sucesivas, y semi-implícito linealizado como casos particulares de este método general.

54

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

.5.1

Métodos de Newton- Raphson: El Método General.

Con el objeto de simplificar la exposición del método, considérese flujo lineal horizontal en la Ec. 4.8. Se define entonces la siguiente función de residuos, Fi en i = 1,2,..., I , en el nivel de tiempo n + 1 . n +1

⎡ p − p − (γΔD ) 1 ⎤ i i + 2 ⎥⎦ i + 12 ⎢⎣ i +1

Fi ( pi −1 , pi , pi +1 )n +1 = T

n +1

+ qi

V − ri Δt

n +1

n +1

⎡ p − p − (γΔD ) 1 ⎤ i −1 i − 2 ⎥⎦ i − 12 ⎢⎣ i

−T

n +1

⎡⎛ φ ⎞ n +1 ⎛ φ ⎞ n ⎤ ⎢⎜ ⎟ −⎜ ⎟ ⎥ = 0 ⎝ B ⎠ i ⎥⎦ ⎢⎣⎝ B ⎠ i

...(4.27)

i = 1,2,..., I n = 1,2,...

La función de residuos del nodo i depende de las incógnitas pi −1 , pi , pi +1 . El método de Newton resuelve las incógnitas pin +1 , i = 1,2,..., I , en forma iterativa. El proceso iterativo se fundamenta en la expansión de la función de residuos, Fin +1 , en

series de Taylor alrededor del nivel de iteración (ν ) . De esta expansión sólo se conservan los términos de menor orden, es decir: (ν ) (ν ) (ν ) ⎛ ∂Fi ⎞ ⎛ ∂F ⎞ ⎛ ∂Fi ⎞ ⎟ δp (ν +1) + ⎜ i ⎟ δp (ν +1) + ⎜ ⎟ δp (ν +1) = 0 ...(4.28) F (ν +1) ( p , p , p ) ≈ F (ν ) + ⎜ i

i −1

i

i +1

i

⎜ ∂p ⎟ ⎝ i −1 ⎠

Por lo tanto ⎛ ∂Fi ⎞ ⎟⎟ ⎜⎜ ⎝ ∂pi −1 ⎠

(ν )

⎛ ∂F ⎞ δpi(ν−1+1) + ⎜⎜ i ⎟⎟ ⎝ ∂p i ⎠

(

(ν )

⎜ ∂p ⎟ ⎝ i⎠ i = 1,2,..., I n = 1,2,...

i −1

⎛ ∂Fi ⎞ ⎟⎟ δp i(ν +1) + ⎜⎜ ⎝ ∂p i +1 ⎠

(ν )

i

⎜ ∂p ⎟ ⎝ i +1 ⎠

δp i(ν+1+1) = − Fi(ν )

i +1

...(4.29)

)

Donde δpiν +1 representa los cambios iterativos de las incógnitas:

δpi(ν +1) = pi(ν +1) − pi(ν )

...(4.30)

Nótese que con el fin de simplificar la escritura, el superíndice correspondiente al nivel de tiempo n + 1 fue eliminado de la Ec. 4.28. La forma particular que adopta la Ec. 4.26 en los nodos i = 1 e i = I , debido al acoplamiento de las condiciones de frontera, es respectivamente:

55

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

(

)

⎛ ∂F ⎞ F1(ν +1) p 1 , p 2 ≈ F1(ν ) + ⎜⎜ 1 ⎟⎟ ⎝ ∂p1 ⎠

(ν )

⎛ ∂F ⎞ δp1(ν +1) + ⎜⎜ 1 ⎟⎟ ⎝ ∂p 2 ⎠

Y,

(

)

(ν )

⎛ ∂FI ⎞ ⎟⎟ FI(ν +1) p I −1 , p I ≈ FI(ν ) + ⎜⎜ ⎝ ∂p I −1 ⎠

(ν )

δp 2(ν +1) = 0

⎛ ∂F ⎞ δp I(ν−+11) + ⎜⎜ I ⎟⎟ ⎝ ∂p I ⎠

(ν )

δp I(ν +1) = 0

...(4.31)

...(4.32)

De las Ecs. 4.28, 4.31 y 4.32, se establece el siguiente sistema algebraico de ecuaciones lineales: Para i = 1 , ⎛ ∂F1 ⎞ ⎜⎜ ⎟⎟ ⎝ ∂p1 ⎠

(ν )

⎛ ∂F ⎞ δp1(ν +1) + ⎜⎜ 1 ⎟⎟ ⎝ ∂p 2 ⎠

Para i = 2,3,..., I − 1 , ⎛ ∂Fi ⎞ ⎜⎜ ⎟⎟ ⎝ ∂pi −1 ⎠

(ν )

(ν )

⎛ ∂F ⎞ δpi(ν−1+1) + ⎜⎜ i ⎟⎟ ⎝ ∂p i ⎠

(ν )

δp 2(ν +1) = − F1(ν )

⎛ ∂Fi ⎞ ⎟⎟ δp i(ν +1) + ⎜⎜

(ν )

⎝ ∂p i +1 ⎠

δp i(ν+1+1) = − Fi(ν )

...(4.33)

...(4.34)

Para i = I ,

⎛ ∂FI ⎞ ⎜⎜ ⎟⎟ ⎝ ∂p I −1 ⎠

(ν )

⎛ ∂F ⎞ δp I(ν−+11) + ⎜⎜ I ⎟⎟ ⎝ ∂p I ⎠

(ν )

δp I(ν +1) = − FI(ν )

El proceso iterativo se inicia generalmente con la siguiente estimación: i = 1,2,3,..., I pi(0 ) = pin ;

...(4.35)

...(4.36)

(ν +1) , se calculan las presiones en la iteración

ν + 1:

Una vez resueltos los cambios iterativos, δpi

pi(ν +1) = pi(ν ) + δp i(ν +1)

...(4.37)

Se considera que el proceso iterativo converge a la solución cuando los cambios iterativos de las incógnitas, en valor absoluto, son menores que una tolerancia estipulada, ε p :

δpi(ν +1) 〈 ε p ; para toda i

...(4.38)

Otro criterio de convergencia es considerar la norma del vector de residuos, F = (F1 , F2 ,..., FI ), a adquirir un valor menor que ε p : Fi(ν +1) 〈 ε p ; para toda i

...(4.39)

56

FI-UNAM

.5.2

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Estructura Matricial del Sistema Lineal de Ecuaciones.

La estructura matricial del sistema de Ecs. 4.33 - 4.35, se muestra a continuación:

⎡ a1 ⎢c ⎢ 2 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

b1 a2 c3

b2 a3 .

b3 . .

. . .

. . cI

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ . ⎥ a I ⎥⎦

(ν ) δp (ν +1) ⎡ 1⎤ ⎢δp ⎥ ⎢ 2⎥ ⎢δp3 ⎥ ⎥ ⎢ ⎢ . ⎥ ⎢ . ⎥ ⎥ ⎢ ⎢ . ⎥ ⎢ δp ⎥ ⎣ i⎦

⎡ F1 ⎤ ⎢F ⎥ ⎢ 2⎥ ⎢ F3 ⎥ ⎢ ⎥ = −⎢ . ⎥ ⎢ . ⎥ ⎢ ⎥ ⎢ . ⎥ ⎢F ⎥ ⎣ i⎦

(ν )

Como puede observarse, el problema de flujo lineal genera un sistema tridiagonal de ecuaciones. En forma resumida se puede escribir el sistema de ecuaciones como sigue: ...(4.40) [J ](ν ) δp (ν +1) = − F (ν ) Donde [J ] es la matriz de derivadas o matriz Jacobiana, cuyos elementos son: ∂Ti − 1 ⎛ ∂Fi ⎞ 2 ⎟⎟ = Ti − 1 − Δ x pi − 1 ci = ⎜⎜ 2 2 ∂p i −1 ⎝ ∂pi −1 ⎠ ⎛ ∂F ai = ⎜⎜ i ⎝ ∂pi

∂Ti + 1 ∂Ti − 1 ⎞ Vr ∂ ⎛ φ ⎞ 2 2 ⎟⎟ = −⎛⎜ Ti 1 + Ti 1 ⎞⎟ − Δ x pi 1 − Δ x pi − 1 − i ⎜ ⎟ + − + Δt ∂pi ⎝ B ⎠ i 2 ∂p i 2 ⎠ 2 ∂p i ⎝ 2 ⎠ ∂Ti + 1 ⎛ ∂Fi ⎞ 2 ⎟⎟ = T 1 + Δ x p 1 bi = ⎜⎜ i+ 2 i + 2 ∂p p ∂ i +1 ⎝ i +1 ⎠

.5.3

...(4.41)

...(4.42)

...(4.43)

Casos Particulares del Método de Newton.

Métodos de linealización directa, de aproximaciones sucesivas y semi-implícito linealizado pueden obtenerse como casos particulares del método de Newton. La relación entre el método de semi-implícito linealizado y el método de Newton es presentada en Aziz y Settari. Primeramente se nota que los elementos de la matriz Jacobiana están constituidos por tres tipos de términos. Estos son: transmisibilidades, derivadas de transmisibilidades y derivadas del término de acumulación. De esta manera se tiene: ...(4.44) ci = (cT )i + (cT ' )i ai = (aT )i + (aT ' )i + (a AC ' )i

y

bi = (bT )i + (bT ' )i

...(4.45) ...(4.46)

Donde:

57

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

(cT )i = Ti − 1

...(4.47)

2

y

(cT ' )i = −Δ x p

∂Ti − 1

...(4.48)

2

i− 1 2

∂pi −1

Los términos (bT )i , (bT ')i , (aT )i , (aT ')i , se definen similarmente. El término (a AC )i esta dado por: Vr ...(4.49) (a AC ' )i = − i ∂ ⎛⎜ φ ⎞⎟ Δt ∂p i ⎝ B ⎠ i Con estos antecedentes, se puede expresar la matriz Jacobiana, [J ], como la suma de las siguientes submatrices: [J ] = [T ] + [T '] + [(φ B )'] ...(4.50) Donde: 1. [T ] : Matriz tridiagonal con elementos (cT )i , (aT )i y (bT )i

2. [T '] : Matriz tridiagonal con elementos (cT ')i , (aT ')i y (bT ')i

3. [(φ B )'] : Matriz diagonal con elementos (a AC )i

De acuerdo con lo anterior, se puede derivar del método de Newton los métodos de linealización citados anteriormente. Estos se obtienen prelinealizando parcialmente la función de residuos, Fi y controlando el número de iteraciones, como se muestra a continuación:

.5.3.1 Método de Linealización Directa. Considérese en la función de residuos, Ec. 4.27, la siguiente aproximación: T n +1 ≈ T n

...(4.51)

La función de residuos del método de linealización directa se expresa entonces como:

(FLD )i = Δ[T n Δp n +1 ]i

+ qin +1 −

Vri ⎛ φ ⎞ Δt ⎜ ⎟ = 0 Δt ⎝ B ⎠i

...(4.52)

Siendo la transmisibilidad evaluada explícitamente al nivel de tiempo n, es evidente que la matriz Jacobiana generada a partir de (FLD )i es la mostrada a continuación:

[J ]LD = [T ] + [(φ B )']

...(4.53)

Se observa que [J ]LD es un caso particular de la matriz Jacobiana del método de Newton, [J ]TI , que se obtiene al eliminar la submatriz de derivadas de transmisibilidad, [T '].

58

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Dados el vector de residuos FLD y la matriz Jacobiana [J ]LD , se puede verificar que el método de linealización directa se obtiene aceptando la solución de la primera iteración, es decir: n [J ]nLD δ p (n +1) = − F LD

...(4.54)

Nótese que:

(FLD )in = Δ[T n Δp n ]i + qin +1 = 0

...(4.55)

.5.3.2 Método Semi-Implícito Linealizado. Aziz y Settari muestran que la solución que se obtiene con el método semi-implícito linealizado es equivalente a la solución generada en la primera iteración del método de Newton, o sea:

[J ]nSL δp (n +1) = − FSLn

...(4.56)

(FSL )in = (FLD )in

...(4.57)

[J ]nSL = [J ]n

...(4.58)

Donde: Y,

5. Simulación Numérica del Flujo Monofásico Multidimensional. En este capítulo se revisara la simulación numérica del flujo monofásico multidimensional en medios porosos. En primer término se trata de un flujo bidimensional y posteriormente el flujo tridimensional de una sola fase.

5.1 Flujo Monofásico Bidimensional. Considérese flujo bidimensional, en coordenadas cartesianas, de un fluido compresible en un yacimiento areal. Considérese que la distribución de presiones al tiempo cero, p ( x, y, t = 0 ), esta dada por el equilibrio gravitacional del fluido. Las fronteras del yacimiento son impermeables y existen fuentes y sumideros, pozos, en el dominio del yacimiento. La ecuación que describe el flujo del fluido en el yacimiento es: ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ∂ ⎛φ ⎞ ⎜ ⎟⎟⎥ + q* = ⎜ ⎟ − + − λ λ γ γ ⎜ ⎟ ⎢ ⎢ ⎥ ⎜ ∂x ⎣ ⎝ ∂x ∂x ⎠⎦ ∂y ⎣ ⎝ ∂y ∂y ⎠⎦ ∂t ⎝ B ⎠

...(5.1)

0 < x < Lx 0 < y < Ly t>0 Las condiciones iniciales y de frontera se expresan como:

59

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

p ( x, y,0 ) = pi ( x, y ) ;

0 < x < Lx

y

0 < y < Ly

...(5.2)

∂D ⎞ ⎛ ∂p = 0; ⎜ −γ ⎟ ∂x ⎠ x = 0, x = Lx ⎝ ∂x

0 < y < Ly ;

t >0

...(5.3)

⎛ ∂p ∂D ⎞ ⎜⎜ − γ ⎟ = 0; ∂y ⎟⎠ y = 0, y = Ly ⎝ ∂y

0 < x < Lx ;

t >0

...(5.4)

5.1.1 Aproximación de las ecuaciones diferenciales en diferencias finitas.

La aproximación en diferencias finitas de la Ec 5.1 en el nodo i, j de una malla de bloques centrados espaciada irregularmente, ver Fig. 5.1, considerando un esquema implícito, parte de la siguiente ecuación: n +1

n +1

n +1 ⎧⎪ ∂ ⎡ ⎛ ∂p ⎧ ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ⎫ ∂D ⎞⎤ ⎫⎪ n +1 ⎧ ∂ ⎛ φ ⎞ ⎫ ⎟⎥ ⎬ + ⎨ ⎢λ ⎜⎜ − γ + {q *}i, j = ⎨ ⎜ ⎟⎬ ⎟ ⎬ ⎨ ⎢λ ⎜ − γ ∂x ⎠⎥⎦ ⎭i , j ∂y ⎟⎠⎦ ⎭⎪ ⎩ ∂t ⎝ B ⎠⎭ i, j ⎩ ∂x ⎣ ⎝ ∂x ⎩⎪ ∂y ⎣ ⎝ ∂y i, j

...(5.5)

i = 1,2,...,I j = 1,2,...,J n = 1,2,... Fig. 5.1 Núcleo Típico de la malla de cálculo bidimensional de un yacimiento rectangular

i, j − 1 Δy j

i − 1, j

i, j

i + 1, j

h

i, j + 1 Δx i Para discretizar la componente del término de flujo en la dirección x, de la Ec. 5.1, se mueve uno a lo largo de y = constante; esto es, se mantiene fijo el subíndice j. De igual forma, para discretizar la componente en la dirección x, se mueve uno a lo largo de x = constante y se mantiene

60

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

fijo el subíndice i. Por similitud con el problema de flujo unidimensional, los términos de flujo de la Ec. 5.1 se escriben como: n +1 ⎧ n +1 ⎪⎛ λ x ⎞ ⎡p − − pi , j − (γΔD )i + 1 , j ⎤ ⎟ ⎨⎜ + i j 1 , ⎢⎣ 2 ⎥ ⎦ x Δ 1 ⎝ ⎠ i+ 2 , j ⎪⎩ n +1 n +1 ⎫ ⎛ λx ⎞ ⎤ ⎪ ⎡p − p ( ) γ D − Δ 1 ⎜ ⎟ ⎬ i, j i −1, j i − 2 , j ⎥⎦ ⎝ Δx ⎠ i − 1 , j ⎢⎣ ⎪⎭ 2 n +1

⎧ ∂ ⎡ ⎛ ∂p 1 ∂D ⎞⎤ ⎫ ≈ ⎟⎥ ⎬ ⎨ ⎢λ x ⎜ − γ Δxi ∂x ⎠⎦ ⎭ i, j ⎩ ∂x ⎣ ⎝ ∂x

...(5.6)

Y n +1 ⎧ n +1 ⎪⎛⎜ λ y ⎞⎟ ⎡p ⎤ ( ) − − Δ − p γ D 1 ⎨⎜ i j + i j , 1 , i j + , ⎟ ⎢ 2⎥ ⎦ ⎪⎝ Δy ⎠ i, j + 1 ⎣ 2 ⎩ n +1 n +1 ⎫ ⎛ λy ⎞ ⎡p − p ⎤ ⎪ ⎜ ⎟ ( ) − Δ γ D 1 ⎬ i, j −1 i , j − 2 ⎥⎦ ⎜ Δy ⎟ ⎢ i, j ⎝ ⎠ i, j − 1 ⎣ ⎪ 2 ⎭ n +1

⎧⎪ ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ⎫⎪ 1 ⎟⎟⎥ ⎬ ≈ ⎨ ⎢λ y ⎜⎜ − γ ∂y ⎠⎦ ⎪⎭ Δy j ⎪⎩ ∂y ⎣ ⎝ ∂y i, j

...(5.7)

La aproximación del término de acumulación es: n +1 n +1 n ⎧ ∂ ⎛ φ ⎞⎫ 1 ⎡⎛ φ ⎞ ⎛φ ⎞ ⎤ ⎢ ⎥ − ≈ ⎟ ⎜ ⎟ ⎟ ⎜ ⎜ ⎨ ⎬ Δt ⎢⎝ B ⎠ i, j ⎝ B ⎠ i, j ⎥ ⎩ ∂t ⎝ B ⎠⎭i, j ⎣ ⎦

...(5.8)

Sustituyendo las Ecs. 5.6 a 5.8 en la Ec. 5.1 y multiplicando por el volumen de roca de la celda i, j , Vri , j = h Δxi Δy j , se obtiene:

Tx n +11 ⎡ pi +1, j − pi, j − (γΔD )i + 1 , j ⎤ i + 2 , j ⎢⎣ 2 ⎥ ⎦

n +1

Ty n +1 1 ⎡ pi, j +1 − pi, j − (γΔD )i , j + 1 ⎤ i , j + 2 ⎢⎣ 2⎥ ⎦ + qin, +j 1 =

− Tx n +11 ⎡ pi , j − pi −1, j − (γΔD )i − 1 , j ⎤ i − 2 , j ⎢⎣ 2 ⎥ ⎦

n +1

n +1

− Ty n +1 1 ⎡ pi , j − pi , j −1 − (γΔD )i , j − 1 ⎤ i , j − 2 ⎢⎣ 2⎥ ⎦

+

n +1

...(5.9)

Vri , j ⎡⎛ φ ⎞ n +1 ⎛ φ ⎞ n ⎤ ⎢⎜ ⎟ −⎜ ⎟ ⎥ Δt ⎢⎝ B ⎠ i, j ⎝ B ⎠ i, j ⎥ ⎣ ⎦

i = 1,2,..., I

j = 1,2,..., J n = 1,2,... Donde T es la transmisibilidad del fluido, dada por:

⎛ h Δy ⎜ j Txi + 1 , j = T pi +1, j , pi , j = ⎜ 2 ⎜ Δx i + 1 2 ⎝

(

)

⎞ ⎟ ⎟λxi + 12 , j ⎟ ⎠

...(5.10)

61

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Y

⎛ ⎞ ⎜ h Δxi ⎟ Tyi , j + 1 = T pi, j +1 , pi, j = ⎜ ...(5.11) ⎟λyi , j + 12 2 ⎜ Δy j + 1 ⎟ 2 ⎠ ⎝ La Ec. 5.9 puede escribirse en términos de operadores, Ecs. 3.29 y 3.32, como sigue:

(

)

Tx n +11

[Δ x p − γΔ x D]in++11, j − Txin−+11, j [Δ x p − γΔ x D]in−+11, j +

Ty n +1 1 i, j + 2

[Δ y p − γΔ y D]

i+ 2, j

2

2

n +1 i , j + 12

qin, +j 1 =

− Ty n +1 1 i, j − 2

[Δ y p − γΔ y D]

2

n +1 i , j − 12

+

...(5.12)

⎛φ ⎞ Δt ⎜ ⎟ Δt ⎝ B ⎠ i, j

Vri, j

Agrupando cada dirección de flujo,

[ (

)]

Vri , j ⎛ φ ⎞ Δ x [Tx (Δ x p − γΔ x D )]in, +j 1 + Δ y Ty Δ y p − γΔ y D n +1 + qin, +j 1 = Δt ⎜ ⎟ i, j Δt ⎝ B ⎠ i, j

...(5.13)

O bien, Δ[T (Δp − γΔD )]in,+j1 + q in, +j 1 =

⎛φ ⎞ Δt ⎜ ⎟ Δt ⎝ B ⎠ i, j

Vri , j

...(5.14)

i = 1,2,..., I j = 1,2,..., J n = 0,1,2,... 5.1.2 Acoplamiento de las Condiciones de Frontera.

El acoplamiento de las condiciones de frontera, Ecs. 5.3 y 5.4 en la Ec. 5.9 se hace de manera similar al problema unidimensional (1D). La Ec. 5.9 entonces adquiere formas particulares en las celdas ubicadas a lo largo de las fronteras. Las Ecs. 5.14 constituyen un sistema algebraico de ecuaciones no lineales. Los métodos de linealización descritos para el problema unidimensional (1D) aplican también en este caso. La diferencia básicamente estriba en las características, o estructura, del sistema de ecuaciones lineales que se genera. Se puede analizar esta diferencia a través del método de Newton. 5.1.3 Método de Linealización de las Ecuaciones de Flujo en Diferencias Finitas. Método de Newton- Raphson. El Método General

Se definen entonces las siguientes funciones de residuos:

62

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

(

)

Fin, j+1 = Fi, j pi, j−1, pi−1, j , pi, j , pi+1, j , pi, j+1 n+1 = Txn+11 ⎡ pi+1, j − pi, j − (γΔD)i+1, j ⎤ i+ 2, j ⎢⎣ 2 ⎥ ⎦ Txn+11 ⎡ pi, j − pi−1, j − (γΔD)i−1, j ⎤ ⎥ i− , j ⎣⎢ 2 ⎦

n+1

2

Tyn+1 1 ⎡ pi, j − pi, j−1 − (γΔD)i, j−1 ⎤ i, j− ⎢⎣ 2⎥ ⎦ 2

+ Tyn+1 1 ⎡ pi, j+1 − pi, j − (γΔD)i, j+1 ⎤ i, j+ ⎢⎣ 2⎥ ⎦ 2

n+1

+ qin,+j1 −

n+1

n+1

−

−

...(5.15)

Vri, j ⎡⎛ φ ⎞ n+1 ⎛ φ ⎞n ⎤ ⎢⎜ ⎟ − ⎜ ⎟ ⎥ = 0 Δt ⎢⎝ B ⎠i, j ⎝ B ⎠i, j ⎥ ⎣ ⎦

El proceso iterativo resulta de expander Fin, j+1 en series de Taylor alrededor de la iteración

(ν ) y conservar solo los términos de menor orden:

(

(ν )

)

(ν ) ⎛ ∂F ⎛ ∂Fi, j ⎞ ⎟ δp(ν +1) + ⎜ i, j +⎜ i, j ⎜ ∂pi, j ⎟ ⎝ ⎠

(ν ) ⎞ ⎛ ∂F ⎟ δp(ν +1) + ⎜ i, j i+1, j

⎜ ∂pi+1, j ⎟ ⎝ ⎠

y

Por lo tanto:

⎛ ∂Fi, j ⎞ ⎜ ⎟ ⎜ ∂pi, j −1 ⎟ ⎝ ⎠

(ν )

⎛ ∂Fi, j +⎜ ⎜ ∂pi +1, j ⎝ Donde:

⎞ ⎟ ⎟ ⎠

(ν ) ⎞ ⎟ δp(ν +1) = 0 i, j +1

...(5.16)

⎜ ∂pi, j+1 ⎟ ⎝ ⎠ i = 1,2,..., I j = 1,2,..., J n = 1,2,... ν = 0,1,2,...

⎛ ∂Fi , j ⎞ ⎟ δpi(ν, j+−11) + ⎜

(ν )

(ν )

⎛ ∂Fi, j ⎞ ⎛ ∂F ⎞ ⎟ δp(ν +1) + ⎜ i, j ⎟ δp(ν +1) i, j −1 ⎜ i−1, j ⎟ ⎟ ⎝ ∂pi, j−1 ⎠ ⎝ ∂pi−1, j ⎠

Fin, j+1 = Fi, j pi, j−1, pi−1, j , pi, j , pi+1, j , pi, j+1 n+1 ≈ Fi(,νj ) + ⎜ ⎜

⎜ ∂pi −1, j ⎟ ⎝ ⎠

(ν )

⎛ ∂Fi , j ⎞ ⎟ δpi(ν−1+,1j) + ⎜

⎛ ∂Fi, j ⎞ ⎟ δpi(ν+1+,1j) + ⎜

(ν )

⎜ ∂pi, j +1 ⎟ ⎝ ⎠

(ν )

⎜ ∂pi, j ⎟ ⎝ ⎠

δpi(ν, j+1) ...(5.17)

δpi(ν, j++11) = − Fi(,νj )

δpi(ν, j+1) = pi(ν, j+1) − pi(ν, j)

...(5.18)

Escrita en forma compacta la Ec. 5.17, se tiene:

f i(,νj )δpi(ν, j+−11) + ci(ν, j)δpi(ν−1+,1j) + ai(ν, j)δpi(ν, j+1) + bi(ν, j)δpi(ν+1+,1j) + ei(ν, j)δpi(ν, j++11) = − Fi(,νj )

(f )

(c )

(a )

(b )

...(5.19)

(e )

La definición de los coeficientes f i, j , ci, j , ai , j , bi, j y ei, j se establece al igualar las Ecs. 5.17 y 5.19. 5.4.4

Estructura Matricial del Sistema Lineal de Ecuaciones.

Para un ordenamiento normal de las incógnitas y de las ecuaciones, cuando se barre la malla de cálculo secuencialmente por renglones, o bien por columnas, el sistema de ecuaciones, Ec. 5.19, genera una matriz pentadiagonal, de la siguiente forma:

63

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

e1,1 ⎡ a1,1 b1,1 ⎢c e2,1 ⎢ 2,1 a2,1 b2,1 ⎢ c3,1 a3,1 b3,1 ⎢ c4,1 a4,1 b4,1 ⎢ ⎢ c5,1 a5,1 ⎢ a1,2 b1,2 ⎢ f1,2 ⎢ f2,2 c2,2 a2,2 ⎢ f3,2 c3,2 ⎢ ⎢ f4,2 ⎢ ⎢ f5,2 ⎢ f1,3 ⎢ ⎢ f2,3 ⎢ ⎢ ⎢ ⎢ ⎣⎢

(ν )

e3,1 e4,1 e5,1 e1,2 b2,2 a3,2 b3,2 c4,2 a4,2 b4,2 c5,2 a5,2

e2,2 e3,2 e4,2 a1,3

f3,3 f4,3

b1,3

c2,3 a2,3 b2,3 c3,3 a3,3 b3,3 c4,3 a4,3 f5,3 c5,3

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ e5,2 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ b4,3 ⎥ ⎥ a5,3 ⎦⎥

(ν +1)

⎡ δp1,1 ⎤ ⎢δp ⎥ ⎢ 2,1 ⎥ ⎢δp3,1 ⎥ ⎢ ⎥ ⎢δp4,1 ⎥ ⎢δp5,1 ⎥ ⎢ ⎥ ⎢δp1,2 ⎥ ⎢δp ⎥ ⎢ 2,2 ⎥ ⎢δp3,2 ⎥ ⎢δp ⎥ ⎢ 4,2 ⎥ ⎢δp5,2 ⎥ ⎢ ⎥ ⎢δp1,3 ⎥ ⎢δp2,3 ⎥ ⎢ ⎥ ⎢δp3,3 ⎥ ⎢δp ⎥ ⎢ 4,3 ⎥ ⎣⎢δp5,3 ⎦⎥

(ν )

⎡ F1,1 ⎤ ⎢F ⎥ ⎢ 2,1 ⎥ ⎢ F3,1 ⎥ ⎢ ⎥ ⎢ F4,1 ⎥ ⎢ F5,1 ⎥ ⎢ ⎥ ⎢ F1,2 ⎥ ⎢F ⎥ ⎢ 2,2 ⎥ = −⎢F3,2 ⎥ ⎢F ⎥ ⎢ 4,2 ⎥ ⎢F5,2 ⎥ ⎢ ⎥ ⎢ F1,3 ⎥ ⎢F2,3 ⎥ ⎢ ⎥ ⎢ F3,3 ⎥ ⎢F ⎥ ⎢ 4,3 ⎥ ⎣⎢ F5,3 ⎦⎥

Como se indicó en la sección anterior, el proceso iterativo comúnmente se inicia con la solución del nivel del tiempo previo, y termina al cumplirse los criterios de convergencia equivalentes con las Ecs. 4.36 y 4.37, respectivamente.

5.5

Ordenamiento Normal.

El esquema de ordenamiento normal o estándar, por ejemplo para un problema bidimensional, consiste en ordenar los nodos de la malla de cálculo, primero en la dirección con el menor número de nodos y por último en la dirección con el mayor número de nodos, esto permite garantizar el menor ancho de banda de la matriz de coeficientes.

5.5 5.4

Ordenamiento D4. Flujo Monofásico Tridimensional.

Considérese el flujo tridimensional, en coordenadas cartesianas, de un fluido compresible en un yacimiento global. Considérese también que la distribución de presiones al tiempo cero, p( x, y, z , t = 0), esta dada por el equilibrio gravitacional del fluido. Las fronteras del yacimiento son Impermeables y existen fuentes y sumideros, pozos, en el dominio del yacimiento. Matemáticamente esto se expresa como: ∂ ⎛φ ⎞ ∂D ⎞⎤ ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ∂ ⎡ ⎛ ∂p ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ⎟⎟⎥ + ⎢λ ⎜ − γ λ⎜ − γ ⎟⎥ + ⎢λ ⎜⎜ − γ ⎟⎥ + q* = ⎜ ⎟ ⎢ ∂t ⎝ B ⎠ ∂x ⎣ ⎝ ∂x ∂x ⎠⎦ ∂y ⎣ ⎝ ∂y ∂y ⎠⎦ ∂z ⎣ ⎝ ∂z ∂z ⎠⎦ 0 < x < Lx

...(5.20)

0 < y < Ly 0 < z < Lz t>0

64

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Las condiciones iniciales y de frontera se expresan como: p(x, y, z,0) = pi(x, y, z ) ; 0 < x < Lx , 0 < y < Ly y 0 < z < Lz

...(5.21)

∂D ⎞ ⎛ ∂p = 0; ⎜ −γ ⎟ ∂x ⎠ x = 0, x = Lx ⎝ ∂x

0 < y < Ly ;

0 < z < Lz ;

t >0

...(5.22)

⎛ ∂p ∂D ⎞ ⎜⎜ − γ ⎟ = 0; ∂y ⎟⎠ y = 0, y = Ly ⎝ ∂y

0 < x < Lx ;

0 < z < Lz ;

t >0

...(5.23)

∂D ⎞ ⎛ ∂p = 0; ⎜ −γ ⎟ ∂z ⎠ z = 0, z = Lz ⎝ ∂z

0 < x < Lx ;

0 < y < Ly ;

t >0

...(5.24)

5.1.1 Aproximación de las ecuaciones diferenciales en diferencias finitas.

La Ec. 5.20 puede aproximarse siguiendo un procedimiento similar al empleado en los problemas previos. Esto es, la ecuación aproximada esta dada por: n +1

n +1

n +1

⎧⎪ ∂ ⎡ ⎛ ∂p ⎧ ∂ ⎡ ⎛ ∂p ⎧ ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ⎫ ∂D ⎞⎤ ⎫ ∂D ⎞⎤ ⎫⎪ ⎟⎟⎥ ⎬ −γ −γ + ⎨ ⎢λ ⎜ −γ + ⎨ ⎢λ ⎜⎜ ⎟⎥ ⎬ ⎟ ⎬ ⎨ ⎢λ ⎜ ∂x ⎠⎦ ⎭ i, j ,k ⎪⎩ ∂y ⎣ ⎝ ∂y ∂y ⎠⎦ ⎪⎭ ∂z ⎝ ∂z ∂z ⎠⎦⎥ ⎭ i , j ,k ⎩ ∂x ⎣ ⎝ ∂x i , j ,k ⎩ ⎣

...(5.25)

n +1

⎧ ∂ ⎛ φ ⎞⎫ + {q *}in,+j1,k = ⎨ ⎜ ⎟⎬ ⎩ ∂t ⎝ B ⎠⎭ i , j ,k i = 1,2,...,I j = 1,2,...,J k = 1,2 ,...,K n = 1,2,...

La Fig. 5.2 muestra el núcleo típico empleado en la discretización del dominio tridimensional. Fig. 5.2 Bloque Típico de celdas para un problema tridimensional.

65

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

i, j , k − 1 i, j − 1, k

i − 1, j , k

i, j, k

Δzk

i + 1, j , k

i, j + 1, k

i, j , k + 1

Δyi

Δx i Aproximación de los términos de flujo: En la dirección x, n +1 n +1 ⎧ n +1 ⎧ ∂ ⎡ ⎛ ∂p 1 ⎪⎛ λ x ⎞ ∂D ⎞⎤ ⎫ ⎡p ⎤ ( ) λ γ p γ D 1 − − Δ − − ≈ ⎜ ⎟ ⎜ ⎟ ⎬ ⎨ ⎨ ⎢ x i +1, j ,k i , j ,k i + 2 , j ,k ⎥⎦ ∂x ⎠⎥⎦ ⎭ i , j ,k Δxi ⎪⎝ Δx ⎠ i + 1 , j ,k ⎢⎣ ⎩ ∂x ⎣ ⎝ ∂x 2 ⎩ n +1 n +1 ⎫ ⎛ λx ⎞ ⎡p ⎤ ⎪ ( ) p γ D − − Δ 1 ⎜ ⎟ ⎬ i , j ,k i −1, j ,k i − 2 , j ,k ⎥⎦ ⎝ Δx ⎠ i − 1 , j ⎢⎣ ⎪⎭ 2

...(5.26)

En la dirección y , n +1 n +1 ⎧ n +1 ⎧⎪ ∂ ⎡ ⎛ ∂p ∂D ⎞⎤ ⎫⎪ 1 ⎪⎛⎜ λ y ⎞⎟ ⎡p ⎤ ⎜ ⎟ ( ) − ≈ p D λ γ γ − − Δ − 1 ⎨ ⎢ y⎜ ⎨ ⎥⎬ i , j ,k i , j + 2 ,k ⎥⎦ ⎢⎣ i , j +1,k ∂y ⎠⎟⎦ ⎪⎭ Δy j ⎪⎜⎝ Δy ⎟⎠ ⎪⎩ ∂y ⎣ ⎝ ∂y 1 i , j + 2 ,k i , j ,k ⎩

...(5.27)

n +1 ⎫ ⎛ λy ⎞ ⎪ ⎡ ⎤ ⎜ ⎟ p i , j ,k − p i , j −1,k − (γΔD )i , j − 1 ,k ⎬ ⎜ Δy ⎟ ⎢ ⎥ 2 ⎣ ⎦ ⎝ ⎠ i , j − 1 ,k ⎪ 2 ⎭ n +1

Y en la dirección z, n +1 n +1 ⎧ n +1 ⎧ ∂ ⎡ ⎛ ∂p 1 ⎪⎛ λ k ⎞ ∂D ⎞⎤ ⎫ ⎡p ⎤ ( ) p γ D ≈ − − Δ − −γ 1 ⎟ ⎟⎥ ⎬ ⎨⎜ ⎨ ⎢λ z ⎜ i , j ,k +1 i , j ,k i , j ,k + 2 ⎥⎦ ∂z ⎠⎦ ⎭ i , j ,k Δz k ⎪⎝ Δz ⎠ i , j ,k + 1 ⎢⎣ ⎩ ∂z ⎣ ⎝ ∂z 2 ⎩ n +1 n +1 ⎫ ⎛ λk ⎞ ⎪ ⎡ ⎤ pi, j ,k − pi, j ,k −1 − (γΔD )i, j ,k − 1 ⎜ ⎟ ⎬ ⎢ ⎥ 2⎦ ⎝ Δz ⎠ i, j ,k − 1 ⎣ ⎪⎭ 2

...(5.28)

66

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Aproximación del término de acumulación: n +1 n +1 n ⎤ ⎧ ∂ ⎛ φ ⎞⎫ 1 ⎡⎛ φ ⎞ ⎛φ ⎞ ⎢⎜ ⎟ ⎥ −⎜ ⎟ ≈ ⎨ ⎜ ⎟⎬ ⎩ ∂t ⎝ B ⎠⎭i, j , k Δt ⎢⎣⎝ B ⎠ i, j , k ⎝ B ⎠ i, j , k ⎥⎦

...(5.29)

Sustituyendo las Ecs. 5.26 a 5.29 en la Ec. 5.25 y multiplicando por el volumen de roca de la celda i, j , k , Vri , j , k = h Δxi Δy j Δz k , se obtiene: ⎡p − pi, j , k − (γΔD )i + 1 , j , k ⎤ ⎥⎦ i + 2 , j , k ⎢⎣ i +1, j , k 2

Tx n +11

⎡p − pi, j , k − (γΔD )i, j + 1 , k ⎤ i, j + 2 , k ⎢⎣ i, j +1, k 2 ⎥ ⎦

Ty n +1 1

Tz n +1 i,

j , k + 12

n +1

n +1

⎡p − pi, j , k − (γΔD )i, j , k + 1 ⎤ ⎢⎣ i, j , k +1 2⎥ ⎦ + qin, +j ,1k =

⎡p − pi −1, j , k − (γΔD )i − 1 , j , k ⎤ ⎥⎦ i − 2 , j , k ⎢⎣ i, j , k 2

− Tx n +11

⎡p − pi , j −1, k − (γΔD )i, j − 1 , k ⎤ ⎥⎦ i, j − 2 , k ⎢⎣ i , j , k 2

− Ty n +1 1

n +1

− Tz n +1 i,

j , k − 12

n +1

n +1

⎡p − pi, j , k −1 − (γΔD )i, j , k − 1 ⎤ ⎢⎣ i, j , k 2⎥ ⎦

+ +

...(5.30)

n +1

Vri , j , k ⎡⎛ φ ⎞ n +1 ⎛ φ ⎞ n ⎤ ⎢⎜ ⎟ ⎥ −⎜ ⎟ Δt ⎢⎝ B ⎠ i , j , k ⎝ B ⎠ i, j , k ⎥ ⎣ ⎦

i = 1,2,..., I j = 1,2,..., J k = 1,2,..., K n = 1,2,...

Donde T es la transmisibilidad del fluido, dada por: En la dirección x,

⎛ Δy Δz ⎜ j k Txi + 1 , j , k = T pi +1, j , k , pi, j , k = ⎜ 2 ⎜ Δx i + 1 2 ⎝ En la dirección y,

⎞ ⎟ ⎟λxi + 12 , j , k ⎟ ⎠

...(5.31)

⎛ ⎜ Δx i Δz k Ty i, j + 1 , k = T pi , j +1, k , pi , j , k = ⎜ 2 ⎜ Δy j + 1 2 ⎝ Y en la dirección z,

⎞ ⎟ ⎟λyi, j + 12 , k ⎟ ⎠

...(5.32)

(

(

)

)

⎛ Δx Δ y ⎜ i j Tz i, j , k + 1 = T pi, j , k +1 , pi, j , k = ⎜ 2 ⎜ Δz k + 1 2 ⎝

(

)

⎞ ⎟ ⎟λz i , j , k + 12 ⎟ ⎠

...(5.33)

La Ec. 5.30 en términos de operadores, Ecs. 3.29 y 3.32, esta dada por:

67

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Tx n +11

[Δ x p − γΔ x D]in++11, j, k − Txin−+11, j, k [Δ x p − γΔ x D]in−+11, j , k +

Ty n +1 1 i, j + 2 , k

[Δ y p − γΔ y D]

Tz n +1

[Δ z p − γΔ z D]in, +j ,1k + 1 − Tz in, +j ,1k − 1 [Δ z p − γΔ z D]in, +j,1k − 1 +

i + 2 , j, k

i,

j , k + 12

2

2

n +1 i , j + 12 , k

− Ty n +1 1 i, j − 2 , k

2

qin, +j ,1k =

Vri, j , k Δt

[Δ y p − γΔ y D]

2

n +1 i , j − 12 , k

2

+ ...(5.34)

2

⎛φ ⎞ Δt ⎜ ⎟ ⎝ B ⎠ i, j , k

Agrupando cada dirección de flujo,

[ (

)]

Δ x [Tx (Δ x p − γΔ x D )]in,+j1,k + Δ y Ty Δ y p − γΔ y D n +1 + i , j ,k Vri , j ,k ⎛φ ⎞ Δ z [Tz (Δ z p − γΔ z D )]in,+j1,k + q in, +j ,1k = Δt ⎜ ⎟ Δt ⎝ B ⎠ i , j ,k

...(5.35)

O bien, Δ[T (Δp − γΔD )]in,+j1,k + q in, +j ,1k =

Vri, j ,k Δt

⎛φ ⎞ Δt ⎜ ⎟ ⎝ B ⎠ i , j ,k

...(5.36)

i = 1,2,..., I j = 1,2,..., J k = 1,2,..., K n = 0,1,2,... 5.4.2 Acoplamiento de las Condiciones de Frontera.

El acoplamiento de las condiciones de frontera, Ecs. 5.22, 5.23 y 5.24 en la Ec. 5.30, es similar al problema unidimensional (1D). Nuevamente se enfatiza que la Ec. 5.30 adquiere formas particulares en todas aquellas celdas ubicadas a lo largo de las fronteras. Las Ecs. 5.36 forman un sistema algebraico de ecuaciones no lineales. Los métodos de linealización descritos para el problema unidimensional (1D) aplican igualmente. Para analizar la estructura de los sistemas de ecuaciones generados por los métodos de linealización, cuando se aplican al sistema no-lineal de ecuaciones 5.36, se aplicará el método de Newton. Se define entonces el método general.

5.4.3 Método de Linealización de las Ecuaciones de Flujo en Diferencias Finitas. Método de Newton- Raphson. El Método General

Se definen entonces las siguientes funciones de residuos:

68

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

(

)

Fin, j+,1k = Fi, j,k pi, j,k−1, pi, j−1,k , pi−1, j,k , pi, j,k , pi+1, j,k , pi, j+1,k , pi, j,k+1 n+1 = ⎡p − pi, j,k − (γΔD)i+1, j,k ⎤ ⎥⎦ i+ , j,k ⎢⎣ i+1, j,k 2 2

Txn+11

Tyn+1 1 ⎡ pi, j+1,k − pi, j,k − (γΔD)i, j+1,k ⎤ i, j+ 2,k ⎢⎣ 2 ⎥ ⎦ ⎡p − pi, j,k − (γΔD)i, j,k+1 ⎤ i, j,k +1 ⎢⎣ i, j,k +1 2⎥ ⎦ 2

Tzn+1

+ qin, +j,1k

n+1

n+1

n+1

⎡p − pi−1, j,k − (γΔD)i−1, j,k ⎤ ⎥⎦ i− , j,k ⎢⎣ i, j,k 2 2

−Txn+11

−Tyn+1 1 ⎡ pi, j,k − pi, j−1,k − (γΔD)i, j−1,k ⎤ i, j−2,k ⎢⎣ 2 ⎥ ⎦ ⎡p − pi, j,k−1 − (γΔD)i, j,k−1 ⎤ i, j,k −1 ⎢⎣ i, j,k 2⎥ ⎦ 2

−Tzn+1

n+1

n+1

n+1

+ +

...(5.37)

+

Vri, j,k ⎡⎛ φ ⎞n+1 ⎛ φ ⎞ n ⎤ ⎥ =0 ⎢⎜ ⎟ −⎜ ⎟ − Δt ⎢⎝ B ⎠i, j,k ⎝ B ⎠i, j,k ⎥ ⎦ ⎣

El proceso iterativo resulta de expander Fin, j+,1k en series de Taylor alrededor de (ν ) y conservar solo los términos de menor orden: (ν ) ( ( ν +1) ν ) ⎛⎜ ∂Fi, j,k ⎞⎟ n+1 δpi(ν, j+,k1)−1 + Fi, j,k = Fi, j,k pi, j,k−1, pi, j−1,k , pi−1, j,k , pi, j,k , pi+1, j,k , pi, j+1,k , pi, j,k +1 ≈ Fi, j,k +

(

)

(ν )

(ν )

(ν )

⎜ ∂pi, j,k−1 ⎟ ⎝ ⎠

⎛ ∂Fi, j,k ⎞ ⎛ ∂F ⎞ ⎛ ∂F ⎞ ⎜ ⎟ δp(ν +1) + ⎜ i, j,k ⎟ δp(ν +1) + ⎜ i, j,k ⎟ δp(ν +1) + , − 1 , − 1 , , i j k i j k i, j,k ⎜ ∂pi, j−1,k ⎟ ⎜ ∂pi−1, j,k ⎟ ⎜ ∂pi, j,k ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠

(ν )

(ν )

...(5.38)

(ν )

⎛ ∂Fi, j,k ⎞ ⎛ ∂F ⎞ ⎛ ∂F ⎞ ⎜ ⎟ δp(ν +1) + ⎜ i, j,k ⎟ δp(ν +1) + ⎜ i, j,k ⎟ δp(ν +1) = 0 i+1, j,k ⎜ ∂p i, j +1,k ⎜ ∂p i, j,k +1 ⎜ ∂pi+1, j,k ⎟ ⎟ ⎟ ⎝ ⎠ ⎝ i, j+1,k ⎠ ⎝ i, j,k+1 ⎠ i = 1,2,..., I j = 1,2,..., J k = 1,2,..., K n = 1,2,... ν = 0,1,2,...

y

Por lo tanto: (ν ) ⎛ ∂Fi, j , k ⎞ ⎜ ⎟ ⎜ ∂pi , j , k −1 ⎟ ⎠ ⎝ ⎛ ∂Fi, j , k ⎜ ⎜ ∂pi , j , k ⎝

⎞ ⎟ ⎟ ⎠

(ν )

⎛ ∂Fi , j , k ⎞ ⎟ δpi(ν, j+, k1)−1 + ⎜

⎜ ∂pi , j −1, k ⎟ ⎝ ⎠

⎛ ∂Fi, j , k ⎞ ⎟ δpi(ν, j+, k1) + ⎜ ⎜ ∂pi +1, j , k ⎟ ⎝ ⎠

⎛ ∂Fi, j , k ⎞ ⎜ ⎟ ⎜ ∂pi , j , k +1 ⎟ ⎝ ⎠

Donde:

(ν )

(ν )

(ν )

⎛ ∂Fi , j , k ⎞ ⎟ δpi(ν, j+−11), k + ⎜

⎜ ∂pi −1, j , k ⎟ ⎝ ⎠

⎛ ∂Fi , j , k ⎞ ⎟ δpi(ν+1+,1j), k + ⎜ ⎜ ∂pi , j +1, k ⎟ ⎝ ⎠

(ν )

(ν )

δpi(ν−1+,1j), k +

δpi(ν, j++11), k +

...(5.39)

δpi(ν, j+, k1)+1 = − Fi(,νj ,)k

δpi(ν, j+, k1) = pi(ν, j+, k1) − pi(ν, j), k

...(5.40)

Escrita en forma compacta la Ec. 5.17, se tiene:

69

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

hi(ν, j), k δpi(ν, j+, k1)−1 + f i(,νj ,)k δpi(ν, j+−11), k + ci(ν, j), k δpi(ν−1+,1j), k + ai(ν, j), k δpi(ν, j+, k1) + bi(ν, j), k δpi(ν+1+,1j), k + ei(ν, j), k δpi(ν, j++11), k + g i(ν, j), k δpi(ν, j+, k1)+1 = − Fi(,νj ,)k

...(5.41)

La definición de hi, j,k , fi, j,k , ci, j,k , ai, j,k , bi, j,k , ei, j,k y gi, j,k se establece al igualar las Ecs. 5.39 y 5.41. 5.4.4

Estructura Matricial del Sistema Lineal de Ecuaciones.

Para un ordenamiento normal de las incógnitas y de las ecuaciones, el sistema de Ecs. 4.41, genera en este caso una matriz de coeficientes heptadiagonal. ⎡ a1,1 b1,1 ⎢c ⎢ 2,1 a2,1 ⎢ c3,1 ⎢ ⎢ ⎢ ⎢ ⎢ f1,2 ⎢ f2,2 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ h1,3 ⎢ h2,3 ⎢ ⎢ ⎢ ⎢ ⎣⎢

e1,1 b2,1

g1,1 e2,1

a3,1 b3,1 c4,1 a4,1 b4,1 c5,1 a5,1

e3,1

f4,2

h3,3

g3,1 e4,1

g4,1 e5,1

a1,2 f3,2

g2,1

b1,2

e1,2

c2,2 a2,2 b2,2 c3,2 a3,2 b3,2 c4,2 a4,2 b4,2 f5,2 c5,2 a5,2 f1,3 a1,3 f2,3 c2,3 f3,3

h4,3

f4,3 h5,3

f5,3

e2,2 e3,2 e4,2 b1,3 a2,3 b2,3 c3,3 a3,3 b3,3 c4,3 a4,3 c5,3

(ν ) δ (ν +1) (ν ) ⎡ p1,1 ⎤ ⎡ F1,1 ⎤ ⎢F ⎥ ⎢δp ⎥ ⎢ 2,1 ⎥ ⎢ 2,1 ⎥ ⎢δp3,1 ⎥ ⎢ F3,1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢δp4,1 ⎥ ⎢ F4,1 ⎥ ⎢δp5,1 ⎥ ⎢ F5,1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢δp1,2 ⎥ ⎢ F1,2 ⎥ ⎢δp ⎥ ⎢F ⎥ ⎢ 2,2 ⎥ ⎢ 2,2 ⎥ = −⎢F3,2 ⎥ ⎢δp3,2 ⎥ ⎢δp ⎥ ⎢F ⎥ ⎢ 4,2 ⎥ ⎢ 4,2 ⎥ ⎢δp5,2 ⎥ ⎢F5,2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢δp1,3 ⎥ ⎢ F1,3 ⎥ ⎢δp2,3 ⎥ ⎢F2,3 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢δp3,3 ⎥ ⎢ F3,3 ⎥ ⎢δp ⎥ ⎢F ⎥ ⎢ 4,3 ⎥ ⎢ 4,3 ⎥ δ p ⎣⎢ 5,3 ⎦⎥ ⎣⎢ F5,3 ⎦⎥

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ g5,1 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ e5,2 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ b4,3 ⎥ ⎥ a5,3 ⎦⎥

El proceso iterativo se inicia con la solución del nivel del tiempo previo, y termina al cumplirse los criterios de convergencia equivalentes con las Ecs. 4.35 y 4.36.

5.5

Anotaciones para concluir.

Es conveniente enfatizar que independientemente del número de dimensiones y del tipo de sistema coordenado, el flujo monofásico de un fluido compresible es descrito por la ecuación diferencial generalizada, Ec. 2.62:

∇ • [λ (∇p − γ∇D )] + q* =

∂ ⎛φ ⎞ ⎜ ⎟ ∂t ⎝ B ⎠

...(5.42)

Y que la aproximación en diferencias finitas de esta ecuación puede expresarse en su forma más general, en términos de operadores de diferencias, como:

70

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Δ[T (Δp − γΔD )]in, +j ,1k + qin, +j ,1k =

⎛φ ⎞ Δt ⎜ ⎟ ⎝ B ⎠ i, j , k

Vri, j , k Δt

...(5.43)

Nótese que existe una gran similitud entre estas dos ecuaciones y que el paso de la Ec. 5.42 a la Ec. 5.43 es directo. Esto es, en el término de flujo se reemplaza el operador diferencial ∇ por el operador de diferencias centrales Δ y la mobilidad λ por la transmisibilidad T. El término fuente, que se expresa en la ecuación diferencial como gasto por unidad de volumen de roca, deberá ser simplemente el gasto del pozo ubicado en la celda. En el término de acumulación se reemplaza la derivada parcial con respecto al tiempo, ∂ ∂t , por el operador de diferencias regresivas en tiempo, Δ t , y a su vez este se multiplica por el volumen de roca de la celda y se divide entre Δt. Se observa que independientemente del número de dimensiones y del sistema de coordenadas empleados, la transmisibilidad Ti + 1 , entre las celdas i e i + 1 por ejemplo, es el 2 producto entre la mobilidad del fluido, evaluada en la frontera de las celdas, λi + 1 , y el área 2

expuesta al flujo entre esas celdas, Ai + 1 , dividida por la distancia que separa a los nodos i e i + 1 , 2

Δη i + 1 , esta es: 2

⎛ A ⎞ ⎟⎟ λi + 1 Ti + 1 = T ( pi +1 , pi ) = ⎜⎜ 2 2 ⎝ Δη ⎠ i + 1

...(5.44)

2

Es evidente que la forma de calcular el área expuesta al flujo y la mobilidad, en la frontera de las celdas, dependerá del sistema de coordenadas, como lo indican las Ecs. 4.23 y 4.24. Nótese que la ecuación diferencial, Ec. 5.42, es no lineal y que el sistema algebraico de ecuaciones aproximadas, 5.43, también lo es. Se mostró que el método iterativo de NewtonRaphson es el método general para resolver el sistema no lineal de ecuaciones. Métodos como linealización directa, aproximaciones sucesivas y semi-implícito linealizado pueden obtenerse como casos particulares del método de Newton, prelinealizando parcialmente las ecuaciones y controlando el número de iteraciones. Se observa también que la estructura matricial del sistema lineal de ecuaciones, generado por los métodos de linealización, depende del número de dimensiones del problema. Así para un ordenamiento normal de las ecuaciones y de las incógnitas se tienen que el problema unidimensional genera un sistema de ecuaciones tridiagonal, el problema bidimensional genera un sistema pentadiagonal y el problema tridimensional genera un sistema heptadiagonal.

6. Modelos de Pozos. 6.1 6.2 6.3

Modelo de pozo de Peaceman para un yacimiento homogéneo e isotrópico. Modelo de pozo de Peaceman para un yacimiento heterogéneo y anisotrópico. Modelo de Hoo-Jeen Su para pozos fuera de centro.

71

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

7. Solución de Sistemas de Ecuaciones Lineales. En este capítulo se presenta una breve discusión de los métodos de solución comúnmente empleados en la simulación numérica de yacimientos. Básicamente los sistemas de ecuaciones, como los que se generan en la simulación de los yacimientos, pueden resolverse mediante métodos directos o iterativos. Los Métodos Directos requieren un número fijo de operaciones en la solución de un sistema dado. En los Métodos Iterativos no es posible predeterminar el número de operaciones, ya que la solución consiste en aplicar algoritmos cíclicos, en espera de que al paso de las iteraciones se obtenga una mejor aproximación de la solución.

7.1 Métodos Directo.

La solución numérica de los problemas del flujo de fluidos a través de medios porosos, generan sistemas lineales de ecuaciones que en forma matricial conforman una matriz dispersa con una estructura bandada, en la cual los elementos diferentes de cero están sobre una banda de ancho predecible. En los métodos directos, el trabajo computacional y el requerimiento de memoria están directamente relacionados con el número de ecuaciones o de incógnitas a resolver. Entre mayor sea este número, el ancho de banda del sistema matricial de ecuaciones es generalmente mayor y la cantidad de memoria de cómputo y de tiempo de procesamiento son también mayores. El número de ecuaciones ó de incógnitas esta en función del número total de bloques o de celdas, en que se discretiza el dominio de interés y del número de componentes considerados. Los métodos directos, consisten de un proceso de eliminación hacia adelante, que transforma a la matriz de coeficientes en una matriz triangular superior, y de un proceso de sustitución hacia atrás, donde se obtiene la solución. El método de eliminación Gaussiana es el método directo básico. Existen una serie de variaciones de este método, motivadas por las peculiaridades de algunos sistemas de ecuaciones y en su solución eficiente. Dentro de estos están, los Algoritmos de Banda y las Técnicas de Matrices Dispersas son de interés en Simulación Numérica de Yacimientos. Los algoritmos de banda toman ventaja de la estructura bandada de la matriz para confinar las operaciones de eliminación dentro de la banda. Las técnicas de matrices dispersas, reconociendo la existencia de un gran número de ceros en la matriz de coeficientes, evitan el almacenamiento de esos ceros y operaciones innecesarias. El trabajo requerido en los algoritmos de banda esta directamente relacionado con el ancho de la banda. El menor ancho de banda se obtiene cuando el ordenamiento de las incógnitas y de las ecuaciones se hace en la dirección del menor número de celdas. La eficiencia de los métodos directos se mejora grandemente cuando se hace uso de técnicas especiales de ordenamiento del sistema lineal de ecuaciones, como es el caso del ordenamiento D4. A continuación se presenta el algoritmo de Thomas empleado en la solución de sistemas tridiagonales de ecuaciones generados típicamente en la solución numérica de un problema unidimensional monofásico, o equivalentes. Posteriormente se presentará el método de

72

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

descomposición LU que resuelve sistemas de ecuaciones con matrices arbitrariamente estructuradas.

7.1.1

Solución de Sistemas Tridiagonales: Algoritmo de Thomas.

Considérese el siguiente sistema tridiagonal de ecuaciones:

a1u1 + b1u 2 = d1 ci ui −1 + ai u i + bi ui +1 = d i ; c I u I −1 + a I u I = d I O bien,

i = 2,3,..., I − 1 ...(7.1)

Au = d

...(7.2)

Donde,

⎡ a1 ⎢c ⎢ 2 ⎢ ⎢ A=⎢ ⎢ ⎢ ⎢ ⎢ ⎣

b1 a2

b2

c3

a3

b3

.

. .

. .

.

.

.

⎤ ⎡ u1 ⎤ ⎡ d1 ⎤ ⎥ ⎢u ⎥ ⎢d ⎥ ⎥⎢ 2 ⎥ ⎢ 2 ⎥ ⎥ ⎢u 3 ⎥ ⎢ d 3 ⎥ ⎥⎢ ⎥ ⎢ ⎥ ⎥⎢ . ⎥ = ⎢ . ⎥ ⎥⎢ . ⎥ ⎢ . ⎥ ⎥⎢ ⎥ ⎢ ⎥ . ⎥⎢ . ⎥ ⎢ . ⎥ a I ⎥⎦ ⎢⎣u I ⎥⎦ ⎢⎣ d I ⎥⎦

cI

Así como u T = (u1 , u 2 ,..., u I ) y d T = (d1 , d 2 ,..., d I ). Se descompone la matriz A en el producto de las matrices W y Q,

A = WQ

...(7.3)

Donde W y Q se define como sigue:

⎡ w1 ⎢c ⎢ 2 ⎢ ⎢ W =⎢ ⎢ ⎢ ⎢ ⎢ ⎣

w2 c3

w3 .

. .

. .

. cI

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ wI ⎥⎦

73

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

⎡ 1 ⎢ ⎢ ⎢ ⎢ Q=⎢ ⎢ ⎢ ⎢ ⎢ ⎣

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ q I −1 ⎥ 1 ⎥⎦

q1 1

q2 1

q3 .

. .

. 1

Efectuando el producto de las matrices WQ e igualando la matriz resultante, elemento a elemento, con la matriz A, se obtiene un sistema de 2I −1 ecuaciones con 2I −1 incógnitas: w1, w2 ,...,wI y q1, q2 ,...,qI −1. Su solución proporciona las siguientes definiciones para los elementos de las matrices Q y W : w1 = a1

b q1 = 1 w1 b qi = i ; wi wi = u i − ci qi −1 ;

i = 2,3,..., I

...(7.4)

i = 2,3,..., I

...(7.5)

La solución del sistema tridiagonal se obtiene como sigue:

Au − WQu = Wg Donde g = Qu . Es decir, se tiene: ⎡ w1 ⎢c ⎢ 2 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

w2 c3

w3

.

. .

. .

. cI

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ wI ⎥⎦

...(7.6)

(ν )

⎡ g1 ⎤ ⎢g ⎥ ⎢ 2⎥ ⎢ g3 ⎥ ⎢ ⎥ ⎢ . ⎥ ⎢ . ⎥ ⎢ ⎥ ⎢ . ⎥ ⎢g ⎥ ⎣ I⎦

(ν +1)

⎡ d1 ⎤ ⎢d ⎥ ⎢ 2⎥ ⎢d3 ⎥ ⎢ ⎥ =⎢ . ⎥ ⎢ . ⎥ ⎢ ⎥ ⎢ . ⎥ ⎢d ⎥ ⎣ I⎦

De donde se obtiene en un barrido hacia delante la solución de g, d g1 = 1 w1

d − ci g i −1 gi = i ; wi

(ν )

...(7.7)

i = 2,3,..., I

74

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Finalmente se puede resolver Qu = g como sigue: ⎡ 1 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

q1 1

q2 1

q3 .

. .

. 1

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ q I −1 ⎥ 1 ⎥⎦

(ν )

⎡ u1 ⎤ ⎢u ⎥ ⎢ 2⎥ ⎢u 3 ⎥ ⎢ ⎥ ⎢ . ⎥ ⎢ . ⎥ ⎢ ⎥ ⎢ . ⎥ ⎢u ⎥ ⎣ I⎦

(ν +1)

⎡ g1 ⎤ ⎢g ⎥ ⎢ 2⎥ ⎢ g3 ⎥ ⎢ ⎥ =⎢ . ⎥ ⎢ . ⎥ ⎢ ⎥ ⎢ . ⎥ ⎢g ⎥ ⎣ I⎦

(ν )

Y en un barrido hacia atrás se obtiene la solución deseada:

uI = g I u i = g i − q i u i +1 ;

i = I − 1, I − 2,...,1

...(7.8)

Resumiendo, para resolver el sistema tridiagonal, Au = d , es necesario primero calcular los elementos qi , wi y gi mediante las Ecs. 7.4, 7.5 y 7.7 y posteriormente obtener la solución de u i mediante la Ec. 7.8. 7.1.2

Algoritmos de matrices dispersas: Descomposición LDU

7.2 Métodos Iterativos.

Los métodos iterativos adquieren interés en el caso de problemas multidimensionales, donde el número de ecuaciones a resolver es relativamente grande. Básicamente existen dos tipos: 1.- Métodos iterativos por punto, (MIP). 2.- Métodos iterativos por línea o bloque, (MIL/B). En los métodos iterativos por punto, MIP, las incógnitas de cada nodo son resueltas explícitamente en cada iteración. Por otro lado, los métodos iterativos por línea/bloque, MIL/B, se resuelven simultáneamente las incógnitas de un grupo de celdas. Entre más implícito sea un método iterativo, esto es entre mayor sea el número de incógnitas que se resuelven simultáneamente, la convergencia a la solución es más rápida. Esto requiere sin embargo una mayor capacidad de memoria de cómputo y un mayor esfuerzo computacional por iteración. Lo anterior indica que en la implementación de un método iterativo por bloques, es necesario hacer un balance entre implicitud y simplicidad para resolver el sistema de ecuaciones generado por los bloques. Los métodos iterativos están ligados a parámetros de iteración, que se introducen en el algoritmo para acelerar la convergencia a la solución. El ritmo de convergencia de un algoritmo dado, depende de varios factores, dentro de los que se pueden mencionar: las características de la malla de cálculo, la anisotropía y las heterogeneidades de la formación, la estimación inicial de la solución y los criterios de convergencia.

75

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Dentro de los métodos iterativos que han recibido atención en la simulación numérica de yacimientos se tienen a los métodos PSOR, LSOR, BSOR y algunas variantes de estos, así como los métodos ADIP y SIP. Ninguno de estos métodos es universal: algunos funcionan bien bajo algunas circunstancias y mal bajo otras. A continuación se presentan algunos de estos métodos. Para ejemplificar los métodos iterativos, considérese el siguiente sistema de ecuaciones:

f i, j ui, j −1 + ci, j ui −1, j + ai, j ui, j + bi, j ui +1, j + ei, j u i, j +1 = d i, j i = 1,2,..., I

...(7.9)

j = 1,2,..., J Como antecedentes a los métodos SOR se presentan primeramente los métodos de Jacobi y de Gauss Seidel. 7.2.1

Métodos de Jacobi.

La solución del sistema de ecuaciones, Ec. 7.9, mediante el método de Jacobi consiste en la aplicación sucesiva del siguiente proceso iterativo:

ui(,mj +1) =

1 ⎧ (m ) (m ) (m ) (m ) ⎫ ⎨d i, j − f i, j u i , j −1 − ci, j u i −1, j − bi, j u i +1, j − ei, j u i , j +1 ⎬ ai, j ⎩ ⎭

...(7.10)

Donde m + 1 y m son los niveles iterativos, desconocido y conocido respectivamente. Nótese que la solución de las incógnitas se obtienen puntualmente, o explícitamente: el orden en que estas se resuelven es irrelevante. Como cualquier proceso iterativo, se requiere de una estimación inicial de la solución y se espera que en la medida en que las iteraciones avancen, el método converja a la solución. Esto ocurre cuando los cambios iterativos de las incógnitas en todos los nodos sean, en valor absoluto, menor que una tolerancia preestablecida: u im, j+ 1 − u im, j ε , en todos los nodos i, j . 7.2.2

Métodos de Gauss-Seidel.

Se observa que si se establece un orden en la solución de las incógnitas mediante el método de Jacobi, el ordenamiento normal por ejemplo, en la medida que se avanza la solución, existe ya en los nodos barridos una mejor estimación de la solución. El método de Gauss-Seidel reconoce este hecho y modifica el algoritmo iterativo de Jacobi de la siguiente manera:

ui(,mj +1) =

1 ⎧ (m +1) (m +1) (m ) (m ) ⎫ ...(7.11) ⎨d i , j − f i, j u i , j −1 − ci, j u i −1, j − bi, j u i +1, j − ei, j u i , j +1 ⎬ ai , j ⎩ ⎭ Donde el barrido se hace por renglones. Esto es, en orden creciente de j y para cada renglón j se barre en orden creciente de i . El método de Gauss-Seidel es computacionalmente más simple que el método de Jacobi y converge más rápido a la solución.

76

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

7.2.3

Método de SOR en Punto o PSOR.

El método iterativo del PSOR se basa en una modificación del método de Gauss-Seidel, donde se introduce el uso de un parámetro de iteración, ω , con el objeto de acelerar el proceso de

(

)

convergencia a la solución. Si se denomina u i*, jm +1 a la solución obtenida mediante el método de Gauss-Seidel en la iteración m + 1, el método PSOR se define como:

ui(,mj +1) = (1 − ω )ui(,mj ) + ωui*,(jm +1)

...(7.12)

1≤ ω ≤ 2

Donde ω es el parámetro de sobrerelajación que, como ya se mencionó, acelera el ritmo de convergencia de la solución. El parámetro ω adquiere valores en el rango de 1 a 2. Más adelante se daran algunos criterios para la estimación de este parámetro.

7.2.4

Método de SOR en Línea o LSOR.

Se establecen, en cada iteración, sistemas de ecuaciones similares a los generados en problemas lineales. De esta manera se resuelve simultáneamente las incógnitas correspondientes a una línea o columna de celdas, que son más fáciles de resolver. En la solución de las incógnitas de una línea determinada, se emplean los nuevos valores de las incógnitas obtenidos en líneas anteriores. Como una extensión del método de Gauss-Seidel, para cada j se establece el siguiente sistema de ecuaciones:

ci, j u *(m +1) + ai, j u *(m +1) + bi, j u *(m +1) = d i, j − f i, j u (m +1) − ei, j u (m ) i −1, j

i, j

i +1, j

i = 1,2,..., I

i , j −1

i , j +1

...(7.13)

para cada j

El sistema de ecuaciones generado con la Ec. 7.13 es tridiagonal. Se puede entonces resolver mediante el algoritmo de Thomas. Una vez que se obtiene la solución de las incógnitas u i*,(jm +1) , i = 1,2,..., I correspondientes a la línea j , esta se sobrerelaja como la Ec. 7.12, esto es:

u i(,mj +1) = (1 − ω )u i(,mj ) + ωui*,(jm +1)

1≤ ω ≤ 2

...(7.14)

El ritmo de convergencia del método iterativo LSOR, depende del valor del parámetro de sobrerelajación, ω. Existe un valor óptimo de este parámetro, el cual se puede obtener por ensaye y error o mediante algún procedimiento o algoritmo, como se verá posteriormente. El proceso iterativo LSOR comienza con la siguiente estimación inicial:

77

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

δui(,0j) = δuin, j

...(7.15)

Y termina cuando los cambios iterativos de las incógnitas son, en valor absoluto, menores que una cierta tolerancia estipulada, esto es:

δui(η, j+1) − δui(η, j) ≤ tolerancia para toda i = 1,2,...,I

...(7.16)

j = 1,2,...,J 7.2.5

Método de SOR en Bloque o BSOR.

La idea principal de este método consiste en rearreglar por planos el sistema de ecuaciones generado por un problema tridimensional: hi, j,kui, j,k−1 + fi, j,kui, j−1,k +ci, j,kui−1, j,k + ai, j,kui, j,k +bi, j,kui+1, j,k + ei, j,kui, j+1,k + gi, j,kui, j,k+1 = di, j

...(7.17)

i = 1,2,..., I j = 1,2,..., J k = 1,2,..., K Es decir, el problema tridimensional, Ec. 7.17, se resuelve mediante una serie de barridos bidimensionales, empleando los nuevos valores de las incógnitas de los planos previamente resueltos. El método de sobrerelajación en Bloques, BSOR, aplicado a la solución del sistema lineal de ecuaciones definido por la Ec. 7.17, consiste en conservar en el lado izquierdo de la ecuación a los términos correspondientes a dos direcciones así como a la diagonal principal y pasar al lado derecho los términos correspondientes a una dirección y a los términos de residuos, con lo que se establece el siguiente proceso iterativo:

(m+1) (m) *(m+1) *(m+1) *(m+1) +1) fi, j,kui*,(jm−+1,1k) + ci, j,kui*−(m 1, j,k + ai, j,kui, j,k +bi, j,kui+1, j,k + ei, j,kui, j+1,k = di, j − hi, j,kui, j,k−1 − gi, j,kui, j,k+1

...(7.18)

El problema bidimensional, Ec. 7.18, genera para un ordenamiento normal, una matriz de coeficientes pentadiagonal que se resuelve más eficientemente que el problema original. Cada uno de estos problemas reducidos, se resuelve empleando un esquema directo de solución, con lo que se genera una solución intermedia que sirve de base para obtener la solución al nivel de iteración desconocida, (m + 1). Una vez que se obtiene la solución intermedia de las incógnitas u i*(, jm, k+ 1) , i = 1,2,..., I y j = 1,2,..., J correspondientes al plano k , ésta se sobrerelaja como la Ec. 7.14.

En problemas donde el método LSOR no converge debido, por ejemplo, que la formación es altamente heterogénea y anisotrópica, el método BSOR tiene un mejor comportamiento. Las propiedades de convergencia de estos métodos iterativos mejoran a medida que el número de ecuaciones resueltas simultáneamente aumenta. Esto; sin embargo, ocasiona un aumento en el trabajo computacional realizado en cada iteración. Al igual que el método anterior, el método BSOR, comienza con una estimación inicial dada por la Ec. 7.15 y termina cuando alcanza la convergencia estipulada por la Ec. 7.16.

78

FI-UNAM

7.2.6

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Algoritmo para el cálculo del parámetro de sobrerelajación.

Una forma de seleccionar el valor óptimo del parámetro de sobrerelajación, ω , de los métodos iterativos LSOR y BSOR es por ensaye y error. Es decir, se realizan corridas de sensibilidad con diferentes valores de este parámetro y se obtiene una gráfica de los valores de ω empleados contra el número de iteraciones requeridas para resolver el problema en cuestión. La Fig. 7.1 representa una gráfica típica del comportamiento de este parámetro ω , donde se observa su valor óptimo, obtenido con el menor número de iteraciones realizado. 130

No. de Iteracione s

120

110

100

90

ωb = 1.15

80

70 0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

Parametro de Sobrerelajación

Fig. 7.1 Parámetro de Sobrerelajación Vs. Número de Iteraciones.

Otra manera de seleccionar el valor óptimo de ω es a través del procedimiento que a continuación se describe: Partiendo de considerar que cualquier esquema iterativo para resolver el sistema de ecuaciones puede expresarse de la siguiente forma: ...(7.19) x (η +1) = Ax (η ) + b Donde η + 1 es el nivel iterativo desconocido y η es el nivel iterativo conocido. El valor óptimo de ω está en función del radio espectral de la matriz A, ρ ( A), de la siguiente manera: 2 ...(7.20) ωb = 1 2 2 1 + 1 − ρ ( A) Donde ωb es el valor óptimo del parámetro de sobrerelajación.

[

]

El radio espectral de la matriz A, ρ ( A), está definido como:

ρ ( A) =

(θ + ω − 1) ωθ

1

...(7.21)

2

79

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Donde θ es el valor promedio del ritmo al cual converge la solución, θ (η ) , el cual se calcula con la siguiente expresión:

θ

(η )

=

δx (η +1)

...(7.22)

δx (η )

Donde δx (η ) es el vector de los cambios iterativos de la solución, definido como:

δx (η ) = x (η +1) − x (η )

...(7.23)

Y δx es la norma del vector solución de los cambios iterativos.

Los valores de θ (η ) deben estar muy próximos a un valor límite de θ y este valor se usa como una aproximación del radio espectral. También se requiere un valor inicial de ω en la Ec. 7.21 para obtener el vector solución δx (η ) . Por tal motivo, se requiere que este proceso iterativo se inicialice con ω = 1 y se cumpla con la siguiente condición:

θ (η +1) − θ (η ) ≤ tolerancia

...(7.24)

De tal manera que este proceso iterativo se detiene hasta lograr la convergencia estipulada, Ec. 7.24. Cabe mencionar, que este proceso converge solo si el valor actual de ω es menor que el nuevo valor de ω calculado. De otra forma, los valores de θ (η ) oscilarán y no se podrá obtener una aproximación correcta del radio espectral, ρ ( A). Este último caso, ω debe reducirse hasta que

θ (η ) logre la convergencia.

7.3

7.2.7

Procedimiento Iterativo Implícito de Dirección Alternada o ADIP.

7.2.8

Procedimiento Fuertemente Implícito o SIP.

Solución de Sistemas de Ecuaciones Estructurados en Bloques.

Los métodos directos e iterativos descritos previamente se extienden en forma natural a la solución de sistemas estructurados en bloques. La diferencia estriba en que en lugar de operar con escalares, como hasta ahora se ha mostrado, debe operarse matricialmente. Obviamente, el orden de las operaciones es importante en este caso. Se toma como ejemplo la extensión del algoritmo de Thomas a la solución de sistemas tridiagonales en bloques, como los generados en la aplicación del método de Newton a la linealización de las ecuaciones de flujo multifásico en diferencias. Se puede representar el sistema como: a1u1 + b1u 2 = d1

ci ui −1 + ai ui + bi ui +1 = d i ; c I u I −1 + a I u I = d I O bien,

Au = a

i = 2,3,..., I - 1 ...(7.25) ...(7.26)

Donde,

80

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

⎡ a1 ⎢c ⎢ 2 ⎢ ⎢ A=⎢ ⎢ ⎢ ⎢ ⎢ ⎣

b1 a2 c3

b2 a3

b3

.

.

. .

.

.

.

. cI

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ . ⎥ a I ⎥⎦

Así como u T = (u1 , u 2 ,..., u I ) y d T = (d1 , d 2 ,..., d I ) . Los elementos ai , bi y ci son para el caso de flujo de los pseudocomponentes aceite, gas y agua, submatrices de (3 x3), como se definen en el problema multifásico. Los elementos ui y d i son equivalentemente subvectores de orden 3. En forma similar con el procedimiento establecido anteriormente, la solución del sistema tridiagonal en bloque consiste en descomponer la matriz A en el producto WQ. Por lo que se establece el sistema g = Qu, como paso intermedio de la solución. En este caso, las matrices W y Q se definen como sigue: ⎡ ⎡ w1 w2 w3 ⎤ ⎢⎢ ⎥ ⎢ ⎢ w4 w5 w6 ⎥ ⎢ ⎢⎣ w7 w8 w9 ⎥⎦ 1 ⎢ ⎢ ⎡ c1 c 2 c3 ⎤ ⎢⎢ ⎥ ⎢ ⎢c 4 c 5 c 6 ⎥ ⎢ ⎢⎣c 7 c8 c9 ⎥⎦ 2 ⎢ ⎢ ⎢ ⎢ ⎢ W =⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣

⎡ w1 ⎢ ⎢ w4 ⎢⎣ w7

w2 w5

⎡ c1 ⎢c ⎢ 4 ⎢⎣c 7

c2 c5

w8

c8

w3 ⎤ w6 ⎥⎥ w9 ⎥⎦ c3 ⎤ c 6 ⎥⎥ c9 ⎥⎦

2

3

⎡ w1 ⎢w ⎢ 4 ⎢⎣ w7

w2 w5 w8

w3 ⎤ w6 ⎥⎥ w9 ⎥⎦

3

⎡ c1 ⎢c ⎢ 4 ⎢⎣c 7

c2 c5 c8

c3 ⎤ c 6 ⎥⎥ c9 ⎥⎦ I

⎡ w1 ⎢w ⎢ 4 ⎢⎣ w7

w2 w5 w8

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ w3 ⎤ ⎥ ⎥ w6 ⎥⎥ ⎥ ⎥ w9 ⎥⎦ I ⎥⎦

Y,

81

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

⎡ 1 ⎢ ⎢ ⎢ ⎢ Q=⎢ ⎢ ⎢ ⎢ ⎢ ⎣

q1 1

q2 1

q3 .

. .

. 1

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ q I −1 ⎥ 1 ⎥⎦

Donde I son matrices identidad de orden (3 x3),

⎡1 0 0⎤ I = ⎢⎢0 1 0⎥⎥ ⎢⎣0 0 1⎥⎦ Las submatrices wi , qi y el subvector g i se obtienen como: w1 = a1

q1 = w1−1b1 qi = wi−1bi ; wi = ai − ci qi −1 ;

i = 2,3,..., I

...(7.27)

i = 2,3,..., I

Y,

g1 = w1−1d1 g i = wi−1 (d i − ci g i −1 ) ;

i = 2,3,..., I

...(7.28)

La solución de u se obtiene en un barrido hacia atrás como sigue:

uI = g I ui = g i − qi ui +1 ; 7.4

i = I - 1, I - 2,..., 1

...(7.29)

Métodos Directos Versus Métodos Iterativos.

En los métodos directos convencionales, eliminación Gaussiana o descomposición LU por ejemplo, el requerimiento de memoria y de tiempo de cómputo aumentan drásticamente en la medida que aumenta el número de ecuaciones. Por esta razón, los métodos directos son generalmente preferidos cuando el número de ecuaciones no es muy grande. ¿Qué es grande y qué es pequeño?: es una pregunta difícil de contestar. Con el desarrollo de las denominadas supercomputadoras el calificativo puede ser demasiado subjetivo. Lo que antes era un sistema grande, puede ahora considerársele un sistema pequeño. Es también indudable que con el desarrollo de los métodos de matrices dispersas y de técnicas especiales de ordenamiento, los métodos directos tienen ahora un mayor potencial de aplicación: el tamaño de los sistemas de ecuaciones que puede resolverse eficientemente mediante estos métodos es cada vez mayor, y lo seguirá siendo en la medida en que los sistemas de computo sean más poderosos. Una desventaja de los métodos iterativos es la falta de metodologías para el cálculo de los parámetros de iteración óptimos. Los

82

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

métodos existentes para el cálculo de estos parámetros han sido desarrollados bajo una serie de simplificaciones. En algunas situaciones los parámetros teóricos no solo no resuelven eficientemente el problema, sino que provocan problemas de convergencia en el método.

8. Introducción a la Simulación Numérica del Flujo Multifásico. En este capítulo se revisaran los fundamentos de la simulación numérica del flujo multifásico en yacimientos. Aprovechando la experiencia obtenida en el capítulo 4, y habiendo ejercitado el uso de operadores de diferencias, seremos concisos en la escritura de las ecuaciones aproximadas. Se verá como el problema multifásico ofrece una serie de posibilidades, a de formulaciones, para el tratamiento del sistema no lineal de las ecuaciones de flujo en diferencias. Como en el problema monofásico, primeramente se revisará el problema de flujo unidimensional y posteriormente se extenderán los métodos al problema multidimensional. 8.1 Flujo Multifásico Unidimensional de Fluidos Tipo Beta.

Considérese el flujo lineal multifásico isotérmico de fluidos tipo beta: aceite, gas y agua en un yacimiento petrolero. Las distribuciones de presiones y de saturaciones al tiempo cero, p p (x, t = 0) y Sp(x,t = 0) , obedecen al equilibrio gravitacional y capilar, como se describe en el capitulo 2. Se considera que no existen fuentes distribuidas a lo largo del dominio del yacimiento. El yacimiento produce a un gasto de aceite qo en la frontera situada en x = 0 y en la frontera en x=L esta cerrada al flujo. Las ecuaciones que describen el flujo de fluidos en el yacimiento son: Para el componente aceite, ∂ ⎡ ⎛ ∂po ∂D ⎞⎤ ∂ ⎛ φS o −γo ⎟⎥ = ⎜ ⎢λo ⎜ ∂x ⎣ ⎝ ∂x ∂x ⎠⎦ ∂t ⎜⎝ Bo

⎞ ⎟⎟ ⎠

...(8.1)

Para el componente gas, ∂p φS o ⎞⎟ ∂D ⎞⎤ ∂ ⎛⎜ φS g ∂ ⎡ ⎛ ∂p g ∂D ⎞⎤ ∂ ⎡ ⎟⎥ + ⎢λo Rs ⎜⎛ o − γ o −γ g + = R ⎟ ⎢λ g ⎜⎜ s ⎥ ∂x ⎢⎣ ⎝ ∂x ∂x ⎟⎠⎥⎦ ∂x ⎣ ∂x ⎠⎦ ∂t ⎜⎝ B g Bo ⎟⎠ ⎝ ∂x

...(8.2)

La ecuación del agua, ∂ ⎡ ⎛ ∂p w ∂D ⎞⎤ ∂ ⎛ φS w ⎞ ⎟ −γw ⎟⎥ = ⎜ ⎢λ w ⎜ ∂x ⎣ ⎝ ∂x ∂x ⎠⎦ ∂t ⎜⎝ B w ⎟⎠

...(8.3)

0< x0

Con las siguientes relaciones adicionales: So + S g + S w = 1

…(8.4)

83

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Pc go = p g − p o

…(8.5)

Pc wo = p o − p w

…(8.6)

y

Las condiciones iniciales, como se mencionó previamente, corresponden a un equilibrio gravitacional y capilar de las fases en el yacimiento. Esto genéricamente se escriben como: Pp ( x,0 ) = p p10 ( x ) ;

0<x
…(8.7)

S p ( x,0 ) = S p10 ( x ) ;

0<x
…(8.8)

Donde p = o, g , w Las condiciones de frontera se expresan, en x = 0, como: ⎛μ B ⎞ ∂D ⎞ ⎛ ∂p o −γo = −⎜⎜ o o ⎟⎟ qo ; ⎜ ⎟ ∂x ⎠ x =0 ⎝ ∂x ⎝ kkro A ⎠ x =0

t〉0

...(8.9)

Los gastos de producción del gas libre y del agua pueden obtenerse, con algunas simplificaciones, mediante las siguientes relaciones: ⎛ λg ⎞ ⎟ q g = ⎜⎜ qo ⎟ λ o ⎝ ⎠ x =0

…(8.10)

⎛λ ⎞ q w = ⎜⎜ w ⎟⎟ qo ⎝ λ o ⎠ x =0

…(8.11)

Y

Nótese que el gasto total de producción de gas, libre más disuelto, esta dado por: ⎛ λg ⎞ + R s ⎟⎟ q gT = q g + q o (Rs ) x =0 = ⎜⎜ qo ⎝ λo ⎠ x =0

…(8.12)

En la frontera en x = L, la condición de cero flujo para todos los fluidos, p = o, g , w, se expresa como: ⎛ ∂p p ∂D ⎞⎟ ⎜ =0 ; ⎜ ∂x − γ p ∂x ⎟ ⎝ ⎠ x= L

t〉0

...(8.13)

8.2 Forma de las Ecuaciones de Fluidos en Diferencias Finitas.

La aproximación en diferencias finitas de las Ec. 8.1 a 8.3 en el nodo i, mediante un esquema implícito se expresa, a través de operadores en diferencias, como:

84

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Para el componente aceite, Δ[To (Δp o − γ o ΔD )]in +1 =

Para el componente gas,

[ (

)]

(

⎛φ 1 − Sg − Sw Vri Δ t ⎜⎜ Δt Bo ⎝

)⎞⎟

...(8.14)

⎟ ⎠i

Δ T g Δp o + ΔPc go − γ g ΔD n +1 + Δ[To R s (Δp o − γ o ΔD )]in +1 = i

(

Vri ⎛⎜ φS g φR s 1 − S g − S w Δt + ⎜ Bg Δt Bo ⎝

)⎞⎟

...(8.15)

⎟ ⎠i

Y para el componente agua, Δ[Tw (Δpo − ΔPc wo − γ w ΔD )]in +1 =

i = 1,2,..., I n = 0,1,2,...

Vri ⎛ φS w ⎞ ⎟ Δt ⎜ Δt ⎝⎜ Bw ⎟⎠ i

...(8.16)

Nótese que las relaciones adicionales, Ecs. 8.4 a 8.6, fueron acopladas en las ecuaciones de flujo en diferencias. Con esto se tiene ahora en cada nodo tres ecuaciones con tres incógnitas: p o , S g , S w , i = 1,2,..., I.

(

)i

8.3 Cálculo de las Transmisibilidades en las Fronteras.

Las transmisibilidades, Tp , p = o, g, w se definen similarmente al caso monofásico, excepto que la permeabilidad absoluta es reemplazada por la permeabilidad efectiva de la fase en cuestión, esto es: ⎛ A ⎞ ⎛ kkr p ⎞⎟ ...(8.17) T p,i + 1 = T p p p,i +1 , p p,i = ⎜ ⎟ ⎜ 2 ⎝ Δx ⎠ i + 1 ⎜⎝ μ p B p ⎟⎠ 1 i+ 2

(

)

2

Observaciones: 1.- Los criterios para calcular las transmisibilidades, T p , mencionados en el problema de flujo monofásico, aplican en este caso en lo que se refiere al grupo k / μ p B p . 2.- La permeabilidad relativa de la fase p, kr p , deben sin embargo evaluarse, según experimentos numéricos publicados en la literatura, en el nodo de mayor potencial de esa fase p, o lo que se denomina corriente arriba. 3.- Algunos autores (Peaceman) recomienda evaluar corriente arriba no solo kr p sino

krp / μ p Bp . De esta manera se tendrá por ejemplo que: T g,i + 1 ⎛⎜ p g ,i + 1 , S g ,i + 1 2⎝ 2 2

( (

⎧T 1 p g ,i +1 , S g ,i +1 ⎞ = ⎪ g ,i + 2 ⎟ ⎨ ⎠ ⎪T g ,i + 1 p g ,i , S g ,i 2 ⎩

)

)

(

si Δp g − λ g ΔD

)i+

1 2

> 0

...(8.18)

caso contrario

85

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

8.4 Acoplamiento de las Condiciones de Frontera a las Ecuaciones Aproximadas.

Las condiciones de frontera se expresan en diferencias finitas como sigue: q o = −T n +1 1 [Δp o − γ o ΔD ]n1 +1 o, 2

Y,

...(8.19)

2

[

]I +

T n+1 1 Δp p − γ p ΔD n+11 = 0 p, I + 2

...(8.20)

2

Por similitus con la Ec. 8.19, los gastos de producción de gas y de agua, q g y q w , pueden escribirse en forma de diferencias como sigue:

[

]

q g = −T n +11 Δp g − γ g ΔD n1 +1 − (To R s )n1 +1 [Δp o − γ o ΔD ]n1 +1 g, 2

Y,

2

2

...(8.21)

2

q w = −T n+11 [Δp w − γ w ΔD ]n1 +1 w, 2

...(8.22)

2

Aunque los gastos de producción de gas y de agua están subordinados al gasto de producción de aceite, según lo muestran las Ecs. 8.10 y 8.11, las ecuaciones anteriores, Ecs. 8.19, 8.21 y 8.22, sirven de enlace para acoplar las condiciones de frontera en las ecuaciones aproximadas de flujo. El acoplamiento de las condiciones de frontera en las ecuaciones aproximadas se hace en forma similar al problema de flujo monofásico. Las ecuaciones de los nodos i =1 e i = I , adquieren entonces formas particulares. Para el nodo i =1, se tiene: Para el componente aceite, Ton,1+.15 [Δp o − γ o ΔD ]1n.+51 + q on +1 =

Para el componente gas,

[

(

Vr1 ⎛ φ 1 − S g − S w Δ t ⎜⎜ Δt Bo ⎝

)⎞⎟

⎟ ⎠1

...(8.23)

]

T gn,+11.5 Δp o + ΔPc og − γ g ΔD n +1 + (To R s )1n.+51 [Δp o − γ o ΔD ]1n.+51 + 1.5

(

n +1

⎡φS φ 1− S g − S w ⎛ Tg ⎞ Vr ⎜ ⎟ q on +1 = 1 Δ t ⎢ g + R + s ⎜T ⎟ Δt Bo ⎢⎣ B g ⎝ o ⎠1 2

)⎤

...(8.24)

⎥ ⎥⎦ 1

Y para el componente agua, Twn,+11.5

[

]

Δp w − γ w ΔD 1n.+51

n +1

⎡φS ⎤ ⎛T ⎞ Vr + ⎜⎜ w ⎟⎟ q on+1 = 1 Δ t ⎢ w ⎥ Δt ⎣ Bw ⎦1 ⎝ To ⎠ 1

...(8.25)

2

86

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Para el nodo i = I , se tiene: Para el componente aceite, − T n +1 1 [Δp o − γ o ΔD ]n+11 = o, I − I− 2

2

(

⎛ φ 1 − S g − Sw VrI Δ t ⎜⎜ Δt Bo ⎝

)⎞⎟

...(8.26)

⎟ ⎠I

Para el componente gas,

[

]

− T n +1 1 Δp o + ΔPc og − γ g ΔD n +11 − (To R s )n +11 [Δp o − γ o ΔD ]n +11 I−2 g,I − I− I− 2

2

(

)

2

...(8.27)

⎡φS g φ 1 − S g − S w ⎤ Vr = I Δt ⎢ + ⎥ Bo Δt ⎥⎦ I ⎣⎢ B g

Y para el componente agua, − T n+1 1 [Δp w − γ w ΔD ]n+11 = I− w, I − 2

2

⎡φS ⎤ VrI Δt ⎢ w ⎥ Δt ⎣ Bw ⎦ I

...(8.28)

8.5 Análisis de no-linealidades y tratamiento.

Las Ecs. 8.14, 8.15 y 8.16 y sus formas particulares para las celdas de fronteras, definen un sistema no lineal de 3I ecuaciones con 3I incógnitas en el nivel de tiempo n + 1 : po , S g , S w (n +1) ;

(

)I

i = 1,2,..., I. Su solución ha originado una serie de formulaciones, o métodos de soluciones, conocidos en la literatura especializada como: IMPES, Solución Simultánea o SS, Solución Secuencial o SEQ, Semi-Implícitos Linealizados o SL y Totalmente implícitos o TI. Cada uno de los métodos mencionados tiene un desempeño diferente en la solución de problemas específicos. Se distinguen por poseer características de estabilidad numérica diferentes, que van desde menos implícito, IMPES, hasta el más implícito: el método totalmente implícito o método de Newton. Entre ellos, y en orden creciente de estabilidad se tiene: IMPES, SS, SEQ, SL y TI. Los métodos menos implícitos, IMPES, requiere el menor esfuerzo computacional para resolver un problema dado: están sin embargo condicionados a usarse en problemas de flujo donde las variaciones de presión y saturación en las celdas, p y S g , durante una etapa de tiempo, Δt , son suaves. Problemas donde esto no se cumple, problemas de conificación, segregación gravitacional y formaciones de frentes, por ejemplo, requieren de métodos más estables, como los métodos semiimplícitos linealizados, SI o el método totalmente implícito, TI. Debe notarse sin embargo que en la medida en que el método es más implícito, los requerimientos computacionales, memoria y tiempo, son mayores. Consecuentemente, los costos también son mayores. 8.5.1

Método IMPES: Stone y Garder, 1961.

El método fue propuesto en 1961 por Stone y Garder. Debido a los requerimientos computacionales relativamente reducidos del método y a las limitaciones de memoria y velocidad de operación de los sistemas de cómputo de la época, el método IMPES resultó ser el método

87

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

viable en estudios globales de comportamiento de yacimientos en los inicios de la simulación numérica. La idea fundamental del IMPES consiste en reducir, en el nivel de tiempo n + 1, el subsistema de ecuaciones de flujo, Ecs. 8.13 a 8.15, que aplica en cada nodo de la malla de cálculo, i = 1,2,..., I, a una sola ecuación de presiones. El sistema de ecuaciones resultantes se resuelve para obtener la distribución de presiones, p i(n + 1 ) , i = 1,2,3,..., I. Después de obtener la distribución de presiones, las saturaciones se calculan explícitamente, o desacopladamente, en cada nodo de la malla de cálculo: S g , i (n + 1) y Sw,i(n+1), i = 1,2,3,..., I. Estrictamente hablando, el desacoplamiento del cálculo de presiones y saturaciones puede lograrse evaluando explícitamente las funciones dependientes de saturación, krp , Pcgo y Pcwo, en los términos de flujo y fuente de las ecuaciones. El método original y versiones posteriores, como las presentadas por Cotas en 1968 y por Breintenbach en 1969, evalúan explícitamente, al nivel de tiempo n, no solo las funciones de n +1 n presión sino también las funciones de saturación en tales términos, esto es: T pn +1 ≈ T pn , Pcgo ≈ Pcgo

n +1 n y Pc wo en los términos de flujo y fuentes; y ≈ Pc wo

γ np +1 ≈ γ np .

Las ecuaciones aproximadas de flujo del método IMPES, de acuerdo con las consideraciones mencionadas son: Para el componente aceite,

[ (

(

)]

⎛φ 1 − S g − Sw Vr Δ Ton Δp on+1 − γ on ΔD i = i Δ t ⎜⎜ Δt Bo ⎝

Para el componente gas,

[ (

)]

[

(

)⎞⎟

...(8.29)

⎟ ⎠i

)]

n Δ T gn Δp on+1 + ΔPc go − γ gn ΔD + Δ (To Rs )n Δp on+1 − γ on ΔD i = i Vri ⎛⎜ φS g φRs 1 − S g − S w ⎞⎟ Δt + ⎜ Bg ⎟ Δt Bo ⎝ ⎠i

(

Y para el componente agua,

[ (

)

)]

⎛ φS ⎞ Vr n +1 Δ Twn Δp on+1 − ΔPc wo − γ wn ΔD i = i Δ t ⎜⎜ w ⎟⎟ Δt ⎝ Bw ⎠ i

...(8.30)

...(8.31)

i = 1,2,..., I

n = 0,1,2,...

Es necesario ahora hacer explícitas las variables primarias, p o ,i (n +1) , S g,i (n+1) y Sw,i(n+1) en los términos de acumulación de las Ecs. 8.29 a 8.31. Para esto se hacen las siguientes operaciones, se toma como ejemplo el término de acumulación de la ecuación del aceite:

(

⎛φ 1 − Sg − Sw Δ t ⎜⎜ Bo ⎝

) ⎞⎟

(

⎛φ 1 − Sg − Sw ⎜ = ⎟ ⎜ Bo ⎠i ⎝

)⎞⎟ n+1 ⎟ ⎠i

(

⎛φ 1 − Sg − Sw − ⎜⎜ Bo ⎝

)⎞⎟ n ⎟ ⎠i

...(8.32)

88

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

(

)

Sumando y restando (φ / Bo )n+1 1 − S g − Sw n en la Ec. 8.31, se le puede llevar a la siguiente forma:

(

⎛ φ 1 − S g − Sw Δt ⎜⎜ Bo ⎝

)⎞⎟

n+1

⎡φ ⎤ ⎟ = ⎢B ⎥ ⎠i ⎣ o ⎦ i

⎡ φ ⎞ n+1 ⎛ φ ⎞ n ⎤ Δt 1 − S g − S w + 1 − S g − S w ⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥ ⎢⎝ Bo ⎠ ⎝ Bo ⎠ ⎥⎦ ⎣ i

(

n ⎢⎛⎜

) (

)

...(8.33)

Nuevamente sumando y restando φ n +1 / Bon dentro de los paréntesis de corchetes, se llega a:

(

⎛ φ 1 − Sg − Sw Δt ⎜⎜ Bo ⎝

)⎞⎟ = ⎡ φ ⎤n+1Δ (1− S ⎟ ⎠i

⎢ ⎥ ⎣ Bo ⎦i

t

n+1 n ⎡ φ n+1 φ n+1 ⎛ φ ⎞ ⎤⎥ n ⎢⎛⎜ φ ⎞⎟ ⎜ ⎟ − + − − − + − S S S 1 g w g w ⎜ ⎜ ⎟ ⎢⎝ Bo ⎟⎠ Bon Bon ⎝ Bo ⎠ ⎥ ⎣ ⎦i

) (

)

...(8.34)

Es decir:

(

⎛ φ 1 − Sg − Sw Δt ⎜⎜ Bo ⎝

)⎟⎞ = ⎡ φ ⎤n+1Δ (1− S ⎟ ⎠i

⎢ ⎥ ⎣ Bo ⎦i

t

g

⎡1 1⎤ − Sw + 1 − Sg − Sw n ⎢ Δtφ + φ n+1Δt ⎥ n Bo ⎦⎥i ⎣⎢ Bo

) (

)

...(8.35)

Ahora bien, si se considera que: n+1

d ⎛⎜ 1 ⎞⎟ dp ⎜⎝ Bo ⎟⎠

≈

(1 / Bo )n+1 − (1 / Bo )n = Δt (1 / Bo ) pon+1 − pon

Δt po

...(8.36)

Se tiene entonces que: n+1

d ⎛ 1 ⎞ Δt (1 / Bo ) = ⎜ ⎟ dp ⎜⎝ Bo ⎟⎠

Δt po

...(8.37)

Por otro lado, a partir de la definición de la compresibilidad de la roca se establece que: cr =

1 dφ 1 Δtφ ≈ φ dp φ n Δt po

...(8.38)

Donde se supone que la porosidad es función exclusiva de la presión de l aceite: φ = φ ( po ). Entonces:

Δtφ = crφ nΔt po

...(8.39)

Sustituyendo las Ecs. 8.37 y 8.39 en la Ec. 8.35 y finalmente sustituyendo ésta en la ecuación del aceite, Ec. 8.29, se obtiene:

89

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Para el componente aceite,

[ (

)] (

)

Δ Ton Δpon+1 − γ on ΔD i = θ oo Δ t po + θ og Δ t S g + θ owΔ t S w i

...(8.40)

Operando en forma similar en los términos de acumulación de las ecuaciones del gas y del agua, se obtienen las siguientes expresiones, Para el componente gas,

[ (

[

)]

(

)]

n Δ T gn Δp on+1 + ΔPc go − γ gn ΔD + Δ (To Rs )n Δp on+1 − γ on ΔD i = i θ go Δ t po + θ gg Δ t S g + θ gw Δ t S w

(

)i

...(8.41)

Y para el componente agua,

[ (

)] (

)

n +1 Δ Twn Δp on+1 − ΔPc wo − γ wn ΔD i = θ wo Δ t p o + θ wg Δ t S g + θ ww Δ t S w i

...(8.42)

Note que el primer subíndice de los coeficientes θ identifica a la ecuación a que corresponde, aceite, gas o agua, y el segundo coeficiente identifica a la variable primaria, po , S g o

S w . Por ejemplo, θ og se refiere al coeficiente de la variable Δ t S g en la ecuación de aceite. Los coeficientes θ para las tres fases se definen como sigue: Para el componente aceite, n Vri n ⎡ ⎛ φ ⎞ d ⎢ S o cr ⎜⎜ ⎟⎟ + φ n+1 θ oo = Δt dpo ⎢ ⎝ Bo ⎠ ⎣

⎛ φ ⎞ ⎟⎟ ⎜⎜ ⎝ Bo ⎠

Vr θ og = − i Δt

θ ow = − Para el componente gas, n ⎡ Vri ⎢ n ⎛⎜ φ ⎞⎟ d S g cr θ go = + S gnφ n+1 ⎜ ⎟ ⎢ dpg Δt ⎝ Bg ⎠ ⎣⎢ Vr θ gg = − i Δt

⎡⎛ ⎢⎜ φ ⎢⎜ B g ⎢⎣⎝

⎞ ⎟ ⎟ ⎠

n+1

⎛ φ ⎜⎜ ⎝ Bo

Vri Δt

⎛ 1 ⎞ ⎜ ⎟ ⎜ Bg ⎟ ⎝ ⎠

⎛ φR − ⎜⎜ s ⎝ Bo

n+1

⎞ ⎟⎟ ⎠

n+1 ⎤

⎥ ⎥ ⎦

...(8.43)

n+1

...(8.44) n +1

...(8.45)

⎛ φR + Soncr ⎜⎜ s ⎝ Bo

n+1

Vr θ gw = − i Δt

⎞ ⎟⎟ ⎠

⎛ 1 ⎞ ⎜⎜ ⎟⎟ ⎝ Bo ⎠

⎧

n

⎞ d ⎛ Rs ⎞ ⎟⎟ + Sonφ n+1 ⎜⎜ ⎟⎟ dp o ⎝ Bo ⎠ ⎠ ⎛

⎪ dPcgo d ⎜ 1 + S gnφ n+1 ⎨ ⎪⎩ dS g dp g ⎜⎝ B g

⎛ φR s ⎜⎜ ⎝ Bo

⎞ ⎟⎟ ⎠

⎞⎫⎪ ⎟⎬ ⎟⎪ ⎠⎭

n+1 ⎤

⎥ ⎥ ⎦⎥

...(8.46)

n+1 ⎤

⎥ ⎥ ⎥⎦

...(8.47)

n +1

...(8.48)

90

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Y para el componente agua, n n+1 ⎛ 1 ⎞ ⎤ Vri n ⎡ ⎛ φ ⎞ d n + 1 ⎟ ⎥ ⎜ ⎟ +φ S w ⎢cr ⎜⎜ θ wo = dpw ⎜⎝ Bw ⎟⎠ ⎥ Δt ⎢ ⎝ Bw ⎟⎠ ⎣ ⎦ θ wg = 0

θ ww = −

Vri Δt

n+1 ⎡⎛ ⎪ dPc wo d n n +1 ⎧ ⎢⎜ φ ⎞⎟ φ ⎨ S − w ⎢⎜⎝ Bw ⎟⎠ ⎪⎩ dS w dp w ⎣

⎛ 1 ⎞⎫⎪ ⎟⎟⎬ ⎜⎜ ⎝ Bw ⎠⎪⎭

...(8.49) ...(8.50)

n+1 ⎤

⎥ ⎥ ⎦

...(8.51)

Nótese que las incógnitas S g , i y S w, i aparecen solo en el lado derecho de las ecuaciones de flujo, Ecs. 8.40 a 8.42. Esta característica es aprovechada por el método IMPES para combinar las tres ecuaciones aproximadas del nodo i y generar una sola ecuación en términos de presión, exclusivamente. Para esto, se multiplica la ecuación del gas, Ec. 8.41, por αg i y la ecuación de l agua, Ec. 8.42 por αwi , y sumando las ecuaciones juntamente con la ecuación del aceite, Ec. 8.40. Esto es:

[ (

)]

[ (

[

)]

(

)]

n Δ Ton Δpon +1 − γ on ΔD i + αgi Δ Tgn Δpon +1 + ΔPcgo − γ gn ΔD + αgi Δ (To Rs )n Δpon +1 − γ on ΔD i i

[ (

)] (

)

n +1 + αwi Δ Twn Δpon +1 − ΔPcwo − γ wn ΔD i = θooΔt po + αgθ go + αwθ wo Δt po,i + i θog + αgθ gg + αwθ wg Δt S g,i + θow + αgθ gw + αwθ ww Δt Sw,i

)i

(

(

)i

...(8.52)

Para eliminar Δ t S g y Δ t S w de la ecuación anterior es necesario que se cumpla lo siguiente: ...(8.53) θ og + αgθ gg + αwθ wg = θog + αgθ gg = 0

(

Y

) (

)i

(θ ow + αgθ gw + αwθ ww )i = 0

...(8.54)

Esto constituye un sistema de dos ecuaciones con dos incógnitas en cada nodo i : αgi y αwi . Resolviendo este sistema se obtiene: ⎛ θ og ⎜ θ gg ⎝

α g i = −⎜

⎞ ⎟ ⎟ ⎠i

...(8.55)

Y ⎛ θ og ⎛ 1 ⎞ ⎡ ⎟⎟ ⎢θ ow − ⎜ αwi = −⎜⎜ ⎜ θ gg ⎝ θ ww ⎠ i ⎢⎣ ⎝

[ (

⎤ ⎞ ⎟θ gw ⎥ ⎟ ⎥⎦ ⎠ i

...(8.56)

Con estas definiciones de αg y αw, la ecuación combinada, Ec. 8.52, se reduce a:

)]

[ (

)]

[

(

)]

n Δ Ton Δpon+1 − γ on ΔD i + αgi Δ Tgn Δpon+1 + ΔPcgo − γ gn ΔD + αgi Δ (To Rs )n Δpon+1 − γ on ΔD i i

[ (

)] (

)

...(8.57)

n+1 + αwi Δ Twn Δpon+1 − ΔPcwo − γ wn ΔD i = θ oo + αgθ go + αwθ wo Δt po,i i

91

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Esta ecuación contiene exclusivamente las incógnitas: po,i −1n+1, po,i −1n+1 y po,i −1n+1. Puede entonces escribirse de la siguiente manera: cin p on,+i −11 + ain p on,+i 1 + bin p on,+i +11 = d in

...(8.58)

i = 1,2,..., I n = 0,1,2,...

La Ec. 8.57 escrita en cada uno de los nodos de las celdas, i = 1,2,...,I al nivel de tiempo n +1, constituye un sistema tridiagonal de I ecuaciones lineales con I incógnitas, p o,i n+1 , i = 1,2,..., I. Resuelta la distribución de presiones, se puede resolver la distribución de saturaciones S g , i n +1 , S w, i n +1 , i = 1,2,..., I , mediante las siguientes expresiones:

S gn,+i1 = S gn,i + Δ t S g ,i

...(8.59)

S wn+,i1 = S wn ,i + Δ t S w,i

...(8.60)

Y

Donde Δ t S g y Δ t S w se calculan a partir de las Ecs. 8.40 y 8.42, como sigue:

{[ (

)]

⎛ 1 ⎞ n+1 ⎟⎟ Δ Twn Δpon+1 − ΔPcwo Δ t S w,i = ⎜⎜ − γ on ΔD i − θ wo,i Δ t po,i ⎝ θ ww ⎠ i

Y ⎛ 1 Δ t S g ,i = ⎜ ⎜ θ og ⎝

{[ (

)] (

}

...(8.61)

)}

...(8.62)

⎞ ⎟ Δ Ton Δpon+1 − γ on ΔD − θ oo Δ t po + θ owΔ t S w i i ⎟ ⎠i

Una vez calculadas ( po, Sg , Sw)in+1 se puede calcular ( pg , pw, So)in+1 y avanzar al siguiente nivel de tiempo. 8.5.1 Estructura matricial del sistema de ecuaciones

Al escribir la Ec. 8.52 en cada uno de los nodos de la malla de cálculo, i = 1,2,...,I se genera un sistema tridiagonal de ecuaciones que se resuelve con el algoritmo de Thomas. ⎡ a1 b1 ⎤ ⎥ ⎢c ⎥ ⎢ 2 a 2 b2 ⎥ ⎢ c3 a3 b3 ⎥ ⎢ A=⎢ . . . ⎥ ⎥ ⎢ . . . ⎥ ⎢ . . . ⎥ ⎢ ⎢ c I a I ⎥⎦ ⎣ Es evidente que en este sistema el vector de incógnitas esta formado por los subvectores u

T

= ( p1 , p 2 , ..., p I ) y el vector de términos conocidos por los subvectores d T = (d1, d 2 ,...,d I ).

92

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

8.5.2 Método de Newton-Raphson.

A partir de las Ecs. 5.15 a 5.17 se definen las siguientes funciones de residuos: Para el componente aceite, Fon,i+1 = Δ[To (Δp o − γ o ΔD )]in+1 −

(

⎛φ 1 − Sg − Sw Vri Δ t ⎜⎜ Δt Bo ⎝

Para el componente gas,

[ (

)]

)⎞⎟

⎟ =0 ⎠i

Fgn,+i 1 = Δ Tg Δ po + Δ Pc go − γ g Δ D n +1 + Δ [To Rs (Δ po − γ o Δ D )]in +1 − i

(

Vri ⎛⎜ φS g φRs 1 − S g − S w Δt + Bo Δ t ⎜⎝ B g

)⎞⎟

⎟ ⎠i

...(8.63)

...(8.64)

=0

Y para el componente agua, Fwn,+i 1 = Δ[Tw (Δp o − ΔPc wo − γ w ΔD )]in+1 −

i = 1,2,..., I

Vri ⎛ φS w ⎞ ⎟⎟ = 0 Δ t ⎜⎜ Δt ⎝ Bw ⎠ i

...(8.65)

n = 0,1,2,... Nótese que las funciones de residuos, Ecs. 8.63 a 8.65, del nodo i, tienen en lo general las siguientes dependencias de las variables primarias:

[(

F p ,i p o , S g , S w

Donde p = o, g , w.

)i −1 , ( p o , S g , S w )i , ( p o , S g , S w )i +1 ]n +1

=0

...(8.66)

A partir de las funciones de residuos se establece el siguiente proceso iterativo para resolver el sistema no lineal de ecuaciones: ⎛ ∂F p, i ⎞ ⎟ ⎜ ⎜ ∂po, i −1 ⎟ ⎠ ⎝ ⎛ ∂F p, i ⎜ ⎜ ∂p , ⎝ oi ⎛ ∂F p, i ⎞ ⎜ ⎟ ⎜ ∂po, i +1 ⎟ ⎝ ⎠

(ν )

(ν )

⎞ ⎟ ⎟ ⎠

⎛ ∂F p, i ⎞ ⎟ δpo(ν, i+−11) + ⎜

(ν )

⎜ ∂S g , i −1 ⎟ ⎝ ⎠

⎛ ∂F p, i ⎞ ⎟ δpo(ν, i+1) + ⎜ ⎜ ∂S g , i ⎟ ⎝ ⎠

⎛ ∂F p, i ⎞ ⎟ δpo(ν, i++11) + ⎜ ⎜ ∂S g , i +1 ⎟ ⎝ ⎠

(ν )

(ν )

(ν )

⎛ ∂F p, i ⎞ ⎟ δS g(ν, i+−11) + ⎜ ⎟ ⎜ ⎝ ∂S w, i −1 ⎠

⎛ ∂F p, i ⎞ ⎟ δS g(ν, i+1) + ⎜ ⎜ ⎟ ⎝ ∂S w, i ⎠

⎛ ∂F p, i ⎞ ⎟ δS g(ν, i++11) + ⎜ ⎟ ⎜ ⎝ ∂S w, i +1 ⎠

(ν )

(ν )

(ν )

δS w(ν, i+−11) +

δS w(ν, i+1) +

...(5.67)

δS w(ν, i++11) + = − F p(ν, i)

Como se mencionó en el problema monofásico, es común iniciar el proceso iterativo en el nivel de tiempo n + 1 partiendo de la solución de nivel de tiempo n , esto es:

( po , S g , S w )in +1 = ( po , S g , S w )in

i = 1,2,..., I

...(5.68)

93

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Los criterios de convergencia son similares a los mostrados en las Ec. 4.35 y 4.36. Notación compacta: Si se define: ∂F p,i a po,i = ∂po,i a pg,i =

Y, a pw,i =

∂F p,i

...(8.69) ...(8.70)

∂S g ,i ∂F p,i

...(8.71)

∂S w,i

Donde p = o, g , w.

b pp, i y c pp, i se definen similarmente, como las derivadas de la función de residuos, F pi , con respecto a las incógnitas en i + 1 e i − 1 , respectivamente. 8.5.2.1 Estructura matricial del sistema de ecuaciones.

De acuerdo con las Ecs. 8.69, 8.70 y 8.71, entonces se puede escribir el subsistema de ecuaciones, Ec. 8.67 para el nodo i como se muestra a continuación:

(ν +1)

(ν )

⎡ coo cog cow⎤ ⎢ ⎥ ⎢ cgo cgg cgw⎥ ⎢cwo cwg cww⎥ ⎦i ⎣

(ν )

⎡ aoo aog aow⎤ ⎢ ⎥ ⎢ ago agg agw⎥ ⎢awo awg aww⎥ ⎣ ⎦i

(ν )

⎡ boo bog bow⎤ ⎢ ⎥ ⎢ bgo bgg bgw⎥ ⎢bwo bwg bww⎥ ⎣ ⎦i

⎡δpo ⎤ ⎢δp ⎥ ⎢ g⎥ ⎢⎣δpw⎥⎦ i−1

(ν +1)

⎡δpo ⎤ ⎢δp ⎥ ⎢ g⎥ ⎢⎣δpw⎥⎦ i

(ν )

⎡ Fo ⎤ = − ⎢⎢Fg ⎥⎥ ⎢⎣Fw⎥⎦ i

...(8.72)

(ν +1)

⎡δpo ⎤ ⎢δp ⎥ ⎢ g⎥ ⎣⎢δpw⎦⎥

i+1

O bien:

ci(ν )δu i(ν−1+1) + ai(ν )δu i(ν +1) + bi(ν )δui(ν+1+1) = − Fi(ν ) i = 1,2,..., I

...(8.73)

n = 0,1,2,... Nótese que c es la submatriz de orden (3x3 ) de la Ec. 8.73 que contiene los elementos c po, c pg , c pw. a y b son también submatrices de orden (3x3) , similarmente definidas . δu es el

subvector de incógnitas: (δpo , δS g , δS w )T y F es el subvector de residuos (Fo , Fg , Fw )T .

94

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

Al escribir la Ec. 8.72 en cada uno de los nodos de la malla de cálculo, i = 1,2,...,I se genera un sistema bloque-tridiagonal de ecuaciones. Esto es, se tiene un sistema tridiagonal donde los elementos son submatrices de (3x3) : ci , ai y bi . ⎡ a1 ⎢c ⎢ 2 ⎢ ⎢ A=⎢ ⎢ ⎢ ⎢ ⎢ ⎣

b1 a2 c3

b2 a3 .

b3 . .

. . .

. . cI

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ . ⎥ a I ⎥⎦

...(8.74)

Es evidente que en el sistema del vector de incógnitas esta formado por los subvectores

δu = (δu1, δu2 ,...,δu I ) y el vector de términos conocidos por los subvectores F T = −(F1, F2 ,..., FI ) . T

8.6 Extensiones al Flujo Multifásico Bidimensional.

En relación con el método IMPES, se tiene que independiente del número de dimensiones, el método combina las ecuaciones de flujo multifásico de manera tal que la ecuación resultante contiene como incógnitas únicamente a la presión del aceite; las incógnitas de saturación, Sg y Sw, son resueltas desacopladamente después de resolver las presiones. Se verifica que la ecuación combinada, escrita en cada uno de los nodos de la mall de cálculo, i = 1,2,..., I, genera un sistema de ecuaciones con la misma estructura del problema monofásico: en la aplicación del método IMPES al problema unidimensional, se encontró una estructura tridiagonal en el sistema de ecuaciones. Consecuentemente, la aplicación del IMPES a problemas bidimensionales y tridimensionales genera respectivamente sistemas pentadiagonales y heptadiagonales. En la aplicación del método de Newton-Raphson al problema multifásico, se nota que las incógnitas: δpo , δS g y δS w son resueltas simultáneamente en cada nodo de la malla de cálculo. Se verifica que el sistema de ecuaciones resultante puede verse, análogamente al caso monofásico, como un sistema tridiagonal. Los elementos de ese sistema tridiagonal son sin embargo submatrices de orden (3x3) . Extendiendo estas ideas, y con los antecedentes del problema monofásico, es evidente que al aplicar el método de Newton a problemas de flujo multifásico bidimensional y tridimensional, se generan respectivamente sistemas de ecuaciones bloque-pentadiagonales y bloque-heptadiagonales.

95

FI-UNAM

FUNDAMENTOS DE SIMULACIÓN NUMÉRICA DE YACIMIENTOS.

BIBLIOGRAFÍA.

1. AZIZ, K. and SETTARI, A.: Petroleum Reservoir Simulation, Applied Science Publishers, Ltd., London, (1979). 2. BEHIE, A. and VINSOME, P.K.W.: “Block Iterative Methods for Fully Implicit reservoir Simulation”, SPEJ (Oct. 1982) 658-668. 3. BREITENBACH, E.A., THURNAU, D.H. and VAN POOLEN, H.K.: “The Fluid Flow Simulation Equations”, SPEJ (1969), 9, No. 2, pp. 155. 4. COATS, K.H.: “Use and Misuse of Reservoir Simulation Models”, JPT, (Nov., 1969), 13911398. 5. COATS, K.H.: “Reservoir Simulation: State of the Art”, JPT, (Aug., 1982), 1633-1642; Trans., AIME, 273. 6. DOUGLAS, J.Jr., PEACEMAN, D.W. and RACHFORD, H.H. Jr.: “A Method for Calculating Multi-Dimensional Immiscible Displacement”, Trans. SPE of AIME, 216, (1959), pp. 297. 7. LETKEMAN, J.P. and RIDINGS, R.L.: “A Numerical Coning Model”, SPEJ (1970), 19, No. 4, pp.418. 8. MacDONALD. R.C. and COATS, K.H.: “Methods for Numerical Simulation of Water and Gas Coning”, Trans. SPE of AIME, 249, (1970), pp. 107. 9. NOLEN, J.S. and BERRY, D.W.: “Tests of the Stability and Time-Step Sensitivity os SemiImplicit Reservoir Simulation Techniques”, Trans. SPE of AIME, 253, (1972), pp. 253. 10. MITCHELL, A.R.: The finite Difference Method in Partial Differential Equations. John Wiley & Sons, (1980). 11. ODHE, A.S.: “Reservoir Simulation… What is it? ”, JPT, (Nov. 1969), 1383-1388. 12. PEACEMAN, D.W.: Fundamentals of Numerical Reservoir Simulation, Elsevier, Amsterdam, (1977). 13. PRICE, H.S. and COATS, K.H.: “Direct Methods in Reservoir Simulation”, SPEJ, (June 1974), 19. 14. RODRIGUEZ de la GARZA, F.: “CONIMP: Un Simulador Numérico del Fenómeno de Conificación de Fluidos”, Revista IMP, (Julio 1985), Vol. 17, No. 3, pp. 17. 15. SETTARI, A. and AZIZ, K.: “Use of Irregular Grid in Reservoir Simulation”, SPEJ (April 1972) 103-114. 16. SPILLETE, A.G., HILLESTAD, J.G. and STONE, H.L.: “A High-Stability Sequential Solution Approach to Reservoir Simulation”, SPE 4542, Presented at the 48th Annual Fall Meeting, Las Vegas, Nev. (1973). 17. STAGGS, H.M. and HERBECK, E.F.: “Reservoir Simulation Models- An Engineering Overview”, JPT, (Dec., 1971), 1428-1436. 18. STONE, H.L. and GARDER, A.O.Jr.: “Analysis of Gas-Cap or Dissolved-Gas Drive Reservoirs”, SPEJ (June 1961) 91-102. 19. STRANG, G.: “Linear Algebra and Its Applications”, Academic Press, New York, 1976. 20. WATTS, J.W.: “A Method for Improving Line Successive Over-relaxation in Anysotropic Problems- A theoretical Analysis”, SPEJ (April 1973) 105-108, Trans., AIME, 255.

96