Universidad Nacional Autónoma de México

Posgrado en Ciencias de la Tierra

Instituto de Geofı́sica
Modelado de sistemas terrestres

Simulación Numérica de Inyección de Agua en Yacimientos

Petroleros empleando el Método de Lı́neas de Corriente

TESIS
QUE PARA OPTAR POR EL GRADO DE:
MAESTRO EN CIENCIAS DE LA TIERRA

P R E S E N T A:
ROBERTO CARLOS M A R T Í N E Z CARRADA

Director de Tesis:
Dr. Luis M. De la Cruz Salas
Instituto de Geofı́sica

México, D.F. Enero 2014

UNAM – Dirección General de Bibliotecas
Tesis Digitales
Restricciones de uso

DERECHOS RESERVADOS ©
PROHIBIDA SU REPRODUCCIÓN TOTAL O PARCIAL

Todo el material contenido en esta tesis esta protegido por la Ley Federal
del Derecho de Autor (LFDA) de los Estados Unidos Mexicanos (México).

El uso de imágenes, fragmentos de videos, y demás material que sea

objeto de protección de los derechos de autor, será exclusivamente para
fines educativos e informativos y deberá citar la fuente donde la obtuvo
mencionando el autor o autores. Cualquier uso distinto como el lucro,
reproducción, edición o modificación, será perseguido y sancionado por el
respectivo titular de los Derechos de Autor.
Agradecimientos

El presente trabajo no habrı́a sido posible sin la generosa ayuda de muchas per-
sonas e intituciones. A mi alma máter, la Universidad Nacional Autónoma de México,
al Instituto de Geofı́sica y al Posgrado en Ciencias de la Tierra, ası́ como a su person-
al académico y administrativo que me apoyaron durante mis estudios de posgrado. Al
programa nacional de Becas de Posgrado del CONACyT que disfruté durante dos años.

A mi director de tesis, el Dr. Luis De la Cruz por la paciencia que me tuvo durante
la elaboración de la tesis, sus consejos, observaciones y al tiempo dedicado en las se-
siones de tutorı́a.

A mis profesores del posgrado que ayudaron en mi formación a lo largo de estos

años, sus conocimientos y experiencia sirvieron también de inspiración para mi trabajo
de investigación.

A mi jurado de examen de grado: Dr. Ismael Herrera Revilla, Dra. Graciela Herrera
Zamarrón, Dr. Ernesto Rubio Acosta, Dr. Carlos Ortiz Alemán y al Dr. Luis M. De la
Cruz Salas. Sus comentarios y correcciones sobre mi trabajo enriquecieron el contenido
de éste.

A todos ellos, ¡Gracias!

Dedicatoria

A mi más dura crı́tica

y mi más ferviente apoyo:
a mis padres Leonor y Prudencio
y a Hugo, mi hermano.

A Érika, un alma libre y creativa.

Por nuestros colores,
la bohemia y nuestra literatura.
... ya viene el fin de semana.

Con personas como ustedes en mi vida,

pase lo que pase, ¡ya he ganado!
¡Éxito!
Índice General

Índice General 4

Índice de Figuras 7

Índice de Tablas 9

1 Introducción 15
1.1. Estructura de la Tesis . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2 Simulación de Yacimientos 19
2.1. Métodos Clásicos en Ingenierı́a de Yacimientos . . . . . . . . . . . . . . 20
2.1.1. Métodos de Balance de Materiales . . . . . . . . . . . . . . . . . 20
2.1.2. Métodos de Curva de Declinación . . . . . . . . . . . . . . . . . 20
2.1.3. Métodos Estadı́sticos . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1.4. Métodos Analı́ticos . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2. Métodos de Simulación de Yacimientos . . . . . . . . . . . . . . . . . . 21
2.2.1. Etapas de la Simulación . . . . . . . . . . . . . . . . . . . . . . 21
2.2.2. Clasificación de Simuladores de Yacimientos . . . . . . . . . . . 22
2.2.3. Aplicaciones de la Simulación de Yacimientos . . . . . . . . . . 23
2.3. Términos Usados en la Simulación Numérica . . . . . . . . . . . . . . . 24

3 Modelo Conceptual 27
3.1. Formación del Aceite y Gas . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1.1. Transformación de la Materia Orgánica . . . . . . . . . . . . . . 28
3.2. Sistema Petrolero . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3. Propiedades de la Roca y los Fluidos en el Yacimiento . . . . . . . . . . 30
3.3.1. Propiedades de la Roca . . . . . . . . . . . . . . . . . . . . . . . 30
3.3.2. Propiedades de los Fluidos . . . . . . . . . . . . . . . . . . . . . 33
3.4. Mojabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4
 Índice General

3.5. Procesos de Desplazamiento de Fluidos . . . . . . . . . . . . . . . . . . 37

3.6. Interacción Roca-Fluidos . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.7. Etapas de la Recuperación de Hidrocarburos . . . . . . . . . . . . . . . 39

4 Modelos Matemáticos 41
4.1. Formulación Axiomática . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2. Forma Conservativa de las Ecuaciones de Balance . . . . . . . . . . . . 42
4.3. Modelo Matemático para Sistemas Multifásicos . . . . . . . . . . . . . 43
4.4. Formulación Presión-Saturación para Flujos Bifásicos . . . . . . . . . . 44
4.5. Método de Lı́neas de Corriente . . . . . . . . . . . . . . . . . . . . . . 45
4.5.1. Conceptos Fundamentales . . . . . . . . . . . . . . . . . . . . . 45
4.5.2. Lı́neas de Corriente . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.5.3. Funciones de Corriente en 2D . . . . . . . . . . . . . . . . . . . 50
4.5.4. Funciones de Corriente y Lı́neas de Corriente en 3D . . . . . . . 51
4.5.5. Lı́neas de Corriente y Tiempo de Vuelo . . . . . . . . . . . . . . 53
[Link]. Tiempo de Vuelo como Coordenada Espacial . . . . . . 55
[Link]. Cálculo del Tiempo de Vuelo y Lı́neas de Corriente en
2D y 3D . . . . . . . . . . . . . . . . . . . . . . . . . . 55

5 Modelo Numérico 61
5.1. Método de Volumen Finito (MVF) . . . . . . . . . . . . . . . . . . . . 61
5.2. Modelo Discreto de Flujo en dos Fases con MVF . . . . . . . . . . . . . 62
5.2.1. Discretización de la Ecuación de Presión . . . . . . . . . . . . . 62
5.2.2. Discretización de la Ecuación de Saturación . . . . . . . . . . . 66
5.2.3. Algoritmo IMPES . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.3. Condiciones Iniciales y de Frontera . . . . . . . . . . . . . . . . . . . . 69
5.4. Cálculo de la Saturación en las Caras . . . . . . . . . . . . . . . . . . . 71
5.5. Cálculo del Transporte a lo largo de las Lı́neas de Corriente . . . . . . 71
5.5.1. Discretización de la Ecuación de Saturación 1D . . . . . . . . . 72

6 Modelo Computacional y Resultados Numéricos 73

6.1. Software TUNAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.2. Implementación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.2.1. Calibración 1D con el modelo de Buckley- Leverett . . . . . . . 75
6.2.2. Estudio de Caso . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.2.3. Caso homogéneo (MVF) . . . . . . . . . . . . . . . . . . . . . . 79
6.2.4. Caso no-homogéneo (MVF) . . . . . . . . . . . . . . . . . . . . 83
6.2.5. Lı́neas de Corriente caso no-homogéneo . . . . . . . . . . . . . . 87

7 Conclusiones 93

A Formulación Presión-Saturación 95
A.1. Ecuación de Presión . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

5
Índice General

A.2. Ecuación de Saturación . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

B Gradiente Conjugado para Sistemas Lineales Dispersos 99

C Aplicación del Método de Volumen Finito 103

Referencias 105

6
Índice de Figuras

2.1. Etapas del proceso de modelado y simulación. . . . . . . . . . . . . . . . . 22

2.2. Malla para un área en 2D. . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.3. Corte transversal para un dominio 2D. . . . . . . . . . . . . . . . . . . . . 25

3.1. Depositación de la materia orgánica en ambientes reductores y su posterior

sepultamiento en la cuenca sedimentaria. . . . . . . . . . . . . . . . . . . . 28
3.2. Trampa tı́pica con sincronı́a de los distintos elementos del Sistema Petrolero. 30
3.3. A la izquierda, observamos poros interconectados y poros aislados a la derecha. 31
3.4. Correlación permeabilidad-porosidad . . . . . . . . . . . . . . . . . . . . . 33
3.5. Relación densidad-presión. . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.6. A la izquierda un fluido no mojador, al centro una fluido de mojabilidad
intermedia y a la derecha un fluido mojador. . . . . . . . . . . . . . . . . . 36
3.7. Procesos de drenajes e imbibición . . . . . . . . . . . . . . . . . . . . . . 37
3.8. En rojo se muestra la curva de drenaje primario e imbibición en negro y
delimitan el comportamiento de la presión capilar. En verde pc para una
saturación intermedia, en amarillo una inversión de los valores para la lı́nea
verde. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.9. Inyección de agua en yacimientos petroleros. . . . . . . . . . . . . . . . . . 40

4.1. Representación esquemática de un medio continuo . . . . . . . . . . . . . . 42

4.2. A la izquierda, el campo de velocidad en un dominio dado. A la derecha
observamos las lı́neas de corriente tangentes al campo de velocidad local.
Imágenes obtenidas con OpenDx. . . . . . . . . . . . . . . . . . . . . . . . 46
4.3. Tubo de corriente. Imágenes obtenidas con OpenDx. . . . . . . . . . . . . 46
4.4. Tiempo de vuelo τ para una partı́cula de prueba. . . . . . . . . . . . . . . 47
4.5. Una función de corriente no depende de la definición de la trayectoria. . . . 48
4.6. Una función de corriente no depende de la definición de la trayectoria. . . . 49
4.7. Lı́neas de corriente, función de corriente y tubo de corriente en 2D. . . . . 51

7
Índice de Figuras

4.8. Dos superficies de corriente en (a) y en (b) y sus intersecciones en (c) definen
un tubo de corriente. Los bordes del tubo de corriente son lı́neas de corriente. 53
4.9. Bloque de volumen finito para el cálculo del tiempo de vuelo. . . . . . . . 56
4.10. Partı́cula en un esquema de celdas 2D para el cálculo del tiempo de vuelo. 59
4.11. Cálculo del tiempo de vuelo para una sola celda. . . . . . . . . . . . . . . . 59
4.12. Representación del algoritmo de Pollock en 2D. . . . . . . . . . . . . . . . 60
4.13. Lı́nea de corriente en un sistema 3D. . . . . . . . . . . . . . . . . . . . . . 60

5.1. Volumen de control en tres dimensiones. Las letras mayúsculas representan

los centros de los volúmenes de control, mientras que las minúsculas las
caras de dichos volúmenes. . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.2. La matriz resultante presenta 7 bandas en el caso tridimensional. . . . . . 66

6.1. (a) Arquitectura general de TUNAM. (b) Los paquetes FVM y Geom se
muestran de manera esquemática: las letras G, S y A, a la derecha de la
figura, significan Generalization, Specialization y Adaptor, respectivamente. 75
6.2. Medio homogéneo de 300 m de longitud, en un principio saturado de aceite. 76
6.3. Saturación obtenida mediante lı́neas de corriente en rojo y con volumen
finito en verde. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4. Geometrı́a del dominio del caso de cinco pozos. . . . . . . . . . . . . . . . 78
6.5. Evolución del frente de saturación para un medio homogéneo. . . . . . . . 82
6.6. Valores de permeabilidad en el dominio de estudio, de izquierda a derecha:
permeabilidad en dirección x, y y z. . . . . . . . . . . . . . . . . . . . . . . 83
6.7. Evolución del frente de saturación para un medio de permeabilidad variable. 86
6.8. Campo de velocidad para el dominio con permeabilidad variable. . . . . . 89
6.9. Comparación de los métodos MVF y SLS en los mismos dı́as de la simulación. 92

8
Índice de Tablas

3.1. Clasificación de permeabilidad de rocas . . . . . . . . . . . . . . . . . . . . 32

3.2. Valores tı́picos de viscosidad en aceites. . . . . . . . . . . . . . . . . . . . . 36

5.1. Condiciones iniciales para la presión y la saturación. . . . . . . . . . . . . 69

6.1. Condiciones iniciales y de frontera para la presión y la saturación. . . . . . 77

6.2. Datos para el caso de estudio. . . . . . . . . . . . . . . . . . . . . . . . . . 78

9
Abstract

In this thesis I study a three-dimensional model of two-phase flow in porous media

to simulate water injection into oil fields. The mathematical models are derived from
axiomatic formulation. This model consists of a set of partial differential equations, non-
linear and coupled, which are solved using the IMPES algorithm. I use the pressures
and saturations of water and oil phases in order to calculate the total velocity field
in the domain of interest. Streamlines are paths whose points are tangent to the
field velocity at given instant of time. To discretize the resulting system of equations
I used the Finite Volume Method (FVM) which is derived from the conservative
version of balance equations. The discrete system satisfies the principle of conservation
of the extensive properties of the model for each control volumes. The relationship
between the numerical algorithm and the principle of conservation is one of the biggest
advantages of FVM.
This thesis also presents the development and application of a three-dimensional, two-
phase streamline simulator applied to field scale multiwell problems. The underlying
idea of the streamline method is to decouple the full 3D problem into multiple 1D
problems along streamlines. Fluids are moved along the natural streamline grid, rather
than between discrete gridblocks as in conventional methods. Permeability effects and
well conditions dictate the paths that the streamlines take in 3D; the geometry and
density of the streamlines reflects the geological impact of preferential flow paths
introducing higher line density in regions of high porosity and permeability. The physics
of the displacement is captured by the 1D solutions mapped along streamlines. The
one-dimensional solution makes this approach extremely fast and effective to model
flows in fields where there are many heterogeneities.

11
Resumen

En este trabajo se estudia un modelo de flujo bifásico en medios porosos en tres

dimensiones a fin de simular numéricamente la inyección de agua en yacimientos
petroleros. El modelo matemático se obtiene mediante la formulación axiomática. Este
modelo consiste de un conjunto de ecuaciones diferenciales parciales, no lineales y
acopladas, las cuales se resuelven utilizando el algoritmo IMPES. La presión y la
saturación de las fases agua y aceite ayudan a calcular el campo de velocidades total
en el dominio de interés. Las lı́neas de corriente son trayectorias que en cada uno
de sus puntos son tangentes a dicho campo de velocidad en un instante de tiempo
dado. La discretización del sistema de ecuaciones resultante se llevó a cabo mediante el
Método de Volumen Finito (MVF), que se deriva a partir de la forma conservativa de
las ecuaciones de balance. El sistema discreto, cumple con el principio de conservación
de las propiedades extensivas del modelo para cada uno de los volúmenes de control.
La relación existente entre el algoritmo numérico y el principio de conservación es una
de las mayores ventajas del MVF.
Esta tesis presenta también el desarrollo y la aplicación de un simulador 3D de dos
fases aplicado a escala de yacimiento para problemas con múltiples pozos. La principal
idea del método de Lı́neas de Corriente es descomponer un problema completamente en
3D en uno con múltiples lı́neas de corriente 1D. Los fluidos se mueven a lo largo de las
lı́neas de corriente discretizadas, en lugar de entre los bloques de la malla de métodos
convencionales. Los efectos de permeabilidad y la distribución de los pozos dictan los
caminos que llevan las lı́neas de corriente en 3D; la geometrı́a y la densidad de las lı́neas
de corriente reflejarán el impacto geológico sobre los caminos preferenciales del flujo,
introduciendo mayor densidad de lı́neas en regiones de alta porosidad y permeabilidad.
La fı́sica del desplazamiento es capturada por las soluciones 1D asignadas a lo largo
de cada lı́neas de corriente. La solución unidimensional hace que este enfoque sea
extremadamente rápido y efectivo para modelar flujos en yacimientos en donde existen
muchas heterogeneidades.

13
Capı́tulo 1

Introducción

El modelado de yacimientos basado en lı́neas de corriente ha sido usado en la

industria del petróleo desde los años 50’s, [15, 16]. Recientemente se ha incrementado
el interés en estas técnicas debido al desarrollo de nuevas formas de caracterizar los
yacimientos, [48]. Hoy en dı́a es posible obtener modelos estáticos que integran datos
geológicos y geofı́sicos en tres dimensiones muy detallados, [43, 50]. Lo anterior se
traduce en modelos de cientos de millones de nodos para una simulación tı́pica [26],
que requieren de la capacidad de cómputo adecuada, [34, 40, 49]. El incremento en
la resolución produce más variables e incertidumbres que no pueden manejarse con
las técnicas estándares (métodos basados en mallas numéricas), ni con los equipos
de cómputo actualmente existentes. En una simulación numérica se busca entender
y cuantificar el impacto de los elementos desconocidos del modelo estático sobre el
flujo de los fluidos y el transporte de hidrocarburos, para realizar un manejo prudente
y eficiente del yacimiento, [25, 28, 36]. Los recientes desarrollos en las técnicas de
simulación con lı́neas de corriente ofrecen un alto potencial para atacar este tipo de
problemas, pues proporciona herramientas para la simulación rápida de yacimientos
a escalas finas, [42]. La evolución de los frentes de inyección y su interacción con las
heterogeneidades del yacimiento pueden ser visualizadas fácil y rápidamente, y por lo
tanto proveen de una manera natural e intuitiva para caracterizar dinámicamente un
yacimiento. Los fundamentos de estas técnicas se remontan al siglo XIX, desde entonces
ha habido un constante desarrollo de este tipo de métodos, [1]. Una lı́nea de corriente
es una trayectoria que en cada uno de sus puntos es tangente al campo de velocidad en
un instante de tiempo dado. Dado que la velocidad es dependiente del tiempo, entonces
las lı́neas se trazan usando campos de velocidades instantáneos.
La velocidad se obtendrá a partir de la solución numérica de la ecuación de presión.
La ecuación de presión satisface una ecuación diferencial (parabólica o elı́ptica) la
cual una vez obtenida sirve de base para definir y resolver una familia de ecuaciones

15
1. Introducción

hiperbólicas (o, casi hiperbólicas). Ası́, la mayor resolución se logrará aumentando la

finura del mallado para obtener la solución de la presión y reduciendo los intervalos de
tiempo de las ecuaciones de transporte para las componentes.
Una vez obtenidas las presiones, se obtienen las velocidades mediante la ley de Darcy.
Usando las velocidades se trazan las lı́neas de corriente mediante el algoritmo de Pollock
y se calcula el Tiempo de Vuelo (TOF), por sus siglas en inglés, a lo largo de las mismas.
Tı́picamente estas lı́neas comienzan en un pozo inyector y se hace el seguimiento
hasta los pozos productores. Las ecuaciones de transporte se resuelven sobre las lı́neas
de corriente, transformando el problema tridimensional a varios problemas en una
dimensión, [24]. Este cálculo se puede realizar en paralelo simplemente repartiendo
grupos de lı́neas de corriente, a procesadores diferentes y en ellos obtener las soluciones
de transporte para cada uno de estos grupos, [8]. En simulaciones dependientes del
tiempo, se deben actualizar estos campos periódicamente, pero los pasos de tiempo
pueden ser largos comparados con los pasos de tiempo usados en la solución de
las ecuaciones de transporte, [17, 32, 45]. Cada vez que se actualiza la presión, las
saturaciones se interpolan de la malla numérica a las lı́neas de corriente y viceversa. Esto
puede introducir errores numéricos en el balance de masa, y posiblemente introducir
dispersión numérica.
Las lı́neas de corriente proveen imágenes instantáneas de los patrones de flujo en el
campo, que los ingenieros pueden utilizar para acortar el ciclo de ajuste histórico al
validar los modelos de los yacimientos, [5]. Las lı́neas de corriente también pueden
ayudar a los ingenieros a desarrollar estrategias de inyección y mejorar la eficiencia
del desplazamiento mediante el análisis de los patrones de flujo y estimando las
relaciones inyector-productor en el campo durante los cálculos de distintos escenarios de
predicción, [3]. Desarrollos recientes en métodos de lı́neas de corriente permiten ahora
para los simuladores de yacimientos basados en estos métodos, [6, 7], ser aplicados a
un conjunto más general de problemas, [9] que antes sólo podı́an ser resueltos usando
métodos convencionales como el de Volumen Finito.

16
 1.1. Estructura de la Tesis

Objetivo
Desarrollar herramientas computacionales que implementen el método de lı́neas
de corriente para simular el desplazamiento de aceite por agua en un yacimiento
petrolero y evaluar el desempeño del Método de Lı́neas de Corriente contra los métodos
tradicionales, particularmente, el Método de Volumen Finito.

1.1. Estructura de la Tesis

En el capı́tulo 2 se aborda el concepto de simulación de yacimientos, las distintas
maneras que existen de predecir su comportamiento, sus principales aplicaciones y la
terminologı́a usada en ingenierı́a de yacimientos, [1, 13, 36].

En el capı́tulo 3 entramos a la primera parte del proceso de modelado y simu-

lación, pues se hace la descripción del modelo conceptual del problema, [38, 25], es
decir, todas aquellas leyes y principios de la fı́sica que gobiernan el flujo de fluidos en
medios porosos, [14] y sus expresiones matemáticas.

Ya en el capı́tulo 4, se hace un desarrollo de los modelos matemáticos que describen

el desplazamiento un flujo de dos fases en tres dimensiones en un medio poroso. Di-
cho modelo matemático consta de ecuaciones diferenciales parciales que se deducen a
partir de la formulación axiomática, [2, 27] y la presentación del método de Lı́neas de
Corriente y su planteamiento matemático.

El capı́tulo 5 describe el Método de Volumen Finito [32, 19, 45] y su aplicación a

la discretización de las ecuaciones de presión y saturación 1D para el método de Lı́neas
de Corriente.

En el capı́tulo 6 encontraremos la implementación [39] de los algoritmos descritos

en el capı́tulo anterior a fin de resolver las ecuaciones discretas que conforman el modelo
matemático. Dicha implementación se llevó a cabo en el lenguaje de programación
C++ bajo la plataforma TUNAM, [18] ası́ como la visualización de los resultados con
OpenDx, [35]. Finalmente, el capı́tulo 7 consta de las conclusiones a las que llegamos
en este trabajo.

17
Capı́tulo 2

Simulación de Yacimientos

Un yacimiento petrolero es un medio poroso que contiene hidrocarburos. Las

dos caracterı́sticas importantes en un yacimiento son la naturaleza de la roca y de
los fluidos que contiene, [10, 30]. Un yacimiento es generalmente heterogéneo; sus
propiedades dependen en gran medida de la ubicación espacial. Por ejemplo, un
yacimiento fracturado es considerado heterogéneo ya que se compone de un conjunto de
bloques de medios porosos (la matriz rocosa) y una red de fracturas. Las propiedades
de las rocas de un reservorio de esta naturaleza suelen cambiar drásticamente; su
permeabilidad puede variar de un milidarcy (md) en la matriz a miles de md en las
fracturas. Mientras que las ecuaciones que gobiernan los yacimientos fracturados son
similares a los de un yacimiento ordinario, presentan dificultades adicionales que deben
ser consideradas, como la incorporación en los modelos de la distribución variable de
porosidad y permeabilidad, véase [36].
Los principales objetivos en la simulación de yacimientos son la determinación de la
reservas y la predicción de las tasas de recuperación de los yacimientos y encontrar
la manera de optimizar la recuperación de los hidrocarburos bajo diversas condiciones
operativas, prediciendo ası́ el desempeño a futuro del yacimiento, [20]. La simulación
de yacimientos es pues, el proceso de inferir el comportamiento real a partir del
comportamiento de un modelo. Los modelos pueden ser fı́sicos, tales como modelos
a escala de laboratorio o matemáticos.
Existen 4 etapas fuertemente ligas de modelado, [44], estableciendo primero
el modelo fı́sico, seguido del desarrollo de los modelos matemáticos,[22, 27],
posteriormente la discretización de dichos modelos [2, 4] y el diseño de algoritmos
computacionales, [21, 23, 47].

19
2. Simulación de Yacimientos

2.1. Métodos Clásicos en Ingenierı́a de Yacimientos

Los métodos clásicos de predicción de comportamiento de los yacimientos incluyen
métodos analógicos, experimentales y matemáticos. Los métodos analógicos utilizan las
caracterı́sticas de yacimientos maduros que son análogos a los del yacimiento objetivo
en un intento de predecir el rendimiento de una zona o de un yacimiento. Los métodos
experimentales miden propiedades fı́sicas, tales como la presión, la saturación y sus
relaciones en los núcleos de laboratorio y luego se busca ampliar su escala para la
acumulación de hidrocarburos en el yacimiento. Por último, los métodos matemáticos
usan sistemas de ecuaciones que suelen ser por lo general sistemas de ecuaciones
diferenciales parciales para pronosticar el desempeño del yacimiento, ver [13]

2.1.1. Métodos de Balance de Materiales

Los métodos matemáticos son los más utilizados en la simulación clásica de
yacimientos en la industria petrolera a fin de predecir el comportamiento del
yacimiento. Estos métodos incluyen el balance de materiales, curva de declinación,
estadı́sticos y métodos analı́ticos.
Los métodos clásicos de balance de materiales usan una representación matemática de
un yacimiento o volumen de drenaje. Su principio básico es la conservación de masa,
es decir, la cantidad de masa de agua, aceite o gas que queda en el depósito después de
un perı́odo de producción es igual a la diferencia de la cantidad de masa originalmente
en su lugar y que fue removido del yacimiento debido a la producción, más la cantidad
de masa añadida debido a la inyección.

2.1.2. Métodos de Curva de Declinación

Los métodos clásicos de curva de declinación usan uno de los tres declines
matemáticos (exponencial, hiperbólico y armónico) para describir la tasa de declinación
en la producción del aceite. Una curva de declinación tiene la forma general
1 dq
Cqb = − , (2.1)
q dt
donde C es un parámetro de declinación, q es la tasa de producción (m3 /d), y t
es el tiempo (dı́as). Los casos para cuando b = 0, 0 < b < 1 y b = 1 corresponden a
declinaciones exponenciales, hiperbólicas y armónicas, respectivamente.
Los métodos de la curva de declinación coinciden con los datos históricos de
producción para seleccionar una forma adecuada de la ecuación (2.1). Después de la
forma que se elija, los datos históricos son igualados por la elección de los parámetros C
y b que minimicen el error (utilizando a menudo el error de mı́nimos cuadrados) entre
los datos y la ecuación. Extrapolando los datos históricos en el futuro, se predice el
desempeño del yacimiento mediante la ecuación igualada. Un supuesto fundamental de
cualquier método de extrapolación es que todos los procesos que ocurren en el pasado
continuarán en el futuro.

20
 2.2. Métodos de Simulación de Yacimientos

2.1.3. Métodos Estadı́sticos

Los métodos estadı́sticos emplean correlaciones empı́ricas que son estadı́sticamente
obtenidos usando los resultados anteriores de algunos yacimientos para pronosticar
el desempeño futuro de los demás. Se trata de una generalización de los métodos
analógicos. Una correlación se desarrolla con datos de yacimientos maduros de la misma
región, con la misma litologı́a y bajo las mismas condiciones de funcionamiento. Para
tener confianza en el uso de un modelo de correlación empı́rica, las propiedades del
yacimiento deben estar dentro del lı́mite de la base de datos de regresión utilizado para
desarrollar un modelo. Los errores de predicción con los métodos estadı́sticos pueden
ser tan altos como del 20 al 50 %.

2.1.4. Métodos Analı́ticos

Los métodos de análisis, como el de Buckley-Leverett, utilizan la solución analı́tica
de un modelo matemático. El modelo consta de un conjunto de ecuaciones diferenciales
que describen el flujo y transporte de fluidos en un medio poroso, junto con un conjunto
adecuado de condiciones iniciales y/o de frontera. Para resolver estas ecuaciones
exactamente, se deban hacer suposiciones que simplifiquen el modelo y reduzcan
la complejidad del mismo. En general, estas suposiciones son muy restrictivas. Por
ejemplo, en el método de Buckley-Leverett para un flujo de dos fases se ignoran las
fuerzas de gravedad y capilares bajo la condición de incompresibilidad. Sin embargo,
dado que gran parte de la fı́sica de un problema se mantiene, los métodos analı́ticos
a menudo se utilizan para determinar cómo los diferentes parámetros influyen en el
rendimiento del yacimiento. Además, estos métodos pueden ser utilizados para validar
simuladores de yacimientos. Existen métodos más sofisticados tales como la separación
de variables, transformada de Laplace, y métodos integrales, [13].

2.2. Métodos de Simulación de Yacimientos

2.2.1. Etapas de la Simulación
La simulación de yacimientos involucra cuatro etapas principales interrelacionadas,
figura 2.1. En la primera, el modelo fı́sico, se incorporan y desarrollan los procesos
relevantes que intervienen en el fenómeno de estudio, en este caso en particular,
el flujo de fluidos en medios porosos. En el modelo matemático, involucramos un
conjunto de ecuaciones diferenciales parciales [33] no lineales y acopladas dependientes
del tiempo, analizando su existencia, convergencia y estabilidad. El modelo numérico
emplea esquemas numéricos adecuados de discretización, tales como diferencias finitas,
[31], volumen finito o elemento finito, integrando las propiedades básicas de las dos
primeras etapas; modelos fı́sico y matemático. En la cuarta etapa, en el modelo
computacional se desarrollan algoritmos computacionales y sus códigos a fin de resolver
eficientemente el sistema lineal y no lineal de ecuaciones algebraicas asociado a la

21
2. Simulación de Yacimientos

discretización numérica. Cada una de estas etapas es esencial para la simulación del
yacimiento, en ocasiones es necesario iterar un determinado número de veces alguna de
las etapas de la simulación a fin de ajustar los modelos fı́sicos, matemáticos, numéricos
y algoritmos computacionales para obtener un pronóstico preciso sobre el rendimiento
del yacimiento.
La creciente aceptación de la simulación de yacimientos puede atribuirse a los avances
en el desarrollo computacional, modelos matemáticos, métodos numéricos, técnicas de
solución y herramientas de visualización cientı́fica.

Figura 2.1: Etapas del proceso de modelado y simulación.

2.2.2. Clasificación de Simuladores de Yacimientos

Los simuladores de yacimientos se pueden clasificar de acuerdo a diferentes enfo-
ques. Los más comunes se basan en el tipo de fluidos contenidos en el yacimiento de
estudio y los procesos de recuperación que están siendo modelados. Otros enfoques
incluyen el número de dimensiones (1D, 2D y 3D), el número de fases (monofásicos,
bifásicos y trifásicos), y el sistema coordenado utilizado en el modelo (rectangulares,
cilı́ndricos y esféricos). Los simuladores también pueden ser determinados por el tipo
de estructura de la roca o su respuesta (ordinarios, doble porosidad/permeabilidad,
acoplamiento hidráulico/térmico y fracturados).
Los simuladores basados en la clasificación del tipo de fluidos en el yacimiento incluyen
gas, petróleo negro, y composicionales. Los simuladores de petróleo negro suelen ser us-
ados para recuperaciones convencionales pues los procesos simulados no son sensibles a
cambios composicionales en los fluidos del yacimiento. Los simuladores composicionales
son usados cuando los procesos de recuperación de hidrocarburos son sensibles al ago-
tamiento primario de aceite volátil y gas condensado, en operaciones de mantenimiento
de presión en los yacimientos y en múltiples procesos miscibles.

Los Simuladores de yacimientos clasificados de acuerdo a los procesos de

recuperación incluyen recuperación convencional (petróleo negro), desplazamiento

22
 2.2. Métodos de Simulación de Yacimientos

miscible, recuperación térmica e incorporación de quı́micos. Los procesos primarios

de recuperación de petróleo, agua, gas en solución, expansión por gas, drenaje por
gravedad e imbibición capilar, pueden ser modelados con simuladores de petróleo negro.
Además, las etapas de recuperación secundaria, tales como la inyección de agua o gas,
pueden también ser modelado con estos simuladores. Existen mecanismos térmicos
de recuperación que implican inyección de vapor o la combustión in situ, y el uso
de ecuaciones tales como las de conservación de masa y energı́a. Los simuladores de
adición de quı́micos incluyen inyección de alcalinos, polı́meros, y/o espuma y se pueden
utilizar para cambiar la relación de movilidad de forma dinámica en el desplazamiento
(polı́mero) o movilizar el aceite residual (surfactantes). Se deben considerar otros
efectos tales como la adsorción en la roca, la reducción de la permeabilidad y los
fluidos no newtonianos.

2.2.3. Aplicaciones de la Simulación de Yacimientos

La simulación de yacimientos es usualmente aplicada en los siguientes pasos:

Establecer los objetivos de estudio de la simulación. El primer paso en cualquier

estudio de simulación de yacimientos es fijar objetivos claros. Estos objetivos
deben ser alcanzables y compatibles con yacimientos disponibles y datos de
producción.

Recopilar y validar los datos de yacimientos. Después de que los objetivos

de la simulación se han establecido, se deben obtener datos del yacimiento
y la producción. Los datos que cumplan los objetivos son incorporados en el
simulador.

Diseñar el simulador. Una vez que los datos son recogidos y validados, se diseña
el simulador. Este paso implica las cuatro etapas principales interrelacionados
descritas arriba: la construcción de un modelo fı́sico conceptual, el desarrollo
de modelos matemáticos y modelos numéricos ası́ como el diseño de códigos
computacionales.

Validación histórica del simulador. Una ves que el simulador es construido, este
debe ser calibrado con datos disponibles de producción ya que gran parte de los
datos en un simulador tı́pico necesitan ser verificados.

Hacer predicciones. En la etapa final de aplicación, se evalúan diversos planes de

desarrollo y producción llevándose a cabo un análisis de sensibilidad de diversos
parámetros de producción en el yacimiento.

Mientras que la simulación de yacimientos es el método más completo, los métodos

clásicos de ingenierı́a de yacimientos se encuentran todavı́a en uso para predecir el
comportamiento de los yacimientos. Estos métodos clásicos se pueden utilizar para
generar datos de entrada para simuladores de yacimientos. Por ejemplo, un análisis de

23
2. Simulación de Yacimientos

acumulación de presión se puede utilizar, en la caracterización del yacimiento, para

obtener la permeabilidad del dominio de estudio, mientras que los métodos de balance
de materiales proporcionan información sobre la intrusión de agua y el tamaño del
acuı́fero durante la validación histórica.

2.3. Términos Usados en la Simulación Numérica

Método numérico.- Un método numérico para resolver problemas de ecuaciones
diferenciales consiste en discretizar el problema, que tiene un número infinito de
grados de libertad, para producir un problema discreto, que tiene un número
finito de grados de libertad y pueda ser resuelto utilizando una computadora.
Existen diferentes métodos numéricos, entre ellos diferencias finitas, volúmenes
finitos, y el método de elementos finitos.

Estructura de malla.- Es la geometrı́a de una malla que se utiliza para

la simulación numérica de un yacimiento. Esta puede ser cartesiana, radial,
logarı́tmica o distorsionada y 1D, 2D o 3D.

Malla 2D.- Es una estructura hecha para mirar hacia abajo en el yacimiento.
Para un sistema de coordenadas cartesianas, es una división del yacimiento en
las direcciones x1 − x2 y utilizando pasos espaciales h1 y h2 , figura 2.2.

Figura 2.2: Malla para un área en 2D.

Modelo en sección transversal 2D.- Es una estructura que se impone a un

corte vertical a través del yacimiento. Para un sistema cartesiano, es una división
del yacimiento en las direcciones x1 y x3 utilizando pasos espaciales h1 y h3 , como
se muestra en la figura 2.3. Este tipo de modelos se emplean para evaluar el
efecto de la estratificación vertical en el yacimiento.

24
 2.3. Términos Usados en la Simulación Numérica

Figura 2.3: Corte transversal para un dominio 2D.

Transmisibilidad.- La transmisibilidad entre dos bloques adyacentes de una

malla mide la facilidad con la que fluye un fluido entre ellos.

Driscretización espacial.- Se refiere al proceso de dividir el dominio de estudio

en pequeños subdominios con pasos espaciales h1 , h2 y h3 para después modelar el
flujo a través de un método numérico. En la simulación numérica de yacimientos,
siempre se divide el yacimiento en bloques de malla y luego se modela el flujo de
fluidos entre los bloques.

Discretización temporal.- Se refiere al proceso de dividir un intervalo de

tiempo de interés en subintervalos con pasos de tiempo ∆t y avanzar la simulación
en tiempo.

Dispersión numérica.- Dispersión numérica es la propagación de un frente

de inyección en un proceso de desplazamiento, como sucede en la inyección
de agua en medios porosos. Este fenómeno es debido a efectos numéricos. En
concreto, es debido a la discretización en tiempo y espacio o a consecuencia
de un error de truncamiento que surge del mallado del dominio. Este frente de
propagación tiende a conducir a la irrupción temprana de agua y otros errores
en la recuperación. La gravedad del error depende del proceso de recuperación
de fluidos que esta siendo simulado (por ejemplo, inyección de agua e inyección
alternada de agua y gas), del los pasos de tiempo y espacio ası́ como de los
métodos numéricos utilizados.

Conservación de masa.- Es un principio general que se utiliza para verificar la

exactitud de un método numérico en la simulación de yacimientos. Se limita a lo
siguiente:

(masa dentro del bloque) - (masa fuera del bloque) = acumulación de masa
dentro del bloque.
Los modelos de simulación de yacimientos se componen básicamente de la
conservación de masa y la ley de Darcy, en relación a la velocidad del fluido a un

25
2. Simulación de Yacimientos

gradiente de presión. En los métodos térmicos también se añade la conservación

de energı́a. El Balance de materiales es un término usado en ingenierı́a para la
conservación de masa en un volumen fijo, que es normalmente un yacimiento de
hidrocarburos.

26
Capı́tulo 3

Modelo Conceptual

3.1. Formación del Aceite y Gas

El Petróleo es una mezcla de hidrocarburos que en forma natural se encuentran en

la corteza terrestre como gas, lı́quido o sólido y puede existir en una o varias fases en el
mismo lugar. Contiene cantidades menores de Nitrógeno, Oxı́geno, Sodio, Azufre, etc.,
como impurezas, ver [30].
Un yacimiento es la acumulación natural en la corteza terrestre de aceite y/o gas de la
misma composición, comprendida en los mismos lı́mites y sometida a un mismo sistema
de presión en una trampa petrolera. La roca generadora debe estar enterrada a una
profundidad suficiente, generalmente mayor a 1000 m, para que la materia orgánica
contenida pueda madurar hasta convertirse en aceite y/o gas. Es necesario que la roca
generadora se encuentre dentro de una Cuenca Sedimentaria que sufra procesos de
subsidencia (hundimiento por su propio peso) y enterramiento, con un aporte suficiente
de sedimentos, figura 3.1.
Dentro de la roca generadora, no toda la materia orgánica se transforma en petróleo, se
estima que el 70 % permanece como residuo orgánico insoluble, por lo que el rendimiento
promedio de las rocas generadoras es de aproximadamente 30 %. Sin embargo, este
porcentaje no es el petróleo que finalmente obtenemos, pues se estima que sólo el 1 %
del petróleo generado es capaz de migrar hacia la roca almacén y acumularse en ella,
mientras que el 99 % no llega a migrar o se pierde debido a que no existe un sello que
impida que el crudo o el gas escape de la roca almacén.
Por otra parte, se tiene el problema de la cantidad de petróleo recuperable con
rendimiento económico de los yacimientos, por lo general menor al 60 %.

27
3. Modelo Conceptual

Figura 3.1: Depositación de la materia orgánica en ambientes reductores y su posterior

sepultamiento en la cuenca sedimentaria.

3.1.1. Transformación de la Materia Orgánica

La materia orgánica acumulada en ambientes reductores, es decir, aquellos que
favorecen su preservación, se cubren por sepultamiento dando lugar a una serie de
cambios junto con los sedimentos que contienen a dicho material orgánico. Estos
procesos se llaman Diagénesis, Catagénesis y Metagénesis, ver [30].

Diagénesis.- Es el proceso mediante el cual los compuestos orgánicos consti-

tuyentes de los seres vivos, tales como carbohidratos, proteı́nas, etc., son someti-
dos a un ataque microbiano, que se realiza a poca profundidad (con presiones
litostáticas entre cero y 300 bares) y bajas temperaturas (entre 0◦ C y 50◦ C). El
hidrocarburo generado durante esta etapa es el metano y compuestos como el CO2
y H2 O. En esta etapa, se presenta generalmente la consolidación del sedimento,
es decir, las fracciones sueltas se convierten en rocas sedimentarias y la mayor
parte de la materia orgánica que se conserva se transforma en kerógeno, que es
la fracción insoluble y en menor proporción se forma betumen que corresponde a
la parte soluble.

Catagénesis.- Los sedimentos con materia orgánica se sepultan rápida o

lentamente en función de las caracterı́sticas propias de la cuenca sedimentaria,
de la taza de sedimentación y de su entorno. Cuando la roca generadora
alcanza profundidades mayores a 1 km inicia la catagénesis, es decir, inicia la
ventana de generación debido al incremento en la presión y la temperatura. Las
temperaturas que se alcanzan en esta etapa son del orden de 50◦ C y hasta 225◦ C
aproximadamente, mientras que la presión varı́a de 300 a 1500 bares. A los 2.6 km
de profundidad y 100◦ C se alcanza el máximo pico de generación de hidrocarburos

28
 3.2. Sistema Petrolero

lı́quidos y termina a los 175◦ C. Entre los 3 y 3.5 km se produce menos aceite y
más gas.

Metagénesis.- Es la tercera y última etapa en la transformación de la materia

orgánica. Ésta se desarrolla a temperaturas mayores a los 225◦ C, siendo aquı́ la
mayor generación de gas. La generación de metano acaba a los 315◦ C, con
profundidades cercanas a los 8 km, es decir, presiones litostáticas mayores a 1500
bares. La porosidad de las rocas en estas condiciones disminuye notablemente, por
lo que es difı́cil que se formen a estas profundidades yacimientos con rendimiento
económico. Cuando el sepultamiento es mayor a los 10 km, el kerógeno residual
se transforma en grafito y es imposible considerar la producción aún mı́nima de
hidrocarburos gaseosos.

3.2. Sistema Petrolero

El llamado Sistema Petrolero constituye un sistema natural que incluye elementos y
procesos geológicos que intervienen en la formación de un yacimiento de hidrocarburos.
El sistema petrolero es un modelo dinámico que estudia las entradas de materia
orgánica a la cuenca sedimentaria, la transformación de dicha materia, la generación
de hidrocarburos y su acumulación en una trampa petrolera. Por otra parte, el sistema
petrolero esta compuesto por los siguientes subsistemas que deben estar concatenados
en tiempo y espacio para que se forme una acumulación natural de petróleo en la
corteza terrestre, susceptible de ser explotada con rendimiento económico, ver [12].

Roca Generadora.- Éste tipo de roca debe ser rica en materia orgánica
y preferentemente de gran espesor, a condiciones adecuadas de presión y
temperatura se forman el aceite y el gas, los cuales se acumulan en la trampa
petrolera.

Migración.- Es el movimiento de los hidrocarburos en los poros o a través de las

discontinuidades de las rocas, tales como fallas y fracturas. Éste proceso describe
el desplazamiento desde la roca generadora hasta su acumulación en la trampa
petrolera.

Roca Almacén.- Éste tipo de roca debe ser porosa, permeable y tener
continuidad lateral y vertical. Las rocas almacenadoras permiten el flujo de fluidos
dentro de ellas debido a la intercomunicación de los poros. Se dice que una roca
tiene una permeabilidad adecuada para permitir el paso de los hidrocarburos si
posee poros interconectados y de tamaño supercapilar.

Roca Sello.- Para que los hidrocarburos puedan quedar confinados en las rocas
almacenadoras, es necesario que las paredes del depósito estén selladas de manera
efectiva. Éste tipo de rocas deben contar con escasa permeabilidad o, por contener

29
3. Modelo Conceptual

poros de tamaño subcapilar, no permiten el paso del petróleo, sirviendo como

cierre en su migración o desplazamiento. Debido a que los yacimientos petroleros
están asociados a zonas de actividad tectónica, las rocas sello, deben tener
comportamiento plástico, de manera que respondan a los esfuerzos mecánicos
deformándose en el campo dúctil, formando pliegues y no fracturándose. El
espesor de la roca sello es variable, siendo reducido si tiene una excelente calidad
o de espesor mediano o grueso si la calidad de la roca es mediana o mala.

Entrampamiento.- La trampa petrolera es una caracterı́stica geológica que

permite que el aceite y el gas se acumulen y conserven de manera natural
durante un cierto periodo de tiempo. Se tratan de receptáculos cerrados en la
corteza terrestre que cuentan con rocas almacenadoras y rocas sello en posición
tal que permiten se acumulen los hidrocarburos. Las trampas petroleras tienen
una determinada forma, tamaño y geometrı́a.

Sincronı́a.- Involucra la sincronización en tiempo y espacio de los elementos

anteriormente citados, la figura 3.2 muestra estos elementos coexistiendo en
sincronı́a. Algunas trampas presentan caracterı́sticas adecuadas para almacenar
hidrocarburos, con buena relación entre la roca almacenadora y la roca sello;
sin embargo se encuentran vacı́as, siendo las principales causas: ausencia de roca
generadora, los hidrocarburos no alcanzaron la trampa, el petróleo se destruyó o
bien, la trampa se formó tardı́amente.

Figura 3.2: Trampa tı́pica con sincronı́a de los distintos elementos del Sistema Petrolero.

3.3. Propiedades de la Roca y los Fluidos en el

Yacimiento
3.3.1. Propiedades de la Roca
Poros.- Son pequeños pasajes interconectados existentes en una roca permeable.
Las conexiones entre los poros se conocen como “garganta de poro” y son éstas

30
 3.3. Propiedades de la Roca y los Fluidos en el Yacimiento

las que controlan la presión capilar de entrada en un proceso de drenaje. Sus

dimensiones van desde 1 a 200 µm.
Porosidad.- La porosidad es una medida de la capacidad de almacenamiento de
fluidos que posee una roca y se define como el porcentaje del volumen poroso
de la roca respecto al volumen total de la misma. De acuerdo a su conectividad,
figura 3.3, la porosidad puede clasificarse como porosidad total, la cual incluye
a los poros interconectados y a los que se encuentran aislados; y la porosidad
efectiva, que incluye solo los poros interconectados y es en realidad la que
interesa para la estimación del hidrocarburo en sitio.
La porosidad total la denotamos como:
V olumen del espacio del poro
φ=
V olumen representativo

y la porosidad efectiva
V olumen del espacio del poro disponible
φe =
V olumen representativo
en general
φe ≤ φ
La porosidad está en función de la presión debido a la compresibilidad de la roca,
la cual comúnmente se asume constante y presenta valores que oscilan entre 10−6
y los 10−7 psi−1 . La compresibilidad de la roca la expresamos como
1 ∂φ
CR =
φ ∂p
integrando es posible expresar a la porosidad de la siguiente manera
0)
φ = φ0 eCR (p−p
donde φ0 es la porosidad a la presión de referencia p0 , la cual comúnmente es la
atmosférica. Podemos hacer una expansión en series de Tylor y despreciando los
términos de alto orden para una roca ligeramente compresible obtenemos
φ ≈ φ0 1 + CR (p − p0 ) .

Figura 3.3: A la izquierda, observamos poros interconectados y poros aislados a la

derecha.

31
3. Modelo Conceptual

Conductividad hidráulica.- Las interacciones entre fluidos y matriz porosa se

relacionan a través de la conductividad hidráulica que expresamos

ρα g
K = krα k .
µα

En esta relación, k es la permeabilidad intrı́nseca de la roca y krα es la

permeabilidad relativa de la fase α. La viscosidad queda denotada como µα .
La densidad la denotamos con ρ, que es el cociente de la masa de un cuerpo entre
su volumen. La gravedad esta expresada por g.

Permeabilidad.- La permeabilidad de una roca es la capacidad que conducir

fluidos a través de sus poros interconectados, se le conoce también como
permeabilidad absoluta y se mide en mili-darcy (md). La permeabilidad se
encuentra muy ligada a la porosidad, figura 3.4. A continuación se muestra una
tabla con distintos valores de permeabilidad para rocas en yacimientos petroleros.

Clasificación Rango de Permeabilidad en (md)

Pobre 1-15
Moderada 15-20
Buena 50-250
Muy buena 250-1000
Excelente mayores a 1000
Tabla 3.1: Clasificación de permeabilidad de rocas

Generalmente los sistemas continuos son anisótropos, para un caso 3D y donde

el sistema coordenado coincide con las direcciones de flujo se tiene entonces que

 
k11 0 0
k =  0 k22 0 
0 0 k33

Si el medio es isotrópico: k11 = k22 = k33 .

32
 3.3. Propiedades de la Roca y los Fluidos en el Yacimiento

Figura 3.4: Correlación permeabilidad-porosidad

3.3.2. Propiedades de los Fluidos

Fase.- Es la región quı́micamente homogénea de un fluido que se separa de otras
regiones por una interfase. Las faces que generalmente se emplean en simulaciones
de yacimientos son aceite (o), agua (w) y gas (g).

Componente.- Es una especie quı́mica que puede estar contenida en una fase.
Por ejemplo, la fase acuosa (w) contiene componentes agua (H2 O), cloruro de
sodio (NaCl), y oxı́geno disuelto (O2 ), mientras que la fase aceite puede contener
cientos de componentes por ejemplo Nitrógeno, Oxı́geno, Sodio, Azufre, etc.

Tipos de fluidos en el yacimiento.- En general, agua, aceite y gas pueden exis-

tir simultáneamente en un yacimiento petrolero. Estos fluidos pueden ser clasifica-
dos como incompresibles, ligeramente compresibles y compresibles, dependiendo
de como responden a la presión, figura 3.5. Un fluido incompresible es aquel
que tiene compresibilidad cero, es decir, su densidad no depende de la presión. El
agua y el aceite sin gas disuelto se consideran incompresibles. Un fluido ligera-
mente compresible es aquel que presenta una pequeña compresibilidad que se
mantiene constante, tı́picamente se encuentra en el rango de 10−5 a 10−6 psi−1 , el
agua y el aceite sin gas, a condiciones de yacimiento, pueden ser considerados lig-
eramente compresibles. Un fluido compresible tiene una compresibilidad tı́pica
en el rango de 10−3 a 10−4 psi−1 ; de modo que la densidad incrementa a medida
que incrementa la presión pero se estabiliza a presiones altas. A condiciones de
yacimiento el gas es compresible.

33
3. Modelo Conceptual

Figura 3.5: Relación densidad-presión.

Compresibilidad.- Se define como el cambio de volumen (V ) o densidad (ρ) en

función de la presión para una temperatura (T ) dada.
1 ∂V 1 ∂ρ
cf = − = −
V ∂p T ρ ∂p T

integrando la ecuación anterior tenemos que la densidad se expresa como

0)
ρ = ρ0 ecf (p−p
donde ρ0 es la densidad a una presión de referencia p0 . Usando una expansión
en series de Taylor y truncando los términos no lineales obtenemos una buena
aproximación para la densidad en función de la presión para fluidos ligeramente
compresibles:
p ≈ ρ0 1 + cf (p − p0 ) .

Factor de solubilidad del gas.- El factor de solubilidad del gas Rso es el

volumen de gas, medido a condiciones estándar, disuelto a presión y temperatura
del yacimiento por unidad de almacenamiento de aceite.
VGs
Rso = (3.1)
VOs
El subı́ndice s denota que el volumen es a condiciones estándar, mientras que
las letras G y O representan las componentes gas y aceite respectivamente.
Usualmente las unidades se expresan en SCF/STB (standard cubic feet / stock
tank barrels). Nótese que
WO WG
VOs = , VGs =
ρOs ρGs

34
 3.3. Propiedades de la Roca y los Fluidos en el Yacimiento

de modo que la ecuación (3.1) se convierte en

WG ρOs
Rso = .
WO ρGs

Factor de formación de volumen.- Describe la razón del volumen V de una

fase medida a condiciones de yacimiento, entre el volumen Vs de la fase medida
a condiciones estándar. Sus unidades están dadas por RB/ST B para lı́quidos
y RB/SCF para gases, donde RB es ”reservoir barrels”. Para una fase (α), el
factor en función de la densidad es:
ρα,s
Bα (p, T ) = .
ρα
Para el aceite que tiene gas disuelto observemos que
WO + WG
Vo =
ρ0
de manera que el factor de formación de volumen esta dado por

(WO + WG )ρOs
Bo = .
WO ρo

Densidad del fluido.- La densidad de una fase (agua, aceite o gas) está dada
por
ρs
ρ= .
B
Las fracciones de masa de aceite y gas disuelto en aceite son respectivamente:
WO ρOs
COo = = ,
WO + WG Bo ρo
WO Rso ρGs
CGo = = .
WG + WG Bo ρo
Tomando COo + CGo = 1 obtenemos la densidad de la fase aceite
Rso ρGs + ρOs
ρo = .
Bo

Viscosidad.-Se denota con la letra µ y es una medida de la energı́a disipada

cuando el fluido está en movimiento resistiendo una fuerza de corte aplicada,
sus dimensiones son (fuerza/área · tiempo) y sus unidades Pa·s = poise. En un
fluido gaseoso, las moléculas están muy separadas y presentan baja resistencia
a fluir a consecuencia del comportamiento aleatorio de las mismas. Por otra
parte, un fluido denso presenta gran resistencia a fluir debido a la naturaleza

35
3. Modelo Conceptual

cercana de sus moléculas. La viscosidad del agua a condiciones estándar es de 1

cp (centipoise). En general, la viscosidad depende de la presión, la temperatura
y de las componentes de cada fase. En la tabla 3.2 se muestran las viscosidades
tı́picas a condiciones de yacimiento (4000-6000 psi y 200◦ F) de los diferentes tipos
de aceite.

Clasificación Rango de Viscosidad en (cp)

Aceite ligero 0.3 - 1
Aceite intermedio 1-6
Aceite moderado 6 - 50
Aceite muy viscoso 50 - 1000
Aceite pesado mayores a 1000
Tabla 3.2: Valores tı́picos de viscosidad en aceites.

3.4. Mojabilidad
El desempeño de un yacimiento se ve afectado por el hecho de que la roca sea
preferencialmente mojable por agua o por aceite, particularmente en las técnicas de
inyección de agua y recuperación mejorada de hidrocarburos. La mojabilidad es la
preferencia de un sólido por estar en contacto con un fluido en lugar de otro. Una gota
de un fluido preferentemente mojante va a desplazar a otro fluido dispersándose por
la superficie, por el contrario un fluido no mojante formará gotas, disminuyendo su
contacto con la superficie del sólido, figura 3.6.
Mojabilidad por agua.- La fase agua es la mojadora y forma una pelı́cula sobre
las paredes del poro aún incluso en aquellos que contienen aceite.

Mojabilidad por aceite.- La fase mojadora es el aceite y forma una capa sobre
la superficie de la roca, aún en los poros que contienen agua.

Mojabilidad intermedia.- Las fases agua y aceite son mojantes hasta cierto
grado en la matriz porosa. El equilibrio de estos casos creará un ángulo de
contacto θ entre los fluidos de la superficie que está determinado por el equilibrio
de fuerzas resultante de la interacción de las tensiones intersticiales.

Figura 3.6: A la izquierda un fluido no mojador, al centro una fluido de mojabilidad

intermedia y a la derecha un fluido mojador.

36
 3.5. Procesos de Desplazamiento de Fluidos

3.5. Procesos de Desplazamiento de Fluidos

Imbibición.- Es el proceso de desplazamiento de un fluido que ocurre cuando la
fase mojadora se incrementa.

Drenaje.- Proceso de desplazamiento de un fluido cuando la fase no mojadora

se incrementa.

Imbibición espontánea.- Este proceso ocurre cuando una fase mojadora invade
un medio poroso en ausencia de fuerzas externas. La fase mojadora invade bajo
la acción de fuerzas superficiales. La figura 3.7 muestra los procesos de drenaje
e imbibición.

Figura 3.7: Procesos de drenajes e imbibición

3.6. Interacción Roca-Fluidos

Saturación.- La saturación (S) de una fase ya sea agua, aceite o gas, se define
como la fracción de el espacio de poro que ésta ocupa, de esta manera, en un
sistema trifásico tenemos que

Sw + So + Sg = 1

La expresión anterior significa que las tres fases llenan completamente el espacio
poroso. La presión capilar, la permeabilidad relativa entre otras, dependen
fuertemente de la saturación.

Saturación residual.- Denotada por Sαr , la saturación residual de una

determinada fase α es la cantidad de dicha fracción que queda atrapada en la
matriz porosa o es irreducible. La fase no mojadora residual es atrapada en los
poros por fuerzas capilares. Sin embargo, la cantidad de fluido atrapado depende
de la permeabilidad y mojabilidad de la roca.

Presión capilar.- Cuando dos fluidos inmiscibles están en contacto dentro de

los poros, una superficie curvada se forma entre los dos. Por ejemplo, para las
fases agua y aceite, la presión en el lado del fluido no-mojante de la interfase (po ),

37
3. Modelo Conceptual

es mayor que la presión para el lado del fluido mojante (pw ). Esta diferencia de
presiones se define como presión capilar (pc ).
pc = po − pw .
La presión capilar depende de la saturación de la fase mojadora y de la dirección
de cambio de ésta (drenaje o imbibición), figura 3.8. Cabe señalar que la pc
depende también de la tensión superficial σ, la porosidad φ, la permeabilidad k
y el ángulo θ de contacto con la superficie de la roca de la fase mojadora, el cual
depende a su vez de la temperatura y de la composición del fluido. En un flujo
de tres fases, se requieren de dos presiones capilares, a saber:
pcow = po − pw , pcgo = pg − po ,
es posible obtener una tercera expresión de presión capilar haciendo
pcgw = pg − pw = pcow + pcgo .

Figura 3.8: En rojo se muestra la curva de drenaje primario e imbibición en negro y

delimitan el comportamiento de la presión capilar. En verde pc para una saturación
intermedia, en amarillo una inversión de los valores para la lı́nea verde.

Permeabilidad relativa.- Denotada por krα , mide la habilidad de una fase para
fluir en una formación porosa y en presencia de otras fases, pues la presencia de
más de una fase en un sistema inhibe el flujo.

Movilidad.- Se define como la razón de la permeabilidad relativa entre su

viscosidad. En un sistema trifásico (w, o, g) tenemos
krw kro krg
λw = , λo = , λg = .
µw µo µg

38
 3.7. Etapas de la Recuperación de Hidrocarburos

y denotaremos a la movilidad total como

λ = λw + λo + λg .

Flujo fraccional.- Determina la razón de flujo volumétrico fraccional de una

fase bajo un gradiente de presiones dado, en presencia de otra fase:
λw λo λg
fw = , fo = , fg = ,
λ λ λ

3.7. Etapas de la Recuperación de Hidrocarburos

La explotación de un yacimiento de petróleo ocurre básicamente en tres etapas.
En las dos primeras etapas se logra recuperar en promedio del 25 % a 30 % del crudo,
con lo cual el yacimiento contiene todavı́a un estimado de 60-80 % de hidrocarburos,
quedando atrapado en los poros de la estructura del reservorio debido a la viscosidad
y efectos de capilaridad.

Recuperación Primaria.- En esta etapa se aprovecha la presión natural del

yacimiento que lleva los hidrocarburos hasta la superficie debido a la diferencia
de presión entre el yacimiento y la presión atmosférica. Cuando la presión del
medio se hace inadecuada, o cuando se están produciendo cantidades importantes
de otros fluidos (agua y gas, por ejemplo) termina la primera etapa. La taza de
recuperación durante la fase primaria es del 12-15 % de los hidrocarburos en el
yacimiento.

Recuperación Secundaria.- Es toda actividad encaminada a la recuperación de

hidrocarburos adicional a la que se obtendrı́a con la energı́a propia del yacimiento
(producción primaria). Consiste en inyectar dentro del yacimiento un fluido menos
costoso que el crudo para mantener un gradiente de presión adecuado de modo que
la producción vuelva a ser económicamente rentable. La recuperación secundaria
básicamente consiste en la inyección de agua (figura 3.9) en el yacimiento o
la inyección de un gas natural en la cima de la estructura, con el propósito
fundamental de mantener la presión, o bien, de desplazar los hidrocarburos de la
zona de aceite, mediante arreglos especı́ficos de pozos inyectores y productores.
En esta etapa se produce alrededor del 15 al 20 % adicional del petróleo.

Recuperación Terciaria o Mejorada.- Después de las recuperaciones primaria

y secundaria, el yacimiento contiene todavı́a un estimado del 60 % del crudo.
Numerosos métodos han sido estudiados para la recuperación, al menos parcial,
de estas grandes cantidades de crudo remanente en los pozos. Los procesos
de Recuperación Mejorada surgen como una alternativa para incrementar la
recuperación de hidrocarburos, modificando las caracterı́sticas de los fluidos
y las fuerzas capilares que actúan sobre ellos. La Recuperación Mejorada se

39
3. Modelo Conceptual

Figura 3.9: Inyección de agua en yacimientos petroleros.

fundamenta principalmente en técnicas sofisticadas en la operación; suele ser

de alto costo, pero muy efectivas, ası́ pues, la Recuperación Mejorada de
hidrocarburos se define como la producción de aceite, mediante la inyección de un
fluido que, además de desplazar el aceite, modifica favorablemente los mecanismos
de recuperación de hidrocarburos. Existen otros métodos pertenecientes a la
tercera fase de recuperación con aditivos quı́micos, sin embargo en ocasiones han
sido desechados principalmente argumentando la baja rentabilidad del proceso,
debido principalmente a los costos de los aditivos. También, bajo condiciones
óptimas, una solución de surfactantes inyectada al reservorio tiene el potencial de
solubilizar el crudo, dispersándolo de manera efectiva en forma de una emulsión.
Durante esta etapa el yacimiento entrega un rendimiento cercano el 70 %.

40
Capı́tulo 4

Modelos Matemáticos

4.1. Formulación Axiomática

En este capı́tulo se obtendrán las ecuaciones que gobiernan el fenómeno de
transporte de fluidos a través de una matriz porosa aplicando la formulación
axiomática, véase [27]. Dicha formulación consiste en identificar las propiedades
intensivas y extensivas del fenómeno de estudio aplicando balances de ellas para
un volumen que denotaremos como B(t) y que representa un sistema continuo. Las
propiedades extensivas son aquellas funciones de variable escalar o vectorial que
podemos representar por medio de una integral de cuerpo,
Z
E(B(t), t) = ψ(x, t) dx,
B(t)

donde x representa la posición y t el tiempo. El integrando de la expresión anterior es

la propiedad intensiva representada por ψ(x, t). Esta ecuación integral establece una
correspondencia biunı́voca, es decir, a cada propiedad intensiva sobre un dominio que
ocupa cualquier cuerpo B(t) le corresponde una y solo una propiedad extensiva. El
cambio temporal en una propiedad extensiva E es debido a su generación dentro del
sistema continuo o a que se importa por la frontera, esta relación se escribe como:
Z Z Z
dE d
= ψ(x, t) dx = q(x, t) dx + τ (x, t) · n dS, (4.1)
dt dt
B(t) B(t) S(t)

donde q(x, t) y τ (x, t) son la generación y el vector de flujo de la propiedad extensiva

respectivamente. Una representación de un medio continuo B(t) se muestra en la
figura 4.1, donde n es el vector normal a la superficie S(t), y ds es una diferencial de
superficie.

41
4. Modelos Matemáticos

Figura 4.1: Representación esquemática de un medio continuo

La ecuación (4.1) se satisface para cada cuerpo B(t) de un sistema. Aplicando varios
resultados matemáticos a la ecuación (4.1) se obtiene:
Z   Z Z
∂ψ
+ ∇ · (vψ) dx = q dx + ∇ · τ dx. (4.2)
∂t
B(t) B(t) B(t)

La expresión anterior trae como consecuencia que podamos escribir una ecuación
diferencial de balance local para la propiedad extensiva ψ,
∂ψ
+ ∇ · (vψ) = q + ∇ · τ . (4.3)
∂t

4.2. Forma Conservativa de las Ecuaciones de

Balance
Las ecuaciones de balance abordadas anteriormente pueden ser escritas en su forma
conservativa. Se requiere definir una función de flujo F = vψ − τ , de modo que ahora
la ecuación (4.2) la expresamos como:
Z Z Z
∂ψ
dx + ∇ · F dx = q dx, (4.4)
∂t
B(t) B(t) B(t)

lo que nos lleva a la siguiente ecuación diferencial

∂ψ
+ ∇ · F = q. (4.5)
∂t
Las ecuaciones (4.4) y (4.5) son la forma conservativa de las ecuaciones (4.2) y
(4.3) respectivamente. El teorema de Gauss puede ser aplicado a la ecuación (4.4),
obteniendo: Z Z Z
∂ψ
dx + F · n dS = q dx, (4.6)
∂t
B(t) S(t) B(t)

42
 4.3. Modelo Matemático para Sistemas Multifásicos

4.3. Modelo Matemático para Sistemas Multifásicos

Al igual que en el modelo de una fase, en los sistemas multifásicos la propiedad
extensiva es la masa de cada una de las fases α, y lo denotamos como:
Z
E= φρα Sα dx ,
B(t)

para una porosidad del medio φ y donde ρα y Sα son la densidad y la saturación de la

fase α (agua (w), aceite (o) o gas (g)) respectivamente. La propiedad intensiva asociada
al sistema es
ψα = φρα Sα .
Para un sistema donde no hay difusión (τ = 0,) y aplicando la ecuación (4.3),
obtenemos la ecuación de balance local para la masa de fluido de la fase α, que
escribimos como:
∂(φρα Sα )
+ ∇ · (vφρα Sα ) = qα .
∂t
La velocidad de Darcy para una fase α se denota por v α = vφSα , de modo que la
ecuación de balance de masa para la fase α, que forma parte de un sistema multifásico
se expresa como:
∂(φρα Sα )
+ ∇ · (ρα v α ) = qα . (4.7)
∂t
La Ley de Darcy para sistemas multifásicos es:
kkrα
vα = − (∇pα − ρα g) , (4.8)
µα
donde krα es la permeabilidad relativa de la fase α. La viscosidad se denota por
µα , la presión es pα y la densidad se expresa por ρα . Sustituyendo la ecuación (4.8) en
(4.7) y definiendo la movilidad de la fase α como:
krα
λα = ,
µα
llegamos a:
∂(φρα Sα )
− ∇ · (ρα λα k(∇pα − ρα g)) = qα . (4.9)
∂t
La ecuación (4.9) es un sı́ un conjunto de ecuaciones diferenciales parciales de
segundo orden y en general no lineales. Se emplean también en los modelos matemáticos
ecuaciones constitutivas, en particular para un sistema multifásico totalmente saturado,
se usan las siguientes ecuaciones constitutivas
N
X
Sα = 1 , (4.10)
α=1

43
4. Modelos Matemáticos

pcα1 α2 = pα1 − pα2 ; α1 6= α2 . (4.11)

Para N número de fases, además α1 y α2 son las fases no-mojadora y mojadora
respectivamente. Las ecuaciones (4.9), (4.10) y (4.11) forman un sistema de ecuaciones
fuertemente acopladas.

4.4. Formulación Presión-Saturación para Flujos

Bifásicos
El modelo de dos fases para un flujo que se desplaza en un medio poroso hace uso
de la siguientes suposiciones :

El sistema consta de las fases agua (w) y aceite (o),

Las fases se consideran incompresibles,
La porosidad φ de la roca es constante.

Bajo las suposiciones anteriores y a partir de la ecuación (A.8), ver apéndice A,

llegamos a:
 
−∇ · λk fw ∇pw + fo ∇po − (fw ρw + fo ρo )g − (Qw + Qo ) = 0. (4.12)
Obtendremos una ecuación en términos de la presión del aceite, por lo que usaremos
pc = pcow = po − pw para eliminar pw en la ecuación (4.12):

−∇ · (λk∇po ) + ∇ · (λw k∇pc ) + ∇ · (k(λw ρw + λo ρo )g) − (Qw + Qo ) = 0.

En problemas de este tipo, es usual que la presión capilar dependa de la saturación

de agua, por lo que resulta conveniente expresar el gradiente de la presión capilar como
sigue:
dpc
∇pc (Sw ) = ∇Sw . (4.13)
dSw
Haciendo uso de la ecuación anterior, obtenemos la siguiente ecuación para la
presión del aceite:

 
dpc
−∇ · (kλ∇po ) + ∇ · kλw ∇Sw + ∇ · (k(λw ρw + λo ρo )g) = Qw + Qo . (4.14)
dSw
A fin de obtener una ecuación de saturación para la fase agua, partimos de la
ecuación (A.9) y tomando en cuenta las suposiciones al inicio de esta sección, obtenemos

44
 4.5. Método de Lı́neas de Corriente

∂Sw
φ + ∇ · uw = Qw . (4.15)
∂t
Sustituimos ahora en la ecuación (4.15) la Ley de Darcy para la fase acuosa, la
ecuación pc = pcow = po − pw con el objeto de eliminar pw y la relación (4.13), de modo
que ahora llegamos a:

 
∂Sw dpc
φ − ∇ · (kλw ∇po ) + ∇ · kλw ∇Sw + ∇ · (kλw ρw g) = Qw . (4.16)
∂t dSw

Las ecuaciones (4.14) y (4.16) se resolverán numéricamente en el capı́tulo 5.

4.5. Método de Lı́neas de Corriente

La idea en SLS es aproximar los cálculos tridimensionales del flujo y transporte
de las diferentes fases, mediante la solución de las ecuaciones de transporte en una
dimensión a lo largo de las lı́neas de corriente. La solución unidimensional hace que
este enfoque sea extremadamente rápido y efectivo para modelar flujos en yacimientos
en donde existen muchas heterogeneidades. La geometrı́a y la densidad de las lı́neas
de corriente reflejarán el impacto geológico sobre los caminos preferenciales del flujo,
introduciendo mayor densidad de lı́neas en regiones de alta porosidad y permeabilidad,
véase [16]. La técnica de SLS ha sido aplicada con éxito en la inyección de agua y
gases, ver [5], tales como CO2, en yacimientos naturalmente fracturados donde los
efectos de las mallas son un problema numérico importante, especialmente cuando
existen múltiples arreglos de pozos inyectores y productores.

4.5.1. Conceptos Fundamentales

Lı́neas de Corriente (Streamline).- Las Lı́neas de Corriente son curvas
tangentes en cada uno de sus puntos al campo de velocidad local, figura 4.2.
Sólo la dirección de la velocidad del fluido es importante, no su magnitud. En un
campo de velocidad variante en el tiempo las lı́neas de corriente se trazan para
un instante de tiempo en particular. Un concepto relacionado es la lı́nea de flujo
y se denomina ası́ a la trayectoria seguida por una partı́cula de un fluido móvil.
En general, a lo largo de la lı́nea de flujo, la velocidad de la partı́cula varı́a tanto
en magnitud como en dirección. Si todo elemento que pasa por un punto dado
sigue la misma trayectoria que los elementos precedentes, se dice entonces que se
trata de un flujo estacionario. Para una velocidad constante, las lı́neas de flujo y
las lı́neas de corriente describen la misma trayectoria. En sistemas con velocidad
variable, las lı́neas de corriente son una representación de un campo de velocidad
instantáneo, no una trayectoria fı́sica. Como consecuencia, las lı́neas de corriente
nunca se cruzan entre ellas a diferencia de las lı́neas de flujo.

45
4. Modelos Matemáticos

Figura 4.2: A la izquierda, el campo de velocidad en un dominio dado. A la derecha

observamos las lı́neas de corriente tangentes al campo de velocidad local. Imágenes
obtenidas con OpenDx.

Tubos de Corriente (Streamtube).- En dos dimensiones, un tubo de corriente

es la región comprendida entre dos lı́neas de corriente, figura 4.3. De la definición
de la lı́nea de corriente se deduce que no pasa fluido a través de las paredes
laterales de un tubo de corriente. Dentro de cada tubo, tenemos una descripción
1D del flujo, ası́ mismo, se cumple la ecuación de continuidad en cualquier sección
normal al tubo. Además, tubos con suficiente espacio corresponden a flujos lentos,
mientras que tubos estrechos presentan flujos rápidos.

Figura 4.3: Tubo de corriente. Imágenes obtenidas con OpenDx.

46
 4.5. Método de Lı́neas de Corriente

Tiempo de Vuelo (Time of Flight).- Introduzca una partı́cula en un pozo de

inyección y deje que la partı́cula se mueva de acuerdo a la velocidad intersticial
instantánea y mida el tiempo que le toma a la partı́cula en llegar a un punto
determinado, esa es la definición de tiempo de vuelo τ (x, y, z) para dicho punto,
(figura 4.4). El tiempo de vuelo se usa como coordenada espacial, en otras
palabras, la distancia desde la entrada de un sistema se mide por este tiempo, no
por la distancia Euclidiana. En simulación con lı́neas de corriente, usar el tiempo
de vuelo como coordenada espacial es fundamental.

Figura 4.4: Tiempo de vuelo τ para una partı́cula de prueba.

Funciones de Corriente (Streamfunction).- Es posible determinar la

velocidad de Darcy de un fluido a partir del gradiente de presión. Para flujo
de fluidos en 2D es posible determinar la velocidad a partir de la derivada de una
función llamada función de corriente. Una construcción gráfica de la función de
corriente se muestra en la figura 4.5, dados una velocidad de Darcy, dos puntos
y una trayectoria que conecta ambos puntos. Se trata de una función escalar,
como la presión, y como ésta, la función de corriente se determina en relación
a un punto de referencia, en este caso, el punto A, ψA = 0. Para determinar
la función de corriente en el punto B, dibujamos una trayectoria arbitraria del
punto A al B y calculamos el flujo total de Darcy que cruza dicha trayectoria.
La función de corriente en el punto B, es definida como el flujo volumétrico
normalizado por unidad de espesor. Del teorema de la divergencia, esta integral,
y por lo tanto la derivada de la función de corriente entre los puntos ψA y ψB es
independiente de la trayectoria, siempre que la trayectoria no cruce alrededor de
un pozo. Por definición, dado que no hay flujo que cruce una lı́nea de corriente,
la función de corriente es constante a lo largo de una lı́nea de corriente. Como
consecuencia, cuando la función de corriente es conocida, los contornos de dicha
función pueden ser usados para determinar las lı́neas de corriente. El concepto de
función de corriente puede ser extendido a 2 y 3 dimensiones usando funciones de
corriente dobles. Esta teorı́a puede ser extendida también a flujos compresibles
usando una densidad efectiva.

47
4. Modelos Matemáticos

Figura 4.5: Una función de corriente no depende de la definición de la trayectoria.

Simulación con Lı́neas de Corriente

En un tubo de corriente, la velocidad del fluido es igual al flujo volumétrico
por unidad de área. El flujo es constante a lo largo del tubo y el área de la
sección transversal del tubo se calcula explı́citamente. En una simulación con
lı́neas de corriente, se calcula la velocidad a través de diferencias finitas y se
trazan las lı́neas de corriente usando dicha velocidad. Al igual que en los tubos
de corriente, se asocia un flujo volumétrico a cada lı́nea de corriente. Es posible
obtener el área transversal efectiva a lo largo de la lı́nea de corriente dividiendo
el flujo volumétrico por la velocidad; por lo tanto, la geometrı́a espacial es
implı́cita. Resumiendo, los tubos de corriente calculan explı́citamente la sección
transversal de un tubo e implı́citamente calculan la velocidad, mientras que, en
contraste, las lı́neas de corriente hacen lo opuesto. Para grandes cambios en un
fluido en movimiento, la velocidad es actualizada periódicamente, nuevas lı́neas
de corriente son determinadas y saturaciones remuestreadas.

4.5.2. Lı́neas de Corriente

Como se mencionó en el apartado anterior, las lı́neas de corriente son curvas que
son localmente tangentes a la dirección de la velocidad. Las componentes del vector
velocidad v son vx , vy , y vz para el caso tridimensional. El vector local de longitud dr
tiene componentes dx, dy y dz. De a cuerdo a la figura 4.6, la pendiente de la lı́nea
de corriente en cualquier punto está dada por el cociente de las componentes de la
velocidad para un instante de tiempo t0 :

dy vy (x, y, z, t0 ) dz vz (x, y, z, t0 )
= , = . (4.17)
dx vx (x, y, z, t0 ) dx vx (x, y, z, t0 )
Estas ecuaciones diferenciales pueden ser integradas analı́tica o numéricamente
desde un punto (x0 , y0 , z0 ) y resolver para y(x) y z(x) a fin de determinar las lı́neas
de corriente que corren a través del punto de partida. Alternativamente, es posible
expresar estas mismas ecuaciones en su forma paramétrica:

48
 4.5. Método de Lı́neas de Corriente

dx dy dz
dt = = = , (4.18)
vx (x, y, z, t0 ) vy (x, y, z, t0 ) vz (x, y, z, t0 )
y determinar x(t), y(t) y z(t).
La ecuación de una trayectoria fı́sica (pathline) es muy similar, excepto que el campo
de velocidad puede ser independiente del tiempo.
dx dy dz
dt = = = . (4.19)
vx (x, y, z, t) vy (x, y, z, t) vz (x, y, z, t)
Esta diferencia nos recuerda que una lı́nea de corriente es definida como una linea en
el espacio obtenida trazando el campo de velocidad instantáneo. Es importante recalcar
algunos aspectos, es posible definir lı́neas de corriente para cualquier velocidad, y si
la velocidad en la ecuación (4.17) varı́a con el tiempo, entonces las lı́neas de corriente
cambiarán con el tiempo. Para las condiciones de estado no estacionario, se emplean
campos de velocidad instantánea para un tiempo de interés. A menudo, se aproximan
problemas de estado no-estacionario como series de campos de velocidad de estado
estacionario. Para definir lı́neas de corriente, la permeabilidad del medio puede ser
homogénea o heterogénea, isótropo o anisótropo y los fluidos pueden ser compresibles
o incompresibles.

Figura 4.6: Una función de corriente no depende de la definición de la trayectoria.

Existe una relación entre el potencial y las lı́neas de corriente que se obtiene a partir
de la Ley de Darcy de una fase:
1
u = − k · ∇Φ, (4.20)
µ
donde Φ denota el potencial del fluido, para nuestro caso en particular se trata de
la presión. Matemáticamente, dr × u = 0. Para un medio poroso isótropo, donde el
tensor de la permeabilidad es un escalar, tenemos que:

dr × ∇Φ = 0.

49
4. Modelos Matemáticos

Por lo tanto, los contornos del potencial son ortogonales a las lı́neas de corriente para
un medio isótropo.

4.5.3. Funciones de Corriente en 2D

Consideremos un flujo multifásico incompresible dentro de un medio permeable y
no-deformable en ausencia de la fuerza de gravedad y los efectos de capilaridad. La
velocidad total se define como la suma de las velocidades de Darcy para cada fase j,
esto es: X
u= uj = −λ∇P. (4.21)
j
P
La movilidad total está definida por λ = λj , y la fuente o sumidero total como
P j
qt = qj . La condición de compresibilidad requiere que:
j

∇ · u = qt , (4.22)

que se puede utilizar para obtener la ecuación diferencial para la presión:

∇ · (−λ∇P ) = qt . (4.23)
Para flujos incompresibles, no hay una dependencia temporal explı́cita en la
ecuación (4.22) ya que u es una velocidad instantánea. A medida que qt varı́a, o se
produzcan cambios en la saturación que provoquen que los coeficientes en la ecuación
(4.23) cambien, u también cambia, de otra manera, permanece fija.
Consideremos ahora la velocidad total lejos de fuentes y sumideros, es decir, qt = 0.
Para fluidos en dos dimensiones, la ecuación (4.22) la resolvemos si representamos la
velocidad en términos de la función de corriente de Lagrange ψ(x, y).
∂ψ ∂ψ
ux = , uy = . (4.24)
∂y ∂x
Si la velocidad u es conocida, puede ser expresada como una integral para ψ(x, y) :

Z Z   Z
∂ψ ∂ψ
ψ(x, y) − ψ0 = dψ = dx + dy = (ux dy − uy dx). (4.25)
∂x ∂y
0 0 0

La diferencia ψ(x, y) − ψ0 es el flujo total que fluye a través de una lı́nea arbitraria
trazada desde (x, y) a un punto inicial. Por ejemplo, en la figura 4.7 observamos una
trayectoria que va de 0 a A y de A a (x, y). La primera porción de la trayectoria es
a lo largo de una lı́nea de corriente donde uy = 0. De la ecuación (4.25), ningún flujo
cruza la lı́nea. Para la segunda porción de la trayectoria, ∂ψ/∂y = ux > 0, la función
de corriente incrementa. En otras palabras, ψ es una coordenada de flujo. La diferencia
entre dos valores de la función de corriente es igual a la tasa de flujo volumétrico entre

50
 4.5. Método de Lı́neas de Corriente

las dos lı́neas de corriente definida por estos valores. Esto define un tubo de corriente
en 2D. De la ecuación (4.24) podemos observar que las unidades para ψ son de flujo
volumétrico por unidad de espesor.

Figura 4.7: Lı́neas de corriente, función de corriente y tubo de corriente en 2D.

La función de corriente tiene un valor constante a lo largo de las lı́neas de

corriente. Por lo tanto, una manera de trazar lı́neas de corriente en dos dimensiones
será calculando la función de corriente y simplemente trazar sus contornos. La ecuación
(4.25) definirá la función de corriente ψ si la velocidad es conocida. Por el contrario,
si la velocidad no es conocida, en su lugar podemos resolver para ψ, y luego obtener
la velocidad de la ecuación (4.24). La ecuación de la función de corriente se obtiene
de la ecuación (4.21) mediante la integración del gradiente de presión alrededor de un
bucle cerrado. Tal integral de lı́nea debe anularse, independientemente de la función de
presión. A esta condición se le conoce como “irrotacionalidad” y es un requisito para
la existencia de un potencial de velocidad. Esto es,
Z
dP = 0,

o en su forma diferencial:

       
∂ ∂P ∂ ∂P ∂ 1 ∂ψ ∂ 1 ∂ψ
0=− + = + = ∇ · (λ−1 ∇ψ). (4.26)
∂x ∂y ∂y ∂x ∂x λ ∂x ∂y λ ∂y
Esta es la ecuación diferencial para la función de corriente. A diferencia de
la solución en lı́neas de fuentes y sumideros, no está limitada a medios porosos
homogéneos, pues la movilidad total (λ) puede depender de la posición.

4.5.4. Funciones de Corriente y Lı́neas de Corriente en 3D

Las funciones de corriente en 3D representan dos familias de superficies cuyas
intersecciones definen las lı́neas de corriente. Es posible representar cualquier campo
de velocidad en 3D en términos de tres funciones; ρ, ψ y χ.

51
4. Modelos Matemáticos

ρu = ∇ψ × ∇χ. (4.27)
Las funciones ψ y χ son conocidas como bifunciones de corriente. La “porosidad
efectiva”ρ es importante en la descripción de los fluidos compresibles. Por el momento
ρ = 1, hasta que se trate a fondo la compresibilidad más adelante.

u = ∇ψ × ∇χ. (4.28)
Las bifunciones de corriente tienen muchas de las ventajas de las funciones de
corriente en dos dimensiones. Del cálculo vectorial llegamos a que ∇ · (∇ψ × ∇χ) = 0.
Es importante hacer notar que cualquier velocidad que se puede expresar de acuerdo
a la ecuación (4.28), debe ser incompresible, esto es, que satisfaga:

∇ · u = 0. (4.29)

Esta formulación es lo suficientemente general para describir el flujo de fluidos tanto

en 2D como en 3D, aunque se requieren dos funciones de corriente para tal propósito.
Puesto que u se trata de una velocidad instantánea, tanto ψ como χ son funciones
de corriente instantáneas. Si u presenta una dependencia temporal, entonces, ψ y χ
variarán también.
A fin de obtener la función de corriente en dos dimensiones de la que hemos estado
hablando, hacemos ψ = ψ(x, y) y χ = z. Expandiendo la ecuación (4.28) encontramos
que la componente z de la velocidad es cero y que las componentes x y y están dados
por la ecuación (4.24).

El campo de velocidad obtenido de las ecuaciones (4.27) o (4.28) tiene una inter-
pretación geométrica. Tanto ψ como χ son funciones en el espacio. Consideremos ahora
las superficies definidas por ψ = constante y por χ = constante, la intersección de am-
bos planos define una lı́nea, que es la lı́nea de corriente trazada por la velocidad u.
Es posible obtener diferentes lı́neas de corriente seleccionando diferentes valores con-
stantes para ψ y χ. En la figura 4.8, se muestran dos superficies definidas por ψ y
otras dos por χ, las cuatro juntas definen un tubo en el espacio, con lı́neas de corriente
a lo largo de sus bordes.
Las trayectorias de las partı́culas nunca cruzan las superficies de corriente; el flujo total
que entra al tubo de corriente permanece en él.

La ecuación (4.28) se puede integrar sobre el área del tubo,

ZZ
∆Q = ∆ψ∆χ = u · n̂ da ≈ u · δa, (4.30)

Esto muestra que la rapidez |u| variará inversamente a la sección transversal del
área del tubo |δa|. La relación inversa entre el área y la velocidad es fundamental en
el modelado con tubos de corriente. Dicha ecuación también demuestra la relación

52
 4.5. Método de Lı́neas de Corriente

Figura 4.8: Dos superficies de corriente en (a) y en (b) y sus intersecciones en (c)
definen un tubo de corriente. Los bordes del tubo de corriente son lı́neas de corriente.

explı́cita entre el flujo total y las diferencias en ψ y χ. En la práctica, nunca

será necesario calcular ψ o χ, pero el flujo a lo largo de un tubo de corriente se
calcula a partir de la ecuación (4.30), de esta misma ecuación es posible observar las
unidades para ψ y χ, por ejemplo, si Q esta dada en bbl/D, entonces ψ tendrá unidades
de bbl/D/ft y χ unidades de ft (pies). Al igual que la función de corriente en dos
dimensiones,R las ecuaciones para ψ y χ pueden ser obtenidas integrando el gradiente
de presión dP = 0 alrededor de una curva cerrada. En su forma diferencial tenemos
que:  
1
0 = ∇ × (−∇P ) = ∇ × (∇ψ × ∇χ) .
λ
En contraste con la ecuación (4.26) en dos dimensiones de una sola incógnita ψ(x, y),
éstas son ahora tres ecuaciones acopladas (uno en cada dirección) para dos funciones
desconocidas, ψ y χ en tres dimensiones. En contraste a la utilidad de la función
de corriente en dos dimensiones, no existe manera obvia de resolver estas ecuaciones
acopladas.
Esta es una construcción formal de las bifunciones de corriente, que se deduce de
las ecuaciones diferenciales para la trayectoria, como se indica en la ecuación (4.17).
Podemos integrar estas dos ecuaciones independientes para derivar dos funciones de
(x, y, z) que no dependen del tiempo t. De nuevo, este es un enfoque complicado
sin beneficios evidentes en comparación con el cálculo numérico de la velocidad.
Generalmente, las bifunciones de corriente no son útiles como un medio de resolver
para la velocidad. Sin embargo, son muy útiles como parte de la formulación de tiempo
de vuelo, que abordaremos en la siguiente sección.

4.5.5. Lı́neas de Corriente y Tiempo de Vuelo

El tiempo de vuelo se refiere a una coordenada en especı́fico que se usa a lo largo
de las lı́neas de corriente. El uso del tiempo de vuelo como una coordenada espacial es
especialmente efectiva representando los efectos de la heterogeneidad del medio en el
flujo.

53
4. Modelos Matemáticos

A partir de un campo instantáneo de velocidad, es posible definir las funciones ψ(x, y, z)

y χ(x, y, z), tal y como se propuso en la sección anterior. El tiempo de vuelo lo
denotamos como τ (x, y, z), podemos colocar una serie de trazadores en cada pozo
inyector y determinar el tiempo que le toma a las partı́culas alcanzar una determinada
posición en el yacimiento. Para que esto sea el tiempo de tránsito real, el trazador debe
moverse a la velocidad intersticial, no la velocidad de Darcy. Por lo tanto, el tiempo
de vuelo se puede representar con la siguiente integral:
Z
φ
τ= dξ, (4.31)
|u|
0

La partı́cula de prueba se mueve con velocidad intersticial u/φ, y ξ es la distancia

espacial a lo largo de la lı́nea de corriente. Durante una simulación con lı́neas de
corriente, el cálculo del tiempo de vuelo no requiere un cálculo explı́cito de las
bifunciones de corriente, ψ(x, y, z) y χ(x, y, z). Una aproximación de la integral en
la ecuación (4.31) es basada en un esquema de diferencias finitas para el campo de
velocidad. La ecuación (4.31) puede ser reescrita como una relación diferencial:

u · ∇τ = φ, (4.32)

o de otra forma,
u ∆ξ
= . (4.33)
φ ∆τ
Aunque el tiempo de vuelo se mide en unidades de tiempo, por ejemplo dı́as, se
usará como una coordenada espacial.
Ahora veamos la transformación espacial del dominio (x, y, z) a (τ, ψ, χ). El
Jacobiano de la transformación relaciona los elementos de volumen en ambos espacios
coordenados,
∂(τ, ψ, χ)
= |(∇ψ × ∇χ) · ∇τ | = |u · ∇τ | = φ, (4.34)
∂(x, y, z)
en términos de volumen tenemos que:

φ dx dy dz = dτ dψ dχ. (4.35)
Una unidad de volumen en las coordenadas de tiempo de vuelo corresponde a
una unidad de volumen de poro en el espacio fı́sico. La ecuación (4.35) muestra un
importante vı́nculo entre la discretización espacial en diferencias finitas y la simulación
con lı́neas de corriente.
Los gradientes espaciales a lo largo de lı́neas de corriente también tienen una forma
muy simple en las coordenadas de tiempo de vuelo. Usando las coordenadas (τ, ψ, χ),
el operador gradiente se expresa ahora como:
∂ ∂ ∂
∇ = (∇τ ) + (∇ψ) + (∇χ) . (4.36)
∂τ ∂ψ ∂χ

54
 4.5. Método de Lı́neas de Corriente

Ya que u es ortogonal tanto a ψ como a χ, podemos escribir,

∂
u·∇=φ . (4.37)
∂τ
La ecuación (4.37) representa una relación importante en la simulación con lı́neas de
corriente, ya que servirá para transformar ecuaciones del espacio fı́sico a coordenadas
de tiempo de vuelo de lı́neas de corriente.

[Link]. Tiempo de Vuelo como Coordenada Espacial

La principal ventaja de la coordenada τ se hace evidente si tenemos en cuenta la
ecuación de conservación para la fase agua en un flujo incompresible en dos fases, lejos
de fuentes y sumideros,

∂Sw
φ + ∇ · (Fw ut ) = 0.
∂t
En la ecuación anterior se han despreciado por ahora los efectos de la aceleración
gravitacional y la capilaridad. Fw representa el flujo fraccional y lo denotamos como
Fw = λw /λ. Usando las ecuaciones (4.29) y (4.37) tenemos que:

∂Fw
∇ · (Fw ut ) = ut · ∇Fw = φ .
∂τ
Nótese el uso del operador identidad de la ecuación (4.29) para transformar del
espacio fı́sico a coordenadas de tiempo de vuelo, por lo tanto;

∂Sw ∂Fw
+ = 0. (4.38)
∂t ∂τ
Esta transformación de coordenadas descompone el flujo del fluido tridimensional
en series de ecuaciones 1D para la saturación a lo largo de las lı́neas de corriente. Esta
ecuación es válida para una, dos y tres dimensiones, ası́ como para medios homogéneos
y heterogéneos. Lo que se requiere para la implementación es el campo de velocidad y
el cálculo de la integral de lı́nea en la ecuación (4.31).

[Link]. Cálculo del Tiempo de Vuelo y Lı́neas de Corriente en 2D y 3D

La única caracterı́stica de los simuladores basados en el método de lı́neas de corri-
ente es el uso del tiempo de vuelo como una coordenada espacial. Para el cálculo del
tiempo de vuelo τ (x, y, z) se sigue una construcción de las lı́neas de corriente que se
conoce como el algoritmo de Pollock, en el cual, el tiempo de tránsito de un punto
inicial a un punto del espacio se calcula una celda a la vez.
En la aproximación de lı́neas de corriente, comenzamos con una solución numérica de
la ecuación de presión. Para flujos incompresibles, por ejemplo, primero se obtiene la

55
4. Modelos Matemáticos

distribución de presión en el yacimiento con un esquema adecuado de discretización

como el de volumen finito. Una vez conocida la distribución de la presión, es posible
calcular los flujos volumétricos usando la Ley de Darcy.

Para ilustrar el algoritmo de trazado de lı́neas de corriente, consideremos una malla

en bloques de volumen finito de la ecuación de presión, como se muestra en la figura
4.9. La solución numérica nos proporciona la presión al centro del volumen de control y
la velocidad del flujo en las caras del volumen. El algoritmo de Pollock usa una submalla
de modelo de velocidad que se deriva de la suposición de que cada componente de la
velocidad varı́a linealmente entre los valores para el par correspondiente de caras del
volumen. Esto quiere decir que la velocidad en el eje x varı́a linealmente sólo en la
dirección del eje x y es independiente de las velocidades en las otras direcciones.

Figura 4.9: Bloque de volumen finito para el cálculo del tiempo de vuelo.

Esto conduce al siguiente modelo de velocidad por celdas:

ux = ux1 + cx (x − x1 ),
uy = uy1 + cy (y − y1 ), (4.39)
uz = uz1 + cz (z − z1 ),
donde los coeficientes dependen de la diferencia de las velocidades de Darcy en las
caras del volumen,

cx = (ux2 − ux1 )/∆x,

cy = (uy2 − uy1 )/∆y, (4.40)
cz = (uz2 − uz1 )/∆z.
Las trayectorias de las lı́neas de corriente serán hipérbolas dentro de los volúmenes,
las ası́ntotas se alcanzan cuando ux = 0, uy = 0 ó uz = 0. De la ecuación (4.40)
obtenemos que:

56
 4.5. Método de Lı́neas de Corriente

3
X
∇·u= cj = cx + cy + cz . (4.41)
j=1

Este es un resultado importante ya que hemos comprobado que si la solución

discreta conserva el flujo, también lo hace la velocidad local dentro de la celda. Por
lo tanto, para flujos incompresibles, lejos de fuentes y sumideros, la solución numérica
proporcionará cx + cy + cz = 0. Sin embargo, para flujos compresibles la suma es
diferente de cero porque el fluido y la compresibilidad de la roca actúan como términos
fuente. El modelo de velocidad es útil para generar lı́neas de corriente para fluidos
compresibles e incompresibles. Las lı́neas de corriente y el tiempo de vuelo dentro de
cada celda pueden ser calculados a través de la integración directa de las velocidades
como se discutió en la ecuación (4.18).
dτ dx dy dz
= = = . (4.42)
φ ux uy uz
Para las velocidades lineales del juego de ecuaciones (4.40), dichas ecuaciones
pueden ser integradas de manera explı́cita e independiente de cada dirección.
Considérese una partı́cula en un punto arbitrario (x0 , y0 , z0 ) dentro de un volumen
de control como se muestra en la figura 4.9. La partı́cula puede salir del volumen a
través de cualquiera de las 6 caras. Integrando la ecuación (4.42) obtenemos el tiempo
de vuelo a cada una de las seis caras,

Zx1  
∆τxi dx 1 uxi
= = ln ,
φ ux0 + cx (x − x0 ) cx ux0
x0
Zy1  
∆τyi dy 1 uyi
= = ln , (4.43)
φ uy0 + cy (y − y0 ) cy uy0
y0
Zz1  
∆τzi dz 1 uzi
= = ln .
φ uz0 + cz (z − z0 ) cz uz0
z0

El ı́ndice i = 1, 2 indica las caras del volumen en cada dirección. el algoritmo de

Pollock especifica la cara de salida correcta la cual requiere de un tiempo de tránsito
mı́nimo y positivo. Por lo tanto, el tiempo de vuelo para una partı́cula estará dado por
el mı́nimo sobre los bordes permisibles, esto es:

∆τ = M inP ositivo (∆τx1 , ∆τx2 , ∆τy1 , ∆τy2 , ∆τz1 , ∆τz2 ), (4.44)

consideramos únicamente los valores positivos para seleccionar el mı́nimo. Cono-
ciendo el tiempo de vuelo de la partı́cula, sus coordenadas de salida pueden ser ahora
calculadas si reordenamos las ecuaciones (4.44):

57
4. Modelos Matemáticos

ecx ∆τ /φ − 1
 
x = x0 + ux0 = x0 + ux0 · ηx ,
cx
 cy ∆τ /φ 
e −1
y = y0 + uy0 = y0 + uy0 · ηy , (4.45)
cy
 cz ∆τ /φ 
e −1
z = z0 + uz0 = z0 + uz0 · ηz .
cz

El término en paréntesis incrementa con el tiempo y puede ser pensado como

pseudotiempos η. Las ecuaciones (4.46) definen la ecuación paramétrica de una lı́nea
de corriente hiperbólica, aunque en términos de η tienen una forma lineal.
Cuando la velocidad es constante a través del volumen en una dirección dada, por
ejemplo en la dirección x, entonces cx = 0. En el lı́mite,

 
1 ux x − x0
ln → , (4.46)
cx ux0 ux0
 cx ∆τ /φ 
e −1 ∆τ
ηx = → , (4.47)
cx φ
x = x0 + ux0 ∆τ /φ. (4.48)

Como se espera, la posición varı́a linealmente para una velocidad constante. Es

fácil identificar las lı́neas y los puntos de estancamiento basándonos en el algoritmo
de Pollock. Cada vez que cualquiera de las velocidades interpoladas cambia de signo
a través de una celda, significa que en algún lugar dentro de dicha celda la velocidad
es cero. Cuando esto ocurre, el tiempo de vuelo calculado a través de una celda en esa
dirección será infinito. El tiempo de vuelo en las otras direcciones aún puede ser finito
todavı́a, sin embargo, si todas las velocidades cambian de signo, entonces tendremos
un punto de estancamiento.

Una vez descrito el algoritmo para el tiempo de vuelo, se ilustrará con un ejemplo
simple en 2D, figura 4.10. Supongamos una partı́cula p situada en el punto (xp , yp ) en
un tiempo t que entra a una celda (i, j) en el plano (x, y). El primer paso es calcular
el tiempo de vuelo de acuerdo a la ecuación (4.44). De acuerdo con la figura 4.11,
el tiempo en en las caras y1 y x2 estará dado por la intersección de la trayectoria con
la extension de dichas caras de las celdas (en lı́neas punteadas). Para y1 el tiempo
será negativo, basándonos en la dirección del flujo. El tiempo para la cara x1 es cero
porque es ahı́ donde la partı́cula entra en la celda. Por lo tanto, el tiempo de vuelo
mı́nimo será en la cara y2 , que es la cara de salida de la partı́cula, una vez conocido el
tiempo de vuelo, podemos calcular la coordenada de salida usando (4.46).
El tiempo al cual la partı́cula abandona la celda está dado por t = tp + ∆τ . Esta
secuencia de cálculos se repite y los tiempos de tránsito serán acumulados a lo largo de

58
 4.5. Método de Lı́neas de Corriente

las múltiples celdas, llevando eventualmente a un productor, o rastrear hacia atrás a

un inyector, figura 4.12. La figura 4.13 muestra una lı́nea de corriente en un sistema
tridimensional.

Figura 4.10: Partı́cula en un esquema de celdas 2D para el cálculo del tiempo de vuelo.

Figura 4.11: Cálculo del tiempo de vuelo para una sola celda.

59
4. Modelos Matemáticos

Figura 4.12: Representación del algoritmo de Pollock en 2D.

Figura 4.13: Lı́nea de corriente en un sistema 3D.

Este es el algoritmo usado para calcular el tiempo de vuelo hacia un pozo productor
para un punto arbitrario dentro del dominio. Para calcular el tiempo de vuelo al pozo
inyector, se utiliza el mismo enfoque, excepto que ahora se hace el seguimiento de la
partı́cula hacia atrás, es decir, para τ negativa.

60
Capı́tulo 5

Modelo Numérico

5.1. Método de Volumen Finito (MVF)

En esta sección se obtiene la solución numérica de las ecuaciones (4.14) y (4.16)

mediante el Método de Volumen Finito (MVF). El MVF se deriva a partir de las forma
conservativa de las ecuaciones de balance. En este método, el dominio de estudio se
divide en un número de volúmenes de control que no se traslapan, de tal manera que
hay un volumen rodeando a cada punto de la malla, figura 5.1. Luego, se integra
la ecuación de balance (4.5) sobre cada volumen, lo cual es equivalente a aplicar la
ecuación de balance (4.4) sobre cada volumen de control. Los flujos a través de las
caras de los volúmenes se aproximan usando esquemas numéricos apropiados, y esto
da como resultado un conjunto de ecuaciones discretas, una para cada volumen de
control, las cuales deben resolverse para obtener una solución numérica aproximada.
Las ecuaciones discretas que resultan usando esta estrategia, expresan el principio
de conservación para la propiedad extensiva correspondiente en cada volumen de
control. De la misma forma, la ecuación diferencial expresa el mismo principio para un
volumen de control infinitesimal. Esta caracterı́stica es válida para cualquier número
de volúmenes sobre la malla y no solamente para un número grande de ellos. Por lo
tanto, aún una solución en una malla gruesa exhibirá un balance exacto o fı́sicamente
realista, aunque para obtener una buena precisión se requiere de un número grande de
volúmenes o esquemas de alto orden. Esta clara relación entre el algoritmo numérico y
el principio fı́sico de conservación es una de las mayores atracciones del MVF. Se dice
que el MVF es un método conservativo. En las secciones siguientes se adoptará por
simplicidad la siguiente notación; po = p y Sw ≡ S.

61
5. Modelo Numérico

Figura 5.1: Volumen de control en tres dimensiones. Las letras mayúsculas representan
los centros de los volúmenes de control, mientras que las minúsculas las caras de dichos
volúmenes.

5.2. Modelo Discreto de Flujo en dos Fases con

MVF
En esta sección se aplica el MVF a las ecuaciones (4.14) y (4.16) descritas en la
sección 3.5 del capı́tulo 3 correspondientes al modelo de flujo en dos fases.

5.2.1. Discretización de la Ecuación de Presión

A fin de discretizar la ecuación (4.14) definimos como función de flujo:

dpc
F = −kλ∇p + kλw ∇S. (5.1)
dS

Tomamos el tensor de la permeabilidad sólo con valores distintos en la diagonal

principal, como se muestra a continuación,
 
k11 0 0
k =  0 k22 0  ,
0 0 k33

con k11 6= k22 6= k33 , ası́ pues, las componentes de la función de flujo son:

62
 5.2. Modelo Discreto de Flujo en dos Fases con MVF

 
∂p dpc ∂S
Fx = −k11 λ − λw ,
∂x dS ∂x
 
∂p dpc ∂S
Fy = −k22 λ − λw ,
∂y dS ∂y
 
∂p dpc ∂S
Fz = −k33 λ − λw .
∂z dS ∂z

Por lo tanto, la ecuación (4.14) se escribe ahora como,

∂Fx ∂Fy ∂Fz

∇·F = + + = Qw + Qo . (5.2)
∂x ∂y ∂z

Aplicamos ahora la ecuación de balance (4.1) a la expresión anterior, de modo que,

Z   Z
∂Fx ∂Fy ∂Fz
+ + dV = (Qw + Qo )dV, (5.3)
∂x ∂y ∂z
∆V ∆V

para dV = dx dy dz, la integración se realiza sobre un volumen de control como

el mostrado en la figura 5.1, por lo que B(t) ≡ ∆V. El volumen de control de dicha
figura representa un hexaedro donde sus ejes son paralelos a los ejes coordenados de un
sistema cartesiano, por lo que favorece la aproximación de las integrales de la ecuación
(5.3).

Con base en el volumen de control mostrado en la figura 5.1, el primer término

del lado izquierdo de la ecuación integra como:

Zf Zn Ze
∂Fx
dx dy dz = ((Fx )e − (Fx )w )Ax , (5.4)
∂x
b s w

donde Ax = ∆y∆z es el área de aquellas caras del volumen de control paralelas

al plano yz. En la notación (Fx )e significa que Fx se debe evaluar en la cara e del
volumen de control. Aplicamos el mismo desarrollo al resto de los miembros de la
ecuación integral (5.3), obteniendo:

((Fx )e − (Fx )w ) Ax + ((Fy )n − (Fy )s )Ay + ((Fz )f − (Fz )b )Az = (Q̄w + Q̄o )∆V, (5.5)

donde Ay = ∆x∆z, Az = ∆x∆y y ∆V = ∆x∆y∆z. Q̄w y Q̄o son promedios de Qw

y Qo dentro del volumen de control, respectivamente. Ahora veamos la evaluación de
Fx en la cara e:

63
5. Modelo Numérico

 
∂p dpc ∂S
(Fx )e = −k11 λ − λw ,
∂x dS ∂x e
      
∂p dpc ∂S
= k11 (−λ)e + λw ,
∂x e dS e ∂x e
      
pE − pP dpc SE − SP
= k11 (−λ)e + λw ,
∆xe dS e ∆xe

donde las derivadas parciales son aproximadas a través del método de diferencias
finitas centrales y los subı́ndices E y P indican que la variable se evalúa en el centro
del volumen de control correspondiente. Cuando evaluamos Fx en la cara w se hace
de forma similar, por lo que el primer término de la ecuación (5.5) quedarı́a expresado
por,

((Fx )e − (Fx )w )Ax = (c)Ax + (v)Ax ,

donde

      
pE − pP dpc SE − SP
c = k11 (−λ)e + λw ,
∆xe dS e ∆xe

      
p P − pW dpc SP − SW
v = k11 (−λ)w + λw ,
∆xw dS w ∆xw

donde ∆xw = xP − xW , ∆xe = xE − xP , esto mismo aplica para las dos

direcciones, y y z. Las componentes Fy y Fz se evalúan de la misma manera en las
caras correspondientes. Sustituyendo las evaluaciones en la ecuación (5.5) y ordenando
términos obtenemos la ecuación discreta de la presión para el volumen de control P .

aP p P = aE p E + aW p W + aN p N + aS p S + aF p F + aB p B + q P , (5.6)

donde los coeficientes se definen como:

64
 5.2. Modelo Discreto de Flujo en dos Fases con MVF

k11 (λ)e Ax
aE = ≡ TE ,
∆xe
k11 (λ)w Ax
aW = ≡ TW ,
∆xw
k22 (λ)n Ay
aN = ≡ TN ,
∆yn
k22 (λ)s Ay
aS = ≡ TS , (5.7)
∆ys
k33 (λ)f Az
aF = ≡ TF ,
∆zf
k33 (λ)b Az
aB = ≡ TB ,
∆zb
donde:

a P = a E + a W + a N + a S + a F + a B ≡ TP ; (5.8)

          !
dpc dpc dpc dpc dpc dpc
qP = Te + Tw + Tn + Ts + Tf + Tb SP −
dSe dS w dS n dS s dS f dS b
       
dpc dpc dpc dpc
Te SE + Tw SW + Tn SN + Ts SS +
dS e dS w dS n dS s
    !
dpc dpc
Tf SF + Tb SB + (Q̄w + Q̄o )∆V ; (5.9)
dS f dS b

donde los coeficientes están definidos como:

k11 (λw )e Ax
Te = ,
∆xe
k11 (λw )w Ax
Tw = ,
∆xw
k22 (λw )n Ay
Tn = , (5.10)
∆yn
k22 (λw )s Ay
Ts = ,
∆ys
k33 (λw )f Az
Tf = ,
∆zf
k33 (λw )b Az
Tb = .
∆zb

65
5. Modelo Numérico

Figura 5.2: La matriz resultante presenta 7 bandas en el caso tridimensional.

Los coeficientes TN B y Tnb se conocen como transmisibilidades de cada volumen de

control, donde N B = P, E, W, N, S, F, B y nb = e, w, n, s, f, b, deben ser calculadas
adecuadamente o de lo contrario se introducen errores numéricos en la solución.
Las transmisibilidades son dependientes de la movilidad total λ ası́ como de la
movilidad del agua λw y del cambio de la presión capilar pc respecto a la saturación
S. Estas cantidades, al igual que la saturación, deben ser evaluadas en las caras
(λ)nb , (λw )nb y (dpc /dS)nb .
Las ecuaciones discretas para todos los volúmenes de control de la malla arroja un
sistema lineal de ecuaciones que consta de 7 bandas distintas de cero para el caso
tridimensional, figura 5.2. Se trata de un sistema diagonalmente dominante dada la
forma de los coeficientes, por lo que se ha empleado el método iterativo de gradiente
conjugado en la solución del sistema, ver apéndice B.
Las ecuaciones diferenciales (4.14) y (4.16) están acopladas, ésta última tiene
dependencia temporal, por lo tanto, la ecuación discreta (5.6) se resolverá para cada
paso de tiempo a fin de obtener pn+1 , usando información de la saturación del paso
anterior S n .

5.2.2. Discretización de la Ecuación de Saturación

Busquemos ahora discretizar la ecuación de saturación (4.16), al igual que en la
ecuación de presión, definimos una función de flujo:

dpc
F = −kλw ∇P + kλw ∇S. (5.11)
dS
A diferencia de la ecuación (5.1), en el primer término del lado derecho aparece la
movilidad del agua λw en lugar de la movilidad total λ. La obtención de las componentes
Fx , Fy y Fz es similar a como se hizo para la ecuación de presión. Tomando la definición
(5.11), la ecuación de la saturación (4.16) la escribimos como:

66
 5.2. Modelo Discreto de Flujo en dos Fases con MVF

∂S ∂S ∂Fx ∂Fy ∂Fz

φ +∇·F =φ + + + = QW . (5.12)
∂t ∂t ∂x ∂y ∂z
En términos de la ecuación de balance (4.1), la ecuación anterior la denotamos
como,
Z Z   Z
∂S ∂Fx ∂Fy ∂Fz
φ dV + + + dV = Qw dV, (5.13)
∂t ∂x ∂y ∂z
∆V ∆V ∆V

donde ∆V es el volumen de control de la figura 5.1, y B(t) ≡ ∆V .

Las derivadas con dependencia temporal se tratan, en el MVF, haciendo una integral
de toda la expresión en el intervalo [t, t + ∆t]. Ası́ pues, integramos la ecuación (5.13)
en dicho intervalo, obteniendo:

n+1 Z 

Z 
∂S ∂Fx ∂Fy ∂Fz
 φ + + + − Qw dV  dt = 0. (5.14)
∂t ∂x ∂y ∂z
n ∆V

El término de la derivada temporal lo aproximamos como se muestra a continuación:

n+1
Z Z  Z
∂S
S n+1 − S n dV.


φ dV dt = φ (5.15)
∂t
∆V n ∆V

Es necesario conocer la forma de S como función de x para poder aproximar la

integral del lado derecho, esto se puede llevar a cabo mediante funciones lineales o
de mayor orden, pero la saturación en el centro del volumen de control representa un
promedio de todo el volumen, SP , por lo tanto tenemos que:
n+1
Z Z 
∂S
φ dV dt = φ(SPn+1 − SPn )∆V . (5.16)
∂t
∆V n

El resto de los términos de la ecuación (5.14) se aproximan de forma parecida a

como se hizo para la ecuación de presión, obteniendo,

n+1 Z
Z n+1
Z 
(∇·F −Qw )dV dt = bP pP −(bE pE +bW pW +bN pN +bS pS +bF pF +bB pB +qP ) dt,
n ∆V n
(5.17)
donde los coeficientes b son similares a los coeficientes a de la ecuación (5.6), sólo
que en esta ocasión usamos la movilidad del agua λw en sustitución de la movilidad
total λ, incluso el término qP es similar al mostrado en la ecuación (5.9), exceptuando

67
5. Modelo Numérico

que en el último término sólo tenemos el promedio de la fuente de agua Q̄w .

En el algoritmo IMPES, la ecuación discreta de la saturación se resuelve de manera

explı́cita, aplicando el esquema θ a la ecuación (5.17) obtenemos lo siguiente:

n+1 Z
Z  
n n n n n n n n n n n n n n n
(∇·F −qw )dV dt = bP pP −(bE pE +bW pW +bN pN +bS pS +bF pF +bB pB +qP ) ∆t.
n ∆V
(5.18)
Sustituyendo (5.16) y (5.18) en (5.14) y reordenando términos obtenemos la relación
explı́cita para la saturación:
 
∆t
SPn+1 = SPn − n n n n n n n n n n n n n n n
bP pP −(bE pE +bW pW +bN pN +bS pS +bF pF +bB pB +qP ) , (5.19)
φ∆V
de modo que los coeficientes quedan expresados como:

k11 (λw )ne ∆t

bnE = ,
∆xe φ∆x
k11 (λw )nw ∆t
bnW = ,
∆xw φ∆x
k22 (λw )nn ∆t
bnN = ,
∆yn φ∆y
k22 (λw )ns ∆t
bnS = , (5.20)
∆ys φ∆y
k33 (λw )nf ∆t
bnF = ,
∆zf φ∆z
k33 (λw )nb ∆t
bnB = ,
∆zb φ∆z

bnP = bnE + bnW + bnN + bnS + bnF + bnB ; (5.21)

 n  n  n  n
∆t n dpc dpc dpc n dpc
q = bnE + bWn n
+ bN + bS +
φ∆V P dS e dS w dS n dS s
 n  n 
n dpc n dpc
bF + bB SPn −
dS f dS b
  n  n  n  n
n dpc n n dpc n n dpc n n dpc
bE SE + bW S W + bN SN + bS SSn +
dS e dS w dS n dS s
n n !
Q̄n
 
dpc dpc
bnF SFn + bnB SBn + w . (5.22)
dS f dS b φ

68
 5.3. Condiciones Iniciales y de Frontera

Las transmisibilidades Te son iguales a la definidas para la ecuación de presión.

Ecuaciones como la (5.19) se derivan para cada volumen de control de la malla, pero a
diferencia de la ecuación de presión, aquı́ no es necesario resolver ningún sistema lineal
ya que todos los términos del lado derecho son conocidos.

5.2.3. Algoritmo IMPES

Las ecuaciones (4.14) y (4.16) representan un problema no lineal y fuertemente
acoplado, de modo que es preciso establecer una estrategia de linealización. Bajo la
formulación presentada aquı́, se empleará el método IMPES (IMplicit Presure Explicit
Saturation) vease [14], el cual se describe en el algoritmo 1.

Algoritmo 1 IMPES
1: Definir condiciones iniciales y de frontera del problema S 0 , p0 , Tmax , ∆t, ∆x, ∆y, ∆z.
2: while t < Tmax do
3: calcular los coeficientes de la ecuación de presión usando (5.8), (5.9) y (5.11).
Dichos coeficientes se calculan usando valores de la saturación del paso anterior.
4: Resolver la ecuación de presión (5.6) de manera explı́cita usando un método
iterativo.
5: Calcular los coeficientes de la ecuación de saturación usando (5.21) y (5.22).
6: Resolver la ecuación de saturación (5.19) de manera explı́cita.
7: t ← t + ∆t
8: end while

5.3. Condiciones Iniciales y de Frontera

Las ecuaciones de presión y transporte requieren, para su solución, la especificación
de condiciones de contorno asociadas a los lı́mites del medio y de condiciones iniciales
que proporcionen información sobre los campos iniciales de presión y saturación.
Imponiendo las condiciones adecuadas nos aseguramos de tener un problema bien
planteado con solución única. En el caso de estudio, se han impuesto condiciones tipo
Neumann de no flujo en las fronteras del dominio y las siguientes condiciones iniciales.

Presión Saturación
Condiciones iniciales p(t0 ) = 1e + 07 Pa S(t0 ) = 0
Tabla 5.1: Condiciones iniciales para la presión y la saturación.

El modelo tiene en cuenta las siguientes consideraciones:

Se desprecian los efectos de la fuerza de gravedad.

69
5. Modelo Numérico

No hay fuentes ni sumideros.

Los fluidos son incompresibles e inmiscibles.

En este problema definimos la saturación efectiva como:

S − Srw
Sef = , (5.23)
1 − Srw − Sro
y las permeabilidades relativas como:

σ
krw = Sef , kro = (1 − Sef )σ . (5.24)
El exponente σ vale 1 en el caso lineal y 2 para el caso cuadrático. Las saturaciones
residuales del agua y del aceite son Srw y Sro respectivamente. Para evaluar los
coeficientes definidos en (5.8) determinamos el valor de λ en las caras nb=e,w,n,s,f,b
del volumen de control. Hacemos uso de las relaciones (5.23) y (5.24), la movilidad
total (λ) se define como la suma de las movilidades de cada fase:

λ = λw + λo ,
o bien
krw kro
λ= + , (5.25)
µw µo
sustituyendo el valor de las permeabilidades relativas en función de la saturación
efectiva tenemos que:
σ
Sef (1 − Sef )σ
λ= + .
µw µo
Haciendo las respectivas sustituciones y reduciendo términos obtenemos:

(S − Srw )σ (1 − Sro − S)σ

 
1
λ= + . (5.26)
(1 − Srw − Sro )σ µw µo
Si suponemos conocidos los valores de S para el instante n, podemos calcular λ
para el mismo instante de tiempo de la siguiente manera:
 n
(Snb − Srw )σ (1 − Sro − Snb
n σ

n 1 )
(λ)nb = + . (5.27)
(1 − Srw − Sro )σ µw µo
Para calcular los coeficientes de la ecuación (5.21) necesitamos el el valor de λw
en las caras del volumen de control, calculamos dicho valor con las relaciones (5.23) y
(5.24) como se muestra a continuación:
σ σ
Sef

krw 1 S − Srw
λw = = = . (5.28)
µo µw µw 1 − Srw − Sro
Ası́ pues, el valor de λw en las caras del volumen de control, para un instante de
tiempo n, lo calculamos como:

70
 5.4. Cálculo de la Saturación en las Caras

 n
σ
1 Snb − Srw
(λw )nnb = . (5.29)
µw 1 − Srw − Sro

5.4. Cálculo de la Saturación en las Caras

En las ecuaciones (5.27) y (5.29) notamos la necesidad de calcular los valores de la
saturación en las caras de los volúmenes de control. El valor de la saturación es conocido
en todos los centros de los volúmenes, mismos que provienen de la solución de la
ecuación (4.16), de modo que requerimos de un esquema que interpole adecuadamente
dicho valor a las caras. Se presenta el esquema Upwind para tales fines.

Upwind.- Se comparan los valores de presión en los puntos vecinos a la cara del
volumen de control donde se desea evaluar la saturación S n tomándose el valor
de la saturación del punto donde la presión es mayor.

if (pnP ≤ pnE ) then

Sen = SEn
else
Sen = SPn
end if

Esta es una aproximación lineal donde se asumen conocidos los valores de presión
para el instante n. Produce dispersión numérica, aunque se trata de un esquema estable,
por lo que es recomendable usar un número considerablemente grande de volúmenes
de control.

5.5. Cálculo del Transporte a lo largo de las Lı́neas

de Corriente
Hemos visto que el principio básico de simulación mediante lı́neas de corriente es
descomponer las ecuaciones de transporte multidimensionales en una serie de ecua-
ciones 1D a lo largo de las lı́neas de corriente.

En secciones previas se revisó la formulación numérica del cálculo de las lı́neas de

corriente y el tiempo de vuelo. Tratemos ahora en este apartado, la solución 1D de la
ecuación de saturación a lo largo de cada lı́nea de corriente usando el tiempo de vuelo
como coordenada espacial dada en la ecuación (4.38).

71
5. Modelo Numérico

5.5.1. Discretización de la Ecuación de Saturación 1D

Tomaremos la discretización 1D en volumen finito de la ecuación (4.38):
∂Sw ∂Fw
+ = 0,
∂t ∂τ
donde Fw se define cómo:
λw
. Fw = (5.30)
λ
En términos de la ecuación de balance tenemos que:
Z Z
∂Sw ∂Fw
dτ + dτ = 0.
∂t ∂τ
∆τ
| {z } ∆τ | {z }
1 2
=⇒ Discretización de 1:
integrando todo en el intervalo [t, t + ∆t];

n+1 Z

Z
 ∂Sw dτ  dt = 0,
∂t
n ∆τ

n+1
Z Z Z
∂Sw
dτ dt = (S n+1 − S n )dτ,
∂t
∆τ n ∆τ
= (Spn+1 − Spn )∆τ.
=⇒ Discretización de 2:

n+1Ze
Z n+1
Z
∂Fw
dtdτ = [(Fw )e − (Fw )w ] dt,
∂τ
n w n
= [(Fw )ne − (Fw )nw ]∆t.
Reacomodando términos:
∆t
Spn+1 = Spn − [(Fw )ne − (Fw )nw ] . (5.31)
∆τ
Finalmente, empleando la relación (5.30) obtenemos la ecuación discreta de la
saturación sobre las lı́neas de corriente:
 n  n 
n+1 n λw λw ∆t
Sp = Sp − − . (5.32)
λ e λ w ∆τ
En el capı́tulo siguiente abordaremos la codificación para la solución computacional
de la ecuación anterior y, en general, de las ecuaciones discretas de presión y saturación.

72
Capı́tulo 6

Modelo Computacional y
Resultados Numéricos

Los modelos matemáticos y numéricos abordados en capı́tulos anteriores necesitan

ser resueltos a través de la codificación de un conjunto de algoritmos especializados a
fin de obtener una solución numérica del problema que sea precisa, estable y eficiente.
En el presente trabajo se utilizará el sofware TUNAM (Templates Units for Numerical
Applications and Modeling), ver [18], el cual contiene algoritmos para el MVF y que
fue adaptado para dar solución al Método de Lı́neas de Corriente.

6.1. Software TUNAM

TUNAM es un software desarrollado en un principio para resolver problemas de
convección natural en dominios rectangulares mediante el método de volumen finito.
En su construcción se hizo uso de los paradigmas de programación orientada a objetos
(POO) y programación genérica, de modo que es posible utilizar las herramientas de
TUNAM para resolver problemas de otras áreas de estudio. Sus componentes estan
hechas en el lenguaje de programación C++ haciendo uso intensivo de templates
que proveen una herramienta efectiva para desarrollar programas genéricos. Hace uso
también de la biblioteca Blitz++, para el manejo de arreglos multidimensionales, véase
[11]. En TUNAM se hace uso de las siguientes definiciones:

Definición 1 Una generalización (Generalization) propone la existencia de un

conjunto de elementos con caracterı́sticas comunes, por ejemplo sus atributos y
operaciones.

73
6. Modelo Computacional y Resultados Numéricos

Definición 2 Una especialización (Specialization) se trata de un caso particular de

una entidad general que adiciona y/o redefine caracterı́sticas tales como atributos y
operaciones especiales.

Definición 3 Un adaptador (Adaptor) es una implementación particular de un

mismo concepto o algoritmo.

Tales definiciones pueden ser empleadas en el uso de templates, pues permiten

implementaciones alternativas de un mismo concepto a fin de generar código
optimizado. En la figura 6.1 se muestra un esquema general de TUNAM. Teniendo
en cuenta, los modelos conceptuales, fı́sicos, matemáticos y numéricos mostrados en
los capı́tulos anteriores, es posible identificar:

Generalizaciones :
GeneralMesh<> : Representa cualquier tipo de malla, por ejemplo la figura 5.1.
GeneralEquation<> : Representa la ecuación general discretizada, por ejemplo
ecuación (5.6), junto con sus respectivos coeficientes.
GeneralMatrix<> : Representa una matriz en general.

Especializaciones :
StructuredMesh<> : Representa mallas estructuradas como que se muestra en
la figura 5.1.
TwoPhaseEquation<> : Son ecuaciones diferenciales particulares para flujo
bifásico.
SparseMatrix<> : Representa matrices ralas o dispersas, [41, 46].

Adaptadores :
Uniform, NonUniform<> : Definen mallas uniformes y no uniformes en dominios
rectangulares.
CDS, Upwind, Quick<> : Definen los esquemas numéricos apropiados que se
emplearán en el cálculos de los coeficientes de volumen finito.
Diagonal<> : Define matrices ralas en formato diagonal.

74
 6.2. Implementación

(a)

(b)
Figura 6.1: (a) Arquitectura general de TUNAM. (b) Los paquetes FVM y Geom
se muestran de manera esquemática: las letras G, S y A, a la derecha de la figura,
significan Generalization, Specialization y Adaptor, respectivamente.

Se ha desarrollado un módulo nuevo para la simulación mediante lı́neas de corriente.

Los objetos en TUNAM que interactúan con el fin de resolver un problema se
definen y crean de la siguiente menara:
S p e c i a l i z a t i o n <Adaptor<p r e c t , dim>> o b j e c t ( arg1 , . . . , argN ) ;

El compilador de C++ analiza la declaración y genera el código de acuerdo

con la implementación de la especialización y el adaptador correspondientes. Donde
object es un objeto perteneciente a la clase Specialization<Adaptor<prec, dim>>
y además arg1, ..., argN son empleados en la construcción de object. El parámetro
prec t define la precisión de los resultados, es decir float, double o long double, el
parámetro dim define la dimensión del problema.

6.2. Implementación
6.2.1. Calibración 1D con el modelo de Buckley- Leverett
Las primeras pruebas para calibrar el método de lineas de corriente se harán sobre
un modelo 1D de Buckley- Leverett, el cual describe el desplazamiento de aceite por
agua en un dominio horizontal como el que se muestra en la figura 6.2. Se han
contemplado las siguientes suposiciones:

75
6. Modelo Computacional y Resultados Numéricos

Se desprecian los efectos de la presión capilar.

Se considera despreciable la acción de la fuerza de gravedad.

No hay fuentes ni sumideros.

Los fluidos son incompresibles e inmiscibles.

El medio poroso es homogéneo.

Figura 6.2: Medio homogéneo de 300 m de longitud, en un principio saturado de aceite.

En este caso, se contempla la inyección de agua a una razón de flujo constante en el

extremo izquierdo del dominio, el agua desplazará al aceite al extremo derecho donde
la presión se mantiene constante. Con las suposiciones anteriores, las ecuaciones (4.14)
y (4.16)se transforman en

−∇ · (kλ∇p) = 0, (6.1)

∂S
φ − ∇ · (kλw ∇p) = 0, (6.2)
∂t
para po = p y Sw = S.
A continuación se muestran los modelos discretos 1D de Buckley-Leverett, (6.3) y
lı́neas de corriente, (6.4) para el cálculo de la saturación, ambas discretizaciones son en
volumen finito.

k[λw ]ne n k[λw ]nw n k[λw ]np k[λw ]nw

   
∆t
Swn+1 = Swn− PE + PW − + PPn , (6.3)
∆x ∆x ∆x ∆x φ∆x
 n  n 
λw λw |ū|∆t
S n+1 = S n
− . (6.4)
λ e λ w φ∆ξ

Con las siguientes condiciones de frontera:

76
 6.2. Implementación

Presión Saturación
Condiciones iniciales p(t0 ) = 1e + 07 Pa S(t0 ) = 0
(S in )A = 0.8
(gpin )A = 3.4722e − 07 m/s
Condiciones de frontera
(S out )B = 0
(pout )B = 1e + 07 Pa
Tabla 6.1: Condiciones iniciales y de frontera para la presión y la saturación.

Las condiciones de frontera para la saturación son de tipo Dirichlet, mientras que
para la presión, tenemos condiciones tipo Dirichlet a la salida y una condición tipo
Neumann a la entrada dada en términos de la velocidad de inyección.
En este ejemplo se ha puesto una lı́nea de corriente que parte del extremo izquierdo,
donde inicia la inyección, y termina en el extremo derecho, coincidente en sus puntos
con los centros de los volúmenes de control, lo que nos permite hacer una comparación
puntual entre ambos métodos. La figura 6.3 muestra una comparación de la solución
para la ecuación de saturación 1D con MVF en verde y SLS en rojo. Nótese que
la solución es muy similar, por lo que ambos métodos resuelven adecuadamente la
ecuación de transporte.

Figura 6.3: Saturación obtenida mediante lı́neas de corriente en rojo y con volumen
finito en verde.

Una manera de cuantificar el error de la aproximación de ambas soluciones es

obteniendo el error cuadrático medio, el cual lo denotamos como:
s
N  2
P (mvf ) (sls)
Si − Si
i=1
, (6.5)
N

77
6. Modelo Computacional y Resultados Numéricos

el valor obtenido para las curvas de la figura anterior es de 387.59e−6, lo que refleja
que la solución a la ecuación de transporte mediante lı́neas de corriente es una buena
aproximación en una simulación de este tipo.

6.2.2. Estudio de Caso

En esta sección se resuelven numéricamente las ecuaciones (5.6) y (5.19) que
corresponden a la discretización en volúmenes finitos de las ecuaciones (4.14) y (4.16)
respectivamente. En esta sección se analiza un problema conocido como Five Spots
Pattern, la geometrı́a del dominio de muestra en la figura 6.4. Se tienen cuatro pozos
productores en las esquinas por uno inyector en el centro del dominio. En las fronteras
del dominio de estudio se impone una condición de no flujo. Los datos del problema se
muestran en la tabla 6.2. Se contempla un dominio en tres dimensiones inicialmente
saturado de aceite y un flujo bisáfico.

Figura 6.4: Geometrı́a del dominio del caso de cinco pozos.

Parámetros Valores Unidades en SI

dimensión x 182.76 [m]
dimensión y 182.76 [m]
dimensión z 9.14 [m]
Razón de inyección de agua (Qw ) 3.86e-04 [m3 /s]
Permeabilidad absoluta (k) 0.9869e-15 [m2 ]
Porosidad (φ) 0.2 -
Viscosidad del agua (µw ) 1.0e-03 [P a · s]
Viscosidad del aceite (µo ) 1.0e-03 [P a · s]
Saturación residual del agua (Srw ) 0 -
Saturación residual del aceite (Sro ) 0.2 -
Tabla 6.2: Datos para el caso de estudio.

78
 6.2. Implementación

6.2.3. Caso homogéneo (MVF)

Para este caso se requiere de una especialización particular de la ecuación
general que tome en cuenta el medio poroso y las dos faces del fluido. La clase
TwoPhaseEquation hereda de GeneralEquation todos los atributos y operaciones para
describir un flujo en dos fases.
Para calcular los coeficientes de la discretización se requieren de varias definiciones
particulares, tanto para la ecuación de presión como para la de saturación. El cálculo
de la saturación en las caras de los volúmenes de control se aproximó mediante el
esquema Upwind. El código a continuación es un fragmento de la implementación del
problema en dos fases. Se definen los arreglos donde se almacenarán la solución e
información durante la simulación. Dichos arreglos son de Blitz++, por lo que TUNA
define una interfaz a estos arreglos.
1 typedef TunaArray<double , 3 > : : huge S c a l a r F i e l d 3 D ;

lo anterior define el tipo ScalarField3D como una clase para construir arreglos
tridimensionales de doble precisión. Enseguida, se define la malla para el problema, se
hace como se muestra a continuación:
2 StructuredMesh<Uniform<double , 3> > mesh ( l e n g t h x , num nodes x ,
3 l e n g t h y , num nodes y ,
4 l e n g t h z , num nodes z ) ;
5 double dx = mesh . g e t D e l t a (X) ;
6 double dy = mesh . g e t D e l t a (Y) ;
7 double dz = mesh . g e t D e l t a ( Z ) ;
8 mesh . p r i n t ( ) ;

El código anterior define una malla estructurada 3D, uniforme y de precisión doble.
El objeto mesh es construido a partir de los argumentos de entrada length x, length y,
length z los cuales representan las longitudes del dominio en las direcciones x, y y z,
los parámetros num nodes x, num nodes y y num nodes z el número de nodos en cada
dirección. Los datos son definidos por el usuario. El objeto mesh nos puede proporcionar
información de la malla. En las lı́neas 5, 6 y 7 obtenemos el tamaño de la malla en
las direcciones x, y y z. Dicho objeto también realiza acciones como la impresión a
salida estándar en la lı́nea 8. Es posible obtener información de la malla en términos
de los nodos y de los volúmenes de control, resulta de utilidad para crear arreglos que
almacenarán la solución del problema y otros que sean definidos sobre la malla. La
siguiente fracción del código implementa estas caracterı́sticas y define las condiciones
iniciales:
9 ScalarField3D p ( mesh . getExtentVolumes ( ) );
10 S c a l a r F i e l d 3 D Sw ( mesh . getExtentVolumes ( ) );
11 S c a l a r F i e l d 3 D p n ( mesh . g e t E x t e n t N o d e s ( ) );
12 S c a l a r F i e l d 3 D Sw n ( mesh . g e t E x t e n t N o d e s ( ) );
13 S c a l a r F i e l d 3 D u1 ( mesh . getExtentVolumes ( ) );
14 S c a l a r F i e l d 3 D u2 ( mesh . getExtentVolumes ( ) );
15 S c a l a r F i e l d 3 D u3 ( mesh . getExtentVolumes ( ) );
16 Sw = 0 ; // I n i t i a l c o n d i t i o n
17 p = p r e s o u t ; // I n i t i a l c o n d i t i o n
18 Range a l l = Range : : a l l ( ) ;
19 NumUtils : : i n t e r p o l a t e T o N o d e s ( p n , p ) ;

79
 6. Modelo Computacional y Resultados Numéricos

20 InOut : : w r i t e T o F i l e D X ( p n , 0 , ” . / DataFS3D/ p r e s . ” , dx , dy , dz ) ;
21 NumUtils : : i n t e r p o l a t e T o N o d e s ( Sw n , Sw) ;
22 InOut : : w r i t e T o F i l e D X ( Sw n , 0 , ” . / DataFS3D/ s a t u . ” , dx , dy , dz ) ;

En las lı́neas 9 y 11 se generan arreglos p y p n para la presión en términos de

la extensión de los volúmenes y los nodos respectivamente, para las lı́neas 10 y 12
los arreglos Sw y Sw n almacenarán la información concerniente a la saturación en
términos de los volúmenes para el primero y en términos de los nodos para el segundo.
Luego, las lı́neas 13 a 15 definen el campo de velocidad del dominio almacenando la
información para cada dirección x, y y z. Por otra parte, las lı́neas 16 y 17 observamos
las condiciones iniciales para la saturación y de la presión. La lı́nea 18 utiliza una
herramienta de BLITZ++ para definir un rango de un arreglo. Las lı́neas 20 y 22
representan la salida a un archivo con formato OpenDX, ver [35], para su posterior
visualización, ambas corresponden de la interpolación a los nodos de la malla para la
presión y la saturación que se hace en las lı́neas 19 y 21 respectivamente. A fin de
definir el sistema lineal que contendrá los coeficientes de la discretización producto del
MVF del sistema de ecuaciones diferenciales hacemos lo siguiente:
23 S p a r s e M a t r i x < Diagonal <double , 3> > A( num nodes x , num nodes y , num nodes z ) ;
24 ScalarField3D b ( num nodes x , num nodes y , num nodes z ) ;

La dicretización en MVF produce una matriz dispersa, en el caso tridimensional

contiene entradas distintas de cero que caen en 7 diagonales como lo muestra la figura
5.2. La clase SparceMatrix es una especialización de GeneralMatrix, mientras que
Diagonal es un adaptador que hace de esta especialización una implementación óptima
para las matrices requeridas por el MVF. En la lı́nea 23 se define un arreglo que
almacenará la matriz diagonal de 7 bandas, la lı́nea 24 define el vector del lado derecho
del sistema lineal, el cual está mapeado a cada punto de la malla. Ya definidos los
objetos p, A, b y mesh construimos las ecuaciones a resolver, comencemos por describir
la ecuación de presión:
25 TwoPhaseEquation< FSIP1<double , 3> > p r e s s u r e ( p , A, b , mesh . g e t D e l t a s ( ) ) ;

donde TwoPhaseEquation es una especialización de GeneralEquation y FSIP1 es

un adaptador para calcular los coeficientes de MVF de la ecuación (5.6). El objeto
pressure representa la ecuación discreta (5.6) y es posible enviarle mensajes para que
ejecute acciones para definir caracterı́sticas del problema, tales como las condiciones
de frontera, como se muestra a continuación:
26 pressure . s e t D e l t a T i m e ( dt ) ;
27 pressure . setPermeability ( permeability ) ;
28 pressure . setPorosity ( porosity ) ;
29 pressure . s e t S r w ( Srw ) ;
30 pressure . s e t S r o ( Sro ) ;
31 pressure . s e t V i s c o s i t y w (mu w) ;
32 pressure . s e t V i s c o s i t y o ( mu o ) ;
33 pressure . setInjection ( injection ) ;
34 pressure . setNeumann (LEFT WALL) ;
35 pressure . setNeumann (RIGHT WALL) ;
36 pressure . setNeumann (TOP WALL) ;
37 pressure . setNeumann (BOTTOM WALL) ;
38 pressure . setNeumann (FRONT WALL) ;

80
 6.2. Implementación

39 p r e s s u r e . setNeumann (BACK WALL) ;

40 p r e s s u r e . s e t S a t u r a t i o n (Sw) ;
41 pressure . print () ;

en la lı́nea 26 se define el paso de tiempo dt que se empleará en la solución de la

evolución temporal del problema. Las siguientes lı́neas de código (27-33) establecen
los valores de permeabilidad, porosidad, saturación residual del agua, del aceite, la
viscosidad del agua, la del aceite y la tasa de inyección de agua. Es importante
mencionar que los valores de permeabilidad y porosidad para este caso son homogéneas
en el dominio de estudio. Las lı́neas 34-40 establecen las condiciones de frontera,
tipo Neumann, para la presión. De manera similar, construimos ahora la ecuación
de saturación con sus respectivos parámetros y condiciones de frontera:
42 TwoPhaseEquation< FSES1<double , 3> > s a t u r a t i o n (Sw , A, b , mesh . g e t D e l t a s ( ) ) ;
43 s a t u r a t i o n . s e t D e l t a T i m e ( dt ) ;
44 saturation . setPermeability ( permeability ) ;
45 saturation . setPorosity ( porosity ) ;
46 s a t u r a t i o n . s e t S r w ( Srw ) ;
47 s a t u r a t i o n . s e t S r o ( Sro ) ;
48 saturation . setInjection ( injection ) ;
49 s a t u r a t i o n . s e t V i s c o s i t y w (mu w) ;
50 s a t u r a t i o n . s e t V i s c o s i t y o ( mu o ) ;
51 s a t u r a t i o n . setNeumann (LEFT WALL) ;
52 s a t u r a t i o n . setNeumann (RIGHT WALL) ;
53 s a t u r a t i o n . setNeumann (TOP WALL) ;
54 s a t u r a t i o n . setNeumann (BOTTOM WALL) ;
55 s a t u r a t i o n . setNeumann (FRONT WALL) ;
56 s a t u r a t i o n . setNeumann (BACK WALL) ;
57 saturation . setPressure (p) ;
58 saturation . print () ;

Ahora, con el problema bien planteado resolvemos con:

59 while ( t <= Tmax) {
60 pressure . calcCoefficients () ;
61 S o l v e r : : TDMA3D( p r e s s u r e , t o l e r a n c e , t d m a i t e r , 1 . 0 ) ;
62 p r e s s u r e . update ( ) ;
63 saturation . calcCoefficients () ;
64 Solver : : solExplicit3D ( saturation ) ;
65 s a t u r a t i o n . update ( ) ;
66 t += dt ;
67 }

Las lı́neas 59-67 implementan el algoritmo IMPES. En el ciclo del código anterior,
se hacen iteraciones desde t=0 hasta Tmax con incrementos de paso de tiempo
dt. La función calcCoefficients llena la matriz A y el vector b, mientras que
Solver::TDMA3D es una función que implementa el algoritmo descrito en el apéndice
B para resolver sistemas lineales. En la lı́nea 62 se actualiza la presión. De manera
análoga, se calculan los coeficientes de la ecuación discreta de la saturación en la
lı́nea 63 y se actualiza el valor de la saturación en la lı́nea 65. Nótese que la presión
se resuelve de manera implı́cita (lı́nea 61) mientras que la saturación es resuelta de
manera explı́cita mediante la función Solver::solExplicit3D (lı́nea 64). La figura
6.5 muestra un cuarto del dominio de estudio y la evolución del frente de saturación
del agua. La simulación se llevó a cabo para un tiempo de 1000 dı́as con incrementos
de 1 hora.

81
6. Modelo Computacional y Resultados Numéricos

Figura 6.5: Evolución del frente de saturación para un medio homogéneo.

82
 6.2. Implementación

En la figura anterior, observamos el frente de saturación para distintos dı́as de

la simulación. Éste evoluciona de manera geométrica, pues el medio es homogéneo e
isótropo, es decir, tiene las mismas caracterı́sticas de porosidad y permeabilidad en
todas las direcciones.

6.2.4. Caso no-homogéneo (MVF)

En esta sección abordaremos la implementación del caso heterogéneo, es decir, con
permeabilidad variable en el dominio de estudio, figura 6.6. Se han introducido val-
oreas aleatorios en las 3 direcciones; x, y y z, ası́ como un cuerpo de baja permeabilidad
al centro del dominio.
La forma de solucionar este problema es igual que el anterior, sólo hay que modificar
la manera en la que se calculan los coeficientes de la presión y la saturación para el
caso 3D y teniendo en cuenta que leeremos un valor aleatorio de permeabilidad en cada
nodo de la malla.

Figura 6.6: Valores de permeabilidad en el dominio de estudio, de izquierda a derecha:

permeabilidad en dirección x, y y z.

Primero hay que modificar la especialización TwoPhaseEquation agregando los

nuevos campos escalares para las tres direcciones distintas como se muestra a
continuación:
i n l i n e void s e t P e r m e a b i l i t y V 1 1 ( S c a l a r F i e l d & p ) { perme 11 . r e f e r e n c e ( p ) ; }
i n l i n e void s e t P e r m e a b i l i t y V 2 2 ( S c a l a r F i e l d & p ) { perme 22 . r e f e r e n c e ( p ) ; }
i n l i n e void s e t P e r m e a b i l i t y V 3 3 ( S c a l a r F i e l d & p ) { perme 33 . r e f e r e n c e ( p ) ; }

y agregando dichos campos en los atributos públicos

public :
S c a l a r F i e l d perme 11 ;
S c a l a r F i e l d perme 22 ;
S c a l a r F i e l d perme 33 ;

Modificamos ahora los adaptadores donde se calculan los coeficientes de la presión

y la saturación. Para la presión modificamos a FSIP1 en la sección donde se calculan
los coeficientes 3D como lo muestra el siguiente extracto del código:

83
 6. Modelo Computacional y Resultados Numéricos

mult o 11 = perme 11 ( i , j , ki ) / ( (1 − Srw − Sro ) ∗ mu o );

mult o 22 = perme 22 ( i , j , ki ) / ( (1 − Srw − Sro ) ∗ mu o );
mult o 33 = perme 33 ( i , j , ki ) / ( (1 − Srw − Sro ) ∗ mu o );
mult w 11 = perme 11 ( i , j , ki ) / ( (1 − Srw − Sro ) ∗ mu w );
mult w 22 = perme 22 ( i , j , ki ) / ( (1 − Srw − Sro ) ∗ mu w );
mult w 33 = perme 33 ( i , j , ki ) / ( (1 − Srw − Sro ) ∗ mu w );

aE (i , j , ki ) = ( (1 − Sro − Sw e ) ∗ mult o 11 + ( Sw e − Srw ) ∗ mult w 11 ) ∗ dydz dx ;

aW (i , j , ki ) = ( (1 − Sro − Sw w ) ∗ mult o 11 + ( Sw w − Srw ) ∗ mult w 11 ) ∗ dydz dx ;
aN (i , j , ki ) = ( (1 − Sro − Sw n ) ∗ mult o 22 + ( Sw n − Srw ) ∗ mult w 22 ) ∗ dxdz dy ;
aS (i , j , ki ) = ( (1 − Sro − Sw s ) ∗ mult o 22 + ( Sw s − Srw ) ∗ mult w 22 ) ∗ dxdz dy ;
aF (i , j , ki ) = ( (1 − Sro − Sw f ) ∗ mult o 33 + ( Sw f − Srw ) ∗ mult w 33 ) ∗ dxdy dz ;
aB (i , j , ki ) = ( (1 − Sro − Sw b ) ∗ mult o 33 + ( Sw b − Srw ) ∗ mult w 33 ) ∗ dxdy dz ;
aP (i , j , ki ) = aE ( i , j , k i ) + aW (i , j, ki ) +
aN ( i , j , k i ) + aS (i , j, ki ) +
aF ( i , j , k i ) + aB (i , j, ki ) ;

Para el cálculo de los coeficientes de la ecuación de saturación modificamos el

adaptador FSES1, al igual que en el código anterior, dichas modificaciones se realizan
en el apartado 3D como se muestra a continuación:
multx = perme 11 ( i , j , k i ) ∗ dt / ( p o r o s i t y ∗ dx ∗ dx ∗ ( 1 − Srw − Sro ) ∗ mu w ) ;
multy = perme 22 ( i , j , k i ) ∗ dt / ( p o r o s i t y ∗ dy ∗ dy ∗ ( 1 − Srw − Sro ) ∗ mu w ) ;
multz = perme 33 ( i , j , k i ) ∗ dt / ( p o r o s i t y ∗ dz ∗ dz ∗ ( 1 − Srw − Sro ) ∗ mu w ) ;
aE ( i , j , k i ) = ( Sw e − Srw ) ∗ multx ;
aW ( i , j , k i ) = ( Sw w − Srw ) ∗ multx ;
aN ( i , j , k i ) = ( Sw n − Srw ) ∗ multy ;
aS ( i , j , k i ) = ( Sw s − Srw ) ∗ multy ;
aF ( i , j , k i ) = ( Sw f − Srw ) ∗ multz ;
aB ( i , j , k i ) = ( Sw b − Srw ) ∗ multz ;
aP ( i , j , k i ) = aE ( i , j , k i ) + aW ( i , j , ki ) +
aN ( i , j , k i ) + aS ( i , j , ki ) +
aF ( i , j , k i ) + aB ( i , j , ki ) ;

Concluidas las adaptaciones al código para llevar a cabo el cálculo de los coeficientes
3D de las ecuaciones de presión y saturación, implementamos el cálculo de la
permeabilidad variable en el dominio de estudio como se muestra en la siguiente fracción
de código:
1 S c a l a r F i e l d 3 D k 1 1 ( mesh . g e t E x t e n t N o d e s ( ) ) ; // v a l o r e s de p e r m e a b i l i d a d
2 S c a l a r F i e l d 3 D k 2 2 ( mesh . g e t E x t e n t N o d e s ( ) ) ; // v a l o r e s de p e r m e a b i l i d a d
3 S c a l a r F i e l d 3 D k 3 3 ( mesh . g e t E x t e n t N o d e s ( ) ) ; // v a l o r e s de p e r m e a b i l i d a d
4 double xmin = num nodes x / 3 ;
5 double xmax = 2 ∗ num nodes x / 3 ;
6 double ymin = num nodes y / 3 ;
7 double ymax = 2∗ num nodes y / 3 ;
8 double zmin = 0 ;
9 double zmax = num nodes z ;
10 f o r ( i n t i = 0 ; i < num nodes x ; i++ ) {
11 f o r ( i n t j = 0 ; j < num nodes y ; j ++){
12 f o r ( i n t k = 0 ; k < num nodes z ; k++){
13 i f ( i < num nodes x / 2 ) { s r a n d ( i ∗ 2 ) ; }
14 else srand ( i ∗2) ;
15 k 1 1 ( i , j , k ) = p e r m e a b i l i t y ∗ rand ( ) ;
16 i f ( i >= num nodes x / 2 ) { s r a n d ( j ∗ 3 ) ; }
17 else srand ( j ∗3) ;
18 k 1 1 ( i , j , k ) = p e r m e a b i l i t y ∗ rand ( ) ;
19 i f ( j < num nodes y / 2 ) { s r a n d ( i +3) ; }
20 e l s e s r a n d ( i +3) ;
21 k 2 2 ( i , j , k ) = p e r m e a b i l i t y ∗ rand ( ) ;
22 i f ( j >= num nodes y / 2 ) { s r a n d ( j ∗ 2 ) ; }

84
 6.2. Implementación

23 else srand ( j ∗2) ;

24 k 2 2 ( i , j , k ) = p e r m e a b i l i t y ∗ rand ( ) ;
25 srand ( k ∗2) ;
26 k 3 3 ( i , j , k ) = p e r m e a b i l i t y ∗ rand ( ) ;
27 i f ( i > xmin && i < xmax && j > ymin && j < ymax && k > zmin && k < zmax ) {
28 k 1 1 ( i , j , k ) ∗= 0 . 1 2 4 5 7 ;
29 k 2 2 ( i , j , k ) ∗= 0 . 1 ;
30 k 3 3 ( i , j , k ) ∗= 0 . 1 ;
31 }
32 }
33 }
34 }

La permeabilidad es un campo que se ha guardado de manera independiente para

cada dirección, pues el tensor de permeabilidad presenta cambios en las 3 direcciones
principales, por ello, en las lı́neas 1 a 3 definimos campos 3D que almacenen los valores
de permeabilidad por nodo para cada componente. Las variables declaradas en las
lı́neas 4-9 nos servirán en la construcción del cuerpo anómalo, figura 6.6. El loop
que inicia en la lı́nea 10 y concluye en la 34, barre las dimensiones del dominio de
estudio, mientras que los condicionales de las lı́neas 13, 16, 19, 22 y 27 establecen las
condiciones necesarias para la construcción de dicho cuerpo y su permeabilidad interna
(lı́neas 27-30) y la permeabilidad aleatoria del resto del dominio.
Una vez que la permeabilidad ya es variable y los coeficientes tienen en cuenta éstos
valores, se incorpora en el código el llamado a los valores de permeabilidad, lı́neas 51-
53 para la presión y 73-75 para la saturación y resolvemos como se hizo en el caso
homogéneo, los resultados se muestran en la figura 6.7.
32 S p a r s e M a t r i x < Diagonal <double , 3> > A( num nodes x , num nodes y , num nodes z ) ;
33 ScalarField3D b ( num nodes x , num nodes y , num nodes z ) ;
34 TwoPhaseEquation< FSIP1<double , 3> > p r e s s u r e ( p , A, b , mesh . g e t D e l t a s ( ) ) ;
35 p r e s s u r e . s e t D e l t a T i m e ( dt ) ;
36 pressure . setPermeability ( permeability ) ;
37 pressure . setPorosity ( porosity ) ;
38 p r e s s u r e . s e t S r w ( Srw ) ;
39 p r e s s u r e . s e t S r o ( Sro ) ;
40 p r e s s u r e . s e t V i s c o s i t y w (mu w) ;
41 p r e s s u r e . s e t V i s c o s i t y o ( mu o ) ;
42 pressure . setInjection ( injection ) ;
43 p r e s s u r e . setNeumann (LEFT WALL) ;
44 p r e s s u r e . setNeumann (RIGHT WALL) ;
45 p r e s s u r e . setNeumann (TOP WALL) ;
46 p r e s s u r e . setNeumann (BOTTOM WALL) ;
47 p r e s s u r e . setNeumann (FRONT WALL) ;
48 p r e s s u r e . setNeumann (BACK WALL) ;
49 p r e s s u r e . s e t S a t u r a t i o n (Sw) ;
50 p r e s s u r e . setPorosityV ( phi ) ;
51 p r e s s u r e . setPermeabilityV 11 ( k 11 ) ;
52 p r e s s u r e . setPermeabilityV 22 ( k 22 ) ;
53 p r e s s u r e . setPermeabilityV 33 ( k 33 ) ;
54 pressure . print () ;

85
 6. Modelo Computacional y Resultados Numéricos

56 TwoPhaseEquation< FSES1<double , 3> > s a t u r a t i o n (Sw , A, b , mesh . g e t D e l t a s ( ) ) ;

57 s a t u r a t i o n . s e t D e l t a T i m e ( dt ) ;
58 saturation . setPermeability ( permeability ) ;
59 saturation . setPorosity ( porosity ) ;
60 s a t u r a t i o n . s e t S r w ( Srw ) ;
61 s a t u r a t i o n . s e t S r o ( Sro ) ;
62 saturation . setInjection ( injection ) ;
63 s a t u r a t i o n . s e t V i s c o s i t y w (mu w) ;
64 s a t u r a t i o n . s e t V i s c o s i t y o ( mu o ) ;
65 s a t u r a t i o n . setNeumann (LEFT WALL) ;
66 s a t u r a t i o n . setNeumann (RIGHT WALL) ;
67 s a t u r a t i o n . setNeumann (TOP WALL) ;
68 s a t u r a t i o n . setNeumann (BOTTOM WALL) ;
69 s a t u r a t i o n . setNeumann (FRONT WALL) ;
70 s a t u r a t i o n . setNeumann (BACK WALL) ;
71 saturation . setPressure (p) ;
72 s a t u r a t i o n . setPorosityV ( phi ) ;
73 s a t u r a t i o n . setPermeabilityV 11 ( k 11 ) ;
74 s a t u r a t i o n . setPermeabilityV 22 ( k 22 ) ;
75 s a t u r a t i o n . setPermeabilityV 33 ( k 33 ) ;
76 saturation . print () ;

A continuación se muestran imágenes a distintos tiempos de la simulación.

Figura 6.7: Evolución del frente de saturación para un medio de permeabilidad variable.

La imagen anterior nos muestra un frente de saturación, que a diferencia del caso
anterior, ya no evoluciona de manera geométrica, esto es por la heterogeneidad del
medio al introducir una permeabilidad aleatoria en el dominio de estudio. Dado que
el agua encuentra más difı́cil su recorrido por el domino, observamos que los dı́as de

86
 6.2. Implementación

simulación son más, del orden de 4800 dı́as, a diferencia de los mil del caso anterior,
para llegar al pozo extractor. Nótese que el cuerpo de baja permeabilidad impacta
directamente sobre el frente de saturación, el agua se ve obligada a rodear el cuerpo y
muy poca penetra en él.

6.2.5. Lı́neas de Corriente caso no-homogéneo

En la simulación con lı́neas de corriente es importante calcular el campo de velocidad
en todo el dominio, pues las lı́neas de corriente son tangentes en cada punto al campo
local. En este apartado se resolverá el problema con los datos de la tabla 6.2 para un
medio no-homogéneo. A continuación la implementación en TUNAM.
1 typedef TunaArray<double , 3 > : : huge S c a l a r F i e l d 3 D ;
2 typedef TunaArray<double , 1 > : : huge S c a l a r F i e l d 1 D ;

A diferencia de la codificación para volumen finito, en la simulación con lı́neas de

corriente, además de declarar un arreglo 3D es necesario definir uno unidimensional,
(lı́nea 2), pues recordemos que la saturación se resolverá sobre cada lı́nea de corriente.
Como en los casos anteriores, creamos una malla numérica de la siguiente manera:
3 StructuredMesh<Uniform<double , 3> > mesh ( l e n g t h x , num nodes x ,
4 l e n g t h y , num nodes y ,
5 l e n g t h z , num nodes z ) ;
6 mesh . p r i n t ( ) ;

La malla sigue siendo del mismo número de nodos, el dominio tiene 80 nodos en
las direcciones x y y y 4 en la dirección z.
13 ScalarField3D p ( mesh . getExtentVolumes ( ) ) ;
14 ScalarField3D p n ( mesh . g e t E x t e n t N o d e s ( ) );
15 ScalarField3D u1 ( num nodes x , num vols y , n u m v o l s z ) ;
16 ScalarField3D u2 ( num vols x , num nodes y , n u m v o l s z ) ;
17 ScalarField3D u3 ( num vols x , num vols y , num nodes z ) ;
18 ScalarField3D u1 n ( mesh . g e t E x t e n t N o d e s ( ) );
19 ScalarField3D u2 n ( mesh . g e t E x t e n t N o d e s ( ) );
20 ScalarField3D u3 n ( mesh . g e t E x t e n t N o d e s ( ) );

En las lı́neas 13-20 vemos la declaración de los campos escalares que almacenarán
la información concernientes tanto a la presión como a la velocidad. Hay que hacer una
distinción especial para el campo de la velocidad, ya que los campos con nomenclatura
u1, u2 y u3 son atribuidos a cada centro de volumen de control de la malla, de modo
que hay que interpolar a los nodos y para eso reservamos un espacio en memoria con
los nombres u1 n, u2 n y u3 n.
21 // I n i c i a l i z a m o s campos para e l c a l c u l o s o b r e l a s l i n e a s de c o r r i e n t e
22 S c a l a r F i e l d 1 D SwSLS [ N l i n e s S L S ] ;

En la lı́nea 22 declaramos el campo unidimensional que almacenará la solución de

la ecuación de saturación sobre las lı́neas de corriente. En la sección anterior, se hizo el
cálculo variable de la permeabilidad, ası́ que incorporamos esos valores en el cálculo de
la velocidad (lı́neas 24-67) como se muestra en el siguiente fragmento de código, una
visualización del cálculo se observa en la figura 6.8.

87
 6. Modelo Computacional y Resultados Numéricos

24 // C a l c u l o de l a v e l o c i d a d
25
26 double Sef = (Sw( i , j , kk )−Srw ) /(1−Srw−Sro ) ;
27 double krw = Sef ;
28 double kro = 1− S e f ;
29 double lam T = ( krw/mu w) + ( k r o /mu o ) ;
30 double lam w = krw / mu w ;
31
32 f o r ( i n t i = b i ; i <= e i ; ++i ) {
33 f o r ( i n t j = b j ; j <= e j ; ++j ) {
34 f o r ( i n t k = bk ; k <= ek ; ++k ) {
35 S e f = (Sw( i , j , kk )−Srw ) /(1−Srw−Sro ) ;
36 krw = S e f ;
37 k r o = 1− S e f ;
38 lam T = ( krw/mu w) + ( k r o /mu o ) ;
39 lam w = krw / mu w ;
40 u1 ( i , j , k ) = −k 1 1 ( i , j , k ) ∗ lam T ∗ ( ( p ( i +1 , j , k )−p ( i −1 , j , k ) ) / ( 2 ∗ dx ) ) ;
41 }
42 }
43 }
44 f o r ( i n t i = b i ; i <= e i ; ++i ) {
45 f o r ( i n t j = b j ; j <= e j ; ++j ) {
46 f o r ( i n t k = bk ; k <= ek ; ++k ) {
47 S e f = (Sw( i , j , kk )−Srw ) /(1−Srw−Sro ) ;
48 krw = S e f ;
49 k r o = 1− S e f ;
50 lam T = ( krw/mu w) + ( k r o /mu o ) ;
51 lam w = krw / mu w ;
52 u2 ( i , j , k ) = −k 2 2 ( i , j , k ) ∗ lam T ∗ ( ( p ( i , j +1 ,k )−p ( i , j −1 ,k ) ) / ( 2 ∗ dy ) ) ;
53 }
54 }
55 }
56 f o r ( i n t i = b i ; i <= e i ; ++i ) {
57 f o r ( i n t j = b j ; j <= e j ; ++j ) {
58 f o r ( i n t k = bk ; k <= ek ; ++k ) {
59 S e f = (Sw( i , j , kk )−Srw ) /(1−Srw−Sro ) ;
60 krw = S e f ;
61 k r o = 1− S e f ;
62 lam T = ( krw/mu w) + ( k r o /mu o ) ;
63 lam w = krw / mu w ;
64 u3 ( i , j , k ) = −k 3 3 ( i , j , k ) ∗ lam T ∗ ( ( p ( i , j , k+1)−p ( i , j , k−1) ) / ( 2 ∗ dz ) ) ;
65 }
66 }
67 }

La velocidad se calcula en tres secciones de bucles, uno para cada componente de

la velocidad, u1, u2 y u3 para las direcciones x, y y z respectivamente. Dichos valores
se encuentran en los volúmenes de control, pero habrán de ser interpolados a los nodos
de la malla para fines de visualización, a continuación las lı́neas de código que llevan a
cabo dicha interpolación.
68 NumUtils : : i n t e r p o l a t e T o N o d e s U s ( u1 n , u1 ) ;
69 NumUtils : : i n t e r p o l a t e T o N o d e s V s ( u2 n , u2 ) ;
70 NumUtils : : i nt e rp o la t eT o No d es W s ( u3 n , u3 ) ;

88
 6.2. Implementación

Figura 6.8: Campo de velocidad para el dominio con permeabilidad variable.

La imagen anterior muestra el resultado del cálculo de la velocidad con la

permeabilidad variable. Al centro del dominio observamos una zona de baja velocidad
que corresponde al cuerpo insertado de baja permeabilidad, dado que la velocidad de
Darcy depende de la permeabilidad, ésta se ve reflejada en los valores que adquiere
sobre el flujo de los fluidos al centro del dominio, donde esperamos poca penetración
del frente de agua pues la velocidad es prácticamente cero. Finalmente impactará en el
cálculo de la saturación sobre las lı́neas de corriente, ver [29]; dicho cálculo corresponde
a la siguiente codificación:
71 double S e f s 1 , krw s1 , k r o s 1 , lam T s1 , l a m w s 1 ;
72 double S e f s 2 , krw s2 , k r o s 2 , lam T s2 , l a m w s 2 ;
73 f o r ( i n t i = 0 ; i < N l i n e s S L S ; i ++) {
74 f o r ( i n t j = 1 ; j < NpointsSLS −1; j ++) {
75 S e f s 1 = ( SwSLS [ i ] ( j )−Srw ) /(1−Srw−Sro ) ;
76 krw s1 = S e f s 1 ;
77 k r o s 1 = (1− S e f s 1 ) ;
78 la m T s 1 = ( k r w s 1 /mu w) + ( k r o s 1 /mu o ) ;
79 la m w s 1 = k r w s 1 / mu w ;
80
81 Sef s2 = ( SwSLS [ i ] ( j −1)−Srw ) /(1−Srw−Sro ) ;
82 krw s2 = Sef s2 ;
83 kro s2 = (1− S e f s 2 ) ;
84 la m T s2 = ( k r w s 2 /mu w) + ( k r o s 2 /mu o ) ;
85 la m w s2 = k r w s 2 / mu w ;
86
87 Dtao = c a l c D t a o ( u1 n , u2 n , u3 n , x s l s , y s l s , z s l s , p o r o s i t y , dx , dy , dz , i , j ) ;
88 SwSLS [ i s ] ( j s ) = SwSLS [ i s ] ( j s ) − ( ( l a m w s 1 / l a m T s 1 ) − ( l a m w s 2 / l a m T s 2 ) )
89 ∗ ( dt / Dtao ) ;
90 i f ( SwSLS [ i s ] ( j s ) < 1 . 0 e −10) SwSLS [ i s ] ( j s ) = 0 . 0 ;
91 i f ( SwSLS [ i s ] ( j s ) > 0 . 8 ) SwSLS [ i s ] = 0 . 8 ;
92 }
93 }

89
6. Modelo Computacional y Resultados Numéricos

En el listado de las lı́neas 71-93 observamos el cálculo de la saturación sobre las

lı́neas de corriente. En particular, la lı́nea 87 hace un llamado a la función, calcDtao,
que es la encargada de calcular el tiempo de vuelo para cada lı́nea de corriente, a
continuación su implementación:
double c a l c D t a o ( S c a l a r F i e l d 3 D& u1 , S c a l a r F i e l d 3 D& u2 , S c a l a r F i e l d 3 D& u3 ,
double ∗∗ x s l s , double ∗∗ y s l s , double ∗∗ z s l s , double p o r o s i t y ,
double dx , double dy , double dz , i n t i , i n t j )
{
double a = x s l s [ i ] [ j ] / dx ;
double b = y s l s [ i ] [ j ] / dy ;
double c = z s l s [ i ] [ j ] / dz ;
i n t I = s t a t i c c a s t <int >(a ) ;
i n t J = s t a t i c c a s t <int >(b ) ;
i n t K = s t a t i c c a s t <int >(c ) ;
double a l p h a = a − I ;
double b e t a = b − J ;
double gamma = c − K;
double UI , VI , WI ;
UI = ( ( u1 ( I ,J ,K ) ∗ ( 1 − a l p h a ) +
u1 ( I +1 ,J ,K ) ∗ a l p h a ) ∗ ( 1 − b e t a ) +
( u1 ( I , J+1 ,K ) ∗ ( 1 − a l p h a ) +
u1 ( I +1 ,J+1 ,K ) ∗ a l p h a ) ∗ b e t a ) ∗ ( 1 − gamma) +
( ( u1 ( I ,J ,K+1) ∗ ( 1 − a l p h a ) +
u1 ( I +1 ,J ,K+1) ∗ a l p h a ) ∗ ( 1 − b e t a ) +
( u1 ( I , J+1 ,K+1) ∗ ( 1 − a l p h a ) +
u1 ( I +1 ,J+1 ,K+1) ∗ a l p h a ) ∗ b e t a ) ∗ gamma ;
VI = ( ( u2 ( I ,J ,K ) ∗ ( 1 − a l p h a ) +
u2 ( I +1 ,J ,K ) ∗ a l p h a ) ∗ ( 1 − b e t a ) +
( u2 ( I , J+1 ,K ) ∗ ( 1 − a l p h a ) +
u2 ( I +1 ,J+1 ,K ) ∗ a l p h a ) ∗ b e t a ) ∗ ( 1 − gamma) +
( ( u2 ( I ,J ,K+1) ∗ ( 1 − a l p h a ) +
u2 ( I +1 ,J ,K+1) ∗ a l p h a ) ∗ ( 1 − b e t a ) +
( u2 ( I , J+1 ,K+1) ∗ ( 1 − a l p h a ) +
u2 ( I +1 ,J+1 ,K+1) ∗ a l p h a ) ∗ b e t a ) ∗ gamma ;
WI = ( ( u3 ( I ,J ,K ) ∗ ( 1 − a l p h a ) +
u3 ( I +1 ,J ,K ) ∗ a l p h a ) ∗ ( 1 − b e t a ) +
( u3 ( I , J+1 ,K ) ∗ ( 1 − a l p h a ) +
u3 ( I +1 ,J+1 ,K ) ∗ a l p h a ) ∗ b e t a ) ∗ ( 1 − gamma) +
( ( u3 ( I ,J ,K+1) ∗ ( 1 − a l p h a ) +
u3 ( I +1 ,J ,K+1) ∗ a l p h a ) ∗ ( 1 − b e t a ) +
( u3 ( I , J+1 ,K+1) ∗ ( 1 − a l p h a ) +
u3 ( I +1 ,J+1 ,K+1) ∗ a l p h a ) ∗ b e t a ) ∗ gamma ;

double Dxi = s q r t ( ( x s l s [ i ] [ j ] − x s l s [ i ] [ j −1]) ∗ ( x s l s [ i ] [ j ] − x s l s [ i ] [ j −1]) +

( y s l s [ i ] [ j ] − y s l s [ i ] [ j −1]) ∗ ( y s l s [ i ] [ j ] − y s l s [ i ] [ j −1]) +
( z s l s [ i ] [ j ] − z s l s [ i ] [ j −1]) ∗ ( z s l s [ i ] [ j ] − z s l s [ i ] [ j −1]) ) ;
double Dtao = ( Dxi ∗ p o r o s i t y ) / s q r t ( UI∗UI + VI∗VI + WI∗WI) ;
return Dtao ;
}

Recibe como argumentos los campos escalares de velocidad para las componentes
x, y y z, las coordenadas de las lı́neas, la porosidad y los deltas espaciales en cada
dirección. Se trata de una interpolación trilineal del campo de velocidad local a la
posición j de la lı́nea i en cuestión. La velocidad interpolada de los nodos de la malla
a los puntos de cada lı́nea de corriente se hace por componentes, de modo que UI, VI
y WI son los valores para la dirección x, y y z respectivamente. Esta función es para
calcular el tiempo de vuelo según la ecuación (4.33). Los resultados que se muestran
en la figura 6.9 son una comparación cualitativa de los métodos MVF y SLS.

90
 6.2. Implementación

91
6. Modelo Computacional y Resultados Numéricos

Figura 6.9: Comparación de los métodos MVF y SLS en los mismos dı́as de la
simulación.

En principio podemos observar que ambos métodos son sensibles al cambio de

permeabilidad, en particular para el cuerpo anómalo al centro del dominio. La baja
velocidad de esta zona repercute en el frente de saturación de los dos métodos,
identificando zonas de drenado pobre de hidrocarburos.

92
Capı́tulo 7

Conclusiones

En el presente trabajo se describieron los modelos matemáticos a partir de la

formulación axiomática y los modelos computacionales para describir un flujo de dos
fases a través de un medio poroso.
Se realizaron adaptaciones al software TUNAM para realizar cálculos de velocidad
sobre medios de permeabilidad variable y resolver la ecuación de saturación sobre
lı́neas de corriente. Se llevaron a cabo simulaciones usando los métodos de volumen
finito y lı́neas de corriente a fin de comparar los resultados de cada uno. Contrastando
los métodos de MVF y SLS, observamos un comportamiento cualitativo similar para
los pasos de tiempo inferiores a los 2000 dı́as de simulación. Aunque ambos métodos
son sensibles al cambio en la permeabilidad del cuerpo al centro del dominio con bajos
valores de esta propiedad, observamos que el frente de saturación de agua en las lı́neas
de corriente evoluciona de manera distinta al método MVF, una posible causa son
las interpolaciones que se hacen de la malla numérica, es decir, de los centros de los
volúmenes de control a los nodos, y de la propia velocidad a las lı́neas de corriente.
Estas interpolaciones acarrean errores de aproximación que se acumulan con el paso
del tiempo en los dı́as de simulación, introduciendo un desajuste en el cálculo de la
saturación del agua. Queda realizar una interpretación cuantitativa de este fenómeno,
llevando el valor de la saturación de cada lı́nea de corriente a la malla numérica y poder
hacer una comparación al respecto contra los resultados mostrados con el método de
MVF.
El método de lı́neas de corriente se muestra como un método para identificar zonas
de drenado pobre de hidrocarburos pues la evolución de los frentes de inyección y
su interacción con las heterogeneidades del yacimiento pueden ser visualizadas fácil y
rápidamente, y por lo tanto proveen de una manera natural e intuitiva para caracterizar
dinámicamente un yacimiento. La optimización de la locación de Pozos basada en el
modelo geológico permite reflejar la geometrı́a y heterogeneidad de los reservorios más

93
7. Conclusiones

detalladamente. La aportación del presente trabajo deja un código computacional en

el lenguaje de programación C++ que es de carácter abierto y queda a disposición
de la comunidad el mejorarlo, adaptarlo y compartirlo, pues cabe señalar que existen
simuladores comerciales de carácter general que incorporan este procedimiento, sin
embargo, no son accesibles a todo público. El problema computacional es de gran escala
y se recomienda paralelizar los cálculos de la saturación sobre las lı́neas de corriente
y el gradiente conjugado para resolver el sistema lineal. Los modelos numéricos que
se utilizan para la simulación de la extracción mejorada de petróleo están concebidos
para que puedan adaptarse a sistemas de cómputo paralelo.

94
Apéndice A

Formulación Presión-Saturación

Un sistema de N fases para un flujo fraccional provoca que el sistema de ecuaciones

diferenciales parciales se desacople. En un sistema desacoplado tenemos una ecuación
para la presión y N − 1 ecuaciones de transporte para las saturaciones. Para resolver
sistemas de ecuaciones donde éstas se encuentran débilmente acopladas, empleamos
métodos iterativos. Bajo esta formulación, la ecuación de presión se resuelve de manera
independiente de la de saturación, de modo que los resultados obtenidos de la ecuación
de presión se usan para resolver las ecuaciones de transporte para la saturación. Como
consecuencia, los resultados de la ecuación de saturación se insertan en la ecuación
de presión a fin de obtener la solución en el paso de tiempo siguiente. Dicho proceso
iterativo se repetirá para un número de pasos de tiempo previamente definido. Las
siguientes secciones muestran la manera de obtener una formulación presión-saturación
para sistemas multifásicos.

A.1. Ecuación de Presión

La ecuación de presión se deriva a partir de la ecuación (4.7), para ello, dividimos
por ρα , obteniendo:
 
1 ∂(φρα Sα )
+ ∇ · (ρα uα ) − qα = 0,
ρα ∂t

sumando ahora para todas las fases del sistema,

X  1  ∂(φρα Sα ) 
+ ∇ · (ρα uα ) − qα = 0,
α
ρα ∂t

95
A. Formulación Presión-Saturación

desarrollando la derivada temporal se tiene que:

X 1  ∂φ ∂ρα ∂Sα

ρα Sα + φSα + ρα φ + ρα ∇ · uα + uα · ∇ρα − qα = 0,
α
ρα ∂t ∂t ∂t

reordenando términos y aplicando la ecuación (4.10) obtenemos:

∂φ X X 1  ∂ρα
 X
qα
+ ∇ · uα + φSα + uα · ∇ρα − = 0. (A.1)
∂t α α
ρα ∂t α
ρα

La velocidad total del sistema se define como:

X
u= uα , (A.2)
α

aplicando el operador divergencia obtenemos:

X X
∇·u=∇· uα = ∇ · uα . (A.3)
α α

Insertando (A.3) en (A.1) llegamos a la siguiente ecuación,

∂φ X 1  ∂ρα
 X
qα
+∇·u+ φSα + uα · ∇ρα − = 0. (A.4)
∂t α
ρ α ∂t α
ρ α

De la definición de la Ley de Darcy para un flujo miltifásico tenemos que:

X Xh i
u= uα = −kλα (∇ · pα − ρα g) . (A.5)
α α

Definamos ahora una función de flujo fraccional para la fase α, de manera que,
λα
fα = =⇒ λα = fα λ, (A.6)
λ
P
P donde λ = α λα representa la movilidad total y como consecuencia se cumple que
α fα = 1. Sustituyendo la ecuación (A.6) en (A.5) llegamos a:
" #
X X
u = −λk fα ∇pα − fα ρ α g . (A.7)
α α

Con esta expresión de la velocidad total, la ecuación (A.4) se transforma en una

ecuación para la presión que denotaremos como:

" #
∂φ X X X 1  ∂ρα
 X
qα
− ∇ · λk fα ∇pα − fα ρ α g + φSα + uα · ∇ρα − = 0.
∂t α α α
ρα ∂t α
ρα
(A.8)

96
 A.2. Ecuación de Saturación

A.2. Ecuación de Saturación

Al dar solución a la ecuación (A.8), obtenemos la presión que ayudará a calcular
la velocidad que emplearemos para la ecuación de saturación. Es posible calcular la
saturación directamente de la ecuación (4.7). En un sistema de dos fases; agua y aceite,
se utiliza la presión del aceite a fin de generar una ecuación para la presión, mientras
que la ecuación de saturación de la fase agua se expresa como

∂(φρw Sw )
+ ∇ · (ρw uw ) = qw , (A.9)
∂t
Hay que hacer notar que tanto la ecuación (A.8) como la ecuación (4.7) fueron
formuladas en términos de una presión de fase (la del aceite por ejemplo).

97
Apéndice B

Gradiente Conjugado para Sistemas

Lineales Dispersos

El método de gradiente conjugado (CGM) por sus siglas en inglés, es un método

efectivo para sistemas simétricos y positivo-definidos, esto es A = AT y uT Au > 0
respectivamente, [37].
Descrito en una frase, el método es una realización de una técnica de proyección
ortogonal sobre el subespacio de Krylov K(r(0) , A), donde r(0) es el residuo inicial.
La idea esencial del método consiste en construir una base de vectores ortogonales
en dicho subespacio y emplearla para realizar la búsqueda iterativa de la solución. El
vector solución puede ser expresado como:

x(j+1) = x(j) + α(j) p(j) .

Con el fin de ajustarse a la notación estándar que se utiliza en la literatura para

describir el algoritmo, el ı́ndice de los vectores p ahora comienza en cero en lugar de
uno, como solı́a acostumbrarse. Ahora bien, en CGM los vectores residuales deben
satisfacer la recurrencia,

r(j+1) = r(j) − α(j) Ap(j) . (B.1)

Si se requiere que los vectores r(j) sean ortogonales, entonces es necesario que
(r − α(j) Ap(j) , r(j) ) = 0 y, como resultado:
(j)

(j) r(j) , r(j)
α = . (B.2)
Ap(j) , r(j)

99
B. Gradiente Conjugado para Sistemas Lineales Dispersos

Además, es conocido que la siguiente dirección de búsqueda p(j+1) es una

combinación lineal de r(j+1) y p(j) , después ajustando la base de los vectores p
apropiadamente, se deduce que,

p(j+1) = r(j+1) − β (j) p(j) . (B.3)

Por lo tanto, una primara consecuencia de la relación anterior es que,
     
Ap(j) , r(j) = Ap(j) , p(j) − β (j−1) p(j−1) = Ap(j) , p(j) ,

porque Ap(j) es ortogonal a p(j−1) . Entonces (B.2) se convierte en:

(j) r(j) , r(j)
α = ,
Ap(j) , p(j)

Por otra parte, escribiendo p(j+1) como se definió en (B.3) y dado que es ortogonal
a Ap(j) , tenemos que:
 
r(j+1) , Ap(j)
β (j) = −   .
p(j) , Ap(j)

Nótese que de (B.1), tenemos que:

1
Ap(j) = − r(j+1) − r(j) ,


α (j)

y, por lo tanto




(j) 1 r(j+1) , r(j+1) − r(j) r(j+1) , r(j+1)
β = (j)   = .
α (j)
Ap , p (j) (r(j) , r(j) )

Poniendo todas estas relaciones juntas da como resultado un algoritmo iterativo,

como se muestra a continuación:

Algoritmo 2 CGM
1: Calcular r (0) = b − Ax(0) , p(0) = r (0)
2: for j = 0, 1, ..., hastaconvergerdo


3: α(j) = r(j) , r(j) / Ap(j) , p(j)
4: x(j+1) = x(j) + α(j) p(j)
5: r(j+1) = r(j) − α(j) Ap(j)



6: β (j) = r(j+1) , r(j+1) / r(j) , r(j)
7: p(j+1) = r(j+1) − β (j) p(j)
8: end for

100
 En este algoritmo, (·, ·) es el producto interno adecuado al sistema lineal en
particular. Se debe considerar, en términos de memoria, almacenar 4 vectores (x, p, Ap
y r), la solución aproximada será x(j+1) y el vector residual es r(j+1) .
La gran ventaja del Método de Gradiente Conjugado radica en que cuando se utiliza
este procedimiento basta con asegurar la ortogonalidad de un nuevo miembro con
respecto al último que se ha construido, para que automáticamente esta condición se
cumpla con respecto a todos los anteriores.

101
Apéndice C

Aplicación del Método de Volumen

Finito

En este apéndice se describe la aplicación del MVF a problemas de una fase, siendo
extensible a más fases en un sistema.

Se construye una malla con volúmenes que no se traslapen como se muestra en

la figura 5.1(a).

Integramos en el espacio del volumen de control de la figura 5.1(b), y en el

tiempo del instante t al tiempo t + ∆t, equivalentes a los pasos de tiempo n y
n + 1 respectivamente.

Zn Ze Z
n+1 n+1Zn Ze
∂ 2p ∂ 2p
Z  
∂p
dt dx dy = Γ + dx dy dt.
∂t ∂x2 ∂y 2
s w n n s w

donde Γ = k/φµcT .

Se aproximan las integrales

103
C. Aplicación del Método de Volumen Finito

n+1
Z e  n 
∂p ∂p
(pn+1
P − pnP )∆x∆y = Γ ∆y + ∆x dt
∂x w ∂y s
n
n+1  
pE − pP pP − pW
Z
= Γ − ∆y
∆xe ∆xw
n
  !
pN − p P pP − pS
+ − ∆x dt
∆yn ∆ys
n+1
Z
= aE pE + aW pW + aN pN + aS pS
n
  !
− aE + aW + aN + aS pP dt.

∆y ∆y ∆x ∆x
donde aE = Γ ∆x e
; aW = Γ ∆x w
; aN = Γ ∆y n
; aS = Γ ∆y s
.
La integral temporal se aproxima usando el esuqema θ:
Dada una función f (x), la integral de n a n + 1 de dicha función, se aproxima de
la siguiente manera
n+1
Z
f dt = [θf n+1 + (1 − θ)f n ]∆t,
n

1
donde para θ = 0 se tiene un esquema explı́cito (f n ∆t) y para θ = 2
se tiene un
esquema conocido como de Crank-Nicolson ([f n+1 + f n ] ∆t2
).
Usando un esquema implı́cito (θ = 1) se tiene que

aP pn+1
P = aE pn+1
E + aW pn+1 n+1
W + aN pN + aS pn+1
S + sP pnP ,

∆x∆y ∆x∆y
donde aP = aE + aW + aN + aS + ∆t
y sP = ∆t
.

Para cada volumen de control de la malla se obtiene una ecuación discreta, lo que
genera un sistema lineal de ecuaciones. Dicho sistema tiene la forma mostrada
en la figura 5.2. La matriz del sistema es rala, simétrica y positivo definida.
Existen varios algoritmos para resolver este tipo de sistemas, en este trabajo se
empleó el método de gradiente conjugado, ver apéndice B.

104
Referencias

[1] Abou-Kassem, J. H., Farouq, S. M., Rafiq, M., “Petroleum Reservoir

Simulation. A Basic Approach.”, Gulf Publishing Company, Houston, Texas, 2006.

[2] Allen, M. B., Herrera, I., Pinder, G., “Numerical Modeling in Science and
Engineering.”, 1988.

[3] Al-Zawawi, A. S., et al, “Using Streamline and Reservoir Simulation to

Improve Waterflood Management”, Saudi Aramco Journal of Technology, 2011.

[4] Ames, W. F., “Numerical Methods for Partial Differential Equations.”, Academic
Press, INC., 1977.

[5] Baker, R., “Streamline Technology: Reservoir History Matching and Forecasting
= Its Succes, Limitations, and Future”, Distinguished Author Series, Journal of
Canadian Petroleum Technology, volumen 40 No. 4.

[6] Batycky, R.P., Thiele, M. R., Blunt, M. J., “A Streamline-Based Reservoir

Simulator of the House Mountain Waterflood.”, SCRF, 1997.

[7] Batycky, R.P., “A three-dimensional two-phase field scale streamline simula-

tor.”, Stanford University, Ph D. Thesis. 1997.

[8] Batycky, R.P., Förster, M., Thiele, M.R., Stüben, K., “Parallelization
of a Commercial Streamline Simulator and Performance on Practical Models.”,
Society of Petroleum Engineers, 2010.

[9] Berre, I., Dahle, H. K., Karlsen, K. H., Nordhaug, H. F., “A Streamline
Front Tracking Method for Two- and Three-Phase Flow Including Capillary
Forces.”, American Mathematical Society, 2000.

[10] Bj∅rlykke, Knut., “Petroleum Geoscience: From Sedimentary Enviroments to

Rock Physics.”, Springer, 2010.

105
Referencias

[11] Blitz++: Object-oriented library for scientific computing.

[Link]

[12] Chapman, R.E., “Petroleum Geology”, Elsevier, 1983.

[13] Chen, Z., “Reservoir Simulation: Mathematical Techniques in Oil Recovery”,

SIAM, 2007.

[14] Chen, Z., Huan, G., y Ma ,Y., “Computational Methods for Multiphase Flows
in Porous Media”, SIAM, 2006.

[15] Datta-Gupta, A. y King, M., “Streamline Simulation, A current perspective”,

Texas A&M University, 1998.

[16] Datta-Gupta, A. y King, M., “Streamline Simulation: Theory and Practice”,

SPE, 2006.

[17] De la Cruz, L., “Flujo en una y dos fases en medios porosos: modelos
matemáticos, numéricos y computacionales”, Reportes Internos 2012-04 Instituto
de Geofı́sica UNAM, 2012.

[18] De la Cruz, L., “Tunam: Template units for numerical applications and
modeling.”
[Link]

[19] Durran, D. R. “Numerical Methods for Fluid Dynamics with Applications to

Geophysics.”, Springer, 2da Ed. 2010.

[20] Economides, M. J., Nolte, K. G. “Reservoir Stimulation”, Wiley, 3er Ed.

2010.

[21] Ferziger, J. H., Perić, M. “Computational Methods for Fluid Dynamics.”,

Springer, 3er Ed. 2002.

[22] Glowinski, R., Neittaanmäki, P. “Partial Differential Equations, Modeling

and Numerical Simulation.”, Springer, 2008.

[23] Golub, H., Van Loan, C. “Matrix Computations”, The Johns Hopkins
University Press, 3er Ed. 1996.

[24] Hægland, HÅkon “Streamline methods with application to flow and transport
in fractured media.”, Ph. D. Thesis, University of Bergen. 2009.

[25] Hantschel, T., Kauerauf, A. I. “Fundamentals of Basin and Petroleum

Systems Modeling”, Springer, 2009.

[26] Hasle, G., Lie, K., Quak, E. “Geometric Modelling, Numerical Simulation,
and Optimization.”, Springer, 2007.

106
 Referencias

[27] Herrera, I. y Pinder, G., “Mathematical Modeling in Science and Engineering:

An Axiomatic Approach”, John Wiley and Sons, 2012.

[28] Iske, A., Randen, T. “Mathematical Methods and Modelling in Hydrocarbon

Exploration and Production.”, Springer, 2000.

[29] Klausen, R. A., Rasmussen, A. F., Stephansen A. F., “Velocity

interpolation and streamline tracing on irregular geometries.”, Computational
Geosciences, Volume 16, 2012.

[30] Lavorsen, A. I., “Geology of Petroleum”, segunda edición, Freeman, 1967.

[31] Leveque, R. J., “Finite-Difference Methods for Differential Equations”, Univer-

sity of Washington, 2005.

[32] Leveque, R. J., “Finite-Volume Methods for Hyperbolic Problems”, segunda

edición, Freeman, 1967.

[33] Mijáilov, V. P., “Ecuaciones diferenciales en Derivadas Parciales”, Editorial

Mir Moscú, 1978.

[34] OpenMP Web Site: [Link]

[35] OpenDX Homepage: [Link]

[36] Peaceman, D. W., “Fundamentals of numerical reservoir simulation.”, ELSE-

VIER, 1977.

[37] Saad, Y., “Iterative Methods for Sparse Linear Systems”, segunda edición, SIAM,
2003.

[38] Selley, R. C., “Elements of Petroleum Geology”, segunda edición, Academic

Press, 1998.

[39] Shapira, Y., “Solving PDEs in C++. Numerical Methods in a Unified Object-
Oriented Approach.”, SIAM, 2006.

[40] Shonkwiler, R. W., Lefton, L., “An introduction to Parallel and Vector
Scientifc Computation.”, Cambridge University Press, 2006.

[41] Stewart, G. W., “Matrix Algorithms. Volume I: Basic Decompositions.”, SIAM,

1998.

[42] Thiele, M. R., Batycky, R. P. y Fenwick, D. H.“Streamline Simulation for

Modern Reservoir-Engineering Workflows”, Distinguished Author Series, SPE,
2010.

107
Referencias

[43] Vasco, D. W, Yoon, S., Datta-Gupta, A., “Integrating Dynamic Data Into
High-Resolution Reservoir Models Using Streamline-Based Analytic Sensitivity
Coefficients ”, SPE, 1999.

[44] Velten, K., “Mathematical Modeling and Simulation. Introduction for Scientists
and Engineers”, Wiley-VCH, 2009.

[45] Versteeg, H. K, Malalasekera, W., “An introduction to fluid dynamics.

The finite volume method.”, Longman Scientific & Technical, 1995.

[46] Watkins, D. S., “Fundamentals of Matrix Computations”, Segunda Edición,

Wiley-Interscience, 2002.

[47] White, R., “Computational Mathematics”, Chapman & Hall/CRC, 2003.

[48] Yoon, S., Malallah, A. H., Datta-Gupta, A., Vasco, D. W., Behrens,
R. A., , “A Multiscale Approach to Production-Data Integration Using Streamline
Models.”, SPE, 2001.

[49] Zhao, C., Hobbs, B. E., Ord, A., “Fundamentals of Computational

Geoscience. Numerical Methods and Algorithms.”, Springer, 2009.

[50] Zhong, H., Yoon, S., Datta-Gupta, A., “Streamline-Based Production Data
Integration With Gravity and Changing Field Conditions.”, SPE, 2002.

108