INTRODUCCIÓN A LA SIMULACIÓN NUMÉRICA DE YACIMIENTOS

CAPÍTULO 1 INTRODUCCIÓN A LA SIMULACIÓN NUMÉRICA DE YACIMIENTOS

El objetivo de la ingeniería de yacimientos es adquirir un mejor conocimiento de las características del yacimiento de tal manera que el ingeniero de yacimientos esté en capacidad de estimar las reservas recuperables, definir el mejor esquema de explotación que permita recuperar la mayor cantidad de hidrocarburos a un bajo costo y predecir el comportamiento futuro del yacimiento. Esto se logra partiendo del hecho de que las condiciones reales del yacimiento se pueden representar a través de un modelo. Para nuestro propósito, debemos considerar un modelo como una entidad que permite el estudio de un fenómeno, bajo condiciones de prueba apropiadas, que tiene la probabilidad de que ocurra en la práctica. En general los modelos pueden clasificarse como físicos y matemáticos. Los físicos son reproducciones de laboratorio tendientes a reproducir lo que ocurre en el yacimiento y las ecuaciones que describen de manera teórica este comportamiento constituyen el modelo matemático. Ambos tipos de modelos han jugado un papel importante en la industria del petróleo. Por ejemplo, las leyes que gobiernan el flujo de fluidos en un medio poroso fueron descubiertas y delineadas empleado modelos físicos. La ley de Darcy, los conceptos de permeabilidad relativa, presión capilar, densidad y correlaciones de viscosidad, entre otras, tienen sus orígenes en experimentos con modelos físicos. Está de más decir que estos modelos han sido y son indispensables en la ingeniería de yacimientos. Sin embargo, tienen sus limitaciones: es poco práctico el modelamiento riguroso de sistemas de gran escala como lo es un yacimiento de petróleo, por lo que podría decirse que los modelos físicos son muy útiles en el estudio de fenómenos a pequeña escala. Cuando se requiere modelar sistemas globales, como un yacimiento, se debe recurrir a un enfoque diferente, usualmente un enfoque matemático. El deseo de tratar en forma adecuada un yacimiento con algún grado de exactitud dió origen a la tecnología conocida como simulación de yacimientos. Esto no quiere decir que las técnicas de simulación de yacimientos están limitadas a soluciones globales. Estas también se usan en el estudio de fenómenos locales alrededor de la cara del pozo y han demostrado ser superiores, en este aspecto, a los modelos físicos.

1

El modelo matemático más familiar es la ecuación de balance de materia. Este es un modelo matemático o simulador de yacimiento en todo sentido. Este modelo se basa en un concepto físico fundamental: el principio de conservación. Este principio cuando se expresa en forma matemática bajo las limitaciones de suposiciones arbitrarias constituye el modelo. Vale la pena anotar que los modelos modernos están basado en los mismos principios. Difieren en la medida en que se re-definen las suposiciones inherentes en la ecuación de balance de materiales y se aproxima más cercanamente a las condiciones reales del yacimiento. Antes de discutir los principios involucrados en la construcción de un simulador numérico de yacimientos, se debe aclarar a lo que se hace referencia por “solución numérica”. Supóngase que se expresa un proceso físico particular por medio de una ecuación matemática (o conjunto de ellas), que el proceso toma lugar en alguna región que tiene dimensiones finitas. Las ecuaciones que describen este proceso siempre harán referencia a un punto dentro de la región. Al realizar un proceso matemático sobre la ecuación, se puede determinar la interacción de este punto con todos los otros puntos en el sistema a todos los tiempos y de ese modo se predice el comportamiento del proceso físico bajo estudio. Este proceso comúnmente se conoce como solucionar la ecuación. Si las ecuaciones matemáticas no son muy complejas, el proceso de solución puede llevar a obtener una fórmula la cual puede ser subsecuentemente manipulada para calcular los parámetros deseados. Esta solución es una solución analítica. Por otra parte, si las ecuaciones son muy complejas y no pueden resolverse analíticamente, debemos conformarnos con algo menos que una fórmula para representar la solución. Lo que se debe hacer en estos casos, es reemplazar la ecuación original por un conjunto más simple que sea fácil de resolver y que esté relacionado, de alguna manera, con la ecuación original. Sin embargo, en lugar de obtener otra(s) fórmula(s), se puede llegar a soluciones de las ecuaciones más simples en forma de tablas de valores numéricos, cada uno de los cuales se refiere a puntos discretos en el espacio y en el tiempo dentro de la región. Esto se conoce como solución numérica, la cual representa una aproximación a la ecuación original que se quería resolver. Si se quiere construir un modelo matemático que supere las limitaciones de la ecuación de balance de materiales, por ejemplo, siempre se van a obtener ecuaciones que caen dentro de la segunda categoría, las cuales requieren soluciones numéricas, de aquí el nombre de simulación numérica de yacimientos.

2

En términos generales, podría decirse que un simulador numérico está integrado por tres modelos: el diferencial, el numérico y el de computador. El modelo diferencial está conformado por el conjunto de ecuaciones diferenciales que describen los procesos físico-químicos que ocurren en el yacimiento en función del espacio y del tiempo. El modelo numérico constituye la forma en que se da solución al modelo diferencial. Por lo menos existen tres formas de hacer este modelo: •

diferencias finitas, que es el más usado

•

elementos finitos: no se usa porque no permite la unión de las ecuaciones de masa y energía

•

Elementos de volumen de control

Para este propósito se divide el yacimiento en una serie de pequeños bloques, de acuerdo con la geometría del yacimiento, a los que se les puede asignar un valor único para las propiedades de la roca, de este modo se pueden tener en cuenta las heterogeneidades y la anisotropía del yacimiento. La variación espacial de las propiedades del fluido también se puede asignar por bloques o zonas a través de todo el sistema. Para reflejar la existencia de pozos, se asignan tasas de producción e inyección a la celda en donde se encuentren. Una vez se han ubicado los pozos y se han definido las propiedades a cada bloque se plantean las ecuaciones fundamentales en cada nodo del sistema a diferentes niveles de tiempo. De esta forma se obtiene el modelo numérico del simulador constituido por un alto número de ecuaciones discretas. Con el fin de dar solución a este gran número de ecuaciones se requiere la elaboración de un programa de computador conocido como el modelo de computador del simulador. La información básica que se requiere para ejecutar el modelo de computador es la siguiente. •

Geometría del yacimiento

•

Ubicación de pozos productores e inyectores con sus respectivas tasas de flujo.

•

Propiedades de la roca y del fluido en cada bloque

•

Distribución inicial de fluidos en el yacimiento 3

•

Presión promedia en cada nodo en el tiempo cero.

Los resultados obtenidos dependen del enfoque que se le haya dado al simulador, del tipo de modelo y de la aplicación para la cual fue construido. En general se busca estudiar el comportamiento del yacimiento bajo diferentes esquemas de producción, de ubicación de pozos y para diferentes tasas de producción e inyección.

HISTORIA DE LA SIMULACION La Simulación de Yacimientos ha sido practicada desde el inicio de la Ingeniería de Petróleos. En la década de los 40, el potencial de la simulación de yacimientos fue reconocido y muchas compañías iniciaron el desarrollo de modelos analógicos y numéricos con la finalidad de mejorar las soluciones analíticas existentes (cálculo de balance de materiales y desplazamiento 1-D de Buckley-Leverett). En la década de los 50, se llevaron a cabo investigaciones en lo que respecta a solución numérica de ecuaciones de flujo. Como resultado, se obtuvieron programas de computador para simulación de yacimientos, aunque sencillos pero útiles. Estos programas representaron el mayor avance y usaron la solución de un conjunto de ecuaciones de diferencias finitas para describir el flujo multifásico 2-D y 3-D en medios porosos heterogéneos. Fue la primera vez que los Ingenieros de Yacimientos lograron resolver problemas complejos. En la década de los 60, el desarrollo de la Simulación de Yacimientos estuvo dirigida a resolver problemas de yacimientos de petróleo en tres fases. Los métodos de recuperación que fueron simulados incluían depletación de presión y varias formas de mantenimiento de presión. Los programas desarrollados operaban en grandes computadores (Mainframe) y usaban tarjetas para el ingreso de datos. Durante la década de los 70, la tendencia cambió bruscamente, debido al creciente número de investigaciones en procesos EOR, avances en técnicas de simulación numérica y la disminución del tamaño e incremento de velocidad de los computadores. Los simuladores matemáticos fueron desarrollados de tal manera que incluían procesos de inyección química, inyección de vapor y combustión in situ. La investigación durante este período resultó en avances significativos en lo que respecta a la caracterización de la física del hidrocarburo, en el desplazamiento bajo la influencia de la temperatura, agentes químicos y comportamiento de fase multicomponente. Durante la década de los 80, el rango de las aplicaciones de la simulación de yacimientos continuó expandiéndose. La descripción de yacimientos avanzó hacia

4

el uso de la GeoEstadística para describir heterogeneidades y proporcionar una mejor definición del yacimiento. Se desarrolló la tecnología para modelar yacimientos naturalmente fracturados, incluyendo efectos composicionales. Así mismo, el fracturamiento hidráulico y pozos horizontales y su aplicación al monitoreo del yacimiento. Al inicio de la década de los 90, las aplicaciones fueron hechas en grandes computadores (Mainframe), al final de la década se empezaron a usar microcomputadores. Actualmente, computadores personales y una gran cantidad de sistemas de simulación de yacimientos, proporcionan al Ingeniero, un medio económico y eficiente para resolver complejos problemas de Ingeniería de Yacimientos.

AVANCES RECIENTES Los avances recientes se han centrado principalmente en los siguientes temas : 1.- Descripción del yacimiento. 2.- Yacimientos Naturalmente Fracturados. 3.- Fracturamiento Hidráulico. 4.- Pozos Horizontales. Referente a descripción del yacimiento, se están aplicando técnicas estocásticas sustentadas en lo siguiente : a).- Información incompleta del yacimiento en todas sus escalas. b).- Propiedades de roca variables. c).- Relación desconocida entre propiedades. d).- Abundancia relativa de muestras con información proveniente de los pozos. Referente a yacimientos naturalmente fracturados, la simulación se ha extendido a aplicaciones composicionales e inyección cíclica de vapor. Respecto a fracturamiento hidráulico, se ha enfatizado en la predicción de la geometría de la fractura. Se dispone de varias técnicas para predecir la distribución de los esfuerzos in situ, mejorando de esta forma la simulación del crecimiento de la fractura en el sentido vertical y lateral.

5

El objetivo de la simulación de pozos horizontales es estudiar los efectos de longitud del pozo, ángulo de inclinación, heterogeneidades locales, permeabilidad direccional, barreras y caída de presión en el pozo. La simulación exacta de los fenómenos cerca al pozo, ha permitido estudiar los efectos que tienen los pozos horizontales sobre la productividad, intersección de fracturas, conificación y recuperación de hidrocarburos.

TIPOS DE SIMULADORES Los simuladores de yacimientos se pueden clasificar de acuerdo al número de dimensiones, al tipo de yacimiento que se quiere simular o al proceso particular que se quiere estudiar. Ver figura 1. De acuerdo al número de dimensiones El modelo más simple es el de dimensión cero o modelo tanque ya que considera el yacimiento como un tanque. En este modelo las propiedades del fluido y de la roca no varían de punto a punto y considera que el disturbio de presión inducido en el yacimiento debido a la producción o inyección de fluidos se transmite en forma instantánea a través de todo el sistema, razón por la cual los cálculos siempre se realizan a condiciones estáticas.

Modelo tanque o dimensión cero Con miras a tener un enfoque más realístico, se pueden unir dos o más de estos tanques, asignando un valor único a las propiedades de la roca y del fluido a cada tanque y permitiendo flujo de uno al otro a través de las caras adyacentes. Este constituye un modelo en una dimensión. Aunque los modelos en una dimensión consideran al yacimiento más detalladamente que el modelo de dimensión cero, solamente dan una idea general del movimiento de fluidos y de la distribución de presiones en el yacimiento en función del tiempo. Estos modelos 1D son usados para estudiar la sensibilidad del comportamiento del yacimiento a las variaciones de los parámetros del mismo. Por ejemplo, la sensibilidad del petróleo recuperable a la relación de movilidad, permeabilidad absoluta o la forma de las curvas de permeabilidad relativa.

6

Dimensión Cero

Una Dimensión

Modelo Lineal Modelo Radial Geometría Horizontal

Número de Dimensiones

Dos Dimensiones

Geometría Vertical Geometría Radial

Tres Dimensiones

Modelo Cartesiano Modelo Cilíndrico

Tipos de Simulador

Yacimientos de Gas Tipo de Yacimiento

Yacimientos de Petróleo Yacimientos de Condensado

Recuperación Primaria

Proceso de Recuperación

Simulación Térmica Inyección de Químicos y Polímeros Desplazamiento Miscible

Figura 1. Clasificación General de los Simuladores

7

Estos modelos son raramente usados en estudios de yacimientos para un campo entero, debido a que no se puede modelar el barrido areal y vertical. Por ejemplo, no se pueden efectuar cálculos confiables de la eficiencia del desplazamiento en regiones invadidas debido a que no se puede representar los efectos gravitacionales que actúan perpendicularmente a la dirección del flujo.

(a) (b)

Modelo de una Dimensión: (a) Lineal; (b) Vertical.

r

h

Modelo radial

MODELOS 2-D AREALES.- Los modelos areales cartesianos 2-D (x,y) son los más usados en los estudios de yacimientos. Se usan principalmente para estudiar el yacimiento entero en casos donde el espesor de la formación es relativamente pequeño o donde no hay una gran variación vertical en las propiedades de los fluidos y la formación.

8

Los modelos areales usan normalmente sistema de coordenadas cartesianas (x,y), sin embargo existen algunas aplicaciones que requieren sistemas de coordenadas radiales (r,θ ) o cilindricas. Estos dos últimos sistemas proporcionan una mejor definición cerca a los pozos.

Modelo en dos dimensiones horizontal y vertical El modelo de geometría horizontal es particularmente útil en la simulación de eficiencias de barrido, efectos de barreras y evaluación de arreglos geométricos de pozos de inyección. Los modelos verticales permiten simular la variación vertical de permeabilidad, efectos de estratificación y efectos de segregación de fluidos. SECCIONES TRANSVERSALES.- Se usan primariamente para desarrollar : a.- Para simular inyección de agua periférica o inyección de gas en la cresta con la finalidad de proporcionar información sobre la uniformidad de la eficiencia de barrido. b.- Se usan también para analizar el efecto de la gravedad, capilaridad y fuerzas viscosas sobre la eficiencia de barrido vertical (conificación de agua). Si la eficiencia de barrido areal, es un aspecto importante a ser tomado en cuenta, no se debe usar este tipo de modelo para estimar el comportamiento total del campo. RADIAL. Es usado para desarrollar funciones de pozo que permitan predecir el comportamiento cuando se usen en modelos 2-D areales y 3-D y permiten evaluar el comportamiento de los pozos cuando los efectos verticales dominan el comportamiento como en el caso de la conificación de agua o gas. Los modelos 2-D radiales son muy usados para simular la convergencia o divergencia del flujo en una región radialmente simétrica del yacimiento. Además se usan estos modelos para estudiar el comportamiento de pozos en yacimientos con empuje de agua de fondo, con capa de gas y yacimientos que

9

tienen una delgada columna de petróleo y se encuentran rodeados por agua o gas.

r

θ

h

Modelo radial en dos dimensiones (r,θ)

MODELOS MULTICAPA.- Usados para modelar yacimientos con varias capas sin flujo cruzado ("crossflow"), sin embargo estas capas tienen las mismas condiciones límites, tales como acuífero común o la producción proveniente de las capas se mezcla en el pozo (commingled).

MODELOS 3-D.- Los modelos en tres dimensiones pueden ser cartesianos o radiales. Son los más versátiles ya que permiten simular la variación de las propiedades de la roca y de las condiciones del fluido en forma areal y vertical. Sin embargo su aplicación es limitada debido a su alto costo y a que requiere una caracterización precisa del yacimiento. Los modelos 3-D son usados donde la geometría del yacimiento es muy compleja como para ser modelado por un 2-D. Los yacimientos en etapa de depleción avanzada tienen una dinámica de fluido muy compleja y requieren ser modelados con un simulador 3-D. También se usan modelos 3-D para simular el desplazamiento de fluidos donde los regímenes de flujo son dominados por el flujo vertical. Un problema que está asociado a los modelos 3-D es el tamaño. Un adecuado modelo puede tener tantos gridblock que consumiría mucho tiempo en proporcionar resultados y retardaría la toma de decisiones.

10

Modelo cartesiano en tres dimensiones (x,y,z)

r

θ

z

Modelo radial en tres dimensiones (r,θ,z)

De acuerdo al tipo de yacimiento Según esta clasificación existen tres grupos: simuladores de yacimientos de gas, de yacimientos de aceite negro y de yacimientos composicionales (aceite volátil y gas condensado). Los simuladores de yacimientos de gas pueden tener modelos monofásicos o bifásicos dependiendo de si hay o no agua móvil.

11

Por otra parte, el modelo black oil, también conocido como modelo β es capaz de simular aquellos sistemas donde hay presencia de agua, aceite y gas teniendo en cuenta la solubilidad del gas en el petróleo, por lo cual la transferencia de fases se lleva a cabo entre el gas y el aceite, aunque no tiene en cuenta cambios en la composición de estas fases. En este modelo el fluido del yacimiento es aproximado por dos componentes: un componente no volátil (crudo) y un componente volátil (gas) soluble en la fase del petróleo. Este modelo asume lo siguiente: •

Se considera que el agua y el crudo son inmiscibles

•

El gas es soluble en el aceite, pero no en el agua

•

A condiciones de yacimiento, el aceite está formado por dos componentes: aceite a condiciones de superficie y gas a condiciones normales

•

Dentro del yacimiento los fluidos se encuentran en equilibrio termodinámico

Los yacimientos de condensado y de aceite volátil usualmente requieren un simulador de propósito especial que tenga en cuenta el comportamiento composicional entre los componentes hidrocarburos individuales en las fases gas y líquido. Este tipo de modelo está catalogado dentro de los más complejos puesto que se enfoca a componentes individuales y no a fases, además presentan problemas tales como la determinación de las propiedades físicas de las fases cerca al punto crítico, la caracterización de los componentes del crudo y los largos tiempos requeridos para correr el simulador. De acuerdo al tipo de proceso Los procesos y fenómenos particulares del yacimiento como la conificación en la cara del pozo, el recobro térmico, la inyección de químicos y los desplazamientos miscibles constituyen otros tipos de modelos de yacimiento. En los modelos de conificación se examina en detalle el comportamiento del pozo con el propósito de determinar los mejores intervalos de completamiento y las mejores tasas de producción necesarios para minimizar la conificación de agua o de gas dentro del pozo. Los procesos de recobro térmico, incluyendo la estimulación con vapor, el desplazamiento con vapor y la combustión in situ han dado origen a modelos de yacimientos muy sofisticados que tratan de tener en cuenta todos los fenómenos físicos involucrados.

12

Los modelos de inyección de químicos tienen en cuenta ecuaciones de conservación adicionales para diferentes especies, por ejemplo, polímeros, surfactantes, etc. Además, deben tener en cuenta el efecto de la adsorción e incluir los efectos de reducción de permeabilidad en la fase acuosa despues del contacto con el polímero. La mayoría de los modelos de desplazamiento miscible se basan en modificaciones de los modelos usados en los procesos de desplazamiento inmiscible.

RAZONES PARA EFECTUAR SIMULACION La simulación puede proporcionar beneficios potenciales en los siguientes aspectos: (1) Estudiar la recuperación final primaria y su comportamiento bajo diferentes modos de operación tales como depleción natural, inyección de agua y/o gas. (2) El tiempo en el cual debe iniciarse un proceso de recuperación mejorada a fin de maximizar la recuperación así como el tipo de patrón que debe ser usado. (3) El tipo de proceso de recuperación mejorada más apropiado y cuál será la recuperación final y el comportamiento con el proceso elegido. (4) Investigar los efectos de nuevas ubicaciones y espaciamientos de pozos. (5) Analizar el efecto de las tasas de producción sobre la recuperación (6) Analizar que tipos de datos tienen el mayor efecto sobre la recuperación y por lo tanto los que deben ser estudiados cuidadosamente con experimentos físicos de laboratorio.

PREPARACION DE LOS DATOS La calidad de los datos de salida no puede ser mejor que la calidad de los datos de entrada. Los datos requeridos para hacer un estudio de simulación proviene de varias fuentes y no están siempre en el formato requerido para ser aplicados directamente al computador. Existen diferentes fuentes que proporcionan la misma información. Se debe diferenciar y seleccionar la mejor data disponible. Si ésta no está disponible para

13

un caso particular, se debe determinar alguna forma alternativa de conseguir la misma información.

DATOS REQUERIDOS Los datos requeridos para efectuar una simulación son los siguientes : (.) Dimensión para el modelo del yacimiento (Grid) (.) Geometría del yacimiento (.) Distribución de Porosidad y Permeabilidad (.) Datos de presión capilar y permeabilidad relativa (.) Datos PVT de los fluidos (.) Distribución dentro del yacimiento de la presión y saturación inicial (.) Método de solución de las matrices (.) Parámetros de diagnóstico y control de la ejecución ("Corrida") (.) Datos de producción y de los pozos.

Fuentes para obtener elevación de la Formación Obtenidos a partir de mapas estructurales del subsuelo. Estos datos son compilados inicialmente de : 1).- Datos de perfiles eléctricos. 2).- Registros de perforación.

Fuentes para obtener el espesor de la Formación Es obtenido de mapas de espesor total o de espesor neto. Muchos simuladores usan mapas de espesor total para calcular las características de flujo del modelo.

14

Los mapas de espesor total proporcionan la dimensión vertical correcta necesaria para evaluar la corrección por potencial; sin embargo para calcular el OOIP (que está basado en espesor de arena neta petrolífera), es costumbre incluir ya sea factores de espesor neto/total para permitir el cálculo de OOIP o un programa separado para calcular el OOIP basado en espesores netos. El espesor de la formación también se puede obtener a partir de datos estructurales restando a los contornos estructurales del fondo de la formación de los del tope.

Fuentes para obtener la Permeabilidad La permeabilidad absoluta debe ser obtenida a partir de varias fuentes : 1).- Datos de presión buildup (DST). 2).- Datos de presión falloff. 3).- Pruebas de interferencia. 4).- Pruebas de potencial inicial. 5).- Análisis de regresión. 6).- Medidas de laboratorio. La fuente más importante de datos de permeabilidad es el análisis de pruebas de presión y el ingeniero debe ser familiar con las técnicas actuales disponibles.

Fuentes para obtener la Porosidad Este parámetro se obtiene usualmente de las fuentes siguientes : 1).- Datos de perfiles eléctricos. 2).- Medidas de laboratorio. 3).- Correlaciones publicadas.

15

Datos de Perfiles Eléctricos En la forma de señales acústicas o sónicas. Se obtienen por medir el tiempo de viaje del sonido a través de la formación. El tiempo de viaje es directamente afectado por los fluidos que se encuentran presentes en el espacio poroso de la roca.

Medidas de Laboratorio Que están basadas en la determinación de parámetros tales como: volumen bruto, volumen de granos y volumen poroso. Los métodos usuales determinan el volumen poroso ya sea por la introducción de un fluido hacia la roca o la remoción de los fluidos de la roca.

Correlaciones Publicadas Que se refieren en su mayoría al estimado de la porosidad a una profundidad determinada en base a la compactación natural.

Fuentes para obtener Permeabilidades Relativas Es a menudo la data más dificultosa de evaluar u obtener. requeridas por los simuladores son las siguientes :

Las relaciones

1).- Permeabilidad relativas gas/petróleo. 2).- Permeabilidad relativa petróleo/agua. 3).- Permeabilidad relativa gas/agua.

Fuentes para obtener datos de Presión Capilar Estos datos son necesarios para evaluar las presiones en las diferentes fases durante los cálculos IMPES y también para plantear las ecuaciones en la solución SIMULTANEA. Las presiones capilares son determinadas de datos de laboratorio.

16

Fuentes para obtener datos de Compresibilidad de la Roca Estos datos son obtenidos a partir del análisis de laboratorio o de correlaciones publicadas.

DISTRIBUCIÓN DENTRO DEL YACIMIENTO DE LA PRESIÓN Y SATURACIÓN INICIAL Fuentes para obtener datos de Saturación de fluidos de Formación En un yacimiento existen dos posibles planos de interés que pueden ser usados para evaluar las saturaciones de los fluidos del yacimiento: El contacto gas/petróleo y el contacto agua/petróleo. La saturación de agua connata puede ser evaluada de : 1).- Datos de núcleos. 2).- Perfiles eléctricos. 3).- Datos de presión capilar.

Fuentes para obtener datos de Tasas de Flujo Son requeridos por el simulador para calcular la capacidad productiva de un pozo dentro de un sistema. Estos datos son generalmente basados en lo siguiente: 1).- Indice de productividad. 2).- Indice de Inyectividad 3).- Tasas de flujo óptimas 4).- Máxima diferencial de presión permisible.

El Indice de Productividad es la relación entre la tasa de producción (STB/dia) para flujo líquido a la diferencial de presión en el punto medio del intervalo productor. J = IP = q / (P - Pwf)

17

El IP es una medida del potencial del pozo, o la facilidad del pozo para producir. El Indice de Inyectividad es usado en pozos inyectores para recuperación mejorada. Es la relación entre la tasa de inyección (STB/día) al exceso de presión sobre la presión del yacimiento. II = q / (Pwf - P) Ambos índices (de productividad e inyectividad) están referidos a las presiones en la cara de la arena y las caídas de presión por fricción en el casing o tubing no están incluidas. En el caso de inyección o producción a altas tasas, estas caídas (pérdidas) de presión pueden ser apreciables.

DISEÑO DEL MODELO A pesar que los resultados de un estudio de simulación son más aceptados debido a la complejidad que se le proporciona al modelo, es mejor diseñar el modelo más simple que permita simular el proceso de desplazamiento con suficiente exactitud. El diseño de un modelo es influenciado por los siguientes factores: a.- Tipo y complejidad del problema b.- Tiempo disponible para completar el estudio c.- Costo del estudio d.- Calidad de los datos disponibles f.- Capacidad del simulador y hardware existente.

PLANIFICACION DE UN ESTUDIO DE SIMULACION A continuación se presenta un posible orden de las actividades más significantes a llevar a cabo durante un estudio de simulación : A.- Definición del Problema. B.- Adquisición y Revisión de datos. C.- Descripción del yacimiento y Diseño del modelo. D.- Ajuste de historia.

18

E.- Predicción. F.- Reporte.

A.- DEFINICION DEL PROBLEMA.- El primer aspecto a tratar cuando se lleva a cabo un estudio de simulación es definir los problemas del comportamiento del yacimiento y problemas operativos asociados. Para efectuar esto se debe reunir la información suficiente acerca del yacimiento y su forma de operación para identificar las alternativas necesarias en lo que respecta a pronósticos. Se debe definir en forma clara y concisa el objetivo práctico del estudio. Así mismo son necesarias evaluaciones rápidas a fin de identificar el mecanismo principal de depletación y reconocer que factores dominarán el comportamiento del yacimiento (gravedad, heterogeneidad, conificación, etc.). Si es posible, determinar el nivel de complejidad del modelo de yacimiento, para iniciar el diseño del mismo e identificar los datos necesarios para su construcción.

B.- ADQUISICION Y REVISION DE LOS DATOS.- Los datos deben ser revisados y reorganizados después que estos hayan sido coleccionados, debido a que estos han sido obtenidos para diferentes razones y normalmente no han sido organizados de tal forma que tengan un uso inmediato. La revisión debe efectuarse cuidadosamente y se debe consumir todo el tiempo necesario a fin de evitar trabajo inútil.

C.- DESCRIPCION DEL YACIMIENTO Y DISEÑO DEL MODELO. El diseño de un modelo de simulación estará influenciado por el tipo de proceso a ser modelado, problemas relacionados con la mecánica de fluidos, los objetivos del estudio, la calidad de los datos del yacimiento y su descripción, restricciones de tiempo y el nivel de credibilidad necesario para asegurar que los resultados del estudio sean aceptados.

D.- AJUSTE DE HISTORIA.- Después que un modelo de yacimiento ha sido construido, debe ser probado a fin de determinar si puede duplicar el comportamiento del mismo. Generalmente la descripción del yacimiento usada en el modelo es validado haciendo "correr" el simulador con datos de producción e

19

inyección histórica y comparar las presiones calculadas y el movimiento de fluido con el comportamiento actual del yacimiento. E.- PREDICCION.- Una vez que se ha obtenido un ajuste de historia aceptable, el modelo puede ser usado para predecir el comportamiento futuro del yacimiento y así alcanzar los objetivos trazados por el estudio. La calidad de las predicciones dependerá de las características del modelo y la exactitud de la descripción del yacimiento.

F.- REPORTE.- El paso final de un estudio de simulación es plasmar los resultados y conclusiones en un reporte claro y conciso. El reporte puede ser un breve memorando para un pequeño estudio o un informe completo de gran volumen para un estudio a nivel yacimiento. En el reporte se debe incluir los objetivos del estudio, descripción del modelo usado y presentar los resultados y conclusiones referentes al estudio específico. Dada la importancia del ajuste histórico a continuación se presentan algunos aspectos a tener en cuenta.

AJUSTE HISTÓRICO (History Matching) El principal objetivo de un estudio de simulación es predecir el comportamiento futuro del yacimiento con mayor exactitud que alguna otra técnica simple de predicción. Es evidente que el comportamiento del modelo numérico debe ser similar al del yacimiento para que los resultados sean aceptables. Debido a la incertidumbre inherente a los datos requeridos para construir el modelo, se debe probar el comportamiento del modelo antes de ser usado para predecir el comportamiento futuro. La única forma de probar el modelo es simular el comportamiento pasado del yacimiento y comparar los resultados con los datos históricos. El proceso de probar el modelo a través de comparar el comportamiento pasado es usado también para identificar las inconsistencias del modelo y corregirlo. El ajuste de historia es, por lo tanto, el proceso de refinar el modelo a través del ajuste de parámetros de geología, roca y fluido, para producir la mínima diferencia entre los datos de campo y los resultados del simulador.

20

PARAMETROS PARA EL AJUSTE DE HISTORIA. a.- Presión b.- Tasas de flujo. c.- GOR d.- WOR e.- Tiempo de irrupción del frente. El objetivo es minimizar la diferencia entre estos parámetros y los obtenidos por el simulador.

PARAMETROS QUE PUEDEN SER MODIFICADOS Existen varios parámetros que pueden ser modificados ya sea solo o en conjunto para lograr un buen ajuste de historia : a.- Permeabilidad y espesor del yacimiento. b.- Permeabilidad y espesor del acuífero. c.- Almacenamiento del acuífero. d.- Datos de permeabilidad relativa. e.- Datos de presión capilar. f.- Datos del pozo (factor skin, ect).

Parámetros adicionales que son conocidos con mayor certeza pero que a veces pueden ser variados : g.- Porosidad y espesor del yacimiento. h.- Definición geológica del yacimiento. i.- Compresibilidad de la roca. j.- Propiedades de los fluidos.

21

k.- Contactos agua/petróleo y gas/petróleo. l.- Presión fluyente de fondo.

MECANICA DEL AJUSTE DE HISTORIA a.- Reunir los datos de historia de producción. b.- Evaluar su calidad c.- Definir los objetivos para el ajuste de historia. d.- Desarrollar un modelo preliminar basado en los mejores datos disponibles. e.- Comparar los resultados del simulador con el comportamiento del yacimiento. f.- Decidir si los resultados del ajuste están dentro de una tolerancia aceptable. g.- Decidir si es necesario un ajuste de historia automático. h.- Efectuar ajustes al modelo y simular otra vez para mejorar el ajuste.

ANALISIS DATOS DE CAMPO Los datos de producción deben ser analizados pozo a pozo para identificar y eliminar inconsistencias. Se puede incluir : a.- Producción de petróleo. b.- Producción e inyección de gas. c.- Producción e inyección de agua. d.- Presiones fluyentes o de cierre corregidas al datum. Resultados inexactos de producción de un pozo o zona productiva deben también ser evaluados. Se debe tener especial cuidado para refinar estos datos, ya que estos pueden representar características anormales si es que no se eliminan. Cabe mencionar que datos de producción e inyección de agua no son medidos tan exactamente como la producción de petróleo.

22

El petróleo in-situ, así como los contactos agua/petróleo y gas/agua deben ser comparados con estimados perfectamente conocidos, y si hubiera diferencia, proceder a revisar a fin de continuar con la predicción.

AJUSTE DE LA HISTORIA DE PRESION Se recomiendan los pasos siguientes para un exitoso ajuste de presión : a.- Identificar los parámetros a ser ajustados. Normalmente la permeabilidad de la roca es la variable menos definida, y es usada para producir un ajuste de presión. La porosidad no debe ser ajustada, a menos que exista incertidumbre en la data. Si la porosidad es obtenida del análisis de perfiles eléctricos o núcleos, no debe ser cambiada. La porosidad, espesor y extensión areal del acuífero son menos conocidos que en el yacimiento de petróleo, y pueden ser ajustados para obtener una buena reproducción de la presión. b.- Estimar el anteriormente.

nivel

de

incertidumbre

para

las

variables

mencionadas

c.- Efectuar una primera corrida de prueba y decidir si la presión volumétrica promedia del yacimiento (entero) es satisfactoriamente reproducida por el modelo. Si no lo es, usar alguna técnica simple, junto con la información geológica disponible para efectuar algunos cambios. En este paso, se deben analizar los diferentes mecanismos de depletación a fin de evaluarlo y ajustarlos para producir un ajuste de presión de todo el yacimiento. d.- Después que se logra un ajuste del yacimiento total, se debe llevar a cabo un ajuste de las mayores regiones del yacimiento. En esta etapa se refinan los parámetros de heterogeneidad del yacimiento, barreras al flujo y acuífero. e.- Dependiendo de los objetivos del estudio, se pueden obtener ajustes de presión para cada pozo.

AJUSTE DE GOR Y WOR La mejor indicación de validez del modelo en la representación del yacimiento, es el ajuste de GOR y WOR. El procedimiento usado para el ajuste puede variar de un yacimiento a otro, sin embargo se puede seguir el siguiente procedimiento: a.- Identificar los parámetros que influyen en el movimiento del agua y gas dentro del yacimiento y acuífero.

23

b.- Estimar los límites superior e inferior para cada parámetro basado en su incertidumbre. c.- Decidir, si es necesario, el uso de función de pozo, para simular ciertas condiciones tales como penetración parcial o conificación. El ajuste del comportamiento de un pozo en el cual el agua o gas rodea la zona de completamiento requerirá el uso de un modelo de conificación. El modelo se ajusta variando la permeabilidad en capas donde la incertidumbre es grande. La permeabilidad vertical es un factor de ajuste muy crítico en el modelo. d.- Examinar las corridas efectuadas en el ajuste de presión. Estas corridas pueden ser usadas para identificar la severidad de la estratificación, la cual requerirá ajuste de la permeabilidad vertical. La permeabilidad vertical no puede ser determinada en forma confiable de medidas de campo o núcleos. Se debe probar la sensibilidad del modelo a la permeabilidad vertical. e.- La distribución areal de la permeabilidad es otro factor importante y puede ser ajustado. f.- Decidir si se efectúa ajuste en los datos de permeabilidad relativa. g.- Determine el efecto de la dimensión del gridblock sobre el comportamiento de un grupo de pozos. Gridblocks grandes generan diferencias aparentes entre el modelo y el comportamiento del campo debido a los errores en el cálculo de las eficiencias de desplazamiento. h.- Cuando lleve a cabo estos cambios, continúe comparando el comportamiento de presión actual y calculada. El comportamiento de presión debe mantenerse mientras se ajusta el GOR y WOR.

AJUSTE DE LA PRESION DE LOS POZOS La dimensión de un bloque que contiene a un pozo productor o inyector en un simulador es normalmente mucho mayor que el radio del pozo. La presión de fondo medida, representa la presión a r = rw y al tiempo de la prueba. Por otro lado, la presión calculada representa la presión promedia dentro del bloque donde se encuentra el pozo al final de cualquier time step. Por lo tanto, la presión de fondo medida en un pozo activo, no puede ser comparada directamente con la presión estimada para el bloque.

24

VENTAJAS Y YACIMIENTOS

DESVENTAJAS

DE

LA

SIMULACIÓN

NUMÉRICA

DE

Algunas de las principales ventajas y/o aplicaciones de los modelos de simulación numérica son: •

A cada bloque se le pueden asignar valores únicos de propiedades de la roca, lo que permite tener en cuenta heterogeneidades y anisotropías del yacimiento.

•

A cada bloque o zona se les puede asignar valores de datos PVT, lo cual permite modelar la variación de las propiedades del fluido en el yacimiento.

•

Se tiene en cuenta el flujo de fluidos entre bloques adyacentes, lo que permite simular el movimiento de los frentes de fluidos en poryectos de inyección, los cambios en la posición del contacto gas-aceite en yacimientos con empuje hidráulico y los cambios en las distribuciones de presión y saturación del fluidos en ambos casos.

•

Se puede tener en cuenta la existencia de pozos inyectores o productores mediante la adición de los términos apropiados a las ecuaciones de flujo.

•

Permite determinar los mejores intervalos y tasas de producción en yacimientos con problemas de conificación.

•

Permite el estudio de variables involucradas en los procesos de recuperación, tales como arreglos y espaciamiento de pozos, intervalos de completamiento, tasas de producción, entre otros.

•

Permite modelar el comportamiento de sistemas pozo-arena productora cuando se produce de varias zonas aisladas.

•

Permite ubicar los pozos y las tasas de producción en lugares donde se explota un yacimiento por parte de varios operadores.

Algunas desventajas son: • Distorsión numérica como consecuencia de la división en bloques del yacimiento y asignación de nodos, dado que a cada nodo le corresponde un volumen de control de tamaño considerable, dentro del cual a las variables dependientes como presión y saturación se les asigna un valor único, lo cual no es consistente con lo que ocurre realmente en el yacimiento. Esta suposición puede conducir a errores considerables en los resultados obtenidos. Igualmente los simuladores consideran cambios abruptos de presión y saturación entre volúmenes de control consecutivos, lo cual tampoco ocurre en la realidad. Esta distorsión puede reducirse incrementado el número de nodos (es decir disminuyendo las dimensiones de los volúmenes de control). Sin embargo esta solución no es práctica, en muchas ocasiones, debido a que se incrementan los costos y los tiempos de ejecución del simulador. 25

• Errores de truncamiento debido a que las ecuaciones diferenciales empleadas son aproximadas por una serie de ecuaciones discretizadas, con lo que la solución del conjunto de ecuaciones numéricas difiere, en cierto grado, de la solución de la ecuación diferencial original • Error de redondeo acumulado debido a la gran cantidad de cálculos que se requieren para dar solución al sistema de ecuaciones discretizadas. • Falta de Infomación. Probablemente esta sea la principal limitación para correr el simulador. Por ejemplo, se necesita tener pleno conocimiento de la distribución de permeabilidad y porosidad en el yacimiento, con la finalidad de asignar a cada bloque valores representativos de la variación de estas propiedades a través del yacimiento. Los análisis de núcleos de formación, las pruebas de presión y los estudios geológicos del yacimiento son fundamentales en la obtención de esta información, pero no siempre se dispone de esta información en suficiente cuantía. Al seleccionar y asignar datos de entrada, se debe recordar que la veracidad de los resultados obtenidos depende de la veracidad de la información de entrada. Las ventajas y aplicaciones que ofrecen los modelos de simulación numérica conducen a pensar que los modelos convencionales, tales como la ecuación de balance de materiales y los modelos electrolíticos, han sido totalmente desplazados. Ciertamente, algunos modelos, como por ejemplo los modelos electrolíticos, han entrado en desuso debido a que son más costosos, de menor aplicabilidad y menos confiables que los modelos de simulación numérica. Sin embargo, es importante tener presente que no siempre la simulación numérica representa la mejor alternativa. Algunas veces es recomendable ó necesario aplicar modelos analíticos mas simples. Esto es particularmente cierto cuando no se tiene certeza acerca de la veracidad de los datos de entrada requeridos para correr el simulador o cuando los costos de correr el simulador no justifica la mejoría de los resultados obtenidos comparados con los resultados de un modelo más simple. Al respecto, Coats anota: "seleccione el modelo más simple que le permita obtener el cálculo deseado acerca del comportamiento del yacimiento". Mattax y Dalton señalan: "Los problemas deben ser resueltos por los métodos menos costosos y más simples, siempre y cuando conlleven a la respuesta adecuada. Por lo tanto, el ingeniero de yacimientos siempre debe determinar el nivel de simplificación primero y luego seleccionar el método de análisis apropiado con la finalidad de evitar 'sobre-trabajo' técnico." Adicionalmente, los modelos de laboratorio son de gran utilidad para investigar los procesos físicos que ocurren en el yacimiento. En este sentido, la

26

experimentación con modelos físicos y el desarrollo de modelos de simulación numérica son complementarios, no excluyentes. Finalmente, es importante anotar que la simulación de yacimientos es una valiosa herramienta de trabajo en el análisis del comportamiento de un yacimiento. Sin embargo, los simuladores numéricos no deben ser considerados como panaceas o "cajas negras" que producen resultados infalibles. El más sofisticado de los simuladores numéricos puede producir resultados totalmente erróneos si se aplica a un yacimiento cuyas características difieren de las suposiciones para las cuales se desarrolló el modelo, o si los datos de entrada no corresponden a la caracterización real del yacimiento, o si simplemente no se posee la destreza suficiente para interpretar los resultados obtenidos. Un simulador puede ser aplicable para describir adecuadamente el comportamiento de un yacimiento, en tanto que puede ser totalmente inapropiado para otro yacimiento, así los dos yacimientos presenten, aparentemente, similitud entre sí. El criterio de ingeniería es imprescindible en este caso. Al respecto, Crichlow señala: "En los procesos de simulación de un yacimiento el ingeniero está al tope de la situación". Nada de lo que el simulador haga puede mejorar la calidad de su trabajo, solamente le ayuda a adquirir un mejor entendimiento de los procesos que ocurren en su proyecto."

27

BIBLIOGRAFIA

1.- "Fundamentals of Numerical Reservoir Simulation" - Donald Peaceman Elsevier Scientific Publishing Company - 1977; 173 pag. 2.- "Reservoir Simulation" - Calvin C. Mattax and Robert L. Dalton - SPE Monograph Volume 13 - 1990; 161 pag. 3.- "Modern Reservoir Engineering - A Simulation Approach" - Henry B. Crichlow Prentice Hall Inc. - 1977 - 354 pag. 4.- “Petroleum Reservoir Simulation”. K. Azis and A. Setari. Elsevier 5.- “Principles of Aplied Reservoir Simularion”. Company, 1997.

J. Fanchi.

Gulf Publishing

6.- “Principles of Hydrocarbon Reservoir Simulation”. G. W. Thomas. International Human Resources Development Corporation, 1982. 7.- “Notas sobre simulación numérica de yacimientos”. Gildardo Osorio Gallego, Ph D. Universidad Nacional, Sede Medellín. Febrero de 2002.

28

1

2. ECUACIONES FUNDAMENTALES DE FLUJO EN MEDIOS POROSOS En su forma más simple, una ecuación fundamental de flujo en un medio poroso es una expresión analítica, diferencial, que expresa la variación de la presión o el potencial en función del espacio y el tiempo. Existen algunas condiciones que deben ser fijadas antes de definir completamente la forma de una ecuación fundamental de flujo. Estas condiciones, denominadas condiciones de flujo, son las siguientes: a. La geometría de flujo (lineal, radial, etc.) b. La naturaleza del fluido en movimiento (fluido incompresible, levemente compresible o compresible). c. Tipo de flujo (continuo, semi-continuo o no continuo). d. El número de fases en movimiento (flujo monofásico, bifásico o trifásico). Cada posible combinación de estas condiciones conlleva a una ecuación fundamental de flujo diferente. A continuación se discuten los principales casos que se pueden presentar dentro de cada una de estas categorías de condiciones de flujo. 2.1 GEOMETRIAS DE FLUJO La geometría de flujo en un medio poroso está determinada por la orientación de las líneas de flujo con respecto a un sistema de referencia. Las clases de geometría de flujo mas comúnmente utilizadas son: lineal, vertical, radial, bidimensional, tridimensional y esférico. 2.1.1 Flujo Lineal. Ocurre cuando las líneas de flujo son paralelas. Por ejemplo, la Figura 2.1 presenta un sistema acuífero lateral - yacimiento. La intrusión de agua desde el acuífero al yacimiento ocurre siguiendo líneas de flujo paralelas al eje x.

Agua

Acuífero

We

Petróleo

Yacimiento

Fig. 2.1 – Sistema Acuífero Lateral – Yacimiento (Ejemplo de Flujo Lineal).

2

2.1.2 Flujo Vertical. Es un caso particular del flujo lineal. Ocurre cuando las líneas de flujo son paralelas entre sí y verticales. Por ejemplo, la Figura 2.2 ilustra la trayectoria de las líneas de flujo que gobiernan la intrusión de agua desde un acuífero de fondo a un yacimiento.

Yacimiento

We Acuífero

Fig. 2.2 – Sistema Yacimiento – Acuífero de Fondo (Ejemplo de Flujo Vertical) 2.1.3 Flujo Radial. Se presenta cuando las líneas de flujo tienen forma de rectas localizadas en un plano horizontal, convergiendo hacia un punto central. Por ejemplo, considérese el flujo en las cercanías de un pozo que produce a través de todo el espesor de la formación (Figura 2.3).

Fig. 2.3 – Esquema que Indica el Flujo Radial de Fluidos hacia un Pozo. 2.1.4 Flujo Bidimensional. Se caracteriza por que las líneas de flujo se encuentran localizadas en un plano sin existir paralelismo entre ellas. Este tipo de flujo se puede representar en un sistema de ejes cartesianos (x, y ) . Por ejemplo, la Figura 2.4 ilustra la discretización areal de un yacimiento, tal como se hace en un simulador areal. Se simula que el flujo de bloque a bloque ocurre siguiendo líneas de flujo perpendiculares entre sí y paralelas a un par de ejes cartesianos (x, y ) .

3 y

Pozo

x

Fig. 2.4 – Discretización Areal de un Yacimiento 2.1.5 Flujo Tridimesional. Ocurre en tres dimensiones. Las líneas de flujo se pueden representar en un sistema de ejes cartesianos (x, y, z ) . Por ejemplo, la Figura 2.5 presenta la discretización de un yacimiento en tres dimensiones. El flujo ocurre siguiendo líneas paralelas a los ejes del sistema ortogonal cartesiano (x, y, z ) . z

y

Pozo

x

Fig. 2.5 – Ejemplo de Discretización de un Yacimiento en Tres Dimensiones 2.1.6 Flujo Esférico. Ocurre cuando las lineas de flujo son rectas localizadas en el espacio tridimensional y convergiendo hacia un punto central. Por ejemplo, considérese el flujo en las cercanías de un pozo en una formación de espesor considerable y el cual produce a través de perforaciones centradas en la formación (Figura 2.6) Pozo

h

Fig. 2.6 – Ejemplo de Flujo Esférico en las Cercanías de un Pozo Productor

4

2.2 TIPOS DE FLUIDOS EN UN YACIMIENTO Los fluidos en un yacimiento se suelen clasificar en una de las siguientes tres categorías: fluido incompresible, fluido levemente compresible y fluido incompresible). 2.2.1 Fluido Incompresible. Un fluido es incompresible cuando la variación en su volumen, debido a la variación en la presión, es insignificante, si se compara con la variación en el volumen de los otros fluidos presentes en la formación. Analíticamente,

dV = 0 ......................................................................................................... (2.1) dp Luego, c=−

1 dV = 0 ............................................................................................ (2.2) V dP

En las Ecuaciones 2.1 y 2.2, c es compresibilidad del fluido, V es volumen y p es presión. La compresibilidad del agua sin gas en solución es muy pequeña si se compara con las compresibilidades del petróleo y el gas natural a condiciones de yacimiento. Por esta razón, en simulación de yacimientos, el agua se suele considerar como un fluido incompresible. 2.2.2 Fluido Levemente Compresible. Se dice que un fluido es levemente compresible cuando el cambio en su volumen, debido a la variación de la presión, es constante. Analíticamente, c=−

1 dV = cons tan te ............................................................................. (2.3) V dP

La compresibilidad del petróleo a condiciones de yacimiento suele variar muy poco. Por lo anterior, en simulación de yacimientos, el petróleo a condiciones de yacimiento se suele tratar como un fluido levemente compresible. 2.2.3 Fluido Compresible. Un fluido se comporta como compresible cuando la variación de su volumen, ante un cambio de presión, es función de la presión. Analíticamente,

c=−

1 dV = f (p ) ......................................................................................... (2.4) V dP

El gas natural se comporta como un fluido altamente compresible.

5

2.3 TIPOS DE FLUJO Los tipos de flujo se definen de acuerdo a la forma como varía la presión con el tiempo en cada punto del yacimiento. En simulación de yacimientos se suelen considerar tres tipos de flujo: flujo no-continuo, flujo semi-continuo o pseudocontinuo y flujo continuo. 2.3.1 Flujo no - Continuo. Ocurre cuando la variación de la presión con respecto al tiempo, en cada punto del yacimiento, es función del tiempo. Analíticamente,

∂p   = f (t ) .................................................................................................. (2.5) ∂t  r En la Ecuación 2.5, t y r denotan tiempo y radio, respectivamente. Físicamente, un yacimiento produce bajo condiciones de flujo no-continuo cuando el yacimiento presenta comportamiento infinito. Es decir, desde que el pozo se abre a producción hasta el momento en el cual el radio de investigación llega a los límites externos del yacimiento. Nótese que durante este intervalo de tiempo, siempre existe una parte del sistema donde el radio de investigación no ha llegado y, en consecuencia, la presión del sistema se conserva uniforme e igual a la presión inicial. Las Figuras 2.7, 2.8 y 2.9 ilustran la posición del radio de investigación, la distribución de presión y la distribución de la rata de flujo, respectivamente, típicas de condiciones de flujo no-continuo. Pi t1

t2

t3

P ri1 ri2

ri3 r

Fig. 2.7 – Posición del radio de investigación en función del tiempo (flujo nocontinuo)

6

Pi t1

t2

t3

P ri1 ri2

ri3

r Fig. 2.8 – Forma típica de la distribución de presión en función del tiempo (flujo no – continuo)

t1

t2 Q t3 Q=0 rw

r

r e

Fig. 2.9 – Forma tipica de la distribución de la tasa de flujo en función del tiempo (flujo no – continuo)

2.3.2 Flujo pseudo-continuo o semi-continuo. Ocurre flujo semi-continuo cuando la variación de presión en cada punto del sistema es una función lineal del tiempo. Analíticamente,

7

∂p   = cons tan te ......................................................................................... (2.6) ∂t  r La Ecuación 2.6 se satisface a partir del momento en el cual el radio de investigación llega a los límites externos de un yacimiento cerrado. Por esta razón, se dice que ocurre flujo semi-continuo en sistemas de comportamiento finito, cerrado. La Figuras 2.10 presenta la forma típica de distribución de presión en función del tiempo en un yacimiento bajo condiciones de flujo semi-continuo.

Pi t1 t2

t3

P

rw

r

re

Fig. 2.10 – Forma típica de la distribución de presión en función del tiempo (flujo semi-continuo)

2.3.3 Flujo Continuo. Se dice que un yacimiento presenta comportamiento de flujo continuo cuando la presión, y por ende todas las propiedades que dependen de la presión, no varían con el tiempo. Este comportamiento es típico en yacimientos que presentan algún mecanismo de mantenimiento de presión; por ejemplo, en yacimientos con empuje hidráulico activo o con inyección de agua. Analíticamente, ∂p   = 0 ...................................................................................................... (2.7) ∂t  r

Cuando se presentan condiciones de flujo continuo, la rata de flujo es constante con la distancia.

8

Físicamente, para que ocurra flujo continuo se requiere que no se presente agotamiento de fluidos en el yacimiento, así que la presión se conserve constante con el tiempo. La cantidad de fluido producido y que tiende a bajar la presión debe ser sustituido con otro fluido que evite la disminución de presión. Esta es la razón por la cual se alcanzan condiciones de flujo continuo en yacimientos con empuje hidráulico o que producen por inyección de agua. En el caso de un yacimiento que produce por empuje hidráulico activo, el petróleo que se produce del yacimiento es reemplazado por la intrusión de agua proveniente del acuífero. En el caso de un sistema que produce por inyección de agua, al extraer determinada cantidad de crudo del yacimiento, se inyecta la misma cantidad de agua desde superficie. Este proceso de desplazamiento asume que los fluidos son incompresibles; sin embargo, desde el punto de vista de aplicación de un modelo, la condición primordial para que existan condiciones de flujo continuo es que la presión permanezca constante con el tiempo en cada punto del sistema. 2.4 FLUJO MONOFASICO Y MULTIFASICO Ocurre flujo monofásico cuando existe flujo de un sólo fluido en el medio poroso (petróleo, agua o gas). En caso de que el medio poroso esté ocupado por varios fluidos, pero solamente uno de ellos sea móvil, se considera que el flujo es monofásico. Ocurre flujo multifásico cuando dos o más fluidos fluyen en el medio poroso. 2.5 ECUACIONES DE ESTADO PARA DIFERENTES TIPOS DE FLUIDOS Una ecuación de estado es una expresión algebraica que describe la relación entre la presión, el volumen y la temperatura (PVT) de una sustancia. Esta relación es diferente dependiendo de la naturaleza del fluido. En esta sección se discute brevemente la forma que toman las ecuaciones de estado para fluidos incompresibles, levemente compresibles y compresibles, en la forma como éstas suelen ser empleadas en simulación de yacimientos. 2.5.1 Ecuación de Estado para un Fluido Incompresible. Tal como se expresa en la Ecuación 2.1, para un fluido incompresible se cumple que dV dp = 0 , por consiguiente: d (V m ) d (1 ρ ) 1 dρ = =− 2 =0 dp dp ρ dp

Lo que se cumple solamente si

ρ = cons tan te .............................................................................................. (2.8) En la Ecuación 2.8, ρ es densidad.

9

En este libro se hará uso de la Ecuación 2.8 cada vez que se requiere utilizar la ecuación de estado para un fluido incompresible. 2.5.2 Ecuación de Estado para un Fluido Levemente Compresible. Ecuación 2.3, para un fluido levemente compresible se tiene:

c=−

De la

1 dV = cons tan te V dP

de donde c=−

c=

 1 dρ  d (1 ρ ) 1 d (V ) m m d (V m ) =− = −ρ = − ρ − 2  =0 V dp m V dp dp  ρ dp 

1 dρ = cons tan te ................................................................................... (2.9) ρ dp

Separando variables e integrando se obtiene:

ρ = ρ o e c ( p− p

o

)

............................................................................................ (2.10)

En la Ecuación 2.10, ρ 0 y p0 son la densidad y presión de referencia, respectivamente. En este libro se utilizará la Ecuación 2.10 cada vez que se requiera aplicar la ecuación de estado de un fluido levemente compresible. Algunas veces, la Ecuación 2.10 suele simplificarse como se indica a continuación. Una función f ( x ) puede expandirse mediante una serie de Maclaurin así: ∞

f (x ) = ∑ n =0

f (n ) (0 ) n x ................................................................................... (2.11) n!

En la Ecuación 2.11, f (n ) (0 ) representa la n − ésima derivada de f con respecto a x evaluada en el punto x = 0 . Si se hace la variable x = p y f = ρ , la Ecuación 2.11 puede ser escrita como: ∞

ρ (n ) (0 )

n =0

n!

∞

ρ (n ) (0 )

n =0

n!

ρ ( p) = ∑

ρ ( p) = ∑

ρ ( p) = ρ 0 e

− cp 0

p n ................................................................................... (2.12) ∞

(cp )n

n =0

n!

p n = ρ 0 e −cp0 = ∑

 (cp )2 + (cp )3 + (cp )4 + K + (cp )n  1 + cp +  2! 3! 4! n!  

10

Dado a que se trata de un fluido levemente compresible, la compresibilidad es 2 3 muy pequeña. Por esta razón, las componentes (cp ) / 2 ! , (cp ) / 3 ! , …,

(cp )n / n !

tienden a cero, lo que simplifica la ecuación anteior:

ρ ( p) = ρ o e

− cp o

(1 + cp ) .................................................................................

(2.13)

En particular, cuando p0 = 0 se tiene:

ρ ( p ) = ρ o (1 + cp ) ......................................................................................... (2.14) La Ecuación 2.14 es una simplificación de la Ecuación 2.10 y en algunas ocasiones se suele tomar como la ecuación de estado para fluidos levemente compresibles. 2.5.3 Ecuación de Estado para un Fluido Compresible. Para un fluido compresible, el volumen depende fuertemente de la presión. La ecuación de estado para fluidos compresibles se suele expresar en términos de la ecuación de estado para gases reales: pV = znRT = z

p=z

ρ=

m RT M

m RT RT = zρ V M M

pM .................................................................................................... (2.15) zRT

En la Ecuación 2.15, z es el factor de compresibilidad del gas, M es peso molecular, R es la constante universal de los gases, y T es temperatura.

2.6 LA LEY DE DARCY La ley de Darcy permite modelar el flujo laminar a través de un medio poroso relacionando la tasa de flujo con el gradiente de presión. Esta relación establece que la velocidad aparente de flujo (rata de flujo por unidad de área total) es proporcional al gradiente de presión e inversamente proporcional a la viscosidad de fluido. En general, la ley de Darcy puede ser escrita de la siguiente forma:

u=−

k dΦ ................................................................................................ (2.16) µ dl

11

En la Ecuación 2.16, u es rata volumétrica de flujo (volumen por unidad de área perpendicular al flujo por unidad de tiempo), k es permeabilidad, µ es viscosidad del fluido y dΦ d l es gradiente del potencial de flujo. El potencial Φ está dado por la siguiente expresión:

Φ = p + ρgz ............................................................................................... (2.17) En la Ecuación 2.17, g es la constante gravitacional y z es elevación. En simulación de yacimientos, la Ecuación 2.16 suele aplicarse en dos sistemas de coordenadas: cartesianas y radiales. 2.6.1 Forma de la Ley de Darcy en Coordenadas Cartesianas. La Ecuación 2.16 puede ser aplicada para expresar la velocidad volumétrica de un fluido monofásico, coordenadas cartesianas, de la siguiente forma: En dirección x :

ux = −

kx ∂p .............................................................................................. (2.18) µ ∂x

En dirección y :

uy = −

k y ∂p ............................................................................................... (2.19) µ ∂y

En dirección z :

uz = −

k z  ∂p  + ρg  .................................................................................... (2.20)  µ  ∂z 

En las Ecuaciones 2.18 a 2.20, las variables x , y y z hacen referencia a las direcciones x , y y z , respectivamente. Si ocurre flujo multifásico, la ley de Darcy puede ser aplicada para expresar la velocidad volumétrica de flujo de cada una de las fases. Por ejemplo, para el petróleo se tiene: En dirección x : uo x = −

k o x ∂p o .......................................................................................... (2.21) µ o ∂x

En dirección y :

12

u oy = −

k o y ∂p o ........................................................................................... (2.22) µ o ∂y

En dirección z : u oz = −

k oz  ∂p o  + ρ o g  ............................................................................... (2.23)  µ o  ∂z 

En las Ecuaciones 2.21 a 2.23, el subíndice " o " hace referencia a la variable evaluada para la fase del petróleo. Las permeabilidades presentes en las Ecuaciones 2.21 a 2.23 hacen referencia a permeabilidades efectivas. En términos de permeabilidades relativas, las Ecuaciones 2.21 a 2.23 pueden ser escritas de la siguiente forma: En dirección x : u ox = −

kk rox ∂p o ......................................................................................... (2.24) µ o ∂x

En dirección y : u oy = −

kk roy ∂p o ....................................................................................... (2.25) µ o ∂y

En dirección z : u oz = −

kk roz  ∂p o  + ρ o g  ............................................................................. (2.26) µ o  ∂z 

En las Ecuaciones 2.24 a 2.26, el subíndice " r " indica que la permeabilidad en cuestión es relativa. Expresiones análogas a las Ecuaciones 2.18 a 2.26 pueden ser obtenidas para el agua y el gas. 2.6.2 Forma de la Ley de Darcy para Flujo Radial. La Ecuación 2.16 puede ser aplicada al flujo radial de un fluido monofásico. En este caso, la Ecuación 2.16 toma la siguiente forma: ur = −

k r ∂p ............................................................................................... (2.27) µ ∂r

En forma análoga al caso de coordenadas cartesianas, la Ecuación 2.16 puede ser aplicada a flujo radial, mutifásico. Por ejemplo, para el petróleo se tiene:

13

uor = −

kk ro ∂p o .......................................................................................... (2.28) µ o ∂r

Se pueden obtener expresiones análogas para flujo radial de gas y agua. 2.7 FORMA NO-LINEAL DE LAS ECUACIONES FUNDAMENTALES DE FLUJO MONOFASICO Una ecuación fundamental de flujo es una forma diferencial de expresar la ley de la conservación de la masa. La forma que toma la ecuación depende de la geometría de flujo, la naturaleza del fluido en movimiento, el número de fases en el fluido y el tipo de flujo. Esta sección discute la forma no-lineal de la ecuación fundamental de flujo monofásico en coordenadas cartesianas y coordenadas radiales. 2.7.1 Ecuación Fundamental de Flujo Monofásico en Coordenadas Cartesianas. Supóngase un elemento infinitesimal de volumen en un medio poroso continuo a través del cual ocurre flujo monofásico en la dirección x , tal como se ilustra en la Figura 2.11.

A( X ) ρu x + ∆( ρu x )

ρu x ∆X X

X + ∆X

Fig. 2.11 – Elemento infinitesimal de volumen, flujo lineal, coordenadas Cartesianas Sea u x la velocidad volumétrica de flujo en la dirección x , definida así: ux =

Volumen ..................................................................................... (2.29) Tiempo ⋅ Area

El flujo másico, definido como el flujo de masa por unidad de tiempo por unidad de área, está dado por:

14

Flujo másico = ρu x ..................................................................................... (2.30)

Si se efectúa un balance de la masa dentro del elemento infinitesimal durante un intervalo de tiempo ∆t , se tiene:  Acumulación (+ ) Cantidad de masa  Cantidad     − de masa  ± que entra o sale por  =  o         ∆t que sale  ∆t  fuentes o sumideros  ∆t  Agotamiento (- )  ∆t ................................................................................................................... (2.31)

Cantidad de masa  que entra

La cantidad de masa que entra al elemento durante el tiempo de observación ∆t será:

Cantidad de masa  que entra

  = ρu A( x )∆t ......................................................................... (2.32) x   ∆ t

La cantidad de masa que sale del elemento durante el mismo intervalo de tiempo de observación está dada por:

Cantidad de masa  que sale

  = [ρu A( x ) + ∆ (ρu A( x ))]∆ t .................................................... (2.33) x x   ∆ t

Notando como q~ a la cantidad de masa que entra o sale por fuentes o sumideros, por unidad de volumen del yacimiento por unidad de tiempo, se tiene que:

Cantidad de masa  que entra o sale por  = ± q~∆V∆t = ± q~A(x )∆x∆t ......................................... (2.34)    fuentes o sumideros  ∆t De otro lado, la cantidad de masa existente en el elemento infinitesimal a un tiempo t está dada por:

[A(x )∆ xφρ ] t

............................................................................................... (2.35)

Similarmente, la cantidad de masa existente en el elemento infinitesimal a un tiempo t + ∆t , será:

[A(x )∆xφρ ]t +∆ t

............................................................................................ (2.36)

La acumulación o agotamiento de masa durante el intervalo de tiempo ∆t será:

15

 Acumulación (+ )   = [ A( x )∆xφρ ] − [A( x )∆xφρ ] ........................................ (2.37) o t + ∆t ∆t    Agotamiento (- ) 

Sustituyendo las Ecuaciones 2.32 a 2.34 y 2.37 en la Ecuación 2.31, asumiendo que el volumen del elemento infinitesimal, A( x )⋅ ∆x , no es función del tiempo, y dividiendo ambos lados de la ecuación resultante por ∆x ⋅ ∆t , se obtiene: − ∆ ( ρu x A( x )) ~ ∆ (ρφ ) ± q A(x ) = A( x ) ∆x ∆t

Simplificando y tomando límites cuando los incrementos infinitesimales tienden a cero se obtiene: −

∂ (ρu x A(x )) ~ ∂ (ρφ ) ............................................................... (2.38) ± q A( x ) = A(x ) ∂x ∂t

En forma general, la Ecuación 2.38 suele escribirse como: −

∂ (αρu x ) ∂ (ρφ ) =α + αq~ ........................................................................... (2.39) ∂x ∂t

En la Ecuación 2.39, α = A( x ) . En la Ecuación 2.39, la variable q~ es positiva para sumideros (por ejemplo, pozos productores) y negativa para fuentes (por ejemplo, pozos inyectores). La Ecuación 2.39 es válida para flujo en una dimensión. Si se considera un elemento infinitesimal a través del cual ocurre flujo monofásico en dos direcciones, tal como se ilustra en la Figura 2.12, y se sigue un procedimiento similar al anterior, se encuentra que la ecuación fundamental de flujo en dos dimensiones puede ser escrita como:

16

ρu y

ρu x

ρu x + ∆( ρu x )

H ( x, y )

∆y

ρu y + ∆(ρu y )

∆x

Fig. 2.12 – Elemento infinitesimal de volumen, flujo bidimensional, coordenadas Cartesianas  ∂ (αρu x ) ∂ (αρu y ) ∂ ( ρφ ) − + + αq~ ......................................................... (2.40)  =α x ∂ y ∂ t ∂  

En la Ecuación 2.39, α = H ( x , y ) , espesor en el punto ( x , y ) . Si se considera un elemento infinitesimal a través del cual ocurre flujo monofásico en tres direcciones, tal como se ilustra en la Figura 2.13, se encuentra que la correspondiente ecuación fundamental de flujo tiene la siguiente forma:  ∂ ( ρu x ) ∂ (ρu y ) ∂ (ρu z )  ∂ (ρφ ) ~ − + + + q ................................................... (2.41) = ∂z  ∂t ∂y  ∂x

17

ρu z + ∆ ( ρu z )

ρu y

ρu x

ρu x + ∆( ρu x )

∆z

∆y

ρu y + ∆(ρu y )

∆x

ρu z

Fig. 2.13 – Elemento infinitesimal de volumen, flujo tridimensional, coordenadas Cartesianas Obsérvese que la ecuación fundamental de flujo monofásico en coordenadas cartesianas puede ser escrita, en forma general, como:  ∂ (αρu x ) ∂ (αρu y ) ∂ (αρu z )  ∂ ( ρφ ) − + + + αq~ .......................................... (2.42)  =α x y z t ∂ ∂ ∂ ∂  

En la Ecuación 2.42, α = A( x ) y u y = u z = 0 , para flujo en una dirección;

α = H ( x , y ) y u z = 0 , para flujo en dos direcciones; α = 1.0 para flujo en tres direcciones. En términos de divergencia y gradiente, la ecuación fundamental de flujo monofásico en coordenadas cartesianas puede ser escrita, en forma general, como:

− ∇ ⋅ (αρu ) = α

∂ ( ρφ ) + αq~ ........................................................................... (2.43) ∂t

En la Ecuación 2.43 ∇ ⋅ es la divergencia en coordenadas cilíndricas y u es el vector de velocidad volumétrica los cuales se definen, respectivamente, de la siguiente forma:

∇⋅ =

∂ ∂ ∂ i + j + k .................................................................................. (2.44) ∂x ∂x ∂x

18

u = u x i + u y j + u z k ....................................................................................... (2.45) En las Ecuaciones 2.44 y 2.45, i , j y k son los vectores unitarios en las direcciones x , y y z , respectivamente. Sustituyendo la Ley de Darcy en coordenadas cartesianas, Ecuaciones 2.18 a 2.20, en la Ecuación 2.45 y la Ecuación resultante en la Ecuación 2.43, se obtiene:  k  ∂ ( ρφ ) − ∇ ⋅ αρ ∇p  = α + αq~ ................................................................... (2.46) ∂ t µ  

En la Ecuación 2.46, ∇ representa el gradiente, definido como:

∇p =

∂p ∂p ∂p i+ j + k ............................................................................... (2.47) ∂x ∂x ∂x

k representa el tensor de permeabilidades. Obsérvese que la Ecuación 2.46 puede ser escrita como: ∂  k x ∂p  ∂  k y ∂p  ∂  k z ∂p  ∂ ( ρφ )  + αρ αρ  + αρ  = α + αq~ ...................... (2.48) ∂x  ∂t µ ∂x  ∂y  µ ∂y  ∂z  µ ∂z 

La Ecuación 2.48 se conoce como la ecuación fundamental de flujo monofásico en coordenadas cartesianas. 2.7.2 Ecuación Fundamental de Flujo Monofásico en Coordenadas Radiales. Considérese un elemento infinitesimal de volumen de un medio poroso continuo a través del cual ocurre flujo monofásico, radial, tal como se ilustra en la Figura 2.14.

19

ρ u r + ∆ (ρ u r ) H (r ,θ )

ρu r

∆r

θ

r

Fig. 2.14 – Elemento infinitesimal de volumen en coordenadas radiales La cantidad de masa que ingresa al elemento durante un intervalo de tiempo ∆t , a un radio r , está dada por: Cantidad de masa  que entra

  = ρu A(r ,θ )∆t = ρu θrH (r ,θ )∆t ............................................. (2.49) r r   ∆t

Similarmente, la cantidad de masa que sale del elemento durante el mismo intervalo de tiempo ∆t , a un radio r + ∆r , está dada por:

Cantidad de masa  que sale

  = ρu H (r ,θ ) + ∆ (ρu H (r ,θ )) (r + ∆r )θ∆ t ................................ (2.50) r r   ∆t

[

]

Notando como q~ a la cantidad de masa que entra o sale por fuentes o sumideros, por unidad de volumen del yacimiento por unidad de tiempo, se tiene que: Cantidad de masa que entra o sale por   fuentes o sumideros

  = ± q~∆V∆t ............................................................. (2.51)   ∆t

En la Ecuación 2.51, el volumen infinitesimal ∆V será:

20

∆V = π (r + ∆r )2 H (r ,θ )

θ θ − πr 2 H (r ,θ ) 2π 2π

θ θ ∆V = 2r∆rH (r ,θ ) − π (∆r )2 H (r ,θ ) 2

2

Debido a que ∆r es un cantidad infinitesimal, (∆r ) es aproximadamente igual a cero comparado con ∆r . En consecuencia: 2

∆V = r∆ rH (r ,θ )θ ....................................................................................... (2.52) Sustituyendo la Ecuación 2.52 en la Ecuación 2.51, se obtiene:

Cantidad de masa  que entra o sale por  = ± q~r∆rH (r ,θ )θ∆t ............................................... (2.53)    fuentes o sumideros  ∆ t El acumulamiento o agotamiento en el elemento infinitesimal durante el intervalo ∆t , será:

 Acumulación    = [∆Vφρ ] − [∆Vφρ ] .................................................... (2.54) o t + ∆t t    Agotamiento  ∆t Llevando la Ecuación 2.53 a la Ecuación 2.54, se obtiene:  Acumulación    = r∆rH (r ,θ )θ [( ρφ ) − ( ρφ ) ] .......................................... (2.55) o t + ∆t t    Agotamiento  ∆t

Sustituyendo las Ecuaciones 2.49, 2.50, 2.53 y 2.55 en la Ecuación 2.31, se obtiene:

ρu rθrH (r ,θ )∆t − [ρu r H (r ,θ ) + ∆ (ρu r H (r ,θ ))](r + ∆r )θ∆t ± q~r∆rH (r ,θ )θ∆t = r∆rH (r ,θ )θ [( ρφ )t + ∆t − ( ρφ )t ] o bien,

− ρu r H (r ,θ )∆rθ∆t − ∆ (H (r ,θ )ρu r )rθ∆t − ∆ (H (r ,θ )ρu r )∆rθ∆t ± q~r∆rH (r ,θ )θ∆t = r∆rH (r ,θ )θ [( ρφ ) − (ρφ ) ] t + ∆t

t

Debido a que cada diferencia finita hace referencia a una cantidad infinitesimal, las diferencias de tercer orden son despreciables en comparación con las diferencias de segundo orden. Por consiguiente, la ecuación anterior se simplifica a:

21

− ρu r H (r ,θ )∆ rθ∆t − ∆ (H (r ,θ )ρu r )rθ∆t ± q~r∆rH (r ,θ )θ∆t =

r∆rH (r ,θ )θ [( ρφ )t + ∆t − ( ρφ )t ]

Dividiendo ambos lados entre r ⋅ ∆r ⋅ ∆t ⋅θ , se obtiene:

[(ρφ )t + ∆t − (ρφ )t ] ∆ (H (r ,θ )ρu r ) 1 −  ρu r H (r ,θ ) + r ± H (r ,θ )q~ = H (r ,θ )  r ∆r ∆t  Tomando límites cuando los incrementos infinitesimales tienden a cero, se obtiene:  1 ∂ ∂ −  H (r ,θ )ρu r + r (H (r ,θ )ρu r ) ± H (r ,θ )q~ = H (r ,θ ) ( ρφ ) r ∂t ∂r 

o bien, −

1 ∂ [H (r ,θ )ρu r r ] = H (r ,θ )q~ + H (r ,θ ) ∂ (ρφ ) .............................................. (2.56) r ∂r ∂t

Donde q~ es la cantidad de masa que entra o sale, a través de fuentes o sumideros, por unidad de volumen del yacimiento, por unidad de tiempo. q~ es positivo para sumideros (por ejemplo, pozos productores) y negativo para fuentes (por ejemplo, pozos inyectores). La Ecuación 2.56 puede ser escrita como: −

1 ∂ [αρu r r ] = αq~ + α ∂ (ρφ ) ................................................................... (2.57) r ∂r ∂t

En la Ecuación 2.57 α = H (r ,θ ) . La Ecuación 2.57 se conoce como la ecuación fundamental de flujo monofásico en coordenadas radiales. 2.7.3 Ecuación Fundamental de Flujo Monofásico en Coordenadas Cilíndricas. Siguiendo una metodología análoga a la seguida en los numerales 2.7.1 y 2.7.2, es posible obtener la siguiente ecuación fundamental de flujo monofásico en coordenadas cilíndricas:  k  ∂ ( ρφ ) − ∇ ⋅ αρ ∇p  = α + αq~ ................................................................... (2.58) ∂t  µ 

En la Ecuación 2.58, ∇ ⋅ W y ∇V representan la divergencia y el gradiente en coordenadas cilíndricas de W y V , respectivamente, definidos de acuerdo a las siguientes expresiones:

22

∇⋅W =

1 ∂ (rWr ) 1 ∂Wθ ∂Wz + + .................................................................. (2.59) r ∂r r ∂θ ∂z

y

∇V =

∂V 1 ∂V ∂V i+ j+ k .......................................................................... (2.60) ∂r ∂z r ∂θ

En las Ecuaciones 2.59 y 2.60 r , θ y z indican dirección radial, tangencial y vertical, respectivamente; Wr , Wθ y Wz son las componentes radial, tangencial y vertical, respectivamente, del vector W ; α es igual a 1.0. Obsérvese que la Ecuación 2.58 puede ser escrita como:

1  ∂  kr ∂p  1 ∂  +   rρ r  ∂r  µ ∂r  r ∂θ

 kθ ∂p  ∂  k z ∂p  ∂ (ρφ )  ρ  +  ρ  = α + αq~ ................. (2.61) ∂t  µ ∂θ  ∂z  µ ∂z 

3. APROXIMACIÓN DE ECUACIONES DIFERENCIALES A DIFERENCIAS FINITAS

Las ecuaciones que gobiernan el flujo de fluidos en medios porosos fueron derivadas en los capítulos anteriores.

Estas ecuaciones son ecuaciones

diferenciales parciales no lineales, las cuales relacionan los cambios de presión y de saturación con el tiempo a lo largo del medio poroso. Estas ecuaciones son extremadamente complejas, y sus aplicaciones son complicadas por la presencia de condiciones límite especializadas. La solución de estas ecuaciones por medios analíticos es, generalmente, imposible. Las soluciones, cuando existen, proporcionan una distribución continua de los parámetros dependientes (presión o saturación). En la mayoría de las aplicaciones la única manera de obtener solución de estas ecuaciones es mediante la solución numérica. Esta solución numérica produce resultados en puntos discretos dentro del sistema.

La transformación de la

ecuación diferencial continua a una forma discreta se logra mediante el uso de diferencias finitas.

En este proceso tanto el espacio como el tiempo, son

discretizados. Este capítulo se dedica a la discusión de los conceptos básicos que rigen la aproximación

de

una

ecuación

diferencial

parcial

a

diferencias

finitas.

Inicialmente, se discuten algunas nociones sobre variables discretas y diferencias finitas.

Luego, se consideran las aproximaciones para la primera y segunda

derivada, así como los esquemas de aproximación implícita y explícita. Posteriormente, se discuten los conceptos de error de truncamiento, estabilidad, convergencia y consistencia. En la parte final del capítulo, se habla de sistemas de malla y condiciones de límite. A través de todo el capítulo se trabaja con una ecuación de flujo simple con la finalidad de introducir los conceptos básicos referentes a los temas antes mencionados.

Simulación de Yacimientos I

Esto facilita considerablemente la

1

aplicación y extensión de conceptos básicos a ecuaciones de flujo más complejas las cuales gobiernan el comportamiento real de un yacimiento.

3.1 DISCRETIZACION, VARIABLES DISCRETAS Y DIFERENCIAS FINITAS La solución a los sistemas de ecuaciones de flujo comúnmente encontrados en el trabajo de ingeniería de yacimientos involucra la determinación de algunos parámetros dependientes del espacio y del tiempo.

Como se mencionó

anteriormente, la solución se obtiene en puntos discretos en el espacio y en el tiempo. El dominio espacial se divide en un número de celdas, grids o bloques superponiendo algún tipo de malla. Este enmallado, generalmente, es rectangular pero no necesariamente. El tiempo también se discretiza en un número de pasos de tiempo, durante cada uno de los cuales se resuelve el problema para obtener nuevos valores de los parámetros dependientes.

El tamaño de estos pasos

depende del problema particular que se está resolviendo, y generalmente, entre más pequeño sea este paso, más exacta es la solución. Para aclarar este proceso considérese un medio poroso lineal a través del cual fluye un fluido compresible o levemente compresible que se encuentra a una presión inicial p i (Figura 3.1). Si el sistema se abre a producción en el punto x= L, la presión disminuye a medida que se incrementa el volumen de fluido extraído del sistema.

Supóngase que para este sistema, la variación de la presión p, en

función de la posición x y el tiempo t, se puede representar por la siguiente ecuación diferencial: ∂ 2 p ∂p = ∂x 2 ∂t

Simulación de Yacimientos I

(3.1)

2

X=0

t4

P

X=0

X=L

X

t3

t2

t1

X

X=L

Fig. 3.1 – Forma típica de la distribución de presión en un sistema lineal La solución analítica a la Ecuación 3.1 conlleva una función continua p = p(x, t ) que permite calcular la presión en cada punto a un tiempo determinado.

Un

gráfico de p en función de x y t se muestra en la Figura 3.1. Para la solución numérica de una ecuación diferencial de flujo, tal como la Ecuación 3.1, se requiere tratar el yacimiento como un conjunto de bloques. Esta solución permite calcular la presión en cada bloque a diferentes intervalos de tiempo. A diferencia de la solución analítica, en lugar de obtenerse una expresión que relacione la presión como una función continua de la posición y el tiempo, se obtienen una serie de valores discretos de presión, cada uno de los cuales corresponde a un bloque determinado del yacimiento.

La acción de dividir el

yacimiento en un número determinado de bloques e intervalos de tiempo se denomina discretización. La Figura 3.2 presenta el ejemplo hipotético de un yacimiento real y su correspondiente discretización. A las longitudes de cada bloque, ∆ x 1 , ∆ x 2 , ... , y

a cada intervalo de tiempo, ∆ t 1 , ∆ t 2 , ... , se les denomina diferencias finitas.

Simulación de Yacimientos I

3

Barrera permeable

φ1

φ2 k P2 2 ∆x 2

k 1 P1

∆x 1 Bloque 1

φ4

P φ3 k 3 3 ∆x 3

Bloque 2

k4

∆x 4

Bloque 3

P4

Bloque 4

Fig. 3.2 – Esquema de un sistema real y su respectiva discretización

X

X=0

t3 t4

P

X=0

t3 t4

t2 t3 t4

X=L

t1

t2

t2

t3

t3

t4

t4

X

t1 t2 t3 t4

t1 t2 t3 t4

t1 t2 t3 t4 X=L

Fig. 3.3 – Forma típica de la distribución de presión en un sistema lineal (solución numérica)

Simulación de Yacimientos I

4

La Figura 3.3 presenta la forma típica de la distribución de presión en un sistema lineal a diferentes tiempos, obtenida mediante la solución numérica de una ecuación fundamental de flujo. A manera de ejemplo, se pudiera considerar que las Figuras 3.1 y 3.3 corresponden a las soluciones analítica y numérica, respectivamente, de una misma ecuación de flujo (por ejemplo de la Ecuación 3.1) para un mismo sistema. Los siguientes comentarios, concernientes a la discretización de un yacimiento, son importantes: a. Los lados de cada bloque, perpendiculares a la dirección de flujo, actúan como barreras de espesor infinitesimal (véase Figura 3.2). b. La rata de entrada y salida de fluidos en cada bloque está determinada por la permeabilidad de las barreras adyacentes a cada bloque y la diferencia de presión entre ellos. c. Las propiedades dentro de cada bloque son las mismas en todos los puntos del bloque.

Por ejemplo: a un tiempo determinado, a cada bloque le

corresponde un único valor de presión y de permeabilidad relativa. d. La variación de las propiedades en el yacimiento se representa mediante la variación de las propiedades en cada bloque. Por ejemplo: a cada bloque se le puede asignar un valor diferente de porosidad. Esto implica que pueden existir variaciones fuertes en las propiedades de bloques adyacentes. e. La precisión con la cual el comportamiento real del yacimiento puede ser descrito por el modelo depende del número de bloques utilizados.

Sin

embargo, el número de bloques está limitado por factores tales como el costo del simulador y la disponibilidad de datos de entrada. f.

La vida del yacimiento se discretiza en intervalos de tiempo. Las condiciones del yacimiento se definen únicamente al principio y al final de cada intervalo de

Simulación de Yacimientos I

5

tiempo, razón por la cual las condiciones en cada bloque pueden cambiar considerablemente de un intervalo de tiempo a otro.

Intervalos de corta

duración (es decir, alto número de intervalos), hacen que estos cambios sean menos fuertes y aumentan la precisión de los resultados. g. La discretización del yacimiento hace que el comportamiento continuo de las condiciones en él sea distorsionado. Por ejemplo: la Figura 3.4 presenta un esquema de la distorsión ocasionada en la saturación de agua en un punto fijo de un sistema lineal debido a su discretización. 1

: Yacimiento : Simulador

Sw

0 Fig. 3.4 – Distorsión ocasionada en la saturación de agua en un punto del yacimiento

3.2 APROXIMACIONES PARA LA PRIMERA DERIVADA

El método más comúnmente utilizado para dar solución numérica a la Ecuación 3.1 es expresar la función p(x, t ) mediante una serie de Taylor. Considérese una función f (x ) continua. La función puede ser expandida alrededor de un punto “a”, localizado a una distancia x-a del punto “x” (véase la Figura 3.5), mediante una serie de Taylor:

Simulación de Yacimientos I

6

x

x–a a

X=0

X

x

X

i

X=L

i+ 1

Fig. 3.5 – Bloques tenidos en cuenta en la aproximación progresiva de la primera derivada f (x) =

f ( n )(a) (x − a ) n ∑ n! n=0 ∞

(3.2)

En la Ecuación 3.2, f (n ) es la n-ésima derivada de f evaluada en el punto x=a. Es decir, la Ecuación 3.2 puede ser escrita como:

(x − a) f (x ) = f (a )⋅ 0!

0

+

(x − a) + ∂ 2 f ∂f   ⋅ ∂x  x = a 1! ∂x 2

3 ( ∂ f ∂nf x − a)  + ⋅ +L+ 3! ∂ x 3  x = a ∂x n 3

 (x − a) 2  ⋅ 2! x=a

n  ( x − a)  ⋅ +L n! x=a

(3.3)

3.2.1 La aproximación progresiva y su error de truncamiento

Supóngase un sistema lineal de bloques de igual longitud, tal como se ilustra en la Figura 3.5. A cada bloque le corresponde una variable discreta x 1 , x 2 , ... , x i − 1 , x i , x i + 1 , ..., x n , de acuerdo a su posición dentro del sistema. El incremento entre dos puntos sucesivos, x i y x i + 1 , está dado por la diferencia finita ∆ x . Si se define:

Simulación de Yacimientos I

7

a = xi

y

x = x i+1

Entonces, x − a = x i+1 − x i

x − a = ∆x

(por ser bloques de igual longitud)

y si además se tiene en cuenta que la función, f (x ) , que se pretende expandir mediante la serie de Taylor es la presión en el sistema, P, entonces, de la Ecuación 3.3 se tiene: P (x i + 1 ) = P (x i ) +

2 3 n ∂P  ∂ 2 P  (∆ x ) ∂ 3 P  (∆ x ) ∂ n P  (∆ x ) (3.5)     ⋅ (∆ x ) + ⋅ + ⋅ + + ⋅ L ∂ x  Xi ∂ x 2  Xi 2 ! ∂ x 3  Xi 3 ! ∂ x n  Xi n !

donde: P (x i + 1 ) : Presión en el bloque i + 1 . P (x i ) : Presión en el bloque i .

Para simplificar la notación, se hace uso de las siguientes expresiones P (x i + 1 ) = Pi + 1 P (x i ) = Pi ∂P  ∂P   =  ∂ x  Xi ∂ x  i

Simulación de Yacimientos I

8

Luego, la Ecuación 3.5 se convierte en:

P i + 1 = Pi +

2 3 n ∂P  (∆ x ) ∂ 2 P  (∆ x ) ∂ 3 P  (∆ x ) ∂ n P  (∆ x )     ⋅ + ⋅ + ⋅ + + ⋅ L ∂ x  i 1! ∂ x 2  i 2! ∂ x 3  i 3 ! ∂ x n  i n!

despejando

(3.6)

∂P   de la ecuación (3.6) se llega a: ∂ x  i

2 n−1 ∂ n P  (∆ x ) ∂P  P i + 1 − Pi ∂ 2 P  (∆ x ) ∂ 3 P  (∆ x )  ⋅  ⋅  ⋅  = +L −L− − − n! ∂ x n  i ∂ x 3  i 3 ! ∂ x  i ∆x ∂ x 2  i 2 !

(3.7)

la cual puede ser agrupada de la siguiente manera: ∂P  Pi + 1 − Pi  = + R if ∂x  i ∆x

(3.8)

donde:  ∂ 2 P  (∆ x ) ∂ 3 P  (∆ x ) 2   R if = −  2  ⋅ + ⋅ + L  ∂ x 3  i 3 !  ∂ x  i 2! 

(3.9)

El valor de R if tiende a cero cuando ∆ x tiende a cero, por lo tanto a partir de la Ecuación 3.8, se puede afirmar que ∂P  Pi + 1 − Pi  ≅ ∂ x  i ∆x

(3.10)

A la Ecuación 3.10 se le conoce como la APROXIMACIÓN PROGRESIVA; esta expresión permite evaluar la primera derivada en el bloque "i" en términos de las presiones en los bloques "i" e "i+1".

Simulación de Yacimientos I

9

A la parte truncada, R if , se le conoce como ERROR DE TRUNCAMIENTO de la aproximación progresiva.

En este caso, se dice que el ERROR DE

TRUNCAMIENTO es de PRIMER ORDEN debido a que la menor potencia a la cual se encuentra elevada la diferencia finita ∆ x es uno (Ecuación 3.9).

En

general, se dice que la aproximación numérica de una ecuación diferencial tiene

[

un error de truncamiento de orden "n" (se denota como O (∆ x )

n

] ) cuando la

mínima potencia de la diferencia finita, ∆ x , en la parte truncada es "n". Obsérvese que a mayor orden del error de truncamiento, mejor es la aproximación numérica de la aproximación diferencial.

Por ejemplo, considérense dos

esquemas de aproximación cuyos errores de truncamiento son, respectivamente,  ∂ 2 P  (∆ x ) ∂ 3 P  (∆ x )2 ∂ 4 P  (∆ x )3     + ⋅ + ⋅ + L R i 1 = −  2  ⋅ ∂ x 3  i 3 ! ∂ x 4  i 4 !  ∂ x  i 2 ! 

]

y

Ri

]

2

 ∂ 3 P  (∆ x ) 2 ∂ 4 P  (∆ x )3    ⋅ + + = −  3  ⋅ L ∂ x 4  i 4 !  ∂ x  i 3 ! 

]

La aproximación correspondiente al error de truncamiento R i 2 , es de mayor exactitud debido a que se está truncando únicamente desde el término de (∆ x ) , 2

en tanto que la aproximación correspondiente al error de truncamiento R i truncando desde el término de (∆ x ) .

[

O (∆ x )

2

] > O [ (∆ x ) ]

Simulación de Yacimientos I

En general, se tiene que:

]

1

se está

O [ (∆ x ) ] >

3

10

3.2.2 La aproximación regresiva y su error de truncamiento x

a – x = - ( x –a) x

X=0

X i-1

a

Xi

X=L

Fig. 3.6 – Bloques tenidos en cuenta en la aproximación regresiva de la primera derivada Considérese el sistema lineal de bloques de la figura 3.6. Si se define a = x i y x = xi − 1 Entonces, x − a = x i − 1 − x i = − (x i − x i − 1 ) = − ∆ x Si se hace f (x ) = P , entonces, de la Ecuación 3.3, se tiene: 2 3 ∂ 3 P  (− ∆ x ) ∂P  ∂ 2 P  (− ∆ x )  ⋅  ⋅  ⋅ (− ∆ x ) + P (x i − 1 ) = P (x i ) + +L + 2! 3! ∂ x 3  Xi ∂ x  Xi ∂ x 2  Xi

∂ n P  (− ∆ x )  ⋅ +L n! ∂ x n  Xi n

+

donde: P (x i − 1 ) : Presión en el bloque i − 1. P (xi ) : Presión en el bloque i .

Simulación de Yacimientos I

11

Para simplificar la notación, se suele escribir Pi − 1 en lugar de P (x i − 1 ) ; luego, de la ecuación anterior, se tiene: Pi−1

2 3 n ∂P  (∆ x ) ∂ 2 P  (∆ x ) ∂ 3 P  (∆ x ) ∂ n P  (− ∆ x )     ⋅ = Pi − + ⋅ − ⋅ +L+ ⋅ + L (3.12) ∂ x  i 1! n! ∂ x 2  i 2 ! ∂ x 3  i 3 ! ∂ x n  i

despejando de esta ecuación

∂P   : ∂ x  i

2 n −1 ∂ n P  (− ∆ x ) ∂P  Pi − Pi − 1 ∂ 2 P  (∆ x ) ∂ 3 P  (∆ x )     = ⋅ +L ⋅ +L+ ⋅ − + n! ∂ x  i ∆x ∂ x n  i ∂ x 3  i 3 ! ∂ x 2  i 2!

(3.13)

o bien: ∂P  Pi − P i − 1  = + R bi ∂x i ∆x

(3.14)

donde: R bi =

n−1 2 ∂ n P  (− ∆ x ) ∂ 2 P  (∆ x ) ∂ 3 P  (∆ x )    L ⋅ +L ⋅ + + ⋅ − n! ∂ x n  i ∂ x 3  i 3 ! ∂ x 2  i 2 !

(3.15)

R bi tiende a cero cuando ∆ x tiende a cero. En este caso, la Ecuación 3.14 toma

la forma: ∂P  Pi − P i − 1  ≅ ∂ x  i ∆x

(3.16)

A la Ecuación 3.16 se le conoce como la APROXIMACIÓN REGRESIVA para la primera derivada; su error de truncamiento está dado por la Ecuación 3.15. Obsérvese que el error de truncamiento de esta aproximación es de primer orden, O [ (∆ x ) ] .

Simulación de Yacimientos I

12

3.2.3 La aproximación central y su error de truncamiento

Si se sustrae la Ecuación 3.12 de la Ecuación 3.6, se obtiene: 2 3   ∂ 3 P  (∆ x ) ∂P  (∆ x ) ∂ 2 P  (∆ x )    ⋅ Pi + 1 − P i − 1 = Pi + ⋅ + L ⋅ + + 2  3  ∂x i 3 ! ∂ x  i 1! ∂x i 2 !   2 3   ∂ 3 P  (∆ x ) ∂P  (∆ x ) ∂ 2 P  (∆ x )    ⋅  - Pi − L ⋅ + ⋅ − + ∂ x 3  i 3 ! ∂ x  i 1! ∂ x 2  i 2 !  

Agrupando términos semejantes, ∂P  (∆ x ) ∂ 3 P  (∆ x ) ∂ n P  (∆ x ) ∂ n P  (− ∆ x )   ⋅  ⋅ Pi + 1 − P i − 1 = 2 ⋅ L + 2 ⋅ 3  ⋅ + ⋅ + n! ∂ x  i 1! ∂x i 3 ! ∂ x n  i n ! ∂ x n  i

Despejando

n

n

3

∂P   , ∂ x  i

n−1 n−1 2 ∂ n P  (− ∆ x ) ∂ n P  (∆ x ) ∂P  Pi + 1 − Pi − 1 ∂ 3 P  (∆ x )     = L ⋅ ⋅ − ⋅ − − − 2 ⋅ (∆ x ) n! n! ∂ x n  i ∂ x n  i ∂x i ∂ x 3  i 3 !

(3.17)

o bien: ∂P  Pi + 1 − Pi − 1  = + R iC ∂x i 2 ⋅ (∆ x )

(3.18)

donde: ∂ n P  (∆ x ) ∂ 3 P  (∆ x )  ⋅  L ⋅ + − n! ∂ x n  i ∂ x 3  i 3 ! 2

R iC = −

n−1

−

∂ n P  (− ∆ x )  ⋅ n! ∂ x n  i

n−1

−L

(3.19)

R iC tiende a cero cuando ∆ x tiende a cero. En este caso, la Ecuación 3.18 se

puede aproximar a:

Simulación de Yacimientos I

13

∂P  Pi + 1 − Pi − 1  ≅ ∂ x  i 2 ⋅ (∆ x )

(3.20)

A la Ecuación 3.20 se le conoce como la APROXIMACIÓN CENTRAL de la primera derivada; su error de truncamiento está dado por la Ecuación 3.19, de

[

donde se observa que éste es de SEGUNDO orden, O (∆ x )

].

2

Por lo tanto, el

error de truncamiento de la aproximación central es menor que el error de truncamiento de las aproximaciones progresiva y regresiva. Las Figuras 3.7, 3.8 y 3.9 presentan la interpretación geométrica de las aproximaciones progresiva, regresiva y central, respectivamente.

De estas

figuras, se observa que la aproximación progresiva estima la primera derivada de la presión en el bloque " i " con base en los valores de presión en este bloque y en el bloque posterior, la aproximación regresiva con base en los valores de presión en este bloque y en el bloque anterior, y la aproximación central con base en los valores de presión en los bloques anterior y posterior. Ésta es justamente la razón por la cual el error de truncamiento de la aproximación central es menor que los errores de truncamiento de las aproximaciones progresiva y regresiva.

Pi + 1 Pi + 1 – Pi Pi

∆X

:

 ∂p   ∂x    progresiva

:

 ∂p   ∂x    verdadera

Xi + 1 Fig. 3.7 – Interpretación geométrica de la aproximación progresiva

Simulación de Yacimientos I

14

Pi Pi – Pi – 1 Pi – 1

∆X

Xi – 1

:

 ∂p   ∂x    regresiva

:

 ∂p   ∂x    verdadera

Xi

Fig. 3.8 – Interpretación geométrica de la aproximación regresiva

Pi + 1 Pi

:

 ∂p   ∂x    centrada

:

 ∂p   ∂x    verdadera

Pi – 1

Xi – 1

Xi

Xi + 1

Fig. 3.9 – Interpretación geométrica de la aproximación central

Simulación de Yacimientos I

15

3.3 APROXIMACIÓN NUMÉRICA PARA LA SEGUNDA DERIVADA

La aproximación numérica para la segunda derivada puede ser obtenida sumando las Ecuaciones 3.6 y 3.12 definidas anteriormente para las aproximaciones progresiva y regresiva, de esta operación se obtiene: ∂ 2P  ∂ 4 P  (∆ x ) ∂ 6 P  (∆ x ) 2    ⋅ Pi + 1 + P i − 1 = 2 ⋅ Pi + x 2 2 ⋅ ( ∆ ) + ⋅ ⋅ + ⋅ +L ∂ x 2  i ∂ x 4  i 4 ! ∂ x 6  i 6 ! 4

6

∂ 2P  , Despejando ∂ x 2  i 4 2 ∂ 6 P  (∆ x ) ∂ 4 P  (∆ x ) ∂ 2 P  Pi − 1 − 2 ⋅ Pi + P i + 1    2 2 ⋅ +L ⋅ − ⋅ = − ⋅ ∂ x 6  i 6 ! ∂ x 4  i 4 ! ∂ x 2  i (∆ x )2

(3.21)

o bien, ∂ 2 P  Pi − 1 − 2 ⋅ Pi + P i + 1  = + R i2 2 2  ∂x i (∆ x )

(3.22)

donde: 4  ∂ 4 P  (∆ x ) 2  ∂ 6 P  (∆ x ) R = − 2 ⋅ 4  ⋅ + 2 ⋅ 6  ⋅ + L ∂x i 6!  ∂ x  i 4 !  2 i

(3.23)

R i2 tiende a cero cuando ∆ x tiende a cero. En este caso, de la Ecuación 3.22 se

tiene que ∂ 2 P  Pi − 1 − 2 ⋅ Pi + P i + 1  ≅ ∂ x 2  i (∆ x )2

Simulación de Yacimientos I

(3.24)

16

A la Ecuación 3.24 se le conoce como la APROXIMACIÓN NUMÉRICA PARA LA SEGUNDA DERIVADA. Esta expresión permite estimar el valor de

∂ 2P en el ∂ x2

bloque "i" en término de las presiones en los bloques "i-1", "i" e "i+1". El error de truncamiento de esta ecuación está dado por la Ecuación 3.23 de donde se puede observar que es de SEGUNDO ORDEN.

3.4 ESQUEMAS DE APROXIMACIÓN EXPLICITO E IMPLICITO

Considérese que la Ecuación 3.1 es válida en el intervalo 0 < x < X y t > 0 . Para dar solución numérica a esta ecuación, se requiere dividir el intervalo [0, X] en diferencias finitas de espacio, ∆ x , el intervalo [0, t ] en diferencias finitas de ∂ 2P ∂P tiempo, ∆ t , y expandir las derivadas y , respectivamente, tal como se 2 ∂x ∂t

ilustra en la Figura 3.10. t

∆t

i–1

i

i+1

X ∆X Fig. 3.10 – Discretización de una sistema lineal en función del espacio y el tiempo 0

Simulación de Yacimientos I

17

APROXIMACIÓN DE LA DERIVADA

La derivada

∂P . ∂t

∂P puede ser expandida mediante una aproximación progresiva ∂t

(Ecuación 3.10), regresiva (Ecuación 3.16) ó central (Ecuación 3.20).

En las

Ecuaciones 3.10, 3.16 y 3.20 se aproximó la variación de la presión con respecto a la posición a un tiempo fijo,

∂P   . En este caso se requiere evaluar la variación ∂ x  t

de la presión con respecto al tiempo en una posición fija,

∂P   ; por lo tanto, en ∂ t  x

lugar de emplear la diferencia finita ∆ x , tal como aparece en las Ecuaciones 3.10, 3.16 y 3.20, se debe emplear la diferencia finita ∆ t . Así mismo, el subíndice "i" (el cual indica posición) en este caso permanece fijo puesto que se desea aproximar la derivada de la presión con respecto al tiempo y, en consecuencia, el índice que varía debe ser aquel que indique el nivel de tiempo. Para este propósito, se denotará por la letra "n" los índices que indiquen nivel de tiempo y serán superíndices de la variable P . Los índices que indiquen posición se denotarán por la letra "i", y serán subíndices de la variable P , de tal manera que Pi indica la presión en el bloque i . P n indica la presión al tiempo n . Pin indica la presión en el bloque i , al tiempo n .

Teniendo

en

cuenta

estos

comentarios,

las

APROXIMACIONES progresiva, regresiva y central de

Simulación de Yacimientos I

ecuaciones

para

las

∂P son: ∂t

18

Aproximación progresiva: n+1 − Pin ∂P  Pi  = + O[ (∆ t ) ] ∆t ∂t 

(3.25)

n n−1 ∂P  Pi − Pi = + O[ (∆ t ) ] ∆t ∂ t 

(3.26)

Aproximación regresiva:

Aproximación central:

[

n+1 − Pin − 1 ∂P Pi 2 = + O (∆ t ) 2 ⋅ ∆t ∂t

]

(3.27)

donde los superíndices "n-1", "n" y "n+1" representan los tiempos t n − 1 , t n

y

t n + 1 , respectivamente. No se hará uso de la aproximación regresiva (Ecuación 3.26), debido a que el principal interés es calcular la presión en cada bloque "i" del sistema a tiempos futuros (es decir, Pin + 1 ), partiendo de valores actuales (es decir, Pin ) y, contrario a esto, la aproximación regresiva relaciona la presión actual ( Pin ) con presiones a tiempos pasados ( Pin − 1 ). La aproximación central tampoco será utilizada debido a problemas de estabilidad, concepto que será discutido más adelante. Por estas razones, se empleará la aproximación progresiva, (Ecuación 3.25), para aproximar el término

∂P . ∂t

Simulación de Yacimientos I

19

APROXIMACIÓN DE LA SEGUNDA DERIVADA

∂ 2P . ∂ x2

De la aproximación (3.24) se tiene que la segunda derivada se puede expresar como: ∂ 2 P  Pi − 1 − 2 ⋅ Pi + P i + 1  ≅ ∂ x 2  i (∆ x )2

(3.24)

Esta ecuación, tal como está escrita, no incluye el nivel de tiempo. La pregunta que surge sería: ¿A qué tiempo hace referencia las presiones en esta ecuación? Dependiendo del nivel de tiempo asignado a cada término de presión en esta ecuación, se suele hablar de diferentes esquemas de aproximación, de los cuales los más comunes son: la aproximación explícita, la aproximación implícita y la aproximación de Crank-Nicholson. Desde el punto de vista numérico, cada uno de estos esquemas posee características diferentes.

3.4.1 Esquema de Aproximación Explícita

Supóngase que las presiones incluidas en la Ecuación 3.24 son evaluadas al tiempo t n . En este caso, se tiene: n n n ∂ 2 P Pi − 1 − 2 ⋅ Pi + P i + 1 ≅ ∂x 2 (∆ x )2

(3.28)

Si se lleva las Ecuaciones 3.25 y 3.28 a la Ecuación 3.1, se obtiene: P in+ 1 − 2 ⋅ Pin + Pin− 1

(∆ x )2

=

Pin + 1 − Pin ∆t

(3.28.a)

De donde,

Simulación de Yacimientos I

20

P in+ 1 ⋅ ∆ t − 2 ⋅ Pin ⋅ ∆ t + Pin− 1 ⋅ ∆ t = Pin + 1 ⋅ (∆ x ) − Pin ⋅ (∆ x ) 2

2

o bien: Pin + 1 = P in+ 1 ⋅

∆t

(∆ x )2

  2 ⋅ ∆t ∆t  + Pin− 1 ⋅ 1 − Pin ⋅  − 2  (∆ x )2   (∆ x )

si se define: λ=

∆t

(∆ x )2

(3.29)

entonces: Pin + 1 = λ ⋅ P in+ 1 − (2 ⋅ λ − 1) ⋅ Pin + λ ⋅ Pin− 1

(3.30)

A la Ecuación 3.30 se le conoce como la APROXIMACIÓN EXPLÍCITA de la Ecuación 3.1.

La aplicación repetitiva de esta aproximación numérica permite

estimar la presión en cada bloque "i" del sistema al tiempo t n + 1 , con base en los valores de presión en los bloques " i − 1 ", "i " e " i + 1" al tiempo t n . Para llevar a cabo esta serie de cálculos, se requiere de una condición inicial para cada bloque " i " y dos condiciones de frontera.

Ejemplo:

Supóngase un sistema lineal de siete bloques similar al presentado en la Figura 3.3 . Sea t 1 el tiempo actual; obtener las expresiones de la aproximación explícita para el cálculo de las presiones en cada bloque a un tiempo futuro t 2 = t 1 + ∆ t y a un tiempo

t 3 = t 2 + ∆ t , asumiendo que en cada caso las condiciones de límite

son:

Simulación de Yacimientos I

21

a. P0 = constante = PCERO b. P8 = 0

Solución: Sea t 1 = t n ; luego, t 2 = t 1 + ∆ t = t n + 1 .

Puesto que t 1 es el tiempo actual,

entonces la presiones a este tiempo, P1 , P2 , ..., son conocidas (Condición inicial a t = t Inicial ). De la Ecuación 3.30, se puede obtener las expresiones para el cálculo de las presiones en cada bloque al tiempo t 2 = t n + 1 :

Para el bloque 1, i = 1: P1(2 ) = λ ⋅ P 2(1) − (2 ⋅ λ − 1) ⋅ P1(1) + λ ⋅ P0(1)

P0 = constante = PCERO (Condición de límite del ejemplo).

Para el bloque 2, i = 2 : P2(2 ) = λ ⋅ P 3(1) − (2 ⋅ λ − 1) ⋅ P2(1) + λ ⋅ P1(1)

Para el bloque 3, i = 3 : P3(2 ) = λ ⋅ P 4(1) − (2 ⋅ λ − 1) ⋅ P3(1) + λ ⋅ P2(1)

Simulación de Yacimientos I

22

M

Para el bloque 7, i = 7 : P7(2 ) = λ ⋅ P 8(1) − (2 ⋅ λ − 1) ⋅ P7(1) + λ ⋅ P6(1)

Pero P8 = 0 (Condición de límite en este ejemplo), luego: P7(2 ) = − (2 ⋅ λ − 1) ⋅ P7(1) + λ ⋅ P6(1)

Al tiempo t 3 = t 2 + ∆ t , se pueden re-definir la variables t n y t n + 1 , así: tn = t2

y

tn + 1 = t3

Teniendo en cuenta lo anterior, de la Ecuación 3.30, se puede escribir: Para el bloque 1, i = 1: P1(3 ) = λ ⋅ P 2(2 ) − (2 ⋅ λ − 1) ⋅ P1(2 ) + λ ⋅ P0(2 )

De nuevo, P0 = constante = PCERO .

Para el bloque 2, i = 2 : P2(3 ) = λ ⋅ P 3(2 ) − (2 ⋅ λ − 1) ⋅ P2(2 ) + λ ⋅ P1(2 )

Simulación de Yacimientos I

23

Para el bloque 3, i = 3 : P3(3 ) = λ ⋅ P 4(2 ) − (2 ⋅ λ − 1) ⋅ P3(2 ) + λ ⋅ P2(2 )

etc. Obsérvese que las presiones incluidas en los miembros derechos de las tres últimas ecuaciones ya han sido previamente calculadas en el paso anterior.

tn+1

Xi – 1

Xi

Xi + 1

tn

Esquema Explícito

tn+1

Xi – 1

X

Xi + 1

tn

Esquema iImplícito

Fig. 3.11 – Forma de los esquemas explícito e implícito esta grafica comparativa debiera aparecer después de explicarse el esquema implícito

En la Figura 3.11 se indica los puntos en el espacio y los niveles de tiempo tenidos en cuenta en la aproximación explícita (Ecuación 3.30). El procedimiento general para la aplicación de este esquema es el siguiente: a. A cada uno de los bloques del sistema se le asigna un valor de presión al tiempo inicial t = 0 . De esta forma se obtienen valores para las variables Pi + 1 , Pi y Pi − 1 de la Ecuación 3.30. b. Se calcula los valores de Pi para cada uno de los bloques de la malla.

Simulación de Yacimientos I

24

c. Una vez barrida toda la malla, se repite el procedimiento para el tiempo t n + 2 , Pi − 1 ,

haciendo que los nuevos valores de

Pi ,

Pi + 1 , sean iguales a los

valores Pi − 1 , Pi , Pi + 1 , apenas calculados. Tal como se discutirá más adelante, para obtener estabilidad en este esquema de solución, se requieren incrementos de tiempo, ∆ t , muy pequeños, lo que limita su aplicación en un gran número de problemas.

3.4.2 Esquema de Aproximación Implícita

Supóngase que las presiones Pi − 1 , Pi y Pi + 1 presentes en la Ecuación 3.24 hacen referencia al tiempo t n + 1 . En este caso, la Ecuación 3.24 puede ser escrita como: n+1 n+1 + P in++11 ∂ 2 P Pi − 1 − 2 ⋅ Pi ≅ ∂x 2 (∆ x )2

(3.31)

Si se llevan las Ecuaciones 3.25 y 3.31 a la Ecuación 3.1, se obtiene: P in++11 − 2 ⋅ Pin + 1 + Pin−+11

(∆ x )2

=

Pin + 1 − Pin

∆t

(3.32)

de donde: P in++11 ⋅

∆t

(∆ x )

2

− 2 ⋅ Pin + 1 ⋅

∆t

(∆ x )

2

+ Pin−+11 ⋅

∆t

(∆ x )

2

− Pin + 1 = − Pin

Teniendo en cuenta que: λ=

∆t

(∆ x )2

Simulación de Yacimientos I

25

entonces, λ ⋅ P in++11 − (2 ⋅ λ + 1) ⋅ Pin + 1 + λ ⋅ Pin−+11 = − Pin

(3.33)

La Ecuación 3.33 representa la APROXIMACIÓN IMPLICITA de la Ecuación 3.1. Si se aplica esta ecuación a cada uno de los bloques de la malla, se genera un sistema de ecuaciones cuya solución simultánea permite obtener las distribución de presiones en el sistema (presión en cada bloque) al tiempo t n + 1 . Ejemplo:

Obtener el sistema de ecuaciones generado por el esquema implícito, Ecuación 3.33, para el ejemplo anterior. Solución: Sea t 1 = t n y t 2 = t n + 1 . La Ecuación 3.33 genera una ecuación para cada bloque del sistema, así: Para el bloque 1, i = 1: λ ⋅ P 0(2 ) − (2 ⋅ λ + 1) ⋅ P1(2 ) + λ ⋅ P2(2 ) = − P1(1)

Para el bloque 2, i = 2 : λ ⋅ P1(2 ) − (2 ⋅ λ + 1) ⋅ P2(2 ) + λ ⋅ P3(2 ) = − P2(1)

Para el bloque 3, i = 3 : λ ⋅ P 2(2 ) − (2 ⋅ λ + 1) ⋅ P3(2 ) + λ ⋅ P4(2 ) = − P3(1)

Simulación de Yacimientos I

26

etc. Este sistema de ecuaciones puede ser escrito en forma más compacta de la siguiente forma:

− (2 ⋅ λ + 1) ⋅ P1(2 ) + λ ⋅ P2(2 ) λ ⋅ P1(2 ) − (2 ⋅ λ + 1) ⋅ P2(2 ) + λ ⋅ P3(2 ) λ ⋅ P2(2 ) − (2 ⋅ λ + 1) ⋅ P3(2 ) O O λ ⋅ P6(2 )

+

−

λ ⋅ P4(2 )

(2 ⋅ λ + 1) ⋅ P7(2 )

= − P1(1) = − P2(1) = − P 3(1) M = − P 7(1)

Las incógnitas en este sistema de ecuaciones son: P1 , P2 , ..., P7 . Su solución permite obtener los valores de presión en cada bloque al tiempo t n + 1 = t 2 . La Figura 3.11 incluye los puntos en el espacio y los niveles de tiempo tenidos en cuenta en el esquema implícito.

3.5 CRITERIOS TENIDOS EN CUENTA PARA LA SELECCION DE UN ESQUEMA DE APROXIMACIÓN NUMÉRICA.

Existen diferentes esquemas de aproximación numérica a las ecuaciones diferenciales de flujo, dos de los cuales han sido discutidos en el numeral 3.4 (la aproximación explícita y la aproximación implícita). La selección de uno u otro de estos esquemas se debe realizar teniendo en cuenta los siguientes criterios: - Error de truncamiento. - Estabilidad. - Convergencia. - Consistencia.

Simulación de Yacimientos I

27

3.5.1 Error de Truncamiento

Tal como se mencionó anteriormente, el error de truncamiento hace referencia al error en que se incurre al utilizar la aproximación en diferencias finitas luego de ser truncada, en lugar de la ecuación diferencial misma. El error de truncamiento se define como:  Ecuación   Ecuación en εT =  −  Diferencial   Diferencias Finitas

  

(3.34)

Es importante resaltar que este error de truncamiento es diferente al error de truncamiento que se comete debido a las operaciones repetitivas del computador. Los computadores tienen una memoria de capacidad finita y, en consecuencia, pueden retener únicamente números de "tamaño finito" (es decir, hasta cierto número de cifras significativas). El error cometido por esta clase de truncamiento, el cual se propaga a medida que se incrementa el número de operaciones, es diferente al error de truncamiento numérico a que hace referencia la Ecuación 3.34. Ejemplo: ERROR DE TRUNCAMIENTO DE LA APROXIMACIÓN EXPLÍCITA E

IMPLÍCITA. Considérese la Ecuación 3.1, cuya aproximación explícita está dada por la Ecuación 3.28.a: P in+ 1 − 2 ⋅ Pin + Pin− 1

(∆ x )2

=

Pin + 1 − Pin ∆t

(3.28.a)

Aplicando la Ecuación 3.34, se tiene:  ∂ 2 p ∂p   P in+ 1 − 2 ⋅ Pin + Pin− 1 Pin + 1 − Pin εT =  − − − 2 ∂ t   ∆t (∆ x )2  ∂x

Simulación de Yacimientos I

  

28

Llevando las Ecuaciones 3.22 y 3.25, en reemplazo de

∂p ∂ 2p , se tiene: y 2 ∂t ∂x

2 4  P in+ 1 − 2 ⋅ Pin + Pin− 1 Pin + 1 − Pin ∂ 6 P  (∆ x ) ∂ 4 P  (∆ x ) +L− − 2 ⋅ 6  ⋅ εT =  − 2 ⋅ 4  ⋅ ∂x  i 6 ! ∆t ∂x i 4 ! (∆ x )2    P n − 2 ⋅ Pin + Pin− 1 Pin + 1 − Pin  ∂ 2P  ∆t − L −  i + 1 − + 2  ⋅  ∂t  i 2 ∆t (∆ x )2   

de donde, εT = −

(∆ x )2 ⋅ ∂ 4 P  12

∆t ∂ 2P   +L=O − ⋅ ∂ x 4  i 2 ∂ t 2  i

[ ( ∆x )

2

+ ( ∆t

)]

(3.35)

Si se procede en forma análoga, se encuentra que el error de truncamiento de la aproximación implícita está dado por una expresión similar a la Ecuación 3.35.

Ejemplo:

LA APROXIMACIÓN DE CRANK NICHOLSON Y SU ERROR DE TRUNCAMIENTO.

La aproximación de Crank-Nicholson puede ser interpretada como un promedio de la aproximación implícita y la aproximación explícita. Considérese la Ecuación 3.1 1  El esquema de Crank-Nicholson aproxima esta ecuación en el punto  i, n +  , 2  tal como se ilustra en la Figura 3.12 . Es decir,

Simulación de Yacimientos I

29

t n+1 ( Xi , t n + ½ ) tn Xi – 1 Xi Fig. 3.12 – Esquema del método de Crack-Nicholson n+1

1 ∂ 2P   ⋅ 2 ∂ x 2  i

n

Xi + 1

Central

1 ∂ 2P  ∂P   + ⋅ 2  = 2 ∂x i ∂ t  i

(3.36)

de donde: n+1 n+1 + Pin−+11  1  P in+ 1 − 2 ⋅ Pin + Pin− 1  Pin + 1 − Pin 1  P i + 1 − 2 ⋅ Pi ⋅  + ⋅ = 2   ∆t  (∆ x )2 (∆ x )2  2   2⋅   2

o bien, n+1 n+1 + Pin−+11  1  P in+ 1 − 2 ⋅ Pin + Pin− 1  Pin + 1 − Pin 1  P i + 1 − 2 ⋅ Pi ⋅  + ⋅ = ∆t 2  (∆ x )2 (∆ x )2  2  

(3.37)

Obsérvese que la expansión numérica de la derivada de la presión con respecto al tiempo, proviene de una aproximación central y, en consecuencia, puede ser escrita como: 2

n+1 i

∂P P = ∂t

n−1 i

−P  ∆t  2⋅   2

1  ∆t   2  ∂ 3P  n + 2  ⋅  − 6 ∂ t 3  i

(3.38)

Llevando las Ecuaciones 3.1 y 3.37 a la definición de error de truncamiento, Ecuación 3.34, se tiene:

Simulación de Yacimientos I

30

 ∂ 2 p ∂p   1 P in++11 − 2 ⋅ Pin + 1 + Pin−+11 1 P in+ 1 − 2 ⋅ Pin + Pin− 1 Pin + 1 − Pin − + ⋅ εT =  − − ⋅ 2 2 ∆t ∂ t   2 (∆ x )2 (∆ x )2  ∂x

  

(3.39)

Reemplazando las Ecuaciones 3.22 y 3.38 en la Ecuación 3.39, se tiene: n n n  1 P in++11 − 2 ⋅ Pin + 1 + Pin−+11 1 (∆ x ) 2 ∂ 4 P  1 P i + 1 − 2 ⋅ Pi + Pi − 1  −L + ⋅ − ⋅ ⋅ εT =  ⋅ 2 12 ∂ x 4  i 2 (∆ x )2 (∆ x )2  2 2  1  ∆t  n+ n + 1 n 2   4 3  2 P P − 2  ∂ P 1 (∆ x ) ∂ P  i i ⋅ 3  − L + − ⋅ ⋅ 4  − L − 2 12 ∂ x  i 6 ∂t  i  ∆t   2⋅   2  

 1 P in++11 − 2 ⋅ Pin + 1 + Pin−+11 1 P in+ 1 − 2 ⋅ Pin + Pin− 1 Pin + 1 − Pin − − ⋅ + ⋅ ∆t 2 (∆ x )2 (∆ x )2  2

  

de donde: 2

εT = −

(∆ x ) 24

n+1

2

⋅

∂ 4P   ∂ x 4  i

1  ∆t  n+   3 4 2 (∆ x ) ⋅ ∂ P  +  2  ⋅ ∂ P  + L − 24 ∂ x 4  i 6 ∂ t 3  i 2

n

o bien: εT = O

[ ( ∆x )

2

+ ( ∆t

)2 ]

(3.40)

Al comparar las Ecuaciones 3.35 y 3.40, se infiere que el esquema de CrankNicholson es más aproximado que los esquemas implícito y explícito.

Simulación de Yacimientos I

31

3.5.2 Estabilidad

El análisis de estabilidad es un procedimiento a través del cual es posible determinar si el error obtenido en un punto (bloque) de la malla, en un momento determinado, disminuye o aumenta al incrementar el tiempo. Si el error disminuye, se dice que la aproximación es ESTABLE; si aumenta, se dice que es INESTABLE. Si el error disminuye únicamente bajo ciertas circunstancias, se dice que es CONDICIONALMENTE ESTABLE. Para que las presiones obtenidas de la aplicación de un esquema numérico tengan aplicación real, se requiere que el esquema sea estable. La Figura 3.13 ilustra la variación de presión en un punto determinado del sistema en función del tiempo para un esquema estable y un esquema inestable. Esquema Estable

Esquema Inestable

∆Pn

∆P = P n + 1 − P n

Nivel de tiempo, tn

Fig. 3.13 – Cambio de presión en un punto del sistema en función del tiempo

A continuación se presentan dos criterios muy útiles para determinar estabilidad: El criterio de Karplus y el Método de Von Neumann.

Simulación de Yacimientos I

32

3.5.2.1 Criterio de Karplus

El método de Karplus no considera el efecto de condiciones de límite en el análisis de estabilidad. El procedimiento del método puede ser resumido de la siguiente forma: a. La ecuación en diferencias finitas es reordenada de manera tal que tome la siguiente forma:

(

)

(

)

(

)

(

)

a ⋅ Pin+ 1 − Pin + b ⋅ Pin− 1 − Pin + c ⋅ Pin + 1 − Pin + L + d ⋅ Pin − 1 − Pin + L = 0

(3.41)

b. Si todos los coeficientes a, b, c,..., d, ... son negativos, la aproximación es estable. c. Si algunos de los coeficientes son negativos, entonces para que la aproximación sea estable la suma de los coeficientes debe ser menor o igual a cero. Debido a que los coeficientes a, b, c, ..., d, ... pueden cambiar de un nodo a otro, el criterio debe ser aplicado a cada nodo. Un esquema determinado puede ser estable en un nodo pero inestable en otros.

Ejemplo: Determinar si la aproximación explícita de la Ecuación 3.1 es estable.

La aproximación explícita tiene la siguiente forma, Ecuación 3.30: Pin + 1 = λ ⋅ P in+ 1 − (2 ⋅ λ − 1) ⋅ Pin + λ ⋅ Pin− 1

Paso 1: La Ecuación 3.30, puede ser escrita como:

Simulación de Yacimientos I

33

(

)

(

) (

)

λ ⋅ P in+ 1 − Pin + λ ⋅ P in− 1 − Pin − Pin + 1 − Pin = 0

Paso 2: No todos los coeficientes son menores que cero.

Paso 3: La condición requerida para estabilidad es: λ + λ − 1= 2⋅λ − 1= 2⋅

( ∆t ) ( ∆ x )2

− 1≤ 0

de donde:

( ∆t ) ( ∆ x )2

≤

1 2

(3.42)

Luego, LA APROXIMACIÓN EXPLÍCITA ES CONDICIONALMENTE ESTABLE; la condición de estabilidad está dada por la Ecuación 3.42 . Siguiendo un procedimiento completamente análogo es posible determinar que LA APROXIMACIÓN IMPLICITA ES INCONDICIONALMENTE ESTABLE

3.5.2.2 Método de Von Neumman ó Análisis Armónico

La relación entre el valor exacto de presión en un punto x i , a un tiempo t n , p (x i, t n ) , y el valor calculado al aplicar un esquema numérico es:

Simulación de Yacimientos I

34

p (x i , t n ) = Pin + ε ni

(3.43)

Donde ε ni es la diferencia (error) entre el valor exacto p (x i, t n ) y el valor aproximado Pin . Si el error disminuye al incrementar el tiempo, los valores de presión obtenidos de la solución exacta satisfacen el esquema numérico; por ejemplo, para el esquema implícito se cumple: p (x i + 1, t n + 1 ) − 2 ⋅ p (x i, t n + 1 ) + p (x i − 1, t n + 1 )

(∆ x )2

p (x i, t n + 1 ) − p (x i , t n )

=

(3.44)

∆t

Llevando la Ecuación 3.43 a la Ecuación 3.44, se tiene:

(P

n+1 i+1

)

(

) (

+ ε ni ++11 − 2 ⋅ Pin + 1 + ε ni + 1 + Pin−+11 + ε ni −+11

(∆ x )

2

) = (P

n+1 i

) (

+ ε ni + 1 − Pin + ε ni ∆t

) (3.45)

de la misma manera para el esquema explícito se obtiene,

(P

n i+1

)

(

) (

+ ε ni + 1 − 2 ⋅ Pin + ε ni + Pin− 1 + ε ni − 1

(∆ x )

2

) = (P

n+1 i

) (

+ ε ni + 1 − Pin + ε ni

)

∆t

restando la Ecuación 3.31 ojo, es la 3.32 de la Ecuación 3.45 se llega a: ε ni ++11 − 2 ⋅ ε ni + 1 + ε ni −+11

(∆ x )2

=

ε ni + 1 − ε ni ∆t

(3.46)

Es decir, el error también satisface el esquema numérico de aproximación.

Una función f ( x ) puede ser expandida en series de Fourier de la siguiente forma:

Simulación de Yacimientos I

35

f (x ) = ∑ γ k ⋅ e

−1⋅β ⋅ x

(3.47)

k

donde β es una constante. En el punto x = x i = i ⋅ ∆ x , se tiene:ojo, deberías decir qué es lambda f (x ) = ∑ γ k ⋅ e

−1⋅β ⋅i⋅ ∆x

(3.48)

k

Si la función expandida es la función error, f (x ) = ε ni , entonces, de la ecuación (3.48), se tiene: ε ni = ∑ γ nk ⋅ e

−1⋅β ⋅i⋅ ∆x

(3.49)

k

El n-‚ésimo término de esta sumatoria, al que le haremos corresponder el tiempo t n , tiene la siguiente forma:

ε ni = γ nk ⋅ e

−1⋅β ⋅i⋅ ∆x

(3.50)

El término n + 1 de la sumatoria, al que haremos corresponder el tiempo t n + 1 , tiene la siguiente forma:

ε ni + 1 = γ nk + 1 ⋅ e

−1⋅β ⋅i⋅ ∆x

(3.51)

Si el error en el punto " i ", ε i , disminuye al pasar de t n a t n + 1 , entonces ε ni + 1 < ε ni , de donde se tiene: ε ni + 1 ε ni

≤1

(3.52)

Llevando las Ecuaciones 3.50 y 3.51 a la Ecuación 3.52, se tiene:

Simulación de Yacimientos I

36

γn + 1 ≤1 γn

La relación γ n + 1 γ n es conocida como el factor de amplificación, γn + 1 ν= γn

(3.53)

:

(3.54)

Para que un determinado esquema sea estable, se hace necesario que ν sea menor o igual a la unidad.

Ejemplo: Aplicando el método de Von Neumman, determinar la estabilidad del

esquema explícito. Solución: Paso 1: Se expresa la ecuación que representa la aproximación numérica bajo análisis en términos del error ε i en lugar de la presión P . Particularmente, para el esquema explícito, de la Ecuación 3.46 se tiene: ε ni + 1 − 2 ⋅ ε ni + ε ni − 1

(∆ x )2

=

ε ni + 1 − ε ni ∆t

Paso 2: Cada término de la ecuación anterior se sustituye por su equivalente de la Ecuación 3.49 . En este caso, se tiene:

Simulación de Yacimientos I

37

γn ⋅ e

− 1 ⋅ β ⋅ (i + 1) ⋅ ∆ x

− 2⋅ γn ⋅ e

−1⋅β ⋅i⋅ ∆x

(∆ x )2

+ γn ⋅ e

− 1 ⋅ β ⋅ (i − 1) ⋅ ∆ x

=

γn + 1 ⋅ e

−1⋅β ⋅i ⋅ ∆x

− γn ⋅ e ∆t

−1⋅β ⋅i ⋅ ∆x

O bien: ∆t

(∆ x )

2

[

− 1⋅ β ⋅ i ⋅ ∆x

⋅ γ n ⋅e

⋅e

− 1⋅ β ⋅ ∆x

− 2⋅ γ n ⋅e

− 1⋅ β ⋅ i ⋅ ∆x

+ γ n ⋅e

− 1⋅ β ⋅ i ⋅ ∆x

⋅e −

γ n + 1 ⋅e

∆t

(∆ x )

2

[

−1⋅β ⋅ ∆x

⋅ γn ⋅e

− 2⋅ γn + γn ⋅e−

−1⋅β ⋅ ∆x

]= γ

n+1

− 1⋅ β ⋅ ∆x

]=

− 1⋅ β ⋅ i ⋅ ∆x

− γ n ⋅e

− 1⋅ β ⋅ i ⋅ ∆x

− γn

De donde,

[

γn + 1 ∆t = 1+ ⋅ e n γ (∆ x )2

−1⋅β ⋅ ∆x

− 2 + e−

−1⋅β ⋅ ∆x

]

De trigonometría: e e−

−1⋅ θ

= cos θ + − 1 ⋅ sen θ

−1⋅ θ

= cos θ − − 1 ⋅ sen θ

Luego, haciendo θ = β ⋅ ∆ x ,

[

γn + 1 ∆t = 1+ ⋅ cos θ + n γ (∆ x )2

=1+

∆t

(∆ x ) 2

= 1 + 2⋅

− 1 ⋅ sen θ − 2 + cos θ −

− 1 ⋅ sen θ

]

⋅ [ 2 ⋅ cosθ − 2 ]

∆t

(∆ x ) 2

⋅ [ cosθ − 1 ]

Simulación de Yacimientos I

38

De acuerdo al criterio de Von Neumman, para que el esquema sea estable, se debe cumplir que: γn + 1 ≤1 ν = γn

es decir, 1− 2⋅

∆t

(∆ x )2

⋅ [ 1 − cos θ ] ≤ 1

Lo que se cumple solamente si: −1≤ 1− 2⋅

∆t

(∆ x )2

⋅ [ 1 − cos θ ] ≤ 1

La solución a esta la desigualdad estará dada por aquellos valores de ∆ t (∆ x )

2

que cumplan ambos lados de la desigualdad. Debido a que − 1 ≤ cos θ ≤ 1 el lado derecho siempre se cumple. En cuanto al lado izquierdo, se tiene: −1≤ 1 − 2⋅

− 2 ≤ − 2⋅

1≥

∆t

(∆ x )2

∆t

(∆ x )

2

≤

∆t

(∆ x )2 ∆t

(∆ x )2

⋅ ( 1 − cos θ )

⋅ ( 1 − cos θ )

⋅ ( 1 − cos θ )

1 1 − cos θ

Simulación de Yacimientos I

39

El valor mínimo de la expresión 1 cuyo caso 1

[ 1 − cos θ ] = 12

[ 1 − cos θ ]

se obtiene para cos θ = − 1 , en

. Por tanto, el esquema será estable siempre y

cuando: ∆t

(∆ x )

2

≤

1 2

(3.55)

Este resultado corrobora el resultado del ejemplo anterior: LA APROXIMACIÓN EXPLÍCITA ES CONDICIONALMENTE ESTABLE; la condición está dada por la Ecuación 3.55. Si se procede en forma análoga para el análisis de estabilidad del esquema implícito, se obtendrá que‚ éste es incondicionalmente estable (se deja al lector como ejercicio). Debido a que el esquema implícito es incondicionalmente estable, éste suele ser preferido en Simulación de Yacimientos.

3.5.3 Convergencia

Supóngase que la función p (x ) es la solución exacta de la Ecuación 3.1 y que P in es la solución en diferencias finitas.

Se dice que la aproximación numérica

utilizada para el cálculo de P in es CONVERGENTE si P in tiende a p (x i, t n ) , cuando ∆ x y ∆ t tienden a cero. Obsérvese que la convergencia concierne al error Pin − p (x i, t n ) al tiempo t n , en un punto determinado de la malla, mientras que la estabilidad tiene que ver con la propagación (incremento) del error en un punto determinado al incrementar el tiempo. La estabilidad es un criterio de mayor peso que la convergencia. En general, estabilidad y consistencia implican convergencia.

Simulación de Yacimientos I

40

En Simulación de Yacimientos, no se suele conocer la solución exacta de la ecuación diferencial de flujo. Por esta razón, la forma de analizar convergencia de un esquema numérico, consiste en dar solución numérica a la ecuación para varios valores de ∆ x y ∆ t . Si las soluciones para valores sucesivamente menores de ∆ x y ∆ t tienden a un cierto valor constante, con cierta tolerancia, se puede concluir convergencia del esquema numérico.

3.5.4 Consistencia o Compactibilidad

Se dice que un esquema numérico es consistente, si la ecuación en diferencias finitas tiende a la ecuación diferencial cuando ∆ x y ∆ t tienden a cero. Dada la definición de error de truncamiento, Ecuación 3.34, un esquema es consistente si su error de truncamiento tiende a cero, cuando sus diferencias finitas, ∆ x y ∆ t , por ejemplo, tienden a cero.

Ejemplo: Determinar si la aproximación explícita es consistente.

El error de truncamiento para la aproximación explícita está dado por la Ecuación 3.35, εT

si

2 ( ∆x) ∂ 4P   =− ⋅

12

∆x

y

∆t ∂ 2P   +L − ⋅ ∂ x 4  i 2 ∂ t 2  i

∆t

tienden a cero, ε T tiende a cero; luego este esquema es

consistente.

Simulación de Yacimientos I

41

BIBLIOGRAFIA

1. Osorio, Gildardo.

“Conceptos básicos sobre aproximación de ecuaciones

diferenciales a diferencias finitas”.

Notas de simulación de yacimientos,

Maestría en Hidrocarburos. Febrero de 2002 2. Azis, K and Settari A.

“Petroleum Reservoir Simulation”.

Elsevier Applied

Science Publishers 3. Peaceman, Donald.

“Fundamentals of Numerical Reservoir Simulation”.

Elsevier Scientific Publishing Company. 4. Thomas, G. W. “Principles of Hydrocarbon Reservoir Simulation”. International Human Resources Development Corporation.

Simulación de Yacimientos I

42

4. FLUJO LINEAL DE UN FLUIDO INCOMPRESIBLE El capítulo 3 ha sido dedicado a la discusión de los fundamentos que permiten aproximar una ecuación diferencial parcial a diferencias finitas. A través de todo el capitulo, se trabajó con una ecuación diferencial de flujo sencilla con la finalidad de ilustrar la aplicación de los conceptos básicos. Este capítulo se dedica a la aplicación de los conceptos presentados en el capítulo 3 al flujo lineal de un fluido incompresible a través de un medio poroso. diferencia

del

capítulo

anterior,

las

ecuaciones

diferenciales

A

utilizadas

corresponden a ecuaciones de flujo que describen en forma más real los procesos que ocurren en el sistema.

4.1 PLANTEAMIENTO DEL MODELO NUMÉRICO BÁSICO La ecuación de continuidad para flujo monofásico en una dirección está dada por la siguiente ecuación:   k ∂ (ρφ ) ∇ ⋅  αρ ∇p  = α + α~ q µ ∂t   la cual , para flujo lineal toma la siguiente forma: k ∂p  ∂  ∂ (ρφ )  A (x ) ρ x  = A (x ) + A (x )~ q ∂x  µ ∂x  ∂t

(4.1)

La ecuación de estado para un fluido incompresible está dada por la Ecuación 2.8:

ρ = constante

(2.8)

Sustituyendo la Ecuación 2.8 en la Ecuación 4.1 y considerando que la porosidad del medio es constante con el tiempo, se obtiene: Simulación de Yacimientos I

1

~ k ∂p  ∂  q  A (x ) x  = A (x ) ∂x  µ ∂x  ρ

(4.2)

Si se define a q v como el volumen de fluido que entra o sale por fuentes o sumideros por unidad de volumen total de yacimiento por unidad de tiempo, entonces la relación entre ~ q y q estará dada por la siguiente ecuación: v

qv =

~ q ρ

(4.3)

Llevando la Ecuación 4.3 a la Ecuación 4.2, se obtiene: ∂  kA ∂p   ⋅  = A ⋅ qv ∂x  µ ∂x 

(4.4)

La Ecuación 4.4 es la ecuación fundamental de flujo que rige el flujo monofásico de un fluido incompresible a través de un medio poroso lineal. Para dar solución numérica a la Ecuación 4.4 es necesario dividir el sistema en n bloques, tal como se ilustra en la Figura 4.1. Para fines prácticos, es necesario considerar que los bloques son de diferente longitud y tanto el área como la permeabilidad y la viscosidad en un bloque pueden diferir de los demás. Debido a que el producto kA µ no es el mismo para todos los bloques, no es posible extraerlo del diferencial de la Ecuación 4.4, tal como se hace con una constante. Debido a esta no linealidad, no se puede aplicar la aproximación para la segunda derivada en la forma como se hace en la Ecuación 3.24.

1

i −1

i

i +1

nx

∆x1

∆x i −1

∆x i

∆x i +1

∆x nx

Fig. 4.1 – Discretización lineal en bloques de diferente longitud Simulación de Yacimientos I

2

Considere la siguiente definición para el bloque i :  kA ∂p  Ui =  ⋅   µ ∂x i

(4.5)

De esta forma, la Ecuación 4.4 puede ser escrita como: ∂Ui = A ⋅ qv ∂x

(4.6)

Considérese el i -ésimo bloque de un sistema lineal como el ilustrado en la Figura 4.2

Figura 2. Esquema del i-ésimo bloque de un sistema lineal.

Con ayuda de la serie de Taylor, la función U(x ) puede expandirse en torno a un punto a mediante la siguiente expresión: U(n) (a)(x − a )n U( x ) = ∑ n! n=0 ∞

(4.7)

Si el punto en torno al cual se expande la función U(x ) es a = xi y la función se evalúa en el punto x = xi − 1 = xi − 2

∆x i ∆x , entonces x − a = − i y de la Ecuación 2 2

4.7 se obtiene:

Simulación de Yacimientos I

3

2

3

U′  ∆x  U′′  ∆x  U′′′ ∆x  = Ui + i  − i  + i  − i  + i  − i  + L 1!  2  2!  2  3!  2 

Ui − 12

2

3

U′  ∆x  U′′  ∆x  U′′′ ∆x  = Ui − i  i  + i  i  − i  i  + L ......................................... (4.8) 1!  2  2!  2  3!  2 

Ui− 1

2

Si la función se evalúa en punto x = xi + 1 = xi + 2

∆x i y de la Ecuación 4.7 se obtiene: 2

entonces x − a =

2

Ui+ 1

2

∆x i , en lugar de x = xi − 1 , 2 2

3

U′  ∆ x  U′′  ∆ x  U′′′ ∆ x  = Ui + i  i  + i  i  + i  i  + L ..................................... (4.9) 1!  2  2 !  2  3!  2 

Restando la Ecuación 4.8 de la Ecuación 4.9, se obtiene: 3

Ui+ 1 − Ui− 1 2

2

U′i′′  ∆ x   ∆x  = 2U′i  i  + 2  i  + K 3!  2   2 

Resolviendo para U′i se tiene:

U′i =

Ui+ 1 − Ui− 1 2

2

∆xi

−

U′i′′ (∆xi )2 3!

4

−K =

Ui+ 1 − Ui− 1 2

∆xi

2

+ Ric

Por consiguiente:

U′i ≅

Ui+ 1 − Ui− 1 2

∆xi

2

..................................................................................... (4.10)

El error de truncamiento de la Ecuación 4.10 es de segundo orden: Ric

U′i′′ (∆xi )2 =− +K 3! 4

Sustituyendo la Ecuación 4.5 en la Ecuación 4.10, se obtiene: Simulación de Yacimientos I

4

∂  kA ∂p   ⋅  ≅ ∂x  µ ∂x i

 k ⋅ A ∂p   k ⋅ A ∂p   ⋅  −  ⋅   µ ∂x i+ 1  µ ∂x i− 1 2

2

∆x i

Por conveniencia, la ecuación anterior suele ser escrita de la siguiente forma:

∂  kA ∂p   ⋅  ≅ ∂x  µ ∂x i

 kA   kA   ∂p   ∂p     −       µ i+ 1  ∂x i+ 12  µ i− 1  ∂x i− 12 2

2

∆x i

.............................. (4.11)

 ∂p  4.1.1 Evaluación del término   . Con ayuda de la serie de Taylor, la  ∂x i+ 1 2

función P(x ) puede ser expandida entorno al punto a mediante la siguiente serie infinita: n

 ∂p   (a ) ∞  ∂x   (x − a) n ...................................................................... (4.12) p( x ) = ∑ n ! n=0

Si el punto a en torno al cual se expande la función p(x ) es a = xi + 1 = x i + 2

se evalúa la función en x = xi , entonces x − a = −

∆xi y 2

∆xi y de la Ecuación 4.12 se 2

obtiene:

Pi = Pi+ 1

2

2  ∂p   ∆x   ∂ p  −   i + 2   ∂x i+ 12  2   ∂x i+ 1

Simulación de Yacimientos I

2

3

3  1  ∆xi   ∂ p   1  ∆xi     −  3     + L (4.13)  2!  2   ∂x i+ 1  3!  2  2 2

5

Si de nuevo se expande la función p(x ) en torno al punto a = x i + evalúa la función en x = xi+1 = xi +

1 2

= xi +

∆x i y se 2

∆xi ∆xi+1 ∆xi+1 + , entonces x − a = y de la 2 2 2

Ecuación 4.12 se obtiene:

Pi+1 = Pi+ 1 + 2

3 2 2 1  ∂ 3p   ∆xi+1   ∆xi+1  1  ∂ p   ∆xi+1  + +       +L 2  2!  ∂x 2  1  2  3!  ∂x 3  1  2   i+ 2 i+ 2 2

1  ∂p    1!  ∂x i+ 1

(4.14) Restando la Ecuación 4.13 de la Ecuación 4.14 se obtiene:  ∆x 2  ∆x 2  ∆x  1  ∂ 2p   ∂p   ∆x  i+1  −  i   + L Pi+1 − Pi =    i +1 + i  +  2    2  2  ∂x  1  2   2    ∂x i+ 12  2  i+ 2 

 ∆xi+1 ∆xi  ∆xi+1 ∆xi  ∆x  1  ∂ 2p   ∂p   ∆x − + Pi+1 − Pi =    i+1 + i  +  2    + L  2  2  ∂x  1  2 2  2 2   ∂x i+ 12  2 i+ 2  ∂p  Resolviendo para   se obtiene:  ∂x i+ 1 2

Pi+1 − Pi 1  ∂ 2p   ∆xi+1 − ∆xi  Pi+1 − Pi  ∂p  = − + Ri + 1     +L = 2   1 2 2  ∂x i+ 12 1 (∆x + ∆x ) 2  ∂x i+ 1   (∆xi+1 + ∆xi ) i +1 i 2 2 2 (4.15) Donde el término Ri + 1

2

 ∂p  es el error de truncamiento de   definido de la  ∂x i+ 1 2

siguiente forma:

Ri + 1

2

1  ∂ 2p  =− 2  ∂x 2 

Simulación de Yacimientos I

i + 12

 ∆xi+1 − ∆xi    +L 2   6

La Ecuación 4.15 conlleva a la siguiente aproximación: Pi+1 − Pi  ∂p  ..................................................................................(4.16) ≅    ∂x i+ 12 1 (∆x + ∆x ) i +1 i 2

 ∂p  4.1.2 Evaluación del término   . Si la serie de Taylor en la Ecuación 4.12  ∂x i− 1 2

se expande en torno al punto a = xi − 1 = x i − 2

punto x = xi−1 = xi −

∆x i y la función p(x ) se evalúa en el 2

∆xi ∆xi−1 − (Figura 4.3), 2 2 i− 1

i+ 1

2

i− 1

∆x

i− 1

i+ 1

i

∆x

2

∆x

i

i+ 1

Figura 4.3 Esquema que indica los limites i ± 1 / 2

entonces x − a = −

∆xi−1 y de la Ecuación 4.12 se puede escribir: 2

2 3 2 1  ∂ 3p   ∆x   ∂p   ∆x  1  ∂ p   ∆x  Pi−1 = Pi− 1 +    − i−1  +  2   − i−1  +  3   − i−1  + L 2 2  2!  ∂x  1  2  3!  ∂x  1  2   ∂x i− 12  i− 2 i− 2

(4.17) Si de nuevo se expande la función p(x ) en torno al punto a = xi − 1 = x i − 2

evalúa en el punto x = xi , entonces x − a =

Simulación de Yacimientos I

∆x i y se 2

∆x i y de la Ecuación se obtiene: 2

7

Pi = Pi− 1

2

2 3 2 1  ∂ 3p   ∆xi   ∂p   ∆xi  1  ∂ p   ∆xi  +    +L  +  3   + 3!  ∂x  1  2   ∂x i − 12  2  2!  ∂x 2 i− 1  2  i− 2 2

(4.18)

Restando la Ecuación 4.17 de la Ecuación 4.18 se tiene: 2 2 2 1  ∂ 2p   ∆xi−1   ∂p   ∆xi   ∂p   ∆xi −1  1  ∂ p   ∆xi  Pi − Pi−1 =    +   +   −  2   +L 2  ∂x  1  2   ∂x i− 12  2   ∂x i− 12  2  2  ∂x 2 i− 1  2  i− 2

2  ∂p   ∆x ∆x  1  ∂ p  Pi − Pi−1 =    i + i−1  +  2  2  2  ∂x  1  ∂x i− 12  2 i−

2

2

 ∆xi ∆xi−1  ∆xi ∆xi−1  + −   + L  2  2 2   2

 ∂p  se obtiene: Resolviendo para    ∂x i− 1 2

Pi − Pi−1 1  ∂ 2p   ∆x − ∆xi−1  Pi − Pi−1  ∂p  +L = + Ri − 1 = −  2  i    1 2 2   ∂x i− 12 1 (∆x + ∆x ) 2  ∂x i− 1  (∆xi + ∆xi−1) i i −1 2 2 2 ................................................................................................................... (4.19) Donde el término Ri − 1

2

es el error de truncamiento para la aproximación de

 ∂p  definido de la siguiente forma:    ∂x i− 1 2

Ri + 1 = − 2

1  ∂ 2p  2  ∂x 2 

i − 12

 ∆xi − ∆xi−1   +L  2  

La Ecuación 4.19 conlleva la siguiente aproximación: Pi − Pi−1  ∂p  ≅ ............................................................................ (4.20)    ∂x i− 12 1 (∆x + ∆x ) i i −1 2 Llevando las Ecuaciones 4.16 y 4.20 a la Ecuación 4.11, se obtiene: Simulación de Yacimientos I

8

∂  kA ∂p   ⋅  ≅ ∂x  µ ∂x i

 kA  2(Pi+1 − Pi )  kA  2(Pi − Pi−1 )    −   µ i+ 12 (∆xi+1 + ∆xi )  µ i− 12 (∆xi + ∆xi−1) ∆xi

..................... (4.21)

 kA   Existen varias formas para evaluar numéricamente los términos   µ i + 1

y

2

 kA    en la Ecuación 4.21. El concepto de transmisibilidad es el método más  µ i − 1 2

ampliamente aplicado.

4.1.3 El Concepto de Transmisibilidad. En general, se conoce como

transmisibilidad al término que multiplica la caída de presión en una ecuación de flujo. Por consiguiente, la transmisibilidad da una idea acerca de la facilidad del medio para permitir el flujo de fluidos a través de él y de la facilidad del fluido para desplazarse a través del medio. Considérese el esquema ilustrado en la Figura 4.4 i− 1 2

i+ 1 2

i −1 ∆xi −1 ∆x i−1 2

∆x i−1 2

∆x i 2

i

i +1

∆x i

∆xi +1 ∆x i 2

∆x i+1 2

∆x i+1 2

Fig. 4.4 – Nomenclatura utilizada en la definición del concepto de transmisibilidad

De la ley de Darcy, la tasa de flujo desde el punto xi al punto x i+1/2 está dado por:

Simulación de Yacimientos I

9

 kA   Pi+ 12 − Pi   qi, i + 1 = − 2  µ i  ∆xi 2  Por lo anterior, Pi − Pi+ 1 = 2

qi, i+ 1 ⋅ ∆xi

.................................................................................. (4.22)

2

 kA   2   µ i

Similarmente, la tasa de flujo desde el punto x i + 1 al punto x i +1 será: 2

qi + 1

2, i +1

 kA   Pi +1 − Pi+ 12   = −  µ i+1 ∆xi +1 2 

Por tanto, Pi+ 1 − Pi+1 =

qi+ 1

2

2, i +1

⋅ ∆xi +1

 kA   2  µ i +1

......................................................................... (4.23)

Sumando las Ecuaciones 4.22 y 4.23, se obtiene:   Pi − Pi+1 = q   2 

∆xi  kA     µ i

  ∆xi+1  + ............................................................. (4.24)  kA     2   µ i+1 

En la Ecuación 4.24 se asume que qi, i+ 1 = qi + 1 2

2, i +1

= q . La Ecuación 4.24 puede

ser escrita de la siguiente forma: Pi − Pi+1 =

q  µi∆xik i+1A i+1 + µi+1∆xi+1k i A i    2 k i A ik i+1A i+1 

Solucionando para q , se obtiene: Simulación de Yacimientos I

10

  2k i A ik i+1A i+1 q=  ⋅ (Pi − Pi+1) ................................................. (4.25)  µi∆xik i+1A i+1 + µi+1∆xi+1k i A i 

De acuerdo a la definición de transmisibilidad, y notando a la transmisibilidad entre el punto i y el punto i + 1 como Ti+ 1 , la rata de flujo del punto i al punto i + 1 2

estará dada por: q = Ti+ 1 ⋅ (Pi − Pi+1) ....................................................................................... (4.26) 2

De las Ecuaciones 4.25 y 4.26 se obtiene:   2k i A ik i+1A i+1 Ti+ 1 =   ........................................................... (4.27) 2  µi∆xik i+1A i+1 + µi+1∆xi+1k i A i 

De otro lado, de la ley de Darcy, el flujo entre los puntos i e i + 1 también puede ser escrito como:  kA  (Pi − Pi+1) ....................................................................... (4.28)  q =  ⋅  µ i+ 12 (∆xi + ∆xi+1 ) 2

Comparando las Ecuaciones 4.26 y 4.28, se concluye que otra expresión para la transmisibilidad entre los puntos i e i + 1 es: Ti + 1 = 2

 kA  2   ....................................................................... (4.29) (∆xi + ∆xi+1)  µ i+ 1 2

Siguiendo un procedimiento análogo al seguido en la deducción de las Ecuaciones 4.27 y 4.29, es posible obtener las siguientes dos expresiones para expresar la transmisibilidad entre los puntos i − 1 e i :   2k i−1A i−1k i A i Ti− 1 =   ........................................................... (4.30) 2  µi−1∆xi−1k i A i + µi∆xik i−1A i−1 

Simulación de Yacimientos I

11

Ti − 1 = 2

 kA  2   ...................................................................... (4.31) (∆x i −1 + ∆x i )  µ i − 1 2

Llevando las Ecuaciones 4.29 y 4.31 a la Ecuación 4.21, se obtiene:

Ti+ 1 (Pi+1 − Pi ) − Ti− 1 (Pi − Pi−1 ) ∂  kA ∂p  2 2  .................... ⋅  ≅ ∂x  µ ∂x  i ∆x i

(4.32)

Finalmente, sustituyendo la Ecuación 4.32 en la Ecuación 4.4, se obtiene: Ti+ 1 (Pi+1 − Pi ) − Ti− 1 (Pi − Pi−1 ) 2

2

∆xi

= A i ⋅ qvi

Ti + 1 (Pi +1 − Pi ) − Ti − 1 (Pi − Pi−1) = Q v i ............................................................. (4.33) 2

2

En la Ecuación 4.33, Q v i = qvi A i∆xi es la tasa de producción (en cuyo caso Q vi es positivo) o inyección (en cuyo caso Q vi es negativo) en el bloque i . La Ecuación 4.33 puede ser reordenada de la siguiente forma:

(

)

Ti + 1 Pi+1 − Ti+ 1 + Ti − 1 Pi + Ti− 1 Pi −1 = Q v i .................................................... (4.34) 2

2

2

2

La Ecuación 4.34 constituye el modelo numérico utilizado para simular el flujo de un fluido monofásico, incompresible, a través de un sistema poroso lineal. La evaluación de esta ecuación en cada bloque i del modelo ( i = 1... nx ), genera un sistema de ecuaciones cuya solución simultánea permite estimar la distribución de presiones del sistema para determinados valores de ratas de producción o inyección, Q v i . Obsérvese que si en el bloque i no hay pozos, el término Q v i se hace igual a cero en la ecuación correspondiente a dicho bloque.

Simulación de Yacimientos I

12

4.2 DISCRETIZACIÓN DE LAS CONDICIONES DE LÍMITE EXTERNO

Un yacimiento permanece en equilibrio a menos que se introduzca algún tipo de disturbio en el sistema. Este disturbio se introduce a través de los límites del yacimiento. En general, existen dos tipos de límites de yacimiento: internos y externos. Los límites internos hacen referencia a la interacción del yacimiento con los pozos y dan origen a las condiciones de límite interno. Al tratamiento de las condiciones de límite interno se le conoce como modelamiento de pozos en simulación numérica de yacimientos. Las condiciones de límite externo hacen referencia a la forma como el yacimiento interactúa con sus alrededores. En general existen dos tipos de condiciones de límite externo: Dirichlet y Neumann.

4.2.1 Discretización de Condiciones de Límite Tipo Dirichlet.

En una

condición tipo Dirichlet se especifica la variable dependiente en los límites externos del sistema. En el caso de la Ecuación 4.4, las condiciones de límite tipo Dirichlet consisten en especificar la presión de poro en los límites x = 0 y x = L , donde L representa la longitud del sistema. Analíticamente, P(0 ) = Pcero (t ) ............................................................................................ (4.35) P(L ) = PL (t ) ................................................................................................ (4.36)

Un caso particular de las Ecuaciones 4.35 y 4.36 ocurre cuando Pcero (t ) y PL (t ) son constantes e iguales a Pcero y PL , respectivamente.

Sin pérdida de

generalidad, este será el caso asumido de aquí en adelante.

Simulación de Yacimientos I

13

P(0) = Pcero

−1

1

0

x=0

P(L ) = PL

i −1

i

i +1

Celda interior

nx

n x + 1 nx + 2

x =L

Celda fantasma Celda ficticia Fig. 4.5 - Malla de bloque centrado en una dirección indicando los diferentes tipos de puntos o celdas. La Figura 4.5 ilustra un sistema lineal discretizado el cual incluye tres tipos de puntos o celdas: (i) puntos o celdas interiores (rotulados como 1 , … i − 1 , i , i + 1 … nx ), (ii) puntos o celdas fantasmas (rotulados como 0 y nx + 1 ), y (iii) puntos o celdas ficticios (rotulados como − 1 y nx + 2 ).

Los puntos interiores hacen

referencia a puntos localizados dentro del dominio físico real del sistema. Los puntos fantasmas son puntos que identifican celdas adyacentes a los límites externos del sistema y son utilizados para la discretización de las condiciones de límite externo.

Los puntos ficticios son puntos que identifican celdas que se

adicionan a los lados de las celdas fantasmas y son utilizadas por conveniencia en programación y eficiencia computacional. Obsérvese que la condición de límite externa izquierda se aplica en las interfase que une el bloque fantasma 0 y el primer bloque interior. Así mismo, la condición de límite externa derecha se aplica en las interfase que une el último bloque interior y el bloque fantasma nx + 1 . Por tratarse de una malla de bloque centrado, el punto de aplicación de las condiciones de límite tipo Dirichlet no coincide con la ubicación de ninguno de los nodos o puntos discretos. Por tal motivo, las condiciones de límite tipo Dirichlet se suelen discretizar aplicando promedios de las variables discretas en los bloques Simulación de Yacimientos I

14

adyacentes a la interfase que representa el límite externo. Por ejemplo, para el sistema ilustrado en la Figura 4.5, la condición en x = 0 se discretiza de la siguiente forma: P0 + P1 = Pcero 2

O bien, P0 + P1 = 2Pcero ......................................................................................... (4.37) En forma análoga, la condición en x = L se discretiza de la siguiente forma: Pnx + Pnx +1 = 2PL ....................................................................................... (4.38)

Obsérvese que la discretización de las condiciones de límite adiciona dos nuevas incógnitas al sistema de ecuaciones representado por la Ecuación 4.34 ( Po y Pnx +1 ) lo que hace necesario el planteamiento de dos ecuaciones adicionales

(Ecuaciones 4.37 y 4.38).

4.2.2 Discretización de Condiciones de Límite Tipo Neumann.

En una

condición tipo Neumann se especifica la rata de flujo a través de la frontera del yacimiento. En el caso de la Ecuación 4.4, una condición Neumann en el límite x = 0 puede expresarse de la siguiente forma:

(q)x =0

 kA ∂p   = q cero ...................................................................... (4.38) =   µ ∂x  x =0

En la Ecuación 4.38, qcero es un valor conocido. La Ecuación 4.38 puede ser escrita de la siguiente forma: q µ  ∂p  = cero ...................................................................................... (4.39)   kA  ∂x  x = 0 Simulación de Yacimientos I

15

En el caso de un sistema cerrado, no existe flujo a través de sus límites externos y, en consecuencia, q cero = 0 . En este caso la Ecuación 4.39 toma la siguiente forma:  ∂p  =0    ∂x  x =0

(4.40)

Análogamente, en el caso general, una condición Neumann en el límite x = L puede expresarse de la siguiente forma: q µ  ∂p  = L .......................................................................................... (4.41)    ∂x  x =L kA

En la Ecuación 4.41, qL es un valor conocido. En particular, para un sistema cerrado, la Ecuación 4.41 toma la siguiente forma:  ∂p  = 0 .............................................................................................. (4.42)    ∂x  x =L Para una malla de bloque centrado, las condiciones de límite tipo Neumann se suelen discretizar calculando el gradiente entre los puntos discretos adyacentes a la interfase que representa el límite externo. Por ejemplo, para el sistema ilustrado en la Figura 4.5, una condición tipo Neumann, Ecuación 4.41, en x = 0 se discretiza de la siguiente forma: P1 − P0 qcero µ = x1 − x 0 kA

O bien, P1 − P0 =

q cero µ (x 1 − x 0 ) .......................................................................... (4.43) kA

En particular, para un sistema cerrado, la Ecuación 4.43 toma la siguiente forma:

Simulación de Yacimientos I

16

P1 − P0 = 0 ................................................................................................ (4.44)

En x = L , la forma discretizada de las Ecuaciones 4.41 y 4.42 toma la siguiente forma: Pnx +1 − Pnx q µ = L x nx +1 − x nx kA

O bien, Pnx +1 − Pnx =

qL µ (x nx +1 − x nx ) ................................................................. (4.45) kA

En particular, para un sistema cerrado, la Ecuación 4.45 toma la siguiente forma: Pnx +1 − Pnx = 0 .......................................................................................... (4.46)

Al igual que el caso de las condiciones de límite tipo Dirichlet, la discretización de las condiciones de límite tipo Neumann adiciona dos nuevas incógnitas al sistema de ecuaciones representado por la Ecuación 4.34. Estas dos nuevas incógnitas son P0 y Pnx +1 . Este hecho hace necesario el planteamiento de dos ecuaciones adicionales (Ecuación 4.43 o 4.44, para x = 0 , y Ecuación 4.45 o 4.46, para x = L ).

4.3 EL CONCEPTO DE STENCILS

El concepto de stencil permite expresar los modelos numéricos en forma generalizada. Con la finalidad de ilustrar el concepto de stencil, considérese un sistema lineal, abierto, conformado por nx celdas, cada una de las cuales tiene sus propias características (permeabilidad, área, longitud, etc.), tal como se ilustra en Simulación de Yacimientos I

17

la Figura 4.5. Al aplicar la Ecuación 4.34 a cada celda se genera el siguiente sistema de nx ecuaciones: •

Para i = 1 :

(

)

T 1 P0 − T 1 + T3 P1 + T3 P2 = Q v 1 ............................................................. (4.47) 2

2

2

2

La Ecuación 4.47 conecta el nodo 1 con los nodos 0 y 2 a través de los lados Oeste y Este, respectivamente. Esta conexión se suele expresar en forma general

expresando la Ecuación 4.47 de la siguiente forma: W1P0 + C1P1 + E1P2 = F1 .............................................................................. (4.48) En la Ecuación 4.48, W , C y E representan las orientaciones Oeste, Central y Este, respectivamente, F indica término independiente y el subíndice 1 de W , C

y E indica que se trata de la ecuación para el nodo o celda 1 . Al comparar las   Ecuaciones 4.47 y 4.48 es evidente que W1 = T 1 , C1 = − T 1 + T3  , E1 = T3 2 2  2 2 y F = Q v1 . •

Para i = 2 :

(

)

T3 P1 − T3 + T5 P2 + T5 P3 = Q v 2 ............................................................. (4.49) 2

2

2

2

O bien, W 2P1 + C 2P2 + E 2P3 = F2 ........................................................................... (4.50) •

Para i = 3 :

(

)

T5 P2 − T5 + T7 P3 + T7 P4 = Q v 3 ........................................................... (4.51) 2

2

2

2

O bien, W3 P2 + C 3 P3 + E 3 P4 = F3 .......................................................................... (4.52) Simulación de Yacimientos I

18

•

En general, para i = i :

(

)

Ti− 1 Pi−1 − Ti− 1 + Ti+ 1 Pi + Ti+ 1 Pi+1 = Q v i .................................................. (4.53) 2

2

2

2

O bien, WiPi−1 + C iPi + E iPi+1 = Fi ............................................................................. (4.54)

•

Finalmente, para i = nx :

(

)

Tn− 1 Pn−1 − Tn− 1 + Tn+ 1 Pn + Tn+ 1 Pn+1 = Q v n .......................................... (4.55) 2

2

2

2

O bien, Wnx Pnx −1 + C nx Pnx + E nx Pnx +1 = Fnx ........................................................... (4.56)

Las Ecuaciones 4.47 a 4.56 representan un sistema de nx ecuaciones con nx + 2 incógnitas. Por consiguiente, se requiere de dos ecuaciones adicionales para que el sistema esté matemáticamente bien planteado.

Estas dos ecuaciones

adicionales se obtienen de las condiciones de límite en x = 0 y x = L .

La

condición en x = 0 puede ser escrita, en forma de stencil, de la siguiente forma: W0P−1 + C 0P0 + E 0 P1 = F0 .......................................................................... (4.57) La Ecuación 4.57 toma la forma general de los stencils de las ecuaciones en las celdas interiores. Obsérvese, sin embargo, que en la Ecuación 4.57 siempre se cumple que W0 = 0 .

Adicionalmente, P−1 hace referencia a la presión en el

bloque ficticio de la izquierda. Los valores que toman las componentes C 0 , E 0 y F0 dependen de las condiciones de límite en x = 0 .

Para condiciones tipo

Dirichlet, Ecuación 4.37, C 0 = 1 , E 0 = 1 y F0 = 2Pcero . Para condiciones de tipo Neumann, Ecuación 4.43, C 0 = −1, E 0 = 1 y F0 =

Simulación de Yacimientos I

q cero µ (x 1 − x 0 ) . kA 19

La condición en x = L puede ser escrita de la siguiente forma: Wnx +1Pnx + C nx +1Pnx +1 + E nx +1Pnx + 2 = Fnx +1| ................................................ (4.58) En la Ecuación 4.58 siempre se cumple que E nx +1 = 0 . Adicionalmente, Pnx + 2 hace referencia a la presión en el bloque ficticio de la derecha. Los valores que toman las componentes Wnx +1 , C nx +1 y Fnx +1 dependen de las condiciones de límite en x = L .

Para condiciones tipo Dirichlet, Ecuación 4.38, Wnx +1 = 1 ,

C nx +1 = 1 y Fnx +1 = 2PL .

Para condiciones de tipo Neumann, Ecuación 4.45,

Wnx +1 = −1 , C nx +1 = 1 , y Fnx +1 =

qL µ (x nx +1 − x nx ) . kA

El sistema de Ecuaciones 4.47 a 4.58 puede ser escrito en forma matricial como se ilustra en la Figura 4.6, donde la matriz de coeficientes está compuesta por las transmisibilidades, Ti− 1 y 2

Ti+ 1 , el vector de incógnitas está formado por las 2

presiones, Pi , y el vector de términos independientes por las tasas de producción o inyección en cada bloque, Q v i .

C0 E 0 W  1 C1  W2       

E1 C2 W3 O

E2 C3 O Wnx

E3 O Cnx Wnx +1

        E nx  Cnx +1 

 P0   F0   P   F   1   1   P2   F2       M = M   M   M       Pnx   Fnx   P    nx +1  F nx +1 

Fig. 4.6 – Matriz de coeficientes de un sistema de ecuaciones tridiagonal. (Simulador de un sistema lineal). Los coeficientes son escritos en términos de las componentes del stencil. Simulación de Yacimientos I

20

Con lo cual se puede concluir que la matriz generada por el modelo numérico de un sistema lineal (1D) es una matriz tridiagonal. Un sistema tridiagonal de n ecuaciones con n incógnitas admite algoritmos especiales que conllevan a una solución más rápida que un sistema no tridiagonal de ecuaciones, como es el caso del Algoritmo de Thomas.

Simulación de Yacimientos I

21

5. FLUJO BIDIMENSIONAL DE UN FLUIDO INCOMPRESIBLE

Este capítulo se dedica a discutir los conceptos relacionados con la simulación del flujo en dos direcciones de un fluido incompresible. Inicialmente se presenta el modelo matemático, partiendo de las ecuaciones de flujo. Luego, se plantean las aproximaciones numéricas a estas ecuaciones de flujo.

Posteriormente, se presentan algunos esquemas de ordenamiento que

facilitan su solución. Por último, se discuten los principales métodos directos e indirectos de solución a las ecuaciones planteadas.

5.1 ECUACIÓN BÁSICA La ecuación de continuidad para flujo en dos direcciones, en coordenadas cartesianas, está dada por la Ecuación 2.56:

(

)

 ∂ (αρ U x ) ∂ αρ U y  ∂ (φρ) − + + α~ q ..............................................................(2.56) =α ∂ ∂ x y ∂ t   donde: α = Espesor del sistema, H (x, y ) .

~ q = Cantidad de masa que entra o sale (fuentes o sumideros) por unidad de volumen del yacimiento por unidad de tiempo. Para un fluido incompresible, se cumple que:

ρ = constante

Simulación de Yacimientos I

1

luego, de la Ecuación 2.56,

(

)

~  ∂ (αU x ) ∂ αU y  ∂φ q − +  =α +α ∂y  ∂t ρ  ∂x Si se considera que la porosidad es constante:

(

)

~  ∂ (αU x ) ∂ αU y  q − +  − α = 0 ............................................................................(5.1) ∂y  ρ  ∂x De la ley de Darcy (Ecuaciones 2.31 y 2.32),

Ux = −

k x ∂p ⋅ ..................................................................................................(2.31) µ ∂x

Uy = −

k y ∂p ⋅ ...................................................................................................(2.32) µ ∂y

Además, si q v se define como el volumen de fluido que entra o sale (fuentes o sumideros) por unidad de volumen del yacimiento por unidad de tiempo, entonces: qv =

~ q ................................................................................................................(5.2) ρ

Llevando las Ecuaciones 2.31, 2.32 y 5.2 a la Ecuación 5.1, se obtiene: ∂  k x ∂p  ∂  k y ∂p  H ⋅  + H ⋅  − Hq v = 0 ............................................................(5.3) ∂x  µ ∂x  ∂y  µ ∂y 

La Figura 5.1 presenta un sistema bidimensional de J filas por I columnas. Supóngase que se extrae la fila j , tal como se ilustra en la figura 5.2.

Simulación de Yacimientos I

2

j=J

L

L

M

M

j= j

L

M

i − 1, j + 1

i, j + 1

i + 1, j + 1

i − 1, j

i, j

i + 1, j

i − 1, j − 1

i, j − 1

i + 1, j − 1

M

L

M

M

L

L

j=1

i=1

i=i

i = I

Fig. – 5.1 – Sistema bidimensional (J filas x I columnas)

j= j

i − 1, j

1, j

∆Xi −1

i, j ∆Xi

I, j

i + 1, j ∆Xi +1

Figura 5.2 Fila “j” de una malla bidimensional

En este caso, de la Ecuación 4.21 se puede escribir:

(

∂  k x H ∂p   ⋅  ≅ ∂x  µ ∂x  i

)

(

)

 2 Pi, j − Pi −1, j   2 Pi +1, j − Pi, j   k x H   k xH         −   µ  i + 12 , j  ∆x i +1 + ∆x i   µ  i − 12 , j  ∆x i + ∆x i −1 

Simulación de Yacimientos I

∆x i

............(5.4)

3

Análogamente, si se fija una columna “ i ”, Figura 5.3, de la Ecuación 4.18 se tiene:

J

j=J

j+1

∆Y

j

∆Y

j

j−1

∆Y

j−1

j= 2

2

j=1

1

j+1

i= i

Figura 5.3 Columna “i” de una malla bidimensional

(

∂  k y H ∂p   ⋅  ≅ ∂y  µ ∂y  i

)

(

)

 2 Pi, j + 1 − Pi, j   k y H   2 Pi, j − Pi, j −1   k yH      −      µ      i, j + 12  ∆y j + 1 + ∆y j   µ  i, j − 12  ∆y j + ∆y j −1  ∆y j

.............(5.5)

EL CONCEPTO DE TRANSMISIBILIDAD EN FLUJO BIDIMENSIONAL

Al igual que en el caso de flujo lineal, en flujo bidimensional la transmisibilidad está representada por el término que multiplica la caída de presión.

Simulación de Yacimientos I

4

Considérese el sistema bidimensional ilustrado en la Figura 5.4.

j+1

i − 1, j + 1

∆ Y j +1

i − 1, j

∆Y j

i − 1, j − 1

∆Y j −1

∆ X i −1

j

∆ Y j +1

∆X i

∆ X i −1 j−1

i, j + 1

i, j

∆ X i +1

∆Y j

i + 1, j

∆X i i, j − 1

i + 1, j + 1

∆ X i +1 ∆ Y j −1

i + 1, j − 1

∆ X i −1

∆X i

∆ X i +1

i −1

i

i +1

Figura 5.4 Nomenclatura utilizada en el concepto de transmisibilidad (Flujo bidimensional)

La tasa de flujo desde el punto

(i, j)

al punto (i +

1 ,j 2

)

está dado por la ley de

Darcy:

qi, j

a i+ 12, j

 kA   Pi, j − Pi+ 12, j   =   µ  i, j  ∆x i 2 

de donde:

Pi, j − Pi+ 1 , j = 2

qi, j

a i+ 12, j

 kA     µ  i, j

⋅

∆x i ...............................................................................(5.6) 2

Análogamente, la tasa de flujo desde punto (i +

Simulación de Yacimientos I

1 ,j 2

) al punto (i + 1, j) será :

5

q i+ 1 , j 2

 kA   Pi+ 12, j − Pi+1, j   =   µ  i+1, j  ∆x i+1 2 

a i+1, j

o bien:

Pi+ 1 , j − Pi+1, j − =

q i+ 1 , j

2

2

a i+1, j

 kA     µ  i+1, j

⋅

∆x i+1 ....................................................................(5.7) 2

Sumando las Ecuaciones 5.6 y 5.7, se tiene:

Pi, j − Pi+1, j

Pi, j − Pi + 1, j

   =q  2 

∆xi  kA     µ i, j

  ∆xi+1  +  kA     2  µ i+1, j  

  ∆x i +1 q  ∆x i = +  2 k i, j A i, j k i +1, j A i +1, j  µ i +1, j  µ i, j

Pi, j − Pi +1, j =

     

q  k i +1, j A i + 1, jµ i, j ∆x i + k i, j A i, jµ i +1, j ∆x i +1    2  k i, j A i, jk i +1, j A i + 1, j 

(

)

  2 ⋅ k i, j A i, jk i +1, j A i + 1, j q=  ⋅ Pi, j − Pi +1, j  µ i, j ∆x ik i +1, j A i + 1, j + µ i + 1, j ∆x i +1k i, j A i, j 

q = Ti + 1

2, j

(

(

)

)

⋅ Pi, j − Pi + 1, j ........................................................................................(5.8)

donde:

Simulación de Yacimientos I

6

(

)

Ti + 1

  2 ⋅ k i, j A i, jk i +1, j A i +1, j =   µ i, j ∆x ik i +1, j A i +1, j + µ i +1, j ∆x i +1k i, j A i, j 

Ti + 1

  2 ⋅ k i, jHi, j ∆y jk i + 1, jHi + 1, j ∆y j =   µ i, j ∆x ik i + 1, jHi + 1, j ∆y j + µ i + 1, j ∆x i +1k i, jHi, j ∆y j 

Ti + 1

  2 ⋅ k i, jHi, j ∆y jk i +1, jHi +1, j =  .....................................................(5.9)  µ i, j ∆x ik i + 1, jHi + 1, j + µ i + 1, j ∆x i + 1k i, jHi, j 

2, j

2, j

2, j

(

(

)

)

Siguiendo un procedimiento análogo, se puede obtener:

Ti − 1

2, j

(

)

(

)

(

)

  2 ⋅ k i, jHi, j ∆y jk i −1, jHi −1, j =  ...................................................(5.10)  k i −1, jHi −1, jµ i, j ∆x i + µ i −1, j ∆x i −1k i, jHi, j 

  2 ⋅ k i, jHi, j ∆x ik i, j +1Hi, j +1 Ti, j + 1 =   ..................................................(5.11) 2  k i, j +1Hi, j +1µ i, j ∆y j + µ i, j +1∆y j +1k i, jHi, j    2 ⋅ k i, jHi, j ∆x ik i, j −1Hi, j −1 Ti, j − 1 =   ...................................................(5.12) 2  k i, j −1Hi, j −1µ i, j ∆y j + µ i, j −1∆y j −1k i, jHi, j 

La tasa de flujo desde el punto (i, j) al punto (i + 1, j) será:  kA   q =   µ  i+ 12, j

(Pi,j − Pi+1,j ) 1 2

(∆x i + ∆x i+1 )

.............................................................................(5.13)

Si se compara las Ecuaciones 5.8 y 5.13, se obtiene:

Simulación de Yacimientos I

7

 kA  2  Ti+ 1 , j =  .........................................................................(5.14) 2  µ  i+ 12, j (∆x i + ∆x i+1 )

o bien:

Ti + 1

2, j

Ti+ 1 , j 2

∆y j

 kH∆y j  2  =    µ  i + 12 , j (∆x i + ∆x i +1 )  kH  2  =  .........................................................................(5.15)  µ  i+ 12, j (∆x i + ∆x i+1 )

Análogamente, Ti− 1 , j 2

∆y j Ti, j+ 1

2

∆x i Ti, j− 1

2

∆x i

 kH  2  =  .........................................................................(5.16)  µ  i− 12, j (∆x i−1 + ∆x i )  kH  2  ........................................................................(5.17) =   µ  i,j+ 12 ∆y j + ∆y j+1

(

)

 kH  2  =  ........................................................................(5.18)  µ  i, j− 12 ∆y j−1 + ∆y j

(

)

Sustituyendo las Ecuaciones 5.15 y 5.16 en la Ecuación 5.4, se obtiene:

(

)

(

)

∂  k x H ∂p  Ti+ 12, j Pi+1, j − Pi, j − Ti− 12, j Pi,j − Pi−1, j  ⋅ ≅ ..........................................(5.19) ∂x  µ ∂x  ∆y j ∆x i

Sustituyendo las Ecuaciones 5.17 y 5.18 en la Ecuación 5.5, se obtiene:

(

)

(

)

∂  k x H ∂p  Ti, j+ 12 Pi, j+1 − Pi, j − Ti, j− 12 Pi, j − Pi, j−1  ⋅ ≅ ..........................................(5.20) ∂y  µ ∂y  ∆y j ∆x i

Simulación de Yacimientos I

8

Finalmente, llevando las Ecuaciones 5.19 y 5.20 a la Ecuación 5.3, se obtiene: Ti + 1

2, j

(Pi +1, j − Pi, j ) − Ti − , j (Pi, j − Pi −1, j ) 1 2

∆y j ∆x i

+

(

)

(

Ti, j + 1 Pi, j +1 − Pi, j − Ti, j − 1 Pi, j − Pi, j −1 2

2

∆y j ∆x i

)

− Hi, j q v i, j = 0

o bien:

(

)

(

)

(

)

(

Ti+ 1 , j Pi+1, j − Pi, j − Ti− 1 , j Pi, j − Pi−1, j + Ti, j+ 1 Pi, j+1 − Pi, j − Ti, j− 1 Pi, j − Pi, j−1 2

2

2

2

) − ∆y j ∆x iHi, j q v i, j = 0

Reordenando términos:

(

)

Ti− 1 , jPi−1, j + Ti, j− 1 Pi, j−1 − Ti− 1 , j + Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j Pi, j + Ti, j+ 1 Pi, j+1 + Ti+ 1 , jPi+1, j 2

2

2

2

2

2

2

2

− Q i, j = 0

.........................................................................................................................(5.21) La Ecuación 5.21 representa el modelo numérico para el flujo de un fluido incompresible en dos dimensiones. Al evaluar esta ecuación en cada punto (i, j) de la malla, se genera un sistema de ecuaciones cuya solución permite estimar la distribución de presiones en el sistema en función de la rata de producción o inyección de cada bloque Qi, j .

5.2 ESQUEMAS DE ORDENAMIENTO

Para simplificar un poco, la Ecuación 5.21 puede ser escrita de la siguiente forma: a i, jPi −1, j + bi, jPi, j −1 + c i, jPi, j + di, jPi, j +1 + e i, jPi +1, j = fi, j ...........................................(5.22) donde:

Simulación de Yacimientos I

9

ai, j = Ti− 1 , j ........................................................................................................(5.23) 2

b i, j = Ti, j− 1 ........................................................................................................(5.24) 2

(

)

c i, j = Ti− 1 , j + Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j .................................................................(5.25) 2

2

2

2

di, j = Ti, j + 1 .......................................................................................................(5.26) 2

e i, j = Ti+ 1 , j ........................................................................................................(5.27) 2

fi, j = Qi, j ............................................................................................................(5.28)

La Ecuación 5.22 genera una ecuación por cada bloque en la malla.

Dicha

ecuación incluye la presión del bloque en cuestión y las de sus cuatro bloques adyacentes en cruz, tal como se ilustra en la Figura 5.5.

j−1

j

j+1

i, j + 1

i − 1, j

i, j

i + 1, j

i, j − 1

Fig 5.5. Bloques cuyas presiones se tienen en cuenta en la aplicación de la ecuación básica para el bloque (i,j)

Simulación de Yacimientos I

10

El conjunto de todas las ecuaciones forman un sistema de ecuaciones cuyas incógnitas son las presiones.

La forma de la matriz de coeficientes de este

sistema de ecuaciones depende del orden seguido para escribir las ecuaciones (es decir, la forma como se recorra la malla). ESQUEMA DE ORDENAMIENTO.

A este orden se le denomina

El trabajo de computador y el número de

operaciones requerido para solucionar el sistema de ecuaciones es función, en cierto grado, del esquema de ordenamiento.

A continuación se presentan los

esquemas de ordenamiento más comúnmente utilizados.

5.2.1 ORDENAMIENTO NORMAL POR FILAS ( O ÍNDICE i ).

Como su nombre lo indica, en el ordenamiento normal por filas, el sistema de ecuaciones se genera recorriendo la malla por filas. Considérese el ejemplo del ordenamiento normal por filas ilustrado en la Figura 5.6. j=4

10

11

12

j=3

7

8

9

j=2

4

5

6

j=1

1

2

3

i =1

i=2

i=3

Figura 5.6 Ordenamiento normal por filas

Partiendo de la ecuación básica, Ecuación 5.22, escribir el conjunto de ecuaciones y la matriz de coeficientes generadas al efectuar el recorrido de la malla.

Simulación de Yacimientos I

11

Ecuación 1: i = 1 , j = 1 (bloque 1 de la Figura 5.6). De la Ecuación 5.22: a1,1 = T1− 1 ,1 = T 1 ,1 = 0 2

2

b1,1 = T1, 1 = 0 2

(

c 1,1 = − T 1 ,1 + T3 2

d1,1 = T1,3 e1,1 = T3

2,1

+ T1, 1 + T1,3 2

2

)

2

2,1

Luego: c 1,1P1,1 + d1,1P1,2 + e1,1P2,1 = f1,1 ...........................................................................(5.29)

Ecuación 2: i = 2 , j = 1 (bloque 2 en la Figura 5.6 ). En este caso: T2− 1 ,1 = Ti, j− 1 = 0 2

2

Luego: b 2,1 = 0 Por tanto, de la Ecuación 5.22: a 2,1P1,1 + c 2,1P2,1 + d 2,1P2,2 + e 2,1P3,1 = f 2,1 ..........................................................(5.30)

Simulación de Yacimientos I

12

Ecuación 3: i = 3 , j = 1 (bloque 3 en la Figura 5.6). En este caso: Ti, j− 1 = Ti+ 1 , j = 0 2

2

luego: b i, j = e i, j = 0 Entonces, de la Ecuación 5.22: a 3,1P2,1 + c 3,1P3,1 + d 3,1P3,2 = f3,1 ........................................................................(5.31)

Ecuación 4: i = 1 , j = 2 (bloque 4 en la Figura 5.6). En este caso: Ti− 1 , j = a i,j = 0 2

Luego, de la Ecuación 5.22: b1,2P1,1 + c 1,2P1,2 + d1,2P1,3 + e1,2P2,2 = f1,2 ..........................................................(5.32)

Ecuación 5: i = 2 , j = 2 (bloque 5 en la Figura 5.6). De la Ecuación (5.22): a 2,2P1,2 + b 2,2P2,1 + c 2,2P2,2 + d 2,2P2,3 + e 2,2P3,2 = f 2,2 ......................................(5.33)

Simulación de Yacimientos I

13

Las anteriores ecuaciones conllevan a un sistema matricial como el indicado en la Figura 5.7.

P11  c11 a  21    b21            

P21

P31

e11 c21 a31

P12

P22

P32

P13

P23

P33

P14

P24

d 11 e21 c31

d 21 d 31 c12 a22

b22 b32

e11 c22

e22

d 12

a32

c32

b13 b23

d 22 d 32 c13

e13

a23

c23 a33

b33

d13 e23 c33

b14

d 23 c14 a24

b24 b34

e14 c24 a34

P34             d 33     e24  c34 

 P11  P   21   P31     P12   P22     P32  P   13   P23  P   33   P14     P24   P34 

=

 f 11  f   21   f 31     f 12   f 22     f 32  f   13   f 23  f   33   f14     f 24   f 34 

Figura 5.7 Sistema matricial obtenido del ordenamiento por filas de una malla 4*3

Observaciones: a. Una matriz como la presentada en la Figura 5.7 se define como una MATRIZ DE BANDA PENTADIAGONAL, pues contiene 5 bandas. b. En una matriz de banda, se define el ANCHO DE BANDA como el número máximo de elementos en una fila abarcados entre dos elementos diferentes de cero. Como se observa de la Figura 5.7, en este ejemplo, el ancho de banda es: 2 × 3 + 1 = 7 . En general, en el ordenamiento por filas, el ancho de banda está dado por: 2 × (Numero de Columnas ) + 1 .....................................................................(5.34)

Simulación de Yacimientos I

14

Es importante anotar que el ancho de banda depende del esquema de ordenamiento. En general, a mayor ancho de banda, mayor es el almacenaje y tiempo de computador requerido para solucionar el sistema de ecuaciones por métodos directos.

5.2.2

ORDENAMIENTO NORMAL POR COLUMNAS ( O ÍNDICE j )

En el ordenamiento normal por columnas, el sistema de ecuaciones se genera recorriendo la malla por columnas, tal como se indica en la Figura 5.8.

j=4

4

8

12

j=3

3

7

11

j=2

2

6

10

j=1

1

5

9

i =1

i=2

i=3

Figura 5.8 Ordenamiento normal por columnas

Para este ordenamiento, el ancho de banda está dado por: 2 × (Numero de Filas ) + 1 ...................................................................................(5.35)

La matriz de coeficientes tiene la forma indicada en la Figura 5.9 (las x indican posiciones diferentes de cero).

Simulación de Yacimientos I

15

x x     x           

x

x

x x

x

x x x

x

x x

x x x

x

x

x x x x

x

x x x x

x

x x x

x x x

x x x x x

x x

x

           x     x x 

Figura 5.9 Forma de la matriz de coeficientes para un ordenamiento normal por columnas.

5.2.3

ORDENAMIENTO D2

En el ordenamiento D 2 la secuencia de las ecuaciones se determina recorriendo la malla diagonalmente en forma continua, como se indica en la Figura 5.10.

j=4

7

10

12

j=3

4

8

11

j=2

2

5

9

j=1

1

3

6

i =1

i=2

i=3

Figura 5.10 Ordenamiento D2 La matriz de coeficientes resultante de este ordenamiento tiene la forma indicada en la Figura 5.11

Simulación de Yacimientos I

16

Figura 5.11 Esquema de la matriz de coeficientes del ordenamiento D2

5.2.4 ORDENAMIENTO D4 O BLANCO - NEGRO

En el ordenamiento D4 se recorre la malla diagonalmente en forma alternada. La Figura 5.12 presenta un ejemplo del recorrido de malla seguido en este tipo de ordenamiento. Se conoce tambien como blanco y negro ya que el ordenamiento se lleva a cabo como si la malla fuera un tablero de ajedrez. j=4

9

5

12

j=3

2

10

6

j= 2

7

3

11

j=1

1

8

4

i =1

i=2

i=3

Figura 5.12 Ordenamiento D4

La matriz de coeficientes presenta la forma indicada en la Figura 5.13.

Simulación de Yacimientos I

17

Figura 5.13 Esquema de la matriz de coeficientes del ordenamiento D4

5.3

MÉTODOS DIRECTOS DE SOLUCIÓN DE ECUACIONES.

Tal como se indicó antes, al aplicar la Ecuación 5.21 a cada bloque de la malla se genera un sistema de ecuaciones cuyas incógnitas son las presiones en los bloques. La forma general de este sistema de ecuaciones es la siguiente: a11x 1 + a12 x 2 a 21x 1 + a 22 x 2 a 31x 1 + a 32 x 2 M a n1x 1 + a n2 x 2

+ a13 x 3 + a 23 x 3 + a 33 x 3 + a n3 x 3

+ K + a1n x n + K + a 2n x n + K + a 3n x n M + K + a nn x n

= c1 = c2 = c3 M = cn

(5.36)

donde a11 , a12 , K , a1n , a 21 , a 22 , K , a 2n , K , a nn son los coeficientes de las incógnitas (en nuestro caso transmisibilidades o combinación de ellas) x 1 , x 2 , x 3 , K , x n representan las incógnitas (en nuestro caso presiones).

Simulación de Yacimientos I

18

La solución de este tipo de sistemas de ecuaciones se puede realizar aplicando métodos directos o métodos iterativos.

Dos de los métodos directos más

comúnmente aplicados son el de eliminación Gaussiana y el método de GaussJordan. Otros métodos también utilizados son: inversión matricial, la regla de Cramer y el método de descomposición matricial. A continuación se discuten los métodos de eliminación Gaussiana y de Gauss-Jordan.

5.3.1 Eliminación Gaussiana

La solución de un sistema de ecuaciones, tal como el representado por la expresión (5.36), mediante el método de eliminación Gaussiana consiste en transformar el sistema inicial de ecuaciones a un sistema equivalente que tenga la siguiente forma: ′ x2 x 1 + a12 x2

′ x3 + a13 + a ′23 x 3 x3 O

+K +K +K

xn

+ a1′ n x n = c 1′ + a ′2n x n = c ′2 + a ′3n x n = c ′3 M M + a ′n−1,n x n−1 = c ′n−1 = c ′n xn

(5.37)

cuya solución es igual a la solución del sistema original.

Considérese el siguiente sistema de cuatro ecuaciones con cuatro incógnitas: a11x 1 + a12 x 2 a 21x 1 + a 22 x 2 a 31x 1 + a 32 x 2 a 41x 1 + a 42 x 2

+ a13 x 3 + a 23 x 3 + a 33 x 3 + a 43 x 3

Simulación de Yacimientos I

+ a14 x 4 + a 24 x 4 + a 34 x 4 + a 44 x 4

= c1 = c2 = c3 = c4

19

La matriz aumentada de este sistema de ecuaciones tiene la siguiente forma:  a11  a 21 a 31  a 41

a12 a 22 a 32 a 42

a13 a 23 a 33 a 43

a14 a 24 a 34 a 44

c1   c2  c3   c 4 

(5.38)

Supóngase que se efectúan las siguientes operaciones elementales sobre esta matriz: a. Se divide la fila 1 por a11 . b. Se multiplica la fila 1 resultante de la operación anterior por (− a 21 ) y se suma a la segunda fila. c. Se multiplica la fila 1 resultante de la operación indicada en el literal “a” por

(− a 31 ) y se suma a la tercera fila. d. Se multiplica la fila 1 resultante de la operación indicada en el literal “a” por

(− a 41 ) y se suma a la cuarta fila. La matriz aumentada resultante tendrá la siguiente forma: 1  0  0 0 

(1) a12

(1) a13

(1) a14

1) a (22

1) a (23

1) a (24

1) a (42

1) a (43

1) a (44

1) a (32

1) a (33

1) a (34

c 1(1)   c (21)   c 3(1)  c (41) 

Considérese que sobre esta matriz se realizan las siguientes operaciones elementales:

Simulación de Yacimientos I

20

(1) . a. Se divide la segunda fila por a 22

(

(1) b. Se multiplica la fila 2 resultante del paso anterior por − a 32

)

y se suma a

la fila 3.

(

1) c. Se multiplica la fila resultante del paso 1 por − a (42

)

y se suma a la fila 4.

La nueva matriz aumentada tendrá la siguiente forma:

 1 a (1) 12  0 1  0 0 0 0 

(1) a13

2) a (23

2) a (33

2) a (43

(1) a14

2) a (24

2) a (34

2) a (44

c 1(1)   c (22 )   c 3(2 )  c (42 ) 

Sobre la matriz anterior se pueden realizar las siguientes dos operaciones:

(2 ) a. Se divide la fila 3 por a 33

(

)

2) b. La nueva fila 3 se multiplica por − a (43 y se suma a la fila 4.

La nueva matriz aumentada tendrá la siguiente forma:  1 a (1) 12  0 1  0 0 0 0 

(1) a13

2) a (23

1 0

(1) a14

2) a (24

3) a (34

3) a (44

c 1(1)   c (22 )   c 3(3 )  c (43 ) 

3) Dividiendo la fila 4 por a (44 , se obtiene la siguiente matriz aumentada:

Simulación de Yacimientos I

21

 1 a (1) 12  0 1  0 0 0 0

(1) a13

(1) a14

2) a (23

2) a (24 3) a (34 1

1 0

c 1(1)   c 2(2 )  c 3(3 )  c (44 ) 

la cual corresponde al siguiente sistema de ecuaciones equivalente al sistema inicial:

(1) x x 1 + a12 2 x2

(1) x + a13 3 ( 2) +a x

(1) x + a14 4 ( 2) +a x 24

4

x3

34

4

23 3

+ a (3 ) x x4

= c 1(1) = c (2 ) 2

= c (3 )

(5.39)

3

= c (4 ) 4

El valor de las incógnitas puede ser estimado mediante un sustitución de atrás hacia adelante: x4 x3 x2 x1

= c (44 ) = c (3 ) 3 = c (22 ) = c 1(1)

(3 ) x − a 34 4 ( 2) −a x 24 4 (1) x − a14 4

(2 ) x − a 23 3 ( 1) −a x 13 3

(1) x − a12 2

Este proceso se puede generalizar mediante el siguiente algoritmo, fácilmente ejecutable mediante un programa de computador: a. Se inicializa la variable k , k = 1 b. Se normaliza la fila k : a k, j =

a k, j a k, k

,

Simulación de Yacimientos I

j = k, k + 1, k + 2, K, n

22

donde n es el número total de filas de la matriz aumentada (e.d., el número de ecuaciones). c. Se reduce la matriz, así: a i, j = a i, j − a i,k ⋅ a k, j ,

j = k + 1, k + 2, K, i = k + 1, k + 2, k + 3, K

d. Se incrementa k en una unidad:

k = k +1

e. Se repite el procedimiento anterior hasta que k sea igual a n .

Al finalizar el anterior ciclo, se obtiene el siguiente sistema de ecuaciones en forma triangular: ′ x2 x 1 + a12 x2

′ x3 + a13 + a ′23 x 3 x3

+ K + a1′,n−1x n−1 a1′,n x n + K + a ′2,n−1x n−1 a ′2,n x n + K + a ′3,n−1x n−1 a ′3,n x n M O + a ′n−1,n x n x n−1 xn

= c 1′ = c ′2 = c ′3 M = c ′n−1 = c ′n

de donde: xn x n−1 x n−2 M

= c ′n = c ′n−1 = c ′n−2

− a ′n−1,n x n − a ′n−2,n−1x n−1 − a ′n−2,n−1x n M

Generalizando: x i = c ′i −

n

∑ ai, j x j ..............................................................................................(5.40)

j = i +1

Simulación de Yacimientos I

23

5.3.2

Método de Gauss-Jordan

Considérese el sistema de ecuaciones presentado en (5.37). La solución de este sistema por eliminación de Gauss-Jordan consiste en reducirlo a un sistema equivalente que tenga la siguiente forma: x1 x2 x3 O x n−1

xn

= c 1′ = c ′2 = c ′3 M = c ′n−1 = c ′n

Es decir, transformar la matriz aumentada: a11  a 21 a 31   a n1

a12 a 22 a 32 M a n2

a13 a 23 a 33 a n3

L a1,n−1 L a 2,n−1 L a 3,n−1 M L a n,n−1

a1,n c 1   a 2,n c 2  a 3,n c 3  M  a n,n c n 

en una matriz aumentada que tenga la siguiente forma: 1  1   1  O   1

c1  c 2  c3   M  c n 

Para lograr este objetivo se debe proceder en forma análoga al caso anterior. La diferencia radica en que en este caso no se trata de obtener una matriz triangular superior, sino una matriz diagonal unitaria.

El procedimiento o algoritmo de

solución es igual al anterior, excepto que en el paso “c” (reducción), la ecuación se aplica para i ≠ k en lugar de i = k + 1, k + 2 , k + 3 , K

Simulación de Yacimientos I

24

5.4

METODOS ITERATIVOS DE SOLUCIÓN DE ECUACIONES.

Los métodos iterativos de solución de un sistema de ecuaciones no requieren de una memoria de computador de tan alta capacidad como los métodos directos; sin embargo, requieren de un alto número de cálculos. En general, en los métodos iterativos se parte de una serie de valores (presiones) asumidos y mediante una serie de iteraciones se obtienen los valores verdaderos. Algunos de los métodos iterativos mas comúnmente aplicados en la solución de un sistema de ecuaciones son: Método de Jacobi, Método de Gauss-Seidel, Método PSOR, Método LSOR y el Método LSORC.

5.4.1 El Método de Jacobi

Considérese el siguiente sistema de ecuaciones: a11x 1 + a12 x 2 + a13 x 3 + K + a1n x n = b1 ............................................................(5.41) a 21x 1 + a 22 x 2 + a 23 x 3 + K + a 2n x n = b 2 ..........................................................(5.42) a 31x 1 + a 32 x 2 + a 33 x 3 + K + a 3n x n = b 3 ..........................................................(5.43) M

M

a n1x 1 + a n2 x 2 + a n3 x 3 + K + a nn x n = b n ...........................................................(5.44)

De la Ecuación 5.41: x1 =

b1 − a12 x 2 − a13 x 3 − K − a1n x n ................................................................(5.45) a11

Simulación de Yacimientos I

25

De la Ecuación 5.42: x2 =

b 2 − a 21x 1 − a 23 x 3 − K − a 2n x n ...............................................................(5.46) a 22

De la Ecuación 5.43: x3 =

b 3 − a 31x 1 − a 32 x 2 − K − a 3n x n ...............................................................(5.47) a 33

M De la Ecuación 5.44:

xn =

b n − a n1x 1 − a n2 x 2 − a n3 x 3 − K − a n,n−1x n−1 a nn

............................................(5.48)

Las Ecuaciones 5.45 a 5.48 pueden ser escritas en forma generalizada de la siguiente forma: n

b i − ∑ a i, j x j xi =

i =1 j≠i

a i,i

..............................................................................................(5.49)

La Ecuación 5.49 constituye la ecuación básica del método de Jacobi.

El

procedimiento de aplicación es el siguiente: a. Se asumen valores iniciales de x 1 , x 2 , K , x n . b. Se evalúa x 1 , x 2 , K , x n de las Ecuaciones 5.45 a 5.48 y teniendo en cuenta los datos asumidos en el literal anterior. c. Si los valores encontrados en el paso “b” coinciden con los valores asumidos en el paso “a”, entonces los valores de x 1 , x 2 , K , x n son los correctos. Si

Simulación de Yacimientos I

26

no coinciden, entonces se asumen valores de x 1 , x 2 , K , x n iguales a los últimos calculados y se repite el proceso iterativo.

El método de Jacobi puede ser aplicado a la solución del sistema de ecuaciones generado por la Ecuación 5.22, reorganizándola de la siguiente forma: Pi,(kj + 1) =

[

]

1 fi, j − a i, jPi(−k1), j − b i, jPi,(kj −)1 − di, jPi,(kj +) 1 − e i, jPi(+k1), j .....................................(5.50) c i, j

donde el índice k

y k + 1 corresponde a las iteraciones k

y

k + 1,

respectivamente.

5.4.2 El Método de Gauss-Seidel

El método de Gauss-Seidel difiere del método de Jacobi en que los valores asignados a cada incógnita corresponde a los últimos valores disponibles y no a los valores de la última iteración.

El procedimiento del método puede resumirse de la siguiente forma: a. Inicialmente, se asumen valores para las incógnitas x 1 , x 2 , x 3 , K , x n . b. Se evalúa x 1 de la Ecuación 5.45. c. Se evalúa x 2 de la Ecuación 5.46, pero el valor de x 1 utilizado será el estimado en el paso anterior.

Simulación de Yacimientos I

27

d. Se evalúa x 3 de la Ecuación 5.47. Los valores de x 1 y x 2 tenidos en cuenta en este cálculo corresponden a los valores calculados en los dos pasos anteriores. e. Se repiten los tres últimos pasos hasta evaluar x n de la Ecuación 5.48. f. Si los valores de x 1 , x 2 , x 3 , K , x n calculados en esta iteración coinciden con los calculados en la iteración anterior, entonces estos serán los verdaderos valores; de lo contrario, se repite el anterior procedimiento hasta encontrar convergencia.

El método de Gauss-Seidel puede ser aplicado a la solución del sistema de ecuaciones generado por la Ecuación 5.22. Para ello, esta ecuación debe ser reorganizada de la siguiente forma:

Pi,(kj + 1) =

[

]

1 fi, j − a i, jPi(−k1+, j1) − b i, jPi,(kj −+11) − di, jPi,(kj +) 1 − e i, jPi(+k1), j ..................................(5.51) c i, j

donde k y k + 1 hacen referencia a las iteraciones k y k + 1, respectivamente. Obsérvese que la Ecuación 5.51 asume que el ordenamiento es normal.

5.4.3

Método PSOR

El método de Gauss-Seidel no hace uso del valor de Pi,(kj ) para el cálculo de Pi,(kj +1) (veáse Ecuación 5.51). El método PSOR ("Point Successive Over Relaxation") ofrece una forma de tener en cuenta el valor de Pi,(kj ) en el cálculo de Pi,(kj +1) .

Simulación de Yacimientos I

28

Supóngase que el valor de Pi,(kj +1) en la Ecuación 5.51 se nota como Pi,*j(k +1) , de tal forma que: Pi,*j(k + 1) =

[

]

1 fi, j − a i, jPi(−k1+, j1) − b i, jPi,(kj −+11) − di, jPi,(kj +) 1 − e i, jPi(+k1), j ................................(5.52) c i, j

El método PSOR se fundamenta en calcular un valor ponderado de Pi,(kj +1) promediando el valor de Pi,*j(k +1) y el valor de Pi,(kj ) , de la siguiente forma: Pi,(kj +1) = (1 − ω) Pi,(kj ) + ω Pi,*j (k +1) ..............................................................................(5.53) donde el factor ω se denomina parámetro de relajación.

Este método se

denomina, en general, SOBRE-RELAJACION; con su aplicación se busca utilizar un valor óptimo de ω de tal forma que el proceso de convergencia sea acelerado. Obsérvese que este método es igual al anterior en el caso particular en que ω = 1.

Para que exista convergencia del método PSOR (como de otros métodos SOR), se debe cumplir las siguientes condiciones: a. Que la matriz de coeficientes sea diagonalmente dominante. b. Que el valor de ω sea menor que 2. Se dice que una matriz es diagonalmente dominante cuando el valor absoluto de la diagonal principal es mayor o igual que la suma de los valores absolutos de los demás coeficientes en la misma fila. En simulación de yacimientos este caso se cumple, siempre y cuando las ecuaciones sean bien planteadas.

Simulación de Yacimientos I

29

Concepto de Factor de Reducción, ρ .

Se define el residuo de una iteración k , ri,(kj ) , como la diferencia entre el valor exacto, Pi, j , y el valor calculado en la iteración k , Pi,(kj ) . Es decir: ri,(kj ) = Pi, j − Pi,(kj ) ..................................................................................................(5.54) Si se grafica el logaritmo del máximo valor de ri,(kj ) , obtenido en un iteración k , como función del número de iteraciones, k ,

se obtendrá un gráfico de

comportamiento similar al ilustrado en la Figura 5.14.

ω = Valor fijo

Log ri,(kj )

Convergencia asintótica máx

K (No de Iteraciones)

En un principio, los residuos se reducen muy rápidamente. Posteriormente, la rata de reducción es lenta, tal como se presenta en la región recta de la Figura 5.14. Esta región se denomina "región de convergencia asintótica". Para esta región se define el FACTOR DE REDUCCION como:

ri,(kj + 1) ρ = (k ) ..........................................................................................................(5.55) r i, j

Si ρ < 1 , el proceso converge.

Simulación de Yacimientos I

30

El valor de ρ depende del valor de ω ; en consecuencia, la velocidad de convergencia depende del valor de “ ω ” utilizado. Para la región de convergencia asintótica, se ha encontrado que la variación de ρ en función de ω es como se ilustra en la Figura 5.15.

ρ1

ρ

0 1

ω

2

Figura 5.15 Factor de reducción en función del parámetro de relajación

El valor de ωopt es el valor de “ ω ” al cual el factor de reducción, ρ , es mínimo. El valor de ωopt está dado por:

ωopt =

2 1 + 1 − ρ1

..............................................................................................(5.56)

donde ρ1 es el valor de ρ cuando ω = ω1 y se puede obtener graficando log ri,j

max

en función de k , tal como se ilustra en la Figura 5.14.

Si 1 ≤ ω ≤ ωopt , ρ puede ser calculado de la siguiente ecuación: 2

ρ=

1  ρ + ω − 1   ...........................................................................................(5.57) ρ1  ω 

Simulación de Yacimientos I

31

Para ω ≥ ωopt ,

ρ = ω − 1 ...........................................................................................................(5.58)

) y que el valor Supóngase que el valor exacto de Pi, j en el bloque (i, j) es Pi,(Exacto j calculado en la iteración k es Pi,(kj ) . Una curva de en función de k presenta un comportamiento como el que se ilustra en la Figura 5.16 Si ω < ωopt , la curva es

). monotónica. Si ω > ωopt , la curva es oscilatoria en torno al eje Pi,(kj ) = Pi,(Exacto j Este comportamiento es útil para encontrar el valor de ωopt mediante un proceso de ensayo y error.

Valor exacto

Valor exacto

Pi (, kj )

Pi (, kj )

ω < ωopt

ω > ωopt

k

Figura 5.16

k

Presión calculada en el bloque (i,j) en función del número de

iteraciones k.

Si como criterio de convergencia se fija un valor de ri, j = ri, j ω es aquel para el cual el valor de ri, j

conv

conv

, el mejor valor de

se alcanza con el menor número de

iteraciones. Este valor puede ser menor, igual o mayor que ωopt , tal como, por ejemplo, se ilustra en la Figura 5.17.

Simulación de Yacimientos I

32

w > w opt Log rij(k )

max

w = w opt

w < w opt

Criterio de Convergencia

k (No. de Iteraciones)

Figura 5.17 Logaritmo del residuo máximo en una iteración en función del número de iteraciones.

El siguiente procedimiento indica como obtener el mejor valor de ω : a. Se fija un valor de ω . b. Se mide el tiempo requerido para resolver el sistema de ecuaciones hasta alcanzar convergencia. c. Se repiten los pasos “a” y “b” para otros valores de ω . d. Se grafica t en función de ω (Figura 5.18). e. Se selecciona el mejor valor de ω como aquel en el cual el tiempo al que se

Tiempo Requerido para Convergencia

hace referencia en el literal “b” es mínimo.

w mejor

w opt

Factor de Relajacion, w

Figura 5.18 Selección del mejor valor de ω de acuerdo al criterio de convergencia

Simulación de Yacimientos I

33

5.4.4 El Método LSOR

El método LSOR ("Line Successive Over Relaxation"), o de sobre-relajación por líneas, consiste en determinar simultáneamente los valores, en la iteración k + 1, para todas las incógnitas de una misma línea. Si se lleva el valor de Pi,*j(k +1) de la Ecuación 5.52 a la Ecuación 5.53, se obtiene:

[

1 Pi,(kj + 1) = (1 − ω)Pi,(kj ) + ω ⋅ fi, j − a i, jPi(−k1+, j1) − b i, jPi,(kj −+11) − di, jPi,(kj +) 1 − e i, jPi(+k1), j c i, j

]

fi, j a i, j (k + 1) b i, j (k +1) di, j (k ) e i, j (k ) Pi,(kj +1) = (1 − ω)Pi,(kj ) + ω ⋅ − ω⋅ Pi −1, j − ω ⋅ Pi, j −1 − ω ⋅ Pi, j +1 − ω ⋅ P c i, j c i, j c i, j c i, j c i, j i +1, j .........................................................................................................................(5.59) Considérese una malla de un sistema bidimensional la cual se recorre por filas, tal como se ilustra en la Figura 5.19. En el momento de realizar los cálculos para la fila j , las presiones de los bloques ubicados en las filas 1 a j − 1, corresponden a los valores obtenidos en la iteración k + 1. Así mismo, las presiones asignadas a las filas j + 1 , K , n , corresponden a los valores de presión obtenidos en la iteración k . n Iteración K j+1

j

Pi −1, j

Pi, j

Pi+1, j

j-1 Iteración K+1 1 1

i-1

i

i+1

Figura 5.19 Esquema del método LSOR

Simulación de Yacimientos I

34

La Ecuación 5.59 puede ser reorganizada de la siguiente forma:

ω⋅

a i, j (k +1) e i, j (k ) fi, j b i, j (k +1) di, j (k ) Pi −1, j + Pi,(kj +1) + ω ⋅ Pi + 1, j = (1 − ω)Pi,(kj ) + ω ⋅ − ω⋅ Pi, j −1 − ω ⋅ P c i, j c i, j c i, j c i, j c i, j i, j +1

.........................................................................................................................(5.60) Como se mencionó, el método LSOR consiste en resolver simultáneamente los valores de el sistema de ecuaciones para todas las incógnitas de una misma fila en la iteración (k + 1) . Por esta razón, Pi(+k1), j

puede ser sustituido por Pi(+k1+, j1) , en la

Ecuación 5.60:

ω⋅

a i, j (k +1) e i, j (k +1) fi, j b i, j (k +1) di, j (k ) Pi −1, j + Pi,(kj +1) + ω ⋅ Pi +1, j = (1 − ω)Pi,(kj ) + ω ⋅ − ω⋅ Pi, j −1 − ω ⋅ P c i, j c i, j c i, j c i, j c i, j i, j +1

.........................................................................................................................(5.61) Los términos del lado derecho en la Ecuación 5.61 son conocidos de las filas ya recorridas en la iteración (k + 1) o de las filas que aún faltan por recorrer en la iteración (k + 1) y que, por tanto, almacenan valores obtenidos en la iteración anterior, k . Las presiones del lado izquierdo de la Ecuación 5.61 corresponden a los bloques de la fila j , y representan las incógnitas cuyos valores deben ser estimados. Si la Ecuación 4.61 se aplica a cada uno de los bloques de la fila j , j fijo, resulta un sistema tridiagonal de ecuaciones. Su solución se lleva a cabo mediante la aplicación del algoritmo de Thomas.

Una vez se resuelve este sistema de

ecuaciones para la línea j , se continúa con la línea j + 1 , y se repite el procedimiento. El proceso se repite hasta barrer toda la malla. El método LSOR puede ser aplicado recorriendo la malla por columnas, en lugar de filas. En este caso, la ecuación resultante toma la siguiente forma:

Simulación de Yacimientos I

35

ω⋅

b i, j (k +1) di, j (k +1) fi, j a i, j (k +1) e i, j (k ) Pi, j −1 + Pi,(kj +1) + ω ⋅ Pi, j +1 = (1 − ω)Pi,(kj ) + ω ⋅ − ω⋅ Pi −1, j − ω ⋅ P c i, j c i, j c i, j c i, j c i, j i +1, j

.........................................................................................................................(5.62) El parámetro de sobre-relajación, ω , se determina de igual manera que en el método PSOR. La convergencia del método se logra cuando se cumple que: Pi,(kj +1) − Pi,(kj ) < ri, j

converg

...................................................................................(5.63)

para todos los puntos de la malla.

5.4.5 Método LSORC o WATTS LSOR

Este método representa una versión mejorada del método LSOR. De la Ecuación 5.22 se tiene: a i, jPi −1, j + bi, jPi, j −1 + c i, jPi, j + di, jPi, j +1 + e i, jPi +1, j = fi, j ...........................................(5.22)

Supóngase que la malla será recorrida por columnas, aplicando el algoritmo del método LSOR. De esta forma, cada línea i se resuelve simultáneamente (Figura 5.18). J

→ → → → → I

Figura 5.20 Esquema del recorrido por columnas en los métodos LSOR y LSORC

Simulación de Yacimientos I

36

Después de la iteración k , la Ecuación 5.22 puede ser escrita como: a i, jPi(−k1), j + b i, jPi,(kj −)1 + c i, jPi,(kj ) + di, jPi,(kj +)1 + e i, jPi(+k1), j − fi, j = ri,(kj ) .................................(5.64) donde el vector ri, j es el residuo obtenido al final de la iteración k . Si la solución alcanzada es la exacta; ri,(kj ) = 0 en particular, para cada línea se tiene: J

∑ ri,(kj ) = 0

(si la solución es exacta) .........................................................(5.65)

j =1

El método LSORC consiste en sumar, después de la iteración k , a las presiones de cada columna de la malla un valor de corrección de tal forma que se cumpla la Ecuación 5.65. De esta forma, una vez introducida la corrección, se obtendrá la siguiente ecuación:

(

)

(

)

(

)

(

)

(

)

a i, j Pi(−k1), j + Pi*−1 + b i, j Pi,(kj −)1 + Pi* + c i, j Pi,(kj ) + Pi* + di, j Pi,(kj +) 1 + Pi* + e i, j Pi(+k1), j + Pi*+ 1 − fi, j = 0 ........................................................................................................................(5.66) Si se escribe la Ecuación 5.64 para cada uno de los bloques de la columna i y luego se suman las ecuaciones resultantes, se tiene: J

J

J

J

J

J

J

j =1

j =1

j =1

j =1

j =1

j =1

j =1

∑ ai, jPi(−k1), j + ∑ bi, jPi,(kj −)1 + ∑ c i, jPi,(kj ) + ∑ di, jPi,(kj +)1 + ∑ ei, jPi(+k1), j − ∑ fi, j = ∑ ri,(kj ) ....(5.67) Si se hace algo análogo con la Ecuación 5.66, se tiene: J

J

J

J

J

J

J

J

j =1

j =1

j =1

j =1

j =1

∑ ai, jPi(−k1), j + ∑ ai, jPi*−1 + ∑ bi, jPi,(kj −)1 + ∑ bi, jPi* + ∑ c i, jPi,(kj ) + ∑ c i, jPi* + ∑ di, jPi,(kj +)1 + ∑ di, jPi* j =1

j =1

j =1

J

J

J

j =1

j =1

j =1

+ ∑ e i, jPi(+k1), j + ∑ e i, jPi*+ 1 − ∑ fi, j = 0 .........................................................................................................................(5.68)

Simulación de Yacimientos I

37

Restando la Ecuación 5.67 de la Ecuación 5.68, se tiene: J

J

J

J

J

J

j =1

j =1

j =1

j =1

j =1

j =1

∑ ai, jPi*−1 + ∑ bi, jPi* + ∑ c i, jPi* + ∑ di, jPi* + ∑ ei, jPi*+1 = −∑ ri,(kj ) o bien: J J J  J   J   J   a P * +  b + c + d P * +  e P * = − r (k ) .....................(5.69) ∑ i, j i, j i, j ∑ i, j ∑ i, j i i, j ∑  i −1  ∑  ∑  i +1 j =1 j =1 j =1  j =1   j =1   j =1 

La Ecuación 5.69 genera un sistema tridiagonal de ecuaciones que permite evaluar la corrección P * para cada fila i , 1 ≤ i ≤ I . El procedimiento para aplicar el método es el siguiente: a. Se aplica el método LSOR para barrer la malla, tal como se discutió en el numeral 5.4.4. b. Se calcula Pi* , 1 ≤ i ≤ I , resolviendo el sistema tridiagonal de ecuaciones generado por la Ecuación 5.69. c. Se obtienen nuevos valores para Pi,(kj ) (es decir, los valores corregidos), así: Pi,(kj ) = Pi,(kj ) + Pi* ...........................................................................................(5.70) Obsérvese que para cada columna i , Pi* es constante. d. Si se cumple el criterio de convergencia para todos los puntos de la malla, se para el proceso. Si no se cumple, se repiten los pasos a, b y c. El método LSORC puede ser aplicado barriendo la malla por filas, en lugar de columnas. El análisis es similar.

Simulación de Yacimientos I

38

1

6. FLUJO LINEAL DE UN FLUIDO LEVEMENTE COMPRESIBLE

La ecuación de continuidad para flujo lineal de un fluido levemente compresible puede ser escrita como: −

∂ (A (x )ρUx ) = A (x ) ∂(φρ) + A (x )~q ∂x ∂t

o bien,  ∂ (A (x )Ux ) + (A (x )Ux ) ∂ρ  = A (x ) ∂(φρ) + A (x )~q ........................................ (6.1) − ρ ∂x  ∂t  ∂x Donde: ~ q : Cantidad de masa que entra o sale, por fuentes o sumideros, por unidad de volumen de yacimiento por unidad de tiempo. De la ley de Darcy,

Ux = −

k x ∂p .................................................................................................. (6.2) ⋅ µ ∂x

La ecuación de estado para un fluido levemente compresible está dada por la siguiente expresión:

ρ = ρ oe c (p −p o ) donde c es la compresibilidad del fluido y su densidad a una presión po . De esta última ecuación, se puede escribir:

Simulación de Yacimientos I

1

2

∂ρ ∂p ∂p ..................................................................... (6.3) = ρoec (P −Po ) ⋅ c ⋅ = ρ⋅c ⋅ ∂x ∂x ∂x ∂ρ ∂p ∂p ..................................................................... (6.4) = ρoec (P −Po ) ⋅ c ⋅ = ρ⋅c ⋅ ∂t ∂t ∂t

Sustituyendo las Ecuaciones 6.2 a 6.4 en la Ecuación 6.1 y asumiendo que la porosidad es constante con el tiempo, se obtiene: 2

k ∂p  k  ∂p  ∂p ∂  + A (x ) ~ q ρ  A (x ) x ⋅  + A (x ) x ⋅   ρ ⋅ c = A (x ) ⋅ φ ⋅ ρ ⋅ c ∂t µ  ∂x  ∂x  µ ∂x 

Debido a que el valor de compresiblidad para un fluido levemente compresible es pequeño

(del orden de 10 −5 lpc −1 ), el segundo término del lado izquierdo se

puede considerar despreciable si se compara con el primer término, lo que simplifica la ecuación de la siguiente forma: ρ

k ∂p  ∂p ∂   A (x ) x ⋅  = A (x ) ⋅ φ ⋅ ρ ⋅ c + A (x )~ q ∂t ∂x  µ ∂x 

Dividiendo por ρ

:

~ ∂  k ∂p  q ∂p  A (x ) x ⋅  = A (x ) ⋅ φ ⋅ c + A (x ) ........................................................ (6.5) ∂x  µ ∂x  ∂t ρ Si se define qv como: qv =

~ q ρ

Simulación de Yacimientos I

2

3

se obtiene: ∂  k ∂p  ∂p  A (x ) x ⋅  = A (x ) ⋅ φ ⋅ c + A (x )qv ........................................................ (6.6) ∂x  µ ∂x  ∂t

De la aproximación numérica para flujo lineal de un fluido incompresible se tiene: ∂  kA ∂p  Ti+ 12 (Pi+1 − Pi ) − Ti− 12 (Pi − Pi−1) .................................................. (6.7)  ⋅  ≅ ∂x  µ ∂x i ∆xi La derivada

∂p puede aproximarse a un esquema progresivo: ∂t

∂p Pin +1 − Pin ................................................................................................ (6.8) = ∂t ∆t

Llevando las Ecuaciones 6.7 y 6.8 a la Ecuación 6.6, se obtiene: Ti+ 1 (Pi+1 − Pi ) − Ti− 1 (Pi − Pi−1) 2

2

∆xi

 Pn +1 − Pin  = A i ⋅ φi ⋅ c i  i  + A i ⋅ qv i .......................... (6.9) ∆ t  

En el lado izquierdo de la Ecuación 6.9 no se especifica el nivel de tiempo al que hacen referencia los términos de presión. Dependiendo de este nivel de tiempo, se suele hablar de diferentes esquemas de aproximación tales como la aproximación explícita y la aproximación implícita.

6.1 LA APROXIMACION EXPLÍCITA

Si las presiones del lado izquierdo de la Ecuación 6.9 se evalúan al tiempo t n se habla de un esquema de aproximación explícita. En este caso:

Simulación de Yacimientos I

3

4

(

)

(

Ti+ 1 Pin+1 − Pin − Ti− 1 Pin − Pin−1 2

2

∆x i

)

 Pin +1 − Pin  = A i ⋅ φi ⋅ c i   + A i ⋅ qv i ∆t  

o bien:

(

)

 Pn +1 − Pin  Ti+ 1 Pin+1 − Ti− 1 + Ti+ 1 Pin + Ti− 1 Pin−1 = Vi ⋅ φi ⋅ c i  i  + Q v i ................... (6.10) 2 2 2 2 ∆t   donde Q v i es la tasa de producción o inyección asignada al bloque i .

Si se define Ti como: Ti = Vi ⋅ φi ⋅ c i ................................................................................................... (6.11)

Entonces:

(

)

[

2

2

2

]

Ti n +1 Pi − Pin + Q v i ............................. (6.12) ∆t

Ti+ 1 Pin+1 − Ti− 1 + Ti+ 1 Pin + Ti− 1 Pin−1 = 2

La Ecuación 6.12 puede ser resuelta explícitamente para Pin +1 , Pin +1 = Pin +

[

(

)

]

∆t T 1 Pin+1 − Ti− 1 + Ti+ 1 Pin + Ti− 1 Pin−1 − Qv i ............................. (6.13) 2 2 2 Ti i+ 2

Es decir, este esquema permite el cálculo de las presiones al tiempo t n +1 conocidas las presiones al tiempo t n .

Simulación de Yacimientos I

4

5

El análisis de error de truncamiento, estabilidad, convergencia y consistencia de este esquema puede ser llevado a cabo teniendo en cuenta los criterios discutidos anteriormente.

Por ejemplo, la estabilidad puede ser analizada de acuerdo al

criterio de Karplus; la Ecuación 6.8 puede ser escrita como:

(

)

(

)

Ti+ 1 Pin+1 − Pin + Ti− 1 Pin−1 − Pin − 2

2

[

]

Ti n +1 Pi − Pin = 0 ∆t

La condición de estabilidad será: Ti + 1 + Ti− 1 − 2

2

Ti ≤0 ∆t

o bien: Ti ≥ Ti + 1 + Ti − 1 ........................................................................................... (6.14) 2 2 ∆t

De donde se infiere que el ESQUEMA EXPLÍCITO es CONDICIONALMENTE ESTABLE y la condición de estabilidad está dada por la Ecuación 6.14.

6.2 APROXIMACION IMPLÍCITA

Si las presiones del lado izquierdo de la Ecuación 6.9 son evaluadas al tiempo t n +1 , se habla de que el esquema es IMPLÍCITO. En este caso se tiene:

(

)

(

)

Ti+ 1 Pin++11 − Pin +1 − Ti− 1 Pin +1 − Pin−+11 = 2

2

[

]

Vi ⋅ φi ⋅ c i n +1 Pi − Pin + Q v i ∆t

de donde:

Simulación de Yacimientos I

5

6

T T  Ti+ 1 Pin++11 −  Ti− 1 + Ti+ 1 + i Pin+1 + Ti− 1 Pin−+11 = − i Pin + Qv i ......................... (6.15) 2 2 2 2 ∆t  ∆t  La Ecuación 6.15 genera un sistema tridiagonal de ecuaciones, en el cual las incógnitas son: Pin++11 , Pin +1 , y Pin−+11 .

El análisis de error de truncamiento. estabilidad, convergencia y consistencia de la Ecuación

6.11

previamente.

puede

ser

realizado

siguiendo los criterios presentados

Por ejemplo, aplicando el criterio de Karplus para analizar

estabilidad, de la Ecuación 6.15 se obtiene:

(

)

(

)

(

)

(

)

Ti+ 1 Pin++11 − Pin − Ti+ 1 Pin +1 − Pin − Ti− 1 Pin +1 − Pin + Ti− 1 Pin−+11 − Pin − 2

2

2

2

[

]

Ti n +1 Pi − Pin = 0 ∆t

La condición de estabilidad es: Ti + 1 − Ti + 1 − Ti − 1 + Ti − 1 − 2

2

2

2

Ti ≤0 ∆t

o bien: −

Ti ≤0 ∆t

Lo que siempre se cumple; luego, el esquema IMPLÍCITO presentado en la Ecuación 6.15 es INCONDICIONALMENTE ESTABLE

Simulación de Yacimientos I

6

7. FLUJO BIDIMENSIONAL DE UN FLUIDO LEVEMENTE COMPRESIBLE

La ecuación de continuidad para flujo bidimensional, coordenadas cartesianas, está dada por la siguiente expresión:

∂  k x ∂p  ∂  k y ∂p  ∂ (ρφ ) =α  αρ  +  αρ + α~ q ............................................. (7.1)  ∂x  µ ∂x  ∂y  µ ∂y  ∂t donde: α : Al espesor del sistema, H(x, y ) . ~

q

: Masa de fluido que entra o sale por fuentes o sumideros por unidad de

volumen del yacimiento por unidad de tiempo. La densidad, ρ , para un fluido levemente compresible está dada por la siguiente expresión: ρ = ρ oe c (p −p o ) ............................................................................................. (7.2) Llevando la Ecuación 7.2 a la Ecuación 7.1, considerando la porosidad constante, se obtiene: ~ ∂  k ∂P  ∂  k ∂P  ∂ρ  Hρ  = Hφ  Hρ  + + Hq ∂x  µ ∂x  ∂y  ∂t µ ∂y 

de donde:

ρ

~ ∂  kH ∂P  kH ∂P ∂ρ ∂  kH ∂P  kH ∂P ∂ρ ∂ρ  +   + + ρ  = Hφ + H q ............ (7.3) ∂x  µ ∂x  µ ∂x ∂x ∂y  µ ∂y  µ ∂y ∂y ∂t

Simulación de Yacimientos I

1

Aplicando la regla de la cadena: ∂ρ ∂ρ ∂p ∂p ∂p ........................................................... (7.4) = = ρ0ec (P −P0 )c = ρc ∂x ∂p ∂x ∂x ∂x Análogamente: ∂P ∂ρ ................................................................................................... (7.5) = ρc ∂y ∂y ∂ρ ∂P ................................................................................................... (7.6) = ρc ∂t ∂t

Llevando las Ecuaciones 7.4 a 7.6 a la Ecuación 7.3, se tiene: 2

2

~ ∂  kH ∂P  kH  ∂P  ∂  kH ∂P  kH  ∂P  ∂P   +  = Hφρc   + ρ ρc ρ c  +Hq  +ρ ∂x  µ ∂x  µ ∂y  µ ∂y  µ ∂t  ∂x   ∂y 

Al igual que para el caso lineal, debido a que para un fluido levemente compresible la compresibilidad es muy pequeña, se puede considerar que los términos que involucran gradientes al cuadrado tienden a cero, en cuyo caso se tiene:

ρ

~ ∂  kH ∂P  ∂  kH ∂P  ∂P   = Hφρc   + ρ +Hq ∂x  µ ∂x  ∂y  µ ∂y  ∂t

o bien: ∂  kH ∂P  ∂  kH ∂P  ∂P   = Hφc   + + Hqv .................................................... (7.7) ∂x  µ ∂x  ∂y  µ ∂y  ∂t

Simulación de Yacimientos I

2

Donde qV representa el volumen de fluido que entra o sale por fuentes o sumideros por unidad de volumen del yacimiento por unidad de tiempo.

De la Ecuaciones 5.18 y 5.19 se tiene que:

(

)

(

)

(

)

(

)

∂  k xH ∂p  Ti+ 12, j Pi+1, j − Pi, j − Ti− 12, j Pi, j − Pi−1, j ........................................ (7.8)  ⋅ ≅ ∂x  µ ∂x  ∆y j∆xi ∂  k xH ∂p  Ti, j+ 12 Pi, j+1 − Pi, j − Ti, j− 12 Pi, j − Pi, j−1 ........................................ (7.9)  ⋅ ≅ ∂y  µ ∂y  ∆y j∆xi La derivada ∂P

∂t

puede ser expresada en forma progresiva, asi:

n+1 − Pin ∂P  Pi ...................................................................................... (7.10) = ∂ t  ∆t

Llevando las Ecuaciones 7.8 a 7.10 Ecuación 7.7 se obtiene:

(

)

(

Ti+ 1 , j Pi+1, j − Pi, j − Ti− 1 , j Pi, j − Pi−1, j 2

2

∆y j∆xi

Hi, jφi, jc i, j

Pin, j + 1 − Pin, j ∆t

)

+

(

)

(

Ti, j+ 1 Pi, j+1 − Pi, j − Ti, j− 1 Pi, j − Pi, j−1 2

∆y j∆xi

2

)

=

+ Hi, jqv i, j

Si se define:

Vbij = ∆xi∆yiHi, j ............................................................................................ (7.11) γ ij = Vbijφijc i, j ................................................................................................ (7.12) Qvij = ∆xi∆yiHi, jqvij ...................................................................................... (7.13)

Simulación de Yacimientos I

3

se obtiene:

(

Ti, j− 1 Pi, j−1 + Ti− 1 , jPi−1, j − Ti− 1 2

(P ∆t γ ij

n +1 − Pijn ij

2

)+ Q

Vij

2, j

)

+ Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j Pi, j + Ti+ 1 , jPi+1, j + Ti, j+ 1 Pi, j+1 = 2

2

2

2

2

(7.14)

7.1 APROXIMACIÓN EXPLÍCITA

La Ecuación 7.14 puede ser escrita en forma de aproximación explícita evaluando las presiones del lado izquierdo de la ecuación al tiempo t n :

(

Ti, j− 1 Pi,nj−1 + Ti− 1 , jPin−1, j − Ti− 1 2

(P ∆t γij

n +1 n i, j − Pi, j

2

)+ Q

2, j

)

+ Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j Pin, j + Ti+ 1 , jPin+1, j + Ti, j+ 1 Pin, j+1 = 2

2

2

2

2

Vi, j

o bien: Pin, j+1 = Pin, j +

[

(

)

∆t Ti, j− 1 Pi,nj−1 + Ti− 1 , jPin−1, j − Ti− 1 , j + Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j Pin, j + Ti+ 1 , jPin+1, j + Ti, j+ 1 Pin, j+1 − Qvij 2 2 2 2 2 2 2 2 γ i, j ................................................................................................................... (7.15) La Euación (7.15 puede ser resuelta para Pin, j+1 , conociendo las presiones al tiempo t n . El análisis de estabilidad, convergencia, consistencia y error de truncamiento para la Ecuación 7.15 puede efectuarse siguiendo los criterios discutidos en el capítulo 3. Por ejemplo, para aplicar el criterio de Karplus, de la Ecuación 7.15, se tiene:

Simulación de Yacimientos I

4

]

( γ − (P ∆t

) − P )− Q

(

)

(

)

(

Ti, j− 1 Pin, j−1 − Pin, j + Ti− 1 , j Pin−1, j − Pin, j + Ti+ 1 , j Pin+1, j − Pin, j + Ti, j+ 1 Pin, j+1 − Pin, j 2

ij

n +1 i, j

n i, j

2

Vi, j

2

2

)

=0

La condición de estabilidad requiere que: Ti, j− 1 + Ti− 1

2, j

2

+ Ti+ 1 , j + Ti, j+ 1 − 2

2

γ ij ∆t

≤ 0 ........................................................ (7.16)

De donde se infiere que la aproximación EXPLICITA de la Ecuación 7.1 es CONDICIONALMENTE ESTABLE y la condición de estabilidad está dada por la Ecuación 7.16.

7.2 APROXIMACIÓN IMPLICITA

La aproximación implícita a la Ecuación 7.1 puede ser obtenida evaluando las presiones de la Ecuación 7.14 al tiempo t n +1 :

(

)

Ti, j− 1 Pi,nj+−11 + Ti − 1 , jPin−+1,1j − Ti − 1 , j + Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j Pin, j+1 + Ti+ 1 , jPin++1,1j + Ti, j+ 1 Pin, j++11 = 2

2

(P ∆t γij

n +1 n i, j − Pi, j

)+ Q

Vi, j

2

2

2

2

2

2

(7.17)

Si se define: Sij = Ti, j− 1 ................................................................................................. (7.18) 2

Wij = Ti−1 2, j ................................................................................................ (7.19) γij   Cij = − Ti+1 2, j + Ti−1 2, j + Ti, j−1 2 + Ti, j+1 2 +  ............................................... (7.20) ∆t  

Simulación de Yacimientos I

5

Eij = Ti+1 2, j .................................................................................................. (7.21) Nij = Ti, j+1 2 ................................................................................................. (7.22) Fij = −

γij ∆t

Pijn + Q Vij ..................................................................................... (7.23)

La Ecuación 7.17 puede ser escrita como:

Si, jPin, j+−11 + Wi, jPin−+1,1j + Ci, jPin, j+1 + Ei, jPin++1,1j + Ni, jPin, j++11 = Fi, j .................................... (7.24)

La Ecuación 7.24 genera un sistema pentadiagonal de ecuaciones que permiten el cálculo de las presiones al tiempo t n +1 , conocidas las presiones al tiempo t n . La solución de este sistema de ecuaciones se puede llevar a cabo mediante métodos directos o iterativos, tales como los discutidos en el cap¡tulo 5. Al igual que en el caso anterior, el análisis de estabilidad, convergencia, consistencia y error de truncamiento de la Ecuación 7.24 puede realizarse aplicando los criterios presentados en el capítulo 4. Se deja como ejercicio al lector demostrar que este esquema es INCONDICIONALMENTE ESTABLE.

7.3 APROXIMACIÓN DE CRANK-NICOLSON

La aproximación de Crank-Nicolson, flujo bidimensional, puede ser escrita como: n

1  ∂  kh ∂P  ∂  kh ∂P  1  ∂  kh ∂P  ∂  kh ∂P    +      + +   2  ∂x  µ ∂x  ∂y  µ ∂y  2  ∂x  µ ∂x  ∂y  µ ∂y 

Simulación de Yacimientos I

n +1

= hijφijc ij

∂P   ∂t 

central

+ hqv

6

En diferencias finitas se tiene:

[ [

( (

) )

]

1 Ti, j− 1 Pi,nj−1 + Ti− 1 , jPin−1, j − Ti− 1 , j + Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j Pin, j + Ti+ 1 , jPin+1, j + Ti, j+ 1 Pin, j+1 + 2 2 2 2 2 2 2 2 2 1 Ti, j− 1 Pi,nj+−11 + Ti− 1 , jPin−+1,1j − Ti− 1 , j + Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j Pin, j+1 + Ti+ 1 , jPin++1,1j + Ti, j+ 1 Pi,nj++11 = 2 2 2 2 2 2 2 2 2 γij n+1 Pi, j − Pin, j + Q Vi, j (7.25) ∆t

(

]

)

De donde: γ i, j  n+1   Pi, j + Ti+ 1 , jPin++1,1j + Ti, j− 1 Pin, j+−11 + Ti− 1 , jPin−+1,1j −  Ti− 1 , j + Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j + 2  2 2 2 2 2 2 2 ∆ t   γ i, j  n   Pi, j − Ti, j+ 1 Pin, j++11 = 2Qvi, j − Ti, j− 1 Pin, j−1 − Ti− 1 , jPin−1, j +  Ti− 1 , j + Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j − 2  2 2 2 2 2 2 2 t ∆   Ti+ 1 , jPin+1, j − Ti, j+ 1 Pin, j+1 2

(7.26)

2

Redefiniendo los términos Ci, j y Fi, j del stencil: Ci, j = Ti− 1

2, j

+ Ti, j− 1 + Ti, j+ 1 + Ti+ 1 2

2

2, j

+2

γ i, j ∆t

γ i, j  n   Pi, j − Fi, j = 2Qvi, j − Ti, j− 1 Pin, j−1 − Ti− 1 , jPin−1, j +  Ti− 1 , j + Ti, j− 1 + Ti, j+ 1 + Ti+ 1 , j − 2 2 2 2 2 2 2 ∆t   −Ti + 1 , jPin+1, j − Ti, j+ 1 Pin, j+1 2

2

Entonces la Ecuación 7.26 toma la forma de la Ecuación 7.24. La Ecuación 7.26 genera un sistema pentadiagonal de ecuaciones cuya solución permite el cálculo de las presiones en todos los bloques de la malla al tiempo t n +1 , conocidas las presiones al tiempo t n .

Simulación de Yacimientos I

7

7.4 EL METODO ADIP

El método ADIP ("Altenating Direction Implicit Pressure") soluciona la Ecuación 7.1 en dos pasos, con la finalidad de obtener la distribución de presiones en la malla al tiempo t n +1 : PASO 1:

En el primer paso se barre la malla en la dirección x, tal como se ilustra en la Figura 7.1,

j = Ny . . . . x

. . . j=2 j =1 i =1

i = 2 . . . . . . . . . . i = N x − 1 i = Nx

Figura 7.1 Esquema que ilustra la dirección de barrido en el paso 1 del método ADIP.

de esta forma, se obtiene la distribución de presiones al tiempo t n +1 2 . Para lograr este objetivo, se evalúa la componente

Simulación de Yacimientos I

∂  kh ∂P   al tiempo tn + 1 y la  ∂x  µ ∂x  2

8

∂  kh ∂P    al tiempo t n . Es decir, la Ecuación 7.1 se expande de ∂y  µ ∂y 

componente

acuerdo al siguiente esquema: n+ 1 2

∂  kh ∂P    ∂x  µ ∂x 

n

progresiva

∂  kh ∂P   ∂p    = φHc  + ∂y  µ ∂y   ∂t 

+ Qv

Numéricamente se tiene: n+ 1

(

Ti− 1 , jPi−1, j 2 − Ti− 1 , j + Ti+ 1 2

=

γ ij ∆t 2

(P

n +1 / 2 i, j

2

)

2, j

)P

n+ 1 2 i, j

n+ 1

(

+ Ti+ 1 , jPi+1, j 2 + Ti, j− 1 Pin, j−1 − Ti, j− 1 + Ti, j+ 1 2

2

2

2

)P

n i, j

+ Ti, j+ 1 Pin, j+1 2

− Pin, j + Q Vi, j

Reagrupando términos: 2Vbi, j φi, jc i, j  n + 1  n+ 1 n+ 1  Pi, j 2 + Ti + 1 , jPi +1, j 2 = Ti − 1 , jPi −1, j 2 −  Ti − 1 , j + Ti + 1 , j +  2 2 2 2 ∆t   2Vbi,j φ i, j c i, j  n   Pi, j − T 1 Pin+1, j + Q vi, j ............ (7.27) − Ti, j− 1 Pin, j−1 +  Ti, j− 1 + Ti, j− 1 − i, j+  ∆t 2 2 2 2  

La Ecuación 7.27 genera un sistema tridiagonal de ecuaciones que permite hallar n+ 1 2,

Pi, j

conocido Pin, j , barriendo la malla por filas.

PASO 2:

En el segundo paso, se barre la malla por columnas (véase la Figura 7.2)

Simulación de Yacimientos I

9

j = Ny . . . . . .

y

. j=2 j =1 i=2 . . . . . . . . . .

i =1

i = Ny

Figura 7.2 -

Esquema que

ilustra la dirección de barrido en el paso 2 del método ADIP. n+1

con el propósito de obtener Pi, j

n+ 1 2

, partiendo de los valores de Pi, j

obtenidos

en el paso anterior. Para ello, se aproxima la Ecuación 7.5 de la siguiente forma: ∂  kh ∂P    ∂x  µ ∂x 

n+ 1 2

∂  kh ∂P    + ∂y  µ ∂y 

n+1

 ∂p  = φHc    ∂t 

Centrak ,n+ 1 2

+ Qv

Numéricamente: n+ 1

(

)

n+ 1 2

Ti− 1 , jPi−1, j 2 − Ti− 1 , j + Ti+ 1 ,j Pi, j 2

Ti, j+ 1 Pi,nj++11 = 2

2

2

n+ 1

(

+ Ti+ 1 ,jPi+1, j 2 + Ti,j− 1 Pin, j−+11 − Ti, j− 1 + Ti,j+ 1 2

2

2

2

)P

n+1 i, j

+

γ ij  n+1 n+ 1   Pi, j − Pi, j 2  + Q Vi, j ∆t 2  

∂p (Obsérvese que la derivada   se ha aproximado utilizando un esquema  ∂t  progresivo desde t 1 a t n+1 ). Reagrupando términos: n+ 2

Simulación de Yacimientos I

10

2Vbi, j φ i, j c i, j  n+1  Pi, j + Ti,j+ 1 Pin, j++11 = Ti, j− 1 Pin, j−+11 −  Ti, j− 1 + Ti, j+ 1 +  2 2 2 ∆ t 2   2Vbi, j φi, jc i, j  n + 1  n+ 1 n+ 1 Pi, j 2 − T 1 Pi, j + 1 2 + Q vi, j − Ti − 1 , jPi −1, j 2 +  Ti + 1 , j + Ti − 1 , j − i+ , j  ∆t 2 2 2 2  

(7.28)

La ecuación (7.28) genera un sistema tridiagonal de ecuaciones que permite el cálculo de las presiones al tiempo t n+1 , de los datos de presión al tiempo t 1 n+

2

calculados en el paso anterior.

Si se define: n+ 1 2

Pij

= Pij*

n+ 1

Pi +1, j 2 = Pi*+1, j n+ 1

Pi −1, j 2 = Pi*−1, j

Entonces, los pasos 1 y 2 pueden ser escritos así:

PASO 1:

Simulación de Yacimientos I

11

2Vbi, j φi, jc i, j  ∗   Pi, j + Ti + 1 , jPi∗+1, j = Ti − 1 , jPi∗−1, j −  Ti − 1 , j + Ti + 1 , j +  2 2 2 2 ∆t   2Vbi,j φi, j c i, j  n  Pi, j − T 1 Pin+1, j + Q vi, j .............. (7.29) − Ti,j− 1 Pin, j−1 +  Ti, j− 1 + Ti, j− 1 − i, j+  ∆t 2 2 2 2  

PASO 2:

2Vbi, j φ i, j c i, j  n+1  Pi, j + Ti,j+ 1 Pin, j++11 = Ti, j− 1 Pin, j−+11 −  Ti, j− 1 + Ti, j+ 1 +  2 2 2 ∆t 2   2Vbi, j φ i, j c i,j  ∗  Pi, j − T 1 Pi,∗j+1 + Q vi, j .............. (7.30) − Ti− 1 , jPi∗−1,j +  Ti+ 1 , j + Ti− 1 , j − i+ , j  ∆t 2 2 2 2  

Procedimiento de Solución: PASO 1:

a. Se fija la primera fila, j = 1, y se aplica la Ecuación 7.29 a cada bloque,

generando un sistema tridiagonal de ecuaciones. Se soluciona este sistema de ecuaciones con la finalidad de estimar los valores de Pij* . b. Se repite el paso anterior para otras filas, j = 2,3,4,...,N y . Para cada valor fijo de j , se resuelve el sistema tridiagonal resultante y se obtiene Pij* para los bloques de dicha fila. La figura 7.3 presenta un esquema donde se indica los bloques en los cuales se ha calculado Pij* un instante después de barrer la fila j y un instante antes de barrer la fila j + 1 .

Simulación de Yacimientos I

12

Figura 7.3

j = Ny . . . . . . j=2 j =1

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n*

n*

n*

n*

n*

n*

n*

n*

n*

n*

n*

n*

n*

n*

i = Nx i =1 i=2 . . . . . . . . . Esquema que indica los bloques barridos (*) en un instante

determinado del PASO 1 del método ADIP.

PASO 2:

Una vez barridas todas las filas de la malla, se dispone de todos los valores de Pij* , tal como se ilustra en el esquema de la Figura 7.4. Con esta información, se continúa el procedimiento: c. Se fija la primera columna, i = 1, y se aplica la Ecuación 7.30 a cada bloque, con la finalidad de obtener un sistema tridiagonal de ecuaciones que permita estimar los valores de Pij* para esta columna. d. Se repite el procedimiento para cada una de las demás columnas hasta obtener los valores de Pijn+1 para toda la malla.

Simulación de Yacimientos I

13

Ny . . . . . . .

N

y

j =1

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

*

i =1

Figura 7.4

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

i = Nx

Esquema que ilustra la información disponible (calculada) una vez

ejecutado el paso 1 del método ADIP.

7.5 EL METODO ADEP

Siguiendo una metodología similar al método ADIP, el método ADEP ("Alternating Direction Explicit Pressure") soluciona la Ecuación 7.1 en dos pasos, así: PASO 1:

El paso 1 consiste en aproximar la Ecuación 7.1 de la siguiente forma: n

∂  kh ∂P  ∂  kh ∂P      + ∂x  µ ∂x  ∂y  µ ∂y 

n,*

∂P  = φHc  ∂t 

n,*

+ Hq V ....................................... (7.31)

Donde:

Simulación de Yacimientos I

14

n

∂  kh ∂P    = ∂x  µ ∂x 

  Ti− 1 , jPin−1, j −  Ti+ 1 , j + Ti− 1 , j Pin, j + Ti+ 1 , jPin+1, j 2 2 2  2  ....................... (7.32) ∆x i ∆y j

y ∂  kh ∂P    ∂y  µ ∂x 

n,∗

= Ti, j − 1 Pi∗, j −1 − Ti, j − 1 Pi∗, j − Ti, j + 1 Pin, j + Ti + 1 , jPin, j +1 .................. (7.33) 2

2

2

2

Además:  ∂P     ∂t 

n,∗

=

(

2 Pi∗, j − Pi, j ∆t

) ................................................................................... (7.34)

Llevando las Ecuaciones 7.32, 7.33 y 7.34 a la Ecuación 7.31, se obtiene:

(

)

Ti− 1 , jPin−1, j − Ti− 1 , j + Ti+ 1 ,j Pin, j + Ti+ 1 , jPin+1,j + Ti, j− 1 Pi∗,j−1 − Ti,j+ 1 Pi∗, j − Ti,j− 1 Pin, j + 2

Ti, j + 1 Pin, j +1 = 2

2

(P ∆t 2 γ ij

∗ i, j

2

)

2

2

2

2

− Pin, j + Q Vij

De donde: 2Vbij φ ij c  n   Pi, j + Ti+ 1 , jPin+1, j + Ti, j− 1 Pi,∗j−1 + Ti, j+ 1 Pi,nj+1 − Ti− 1 , jPin−1, j −  Ti− 1 , j + Ti+ 1 , j + Ti, j− 1 −  2 2 2 2 2 2 2 ∆ t  

 2Vbij φijc  + Ti, j + 1  Pi∗, j Q Vi, j =  2  ∆t

o bien: 2Vbij φ ij c  n  Pi, j + Ti+ 1 , jPin+1,j + Ti, j− 1 Pi,∗j−1 Ti− 1 , jPin−1, j −  Ti− 1 ,j + Ti+ 1 , j + Ti,j− 1 − 2 2 2 2 2 2 ∆t   Pi,∗j = +   2Vbij φij c  + Ti,j+ 1   ∆t 2 

Simulación de Yacimientos I

15

Ti, j + 1 Pin, j +1 − Q Vi, j 2

 2Vbij φijc    T + 1 i, j +  ∆t  2  

..................................................................................... (7.35)

PASO 2:

EL paso 2 consiste en aproximar la Ecuación (7.1) de la siguiente forma: ∂  kh ∂P    ∂x  µ ∂x 

n+1,*

∂  kh ∂P    + ∂y  µ ∂y 

n+1,*

= φhc

∂P   ∂t 

n+1,*

+ Q Vij ............................... (7.36)

Donde:

∂  kh ∂P    ∂x  µ ∂x 

n+1,∗

=

Ti− 1 , jPi∗−1, j − Ti− 1 , jPi∗, j − Ti+ 1 , jPin,j+1 + Ti+ 1 , jPin++1,1j 2

2

2

2

∆x i ∆y j

................ (7.37)

y

∂  kh ∂P    ∂y  µ ∂y 

n+1,∗

=

Ti, j− 1 Pi,∗j−1 − Ti, j− 1 Pi,∗j − Ti, j+ 1 Pin, j+1 + Ti, j+ 1 Pi,nj++11 2

2

2

∆x i ∆y j

2

................ (7.38)

Además:  ∂P     ∂t 

n+1,∗

=

(

2 Pin, j+1 − Pi∗, j ∆t

) .............................................................................. (7.39)

Llevando las Ecuaciones 7.37, 7.38 y 7.39 a la Ecuación 7.36, se tiene:

Simulación de Yacimientos I

16

Ti− 1 , jPi∗−1, j − Ti− 1 ,jPi,∗j − Ti+ 1 , jPin,j+1 + Ti+ 1 ,jPin+1+,1j + Ti, j− 1 Pi,∗j−1 − Ti,j− 1 Pi,∗j − Ti,j+ 1 Pin,j+1 + 2

2

Ti, j + 1 Pin, j++11 =

2Vbij φijc ∆t

2

(P

2

n +1 * i, j − Pi, j

)+ Q

2

2

2

2

Vi, j

De donde: 2 Vbij φ ij c  n + 1 2V φ c     Pi, j + T 1 Pi,nj++11 =  Ti − 1 , j − bij ij + Ti, j − 1  Pi∗, j − Ti + 1 , jPin++1,1j −  Ti + 1 , j + Ti, j + 1 + + i , j   2 2 2 2 2 ∆t ∆t  2    T 1 Pi∗−1, j − T 1 Pi∗, j−1 + Q Vi, j ...................................................................... (7.40) i−

,j 2

i, j−

2

El procedimiento para aplicar el método consiste en ejecutar el paso 1, Ecuación 7.35, para estimar los valores de Pi∗, j ; luego, ejecutar el paso 2, Ecuación 7.40, el cual conduce a un sistema tridiagonal de ecuaciones cuya solución permite estimar los valores de presión al tiempo t n+1 .

7.6 EL CRITERIO DE BALANCE DE MATERIALES

Las ecuaciones de flujo han sido deducidas partiendo del principio de conservación de la masa. Por esta razón, la aproximación numérica de estas ecuaciones también debe satisfacer este principio. Por definición, la compresibilidad es c=−

1 ∆V .................................................................................................. (7.41) V ∆P

Si el incremento de tiempo es:

Simulación de Yacimientos I

17

(7.42)

∆t = t n+1 − t n

Entonces el volumen de fluido removido a través del bloque (i, j) durante este intervalo de tiempo será: ∆V = Q ij ∆t ................................................................................................... (7.43) El volumen inicial de fluido en el bloque (i, j) es: V = Vbij φ ij .................................................................................................... (7.44) La caída de presión requerida en el bloque (i, j) para producir un volumen de fluido

∆V durante el intervalo de tiempo ∆t , será: Pin, j+1 − Pin, j .................................................................................................... (7.45) Llevando las Ecuaciones 7.45 y 7.44 a la Ecuación 7.41 se tiene:

C=−

∆Vij 1 Vbij φij Pin, j+1 − Pin, j

o bien:

(

)

∆Vij = CVbij φij Pin, j − Pin, j+1 ........................................................................... (7.46) Sustituyendo la Ecuación 7.12 en la Ecuación 7.46: ∆V j = γ i, j (Pin, j − Pin, j+1 ) Luego, el volumen total será: nx

ny

nx

ny

∑ ∑ ∆Vij = ∑ ∑ γ i, j (Pin, j − Pin, j+1 ) i =1 j =1

............................................................. (7.47)

i =1 j =1

Simulación de Yacimientos I

18

Llevando la Ecuación 7.43 a la Ecuación 7.47: nx

∑

ny

nx

∑ Qi, j ∆t = ∑

i =1 j =1

ny

∑ γ i, j (Pin, j − Pin, j+1 )

i =1 j =1

o bien: nx

ny

∑ ∑ γ i, j (Pin, j − Pin, j+1 ) 1=

i =1 j =1

............................................................................. (7.48)

ny

nx

∆t ∑

∑ Qi, j

i =1 j =1

La Ecuación 7.48 es muy útil como un criterio de chequeo. Después de realizar los cálculos para los tiempos t n y t n+1 , en todos los bloques de la malla, la Ecuación 7.48 debe cumplirse. Es decir, este criterio se aplica una vez se han cumplido los criterios de convergencia: nx

ny

∑ ∑ γ i, j (Pin, j − Pin, j+1 ) MB =

i = 1 j =1

nx

∆t ∑

ny

.......................................................................... (7.49)

∑ Qi, j

i =1 j =1

MB debe ser tan cercano a la unidad como la tolerancia previamente definida lo permita.

Simulación de Yacimientos I

19

1

8. FLUJO LINEAL DE UN FLUIDO COMPRESIBLE Considérese la ecuación de continuidad para el flujo de un fluido en una sola dirección: ∂ (ρ g ) ∂  kA ∂P    = Aφ ρ + A~ q ........................................................................ (8.1) ∂x  µ ∂x  ∂t Por definición Bg =

(V ) (V )

g CY

=

g CN

(ρ ) (ρ )

g CN

g CY

De donde:

(ρ )

g CY

(ρ )

g CN

=

De otro lado Bg

.............................................................................................. (8.2)

Bg

(ρ ) = (ρ )

g CN

g CY

(PM zRT) = (PM zRT)

CN

CY

=

PCN zT ............................................................... (8.3) TCNP

Llevando la Ecuación 8.3 a la Ecuación 8.2

(ρ )

g CY

(ρ )

g CN

=

PCN zT

es decir,

∂ (ρ g )CY ∂t

=

=

TCNP

TCN (ρ g )CN PCN T

⋅

TCN (ρ g )CN P PCN T z

∂ P   ......................................................................... (8.4) ∂t  z 

Sustituyendo las Ecuaciones 8.2 y 8.4 en la Ecuación 8.1 se tiene:

2

TCN (ρ g )CN ∂  P  ∂  kA (ρ g )CN ∂P  ⋅ = Aφ ⋅   + A~ q   ∂x  µ B g ∂x  PCN T ∂t  z  ~ Dividiendo por (ρ g )CN y teniendo en cuenta que (q v )CN = q

(ρ )

, se tiene:

g CN

∂  kA ∂P  AφTCN ∂  P  ⋅ = ⋅   + A (q v )CN ..................................................... (8.5) ∂x  µB g ∂x  PCN T ∂t  z  Si se define la función:  kA ∂P   .............................................................................................. (8.6) ui =   µB ∂x  g  i

Siguiendo un procedimiento completamente análogo al seguido para discretizar la ecuación para flujo lineal de un fluido incompresible se obtiene:

∂  kA ∂P  ≅ ∂x  µB g ∂x 

 kA   µB  g

 kA   ∂P    −    ∂x  1  µB i+  g  i+ 12 2

  ∂P      ∂x  1 i−  i− 12 2

∆x i

............................... (8.7)

O bien:

∂  kA ∂P  ≅ ∂x  µB g ∂x 

Ti+ 1 = 2

 kA   µB  g

  2(Pi+1 − Pi )   kA   −   ∆x + ∆x   µB i 1 i +    g  i+ 12

 kA 2  [∆x i + ∆x i+1 ]  µB g

∆x i

  2(Pi − Pi−1 )       ∆x + ∆x  + i 1 i    i− 12

................. (8.8)

  ........................................................................ (8.9)   i+ 12

3

Ti− 1 = 2

 kA 2  [∆x i−1 + ∆x i ]  µB g

  ...................................................................... (8.10)   i− 12

Así mismo, las transmisibilidades estarán dadas por las siguientes expresiones: 2k i A i k i+1 A i+1 ................................................. (8.11) k i+1 A i+1∆x i µ iB g i + k i A i ∆x i+1µ i+1B g i+1

Ti+ 1 =

[

Ti− 1 =

2k i A ik i−1A i−1 ................................................... (8.12) k i−1A i−1∆x i µ iB g i + k i A i ∆x i+1µ i−1B g i−1

2

2

]

[

]

Las transmisibilades en las ecuaciones numéricas que describen el flujo de un fluido

incompresible

o

levemente

compresible

se

suelen

tomar

independientemente del tiempo, pues en estos casos la viscosidad se considera independiente o poco dependiente del tiempo. Sin embargo, tal como se observa en las Ecuaciones 8.9 a 8.12, para el caso de un fluido compresible, la transmisibilidad depende de la viscosidad y el factor volumétrico, los cuales a su vez dependen de la presión y, por ende, del tiempo.

Por esta razón, las

ecuaciones (8.9) a (8.12) se suelen escribir y aproximar de la siguiente forma:

 kA  2   = [∆x i + ∆x i+1 ]  µB g 

*

∗ i+ 1 2

 kA  2   = [∆x i−1 + ∆x i ]  µB g 

*

∗ i− 1 2

T

T

T * i+ 12 =

i+ 1 2

...................................................................... (8.13)

i− 1 2

...................................................................... (8.14)

2k i A ik i+1A i+1 1 .................................................... (8.15) ∗ [k i+1A i+1∆x i + k i A i ∆x i+1 ] µ i+ 1 B ∗g i+ 1 2

2

4

T * i− 12 =

2k i A ik i−1A i−1 1 ................................................... (8.16) ∗ [k i−1A i−1∆x i + k i A i ∆x i−1 ] µ i− 1 B ∗g i− 1 2

2

Donde el asterisco, ∗ , indica el nivel de tiempo al cual se evalúa la viscosidad y el Si ∗ = n , T i∗+ 1 = T in+ 1

factor volumétrico del gas.

2

transmisibilidad explícita.

2

y T i∗− 1 = T in− 1 , se habla de 2

2

Si ∗ = n + 1, T i∗+ 1 = T in+ +11 y T i∗− 1 = T in− +11 , se habla de 2

2

2

2

transmisibilidad implícita. Una vez definido el nivel del tiempo, ∗ = n o ∗ = n + 1, los valores de µ ∗gi + 1 y B ∗gi + 1 2

son evaluados a una presión promedia; por ejemplo a Pi∗+ 1 , donde: 2

∗ i+ 1 2

P

ViPi∗ + Vi+1Pi∗+1 = ................................................................................... (8.17) Vi + Vi+1

Llevando las Ecuaciones 8.13 y 8.14 a la Ecuación 8.8, se tiene: ∗ ( ) ∗ ( ) ∂  kA ∂P  Ti+ 12 Pi+1 − Pi − Ti− 12 Pi − Pi−1 ≅ ............................................. (8.18) ∂x  µB g ∂x  ∆x i

Y llevando la Ecuación 8.18 a la Ecuación 8.5: Ti∗+ 1 (Pi+1 − Pi ) − Ti∗− 1 (Pi − Pi−1 ) 2

2

∆x i

=

AφTCN ∂  P    + A (q vij )CN PCN T ∂t  z 

Ti∗+ 1 (Pi+1 − Pi ) − Ti∗− 1 (Pi − Pi−1 ) = Vp ij 2

De donde:

2

TCN ∂  P    + (Q vij )CN PCN T ∂t  z 

2

5

T ∂ P Ti∗+ 1 Pi+1 −  Ti∗+ 1 + Ti∗− 1  Pi + Ti∗− 1 Pi−1 = Vp ij CN   + Q vij ......................... (8.19) 2 2 2 2  PCN T ∂t  z 

Si se aproxima: n +1 n ∂  P  1  P  P   ≅   −    ......................................................................... (8.20) ∂t  z  ∆t  z   z  

Se obtiene: ∗ i+ 1 2

T

Pi+1

Vp ij TCN  Pi  −  Ti∗+ 1 + Ti∗− 1 Pi + Ti∗− 1 Pi−1 = 2 2  2  PCN T∆t  z i 

  

n +1

P −  i  zi

  

n

  + (Q vij )CN ... (8.21) 

Sea γ ij =

Vp ij TCN PCN T

................................................................................................. (8.22)

Luego: ∗ i+ 1 2

T

γ ij  Pi    Pi+1 −  Ti∗+ 1 + Ti∗− 1  Pi + Ti∗− 1 Pi−1 = 2 2 2  ∆t  z i  

n +1

P  −  i   zi 

n

  + (Q vij )CN ............ (8.23) 

Las Ecuaciones 8.21 y 8.23 pueden ser aplicadas para generar un esquema de presión implícita y transmisibilidad explícita, haciendo ∗ = n y evaluando las presiones del lado izquierdo al tiempo t = t n+1 :

n i+ 1 2

T

n +1 i+1

P

γ ij  P  −  Tin+ 1 + Tin− 1  Pin+1 + Tin− 1 Pin−1+1 =  i  2 2 2  ∆t  z i 

De donde:

n +1

n

γ ij  P  −  i  + (Q vij )CN ∆t  z i 

6

n

n i+ 1 2

T

n +1 i+1

P

γ P  γ ij  n+1   Pi + Tin− 1 Pin−1+1 = − ij  i  + (Q vij ) ......... (8.24) −  Tin+ 1 + Tin− 1 + n +1  CN 2 2 2 ∆t  z i  ∆tz i  

La ecuación (8.24) representa el modelo numérico a resolver. En esta ecuación se tienen cuatro incógnitas: Pin−1+1 , Pin+1 , Pin+1+1 y Z ni +1 . El hecho de que la variable Z ni +1 sea una incógnita más hace que la solución del sistema de ecuaciones

generado por la ecuación (8.24) sea un proceso de ensayo y error. A continuación se presentan algunos algoritmos

basados en los métodos de solución de

ecuaciones discutidas en los capítulos 4 y 5.

8.1. SOLUCIÓN MEDIANTE EL ALGORITMO DE LA TRIDIAGONAL

El procedimiento de solución incluye los siguientes pasos: a.

Se asume un valor para z ni +1 en la iteración cero: (0 )

z ni+1 = z ni (0 )

z ni+1 representa el factor de compresibilidad z correspondiente al bloque i ,

evaluado al tiempo t n+1 , iteración 0 (cero). b.

Se resuelve el sistema tridiagonal generado por la ecuación (8.24). De esta forma se obtiene: (0)

(0)

(0)

Pin+1 , Pin+1+1 y Pin−1+1

7

(0)

c.

(0)

(0)

Con los valores de Pin+1 , Pin+1+1 y Pin−1+1 obtenidos en el paso anterior, se estima (1)

z ni+1 de datos PVT. (1)

d.

Una vez conocido el valor de z ni+1 , se resuelve de nuevo el sistema el sistema tridiagonal , ecuación (8.24), con la finalidad de obtener nuevos valores de presión: (1)

(1)

(1)

Pin+1 , Pin+1+1 y Pin−1+1 e.

Se repite los pasos b, c y d hasta que las presiones o los factores de compresibilidad, de cada bloque, en dos iteraciones consecutivas sean iguales dentro de cierto grado de tolerancia.

8.2. SOLUCIÓN POR EL MÉTODO DE GAUSS - SEIDEL

La Ecuación 8.24 puede ser escrita como:

Pin+1 =

Tin+ 1 Pin+1+1 + Tin− 1 Pin−1+1 + 2

2

Tin+ 1 + Tin− 1 2

2

γ i Pin ⋅ − (Q vij )CN ∆t z ni ............................................... (8.25) γi + ∆tz ni+1

Para resolver esta ecuación por el método de Gauss – Seidel, se hace necesario escribirla de la siguiente forma:

8

(k )

(k +1)

Pin+1 =

(k +1)

γ i Pin ⋅ − (Q vij )CN ∆t z ni ................................................ (8.26) γi + (k ) ∆t z ni+1

Tin+ 1 Pin+1+1 + Tin− 1 Pin−1+1 + 2

2

Tin++11 + Tin− 1 2

2

El procedimiento de solución puede ser resumido así: a.

Para el tiempo t n , se conocen los valores de Pin y z ni , para cada uno de los bloques de la malla.

b.

Para la primera iteración, k = 0 , se asume: (0 )

Pin+1 = Pin (0 )

z ni+1 = z ni (1)

c.

Se calcula Pin+1 de la ecuación (8.26), para cada una de los bloques de la malla. (1)

d.

Se calcula z ni+1 para cada bloque de la malla.

e.

Se repite el anterior procedimiento hasta que se cumplan los criterios de convergencia: (k +1)

I.

(k )

Pin+1 − Pin+1 (k +1)

n +1 i

P

< ε1

9

II.

z ki +1 − z ki < ε2 z ki +1

III.

MB − 1 < ε 3

8.3. SOLUCIÓN POR EL MÉTODO DE PSOR (k +1)

∗(k +1)

Supóngase que la presión Pin+1 en la ecuación (8.26) se nota como Pin+1 . Es decir: (k )

∗(k +1) n +1 i

P

=

(k +1)

Tin+ 1 Pin+1+1 + Tin− 1 Pin−1+1 + 2

2

Tin+ 1 + Tin− 1 2

2

γ i Pin ⋅ − (Q vij )CN ∆t z ni ............................................... (8.27) γi + (k ) ∆t z ni+1 (k +1)

El método PSOR se puede aplicar a este caso calculando el valor de Pin+1 de la siguiente forma: (k +1)

(k )

∗(k +1)

Pin+1 = (1 − ω)Pin+1 + ωPin+1 .............................................................................. (8.28) ∗(k +1) n+1 i

El valor de P

está dado por la ecuación (8.27).

El procedimiento es el

siguiente: (k +1)

a.

En el momento de calcular Pin+1 se conoce la siguiente información: Pin y z ni .

b.

Para la iteración cero, k = 0 , se asume:

10

(0 )

Pin+1 = Pin (0 )

z ni+1 = z ni (1)

c.

Se calcula Pin+1+1 , para todos los bloques de la malla.

d.

Se calcula z ni+1 de datos PVT para cada uno de los bloques de la malla.

e.

Se repite el anterior procedimiento hasta que se cumplan los criterios de

(1)

convergencia (véase numeral 8.2).

1

9. FLUJO DE UN FLUIDO COMPRESIBLE (2D) Considérese la ecuación de flujo, 2D, coordenadas cartesianas: ∂ (ρ g ) ∂  kH ∂P  ∂  kH ∂P    = Hφ   + ρg ρg + H~ q ............................(9.1) ∂x  µ ∂x  ∂y  µ ∂y  ∂t Siguiendo un procedimiento similar al seguido para flujo lineal de un fluido compresible (capitulo 8), la Ecuación 9.1 puede ser escrita como:

T ∂  kH ∂P  ∂  kH ∂P  ∂ P + = Hφ CN   + H(q V )CN ...............(9.2)     ∂x  µβ g ∂x  ∂y  µβ g ∂y  PCN T ∂t  Z  Siguiendo un procedimiento análogo al seguido para obtener la ecuación para el caso lineal, se puede obtener:

∂  kH ∂P  = ∂x  µβ g ∂x 

 kH   2(Pi+1, j − Pi, j )  kH   2(Pi,j − Pi−1, j )    −      µβ     g  i+ 12, j  ∆x i+1 + ∆x i   µβ g  i− 12,j  ∆x i + ∆x i−1 

∂  kH ∂P  = ∂y  µβ g ∂y 

∆x i  kH   2(Pi, j+1 − Pi, j )  kH   2(Pi, j − Pi−1, j )      −     µβ    g  i, j+ 12  ∆y i+1 + ∆y i   µβ g  i, j− 12  ∆y i + ∆y i−1  ∆y j

Adicionalmente,

Ti+* 1 , j 2

∆y j

*

 kH  2  .....................................................(9.5) =  µβ  [ ] ∆ x + ∆ x i i + 1 g   i+ 12, j

(9.3)

(9.4)

2

Ti*− 1 , j 2

∆y j

Ti,*j− 1

*

 kH  2  = .....................................................(9.6)  µβ   g  i− 12, j [∆x i + ∆x i−1 ]

*

2

∆x i

Ti,*j+ 1

 kH  2  = .....................................................(9.7)  µβ   g  i, j− 12 ∆y j + ∆y j−1

[

]

*

2

∆x i

 kH  2  = .....................................................(9.8)  µβ   g  i, j+ 12 ∆y j + ∆y j+1

[

]

De esta forma las ecuaciones (9.3) y (9.4) pueden ser escritas como: * * ∂  kH ∂P  Ti+ 12,j (Pi+1, j − Pi,j ) − Ti− 12, j (Pi, j − Pi−1, j ) = ..........................(9.9) ∂x  µβ g ∂x  ∆y j ∆x i

* * ∂  kH ∂P  Ti, j+ 12 (Pi,j+1 − Pi,j ) − Ti, j− 12 (Pi, j − Pi,j−1 ) = ..........................(9.10) ∂y  µβ g ∂y  ∆y j ∆x i

Llevando las ecuaciones (9.9) y (9.10) a (9.2), se obtiene: Ti+* 1 ,j (Pi+1,j − Pi,j ) − Ti−* 1 ,j (Pi,j − Pi−1,j ) + Ti,*j+ 1 (Pi,j+1 − Pi,j ) − Ti,*j− 1 (Pi,j − Pi,j−1 ) 2

2

2

2

∆y j ∆x i

=

n+1 n  Pi, j   TCN 1  Pi, j    −    + H(q V )CN .......................................(9.11) Hφ PCNT ∆t  Z   Z   

Si se define:

Vpij = ∆x i ∆y iHφ ............................................................................

(9.12)

3

γ ij = Vpij

TCN ............................................................................... PCN T

(9.13)

Reorganizando términos la ecuación (9.23) puede ser escrita como:

Ti*+ 1 , jPi+1, j + Ti*− 1 , jPi−1, j −  Ti*+ 1 , j + Ti*− 1 , j + Ti,*j+ 1 + Ti,*j− 1  Pi, j + Ti,*j− 1 Pi, j−1 + 2 2 2 2 2 2 2  n+1 n P   γ ij  Pij    −  ij   + (Q vij ) .................................... Ti,*j+ 1 Pi,j+1 = CN  Z  2 ∆t  Z     

(9.14)

Donde: Ti+* 1 ,j = 2

2k i,jHi,j ∆y ik i+1,jHi+1,j k i+1,jHi+1,j ∆x i + k i,jHi,j ∆x i+1

1 ................................. µ gi+ 1 ,jB gi+ 1 ,j

(9.15)

2

2

Expresiones análogas se pueden obtener para Ti−* 1 ,j , Ti,*j+ 1 y Ti,*j− 1 . 2

2

2

La ecuación (9.14) genera el sistema de ecuaciones cuya solución conlleva al cálculo de distribución de presiones en el sistema.

Para transmisibilidades

explícitas y presiones implícitas, se tiene: Tin+ 1 ,jPin+1+,1j + Tin− 1 , jPin−1+,1j −  Tin+ 1 , j + Tin− 1 , j + Ti,nj+ 1 + Ti,nj− 1 Pi,nj+1 + Ti,nj+ 1 Pi.nj++11 + 2 2 2 2 2 2 2  n i, j− 1 2

T

n+1 i, j−1

P

n+1 P γ ij  Pij    −  ij =  Z ij ∆t  Z ij   

   

n

  + (Q vij ) CN  

De donde:  γ ij  n+1 P + Tin+ 1 , jPin+1+,1j + Tin− 1 , jPin−1+,1j −  Tin+ 1 , j + Tin− 1 , j + Ti,nj+ 1 + Ti,nj− 1 + n +1  i, j  2 2 2 2 2 2 tz ∆ i, j   n i, j+ 1 2

T

n +1 i. j+1

P

+T

n i, j− 1 2

n +1 i, j−1

P

= (Q vij )CN

γ ij  Pij −  ∆t  Z ij

n

  ....................................  

(9.16)

4

La ecuación (9.16) puede ser escrita como: S ijPi.nj+−11 + WijPin−1+,1j + CijPi,nj+1 + E ijPin+1+,1j + NijPi.nj++11 + = Fij ...........................

(9.17)

Donde: S ij = Ti,nj− 1 ....................................................................................

(9.18)

Wij = Tin− 1 , j ...................................................................................

(9.19)

   n  γ ij  γ ij  n n n  = − c ij + C ij = − Ti+ 1 , j + Ti− 1 , j + Ti, j+ 1 + Ti, j− 1 + (k )  .... n +1   2 2 2 2 ∆ tZ i, j    ∆t Z nij+1  

(9.20)

E ij = Tin+ 1 , j ....................................................................................

(9.21)

Nij = Ti,nj+ 1 ....................................................................................

(9.22)

2

2

2

2

Fij = (Q vij )CN

γ ij  Pij −  ∆t  Z ij

n

  ................................................................  

(9.23)

La ecuación (9.17) genera un sistema penta-diagonal de ecuaciones el cual puede ser resuelto por algún método directo o iterativo.

Tal como se discutió en el

capítulo 5. A continuación se discute la aplicación de algunos de estos métodos. SOLUCIÓN POR MÉTODOS DIRECTOS

El procedimiento de solución mediante la aplicación de métodos directos puede ser resumido de la siguiente forma:

5

1. Se conoce Z nij .

2. Se asume Z nij+1 = Z nij .

3. Se resuelve el sistema se ecuaciones generado por la

Ecuación 9.17

mediante algún método directo (eliminación Gaussiana, Gauss – Jordan, etc.). 4. Se evalúa Z nij+1 conociendo Pijn+1 calculado en el paso 3 y de datos PVT, ∀ij .

5. Si el valor Z nij+1 obtenido del paso 4 es igual al valor asumido en el paso 2, se ha obtenido el verdadero valor; de lo contrario se asume un nuevo Z nij+1 y se repiten lo cálculos. SOLUCIÓN POR MÉTODOS ITERATIVOS

Esta sección discute algunos de los métodos iterativos discutidos en capítulos anteriores. Solución Mediante Aplicación del Método de Gauss-Seidel

La ecuación (9.17) puede ser escrita como:

Pijn+1 =

[

Fij − S ijPi,nj−+11 + WijPin−1+,1j + E ijPin+1+,1j + NijPi,nj++11    γ ij  Ci,j = − c ij +   ∆t Z ij  

    

] ...............................

Para los cálculos de la iteración k, se puede escribir:

(9.24)

6

(k +1)

Pijn+1

(k ) (k +1) (k ) (k +1)   Fij − a ij Pin+1+,1j + b ij Pin−1+,1j + dij Pi,nj++11 + e ij Pi,nj−+11    ............................... =     γ ij   (k ) (k )  Ci, j = − c ij +   ∆t Z nij+1   

(9.25)

El procedimiento puede ser resumido de la siguiente forma: 1. Se conoce Pijn , Znij , ∀i, j .

2. Para k = 0 , se asume:

(0 )

n+1 i+1, j

P

=P

n i+1, j

(0 )

n+1 i, j+1

,P

=P

n i, j+1

(0 )

y Z nij+1 = Z nij .

(1)

3. Se calcula Pijn+1 , ∀i, j . (1)

4. Se calcula Z nij+1 de datos PVT.

5. Se repite el procedimiento hasta que: (k +1)

a.

(k )

Pijn+1 − Pijn+1 (k +1)

≤ ε1 .

n +1 ij

P

(k +1)

b.

(k )

Z ij − Z ij (k +1)

≤ ε2

Z ij

c. MB − 1.0 ≤ ε 3

Solución Mediante Aplicación del Método PSOR (k +1)

* (k +1)

Supóngase que la presión Pijn+1 de la ecuación (9.25) se nota como: Pijn+1 ; luego, se tiene:

7

* (k +1)

Pijn+1

(k +1) (k ) (k )   (k +1) Fij − S ij Pi,nj−+11 + Wij Pin−1+,1j + E ij Pin+1+,1j + Nij Pi,nj++11    ............................ =     γ ij   (k ) (k )  Ci, j = − c ij +   ∆t Z nij+1   

(9.26)

(k +1)

El método PSOR en este caso consiste en calcular un valor ponderado de Pijn+1 * (k +1)

(k )

utilizando el valor de Pijn+1 y el valor de Pijn+1 en la siguiente forma: (k +1)

* (k +1)

(k )

Pijn+1 = w Pijn+1 + (1 − w )Pijn+1 .............................................................

(9.27)

* (k +1)

El valor de Pijn+1 está dado por la ecuación (9.25). El procedimiento de solución es como se describe a continuación: (k )

(k +1)

(k )

(k +1)

(k )

1. Se conoce Pijn+1 , Pin−1+,1j , Pi,nj++11 , Pi,nj−+11 , Z nij+1 , Znij y Pijn ∀i, j . Para k = 0 se asume: (0 )

(0 )

Pijn+1 = Pijn y Z nij+1 = Z nij (1)

2. Se calcula Pijn+1 ∀i, j . (1)

3. Se calcula Z nij+1 ∀i, j , en base a datos PVT. 4. Se verifican los criterios de convergencia ∀i, j .

8

Solución Mediante Aplicación del Método LSOR

Si se lleva la ecuación (9.26) a la ecuación (9.27), se obtiene:

(k +1)

Pijn+1

    (k +1) (k ) (k )   (kn++11)   n+1 n +1 n +1  Fij − S ij Pi, j−1 + Wij Pi−1, j + E ij Pi+1, j + Nij Pi, j+1 +   (k )  .  = (1 − w )Pijn+1 + w         γ ij     (k )   C i, j = − c ij + (k )  n +1       ∆ t Z ij     

(9.28)

Supóngase que la malla se recorre por filas tal como se ilustra en la Figura 9.1:

}

M . . .

Iteración k

j+1

j

(k +1)

Pin−1+,1j

(k +1)

(k +1)

Pi,nj+1 Pin+1+,1j

}

j −1

. .

Iteración k+1

.

1 1. . . . . . i − 1

i

i +1 . . . .

Figura 9.1 – Recorrido de Malla por Filas.

N

9

El método consiste en resolver el sistema tridiagonal de ecuaciones generado al reorganizar la ecuación (9.28) de tal forma que las incógnitas sean las (k +1)

(k +1)

(k +1)

presiones de la fila j en la iteración k + 1: Pin−1+,1j , Pijn+1 y Pin+1+,1j . De la ecuación (9.28):

(k +1)

Pijn+1

(k +1)

(k )

(k )

wWij Pin−1+,1j wE ij Pin+1+,1j = (1 − w )Pijn+1 − − γ ij γ ij c ij + c + ij (k ) (k ) ∆t Z nij+1 ∆t Z nij+1

   (k +1) (k )   F − S P n+1 − N P n+1  ij ij i, j−1 ij i, j+1  + w γ   ij c ij + (k )     ∆t Z nij+1

De donde: (k +1)

(k )

(k +1) wE ij Pin+1+,1j wWij Pin−1+,1j + Pijn+1 + γ ij γ ij c + c ij + ij (k ) (k ) ∆t Z nij+1 ∆t Z nij+1

   (k +1) (k )  (k )  Fij − S ij Pi,nj−+11 − Nij Pi,nj++11   = (1 − w )Pijn+1 + w  γ ij   c ij + (k )     ∆t Z nij+1

(9.29)

Las presiones del lado derecho de la ecuación (9.29) son conocidas, ó bien de iteración anterior o bien de la iteración actual, de filas ya barridas.

El

procedimiento de solución puede ser resumido de la siguiente forma: 1. Se conoce Pijn , Znij ∀i, j . (0 )

n+1 ij

2. Para k = 0 , se asume: P

(0 )

= P , Z nij+1 = Z nij . n ij

3. Se resuelve la malla barriéndola por filas.

En cada fila se resuelve el

sistema de ecuaciones generado por (9.29). (1)

4. Se obtiene Z nij+1 ∀i, j .

5. Se prueban criterios de convergencia. Si se cumplen se continua con el siguiente nivel de tiempo; si no se cumplen se repite el proceso.