• Author / Uploaded
• elliotarmo

• Categories
• Mecánica de Medios Continuos
• Objetos matemáticos
• Física y matemáticas
• Matemáticas
• Conceptos matemáticos

Citation preview

FACULTAD DE I NGENIERÍA -UNAM

Ecuación de Difusión en Coordenadas Radiales: Solución Numérica

Simulación Numérica de Yacimientos Borrador Oct, 2014

C ONTENTS 1

Principios Físicos

1.1 1.2 1.3 1.4

Ecuación de Continuidad Ecuación de Estado . . . . Ecuación de Transporte . Ecuación de Energía . . .

3

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

2

Derivación de la Ecuación de Difusión

3

Solución Analítica de la Ecuación de Difusión para Yacimientos No Fracturados

3.1 Solución Fuente Lineal . . . . . . . . . . . . . . . . . 3.1.1 Integral Exponencial . . . . . . . . . . . . . . 3.1.2 Ejemplo de Aplicación . . . . . . . . . . . . . 3.2 Cálculo de presión en la cara de la formación, p w f 3.2.1 Factor de Daño . . . . . . . . . . . . . . . . . 3.2.2 Caída de Presión Total en la cara del pozo . 3.2.3 Ejemplo de Aplicación . . . . . . . . . . . . . 4

3 3 3 4 5

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

9

10 11 12 15 17 18 20

Solución Numérica de la Ecuación de Difusión para Yacimientos no Fracturados 21

4.1 Contribuciones al Gasto Total en Superficie . . . . . . . . . . . . . . . . . . . . . 4.1.1 Contribución al Gasto por almacenamiento en el Pozo, q w s . . . . . . .

22 22

1

4.2 4.3

4.4

4.5 5

6

4.1.2 Contribución al Gasto por Aporte del yacimiento, q s f . . . . . Discretización de la Ecuación de Difusión . . . . . . . . . . . . . . . . . Malla Radial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Bloques Distribuidos . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Bloques Centrados . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.3 Ejercicio: Generación de Mallas Radiales . . . . . . . . . . . . . . Aproximación de la Ecuación de Difusión mediante Diferencias Finitas 4.4.1 Discretizacíon del Término de Flujo . . . . . . . . . . . . . . . . . 4.4.2 Discretizacíon Término de Acumulación . . . . . . . . . . . . . . 4.4.3 Acoplamiento de Condiciones de Frontera . . . . . . . . . . . . . 4.4.4 Construcción del Sistema de Ecuaciones Ax = b . . . . . . . . . . Algoritmo de Solución . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

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

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

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

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

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

24 28 28 29 29 30 31 31 33 36 40 43

Solución Numérica para Yacimientos Naturalmente Fracturados

47

5.1 Conceptualización de un Yacimientos Naturalmente Fracturado . . . . . . . . . 5.2 Flujo de Monofásico de aceite en Yacimientos Naturalmente Fracturados . . . . 5.2.1 Modelo de Yacimiento No Fracturado. 1φ − 1k . . . . . . . . . . . . . . . . 5.2.2 Modelo de Yacimiento Fracturado. Doble Porosidad - Doble Permeabilidad. 2φ − 2k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.3 Modelo de Yacimiento Fracturado. Modelo de Doble Porosidad - Una Permeabilidad. 2φ − 1k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Discretización del Modelo de Doble Porosidad - Una Permeabilidad . . . . . . . 5.3.1 Término de Transferencia de Fluidos Matriz-Fractura . . . . . . . . . . . 5.3.2 Discretización de medio Continuo - Fracturas . . . . . . . . . . . . . . . . 5.3.3 Acoplamiento de la Ecuación de la Matriz en la Ecuación de la Fractura . 5.3.4 Construcción del Sistema de Ecuaciones Ax = b . . . . . . . . . . . . . . . 5.4 Algoritmo de Solución para el Modelo de 2φ − 1k . . . . . . . . . . . . . . . . . .

47 49 49

52 53 53 54 55 59 62

Solución Analítica para el Modelo de Doble Porosidad - Una permeabilidad

66

49

2

1 P RINCIPIOS F ÍSICOS La ecuación de difusión describe el comportamiento del flujo de fluidos a través de medios porosos. Esta ecuación, es una de las más importantes en el área de la ingeniería de yacimientos y es la base de la teoría de pruebas de presión. La formulación matemática de esta ecuación es basada en la combinación de la ecuación de continuidad, una ecuación de estado, una ecuación de transporte y la ecuación de energía. Estas ecuaciones se discuten acontinuación:

1.1 E CUACIÓN DE C ONTINUIDAD Esta ecuación es simplemente la conservación de masa en el sistema. Este principio de continuidad establece, considerando un volumen de control, lo siguiente: [ Masa que entra al VC] − [ Masa que sale del VC] = [ Masa que acumula en el VC]

1.2 E CUACIÓN DE E STADO Cualquier ecuación que relaciones la presión, temperatura y volumen, puede ser considera como una ecuación de estado y ser considerada. En el caso de flujo de un fluido ligeramente compresible a través de medios porosos, , por ejemplo, aceite, la compresibilidad a temperatura, T , constante, esta definida como: ·

cf =

1 ∂ρ ρ ∂p

¸

(1.1) T

1.3 E CUACIÓN DE T RANSPORTE La ecuación de transporte relaciona la velocidad del fluido con el potencial del fluido dentro del volumen de control, VC. La lay de Darcy describe la velocidad de un fluido a través de un medio poroso. Esta establece lo siguiente: · ¸ k ∂p ∂D v =− −γ µ ∂r ∂r

(1.2)

La profundidad D es positiva hacia abajo. Cuando el flujo es horizontal, la ecuación de Darcy, se simplifica: · ¸ k ∂p v =− (1.3) µ ∂r

3

1.4 E CUACIÓN DE E NERGÍA Debido a que el flujo de fluidos a través de medios porosos es considerado a temperatura constante, la ecuación de conservación de energía no es considerada en el análisis. El flujo de calor es nulo.

4

2 D ERIVACIÓN DE LA E CUACIÓN DE D IFUSIÓN El proceso para la derivación de las ecuacion de difusión es el siguiente: • Elegir el volumen de control representativo del sistema, VC. • Identificar los flujos másicos que entran y salen del VC en un intervalo de tiempo. • Realizar un balance de masa, i.e. verificar la conservación de masa • Tomar el límite cuando el VC tiende a cero y el intervalo de tiempo tiende a cero. • Incluir la ecuación de Estado • Inlcuir la ecuación de Transporte Se procederá con la derivación. Considere el volumen de control mostrado en la Figura 2.1.Nótese que el flujo tiene la dirección opuesta a la dirección radial.

Dirección radial (+) Dirección de Flujo

Figure 2.1: Volumen de control

5

˜ , que entra al volumen de control es, El flujo másico m

˜ ent r a = −ρu r 2πh(r + ∆r ) m

(2.1)

˜ , que sale del volumen de control es, El flujo másico m

£ ¤ £ ¤ ˜ sal e = −ρu r − ∆(−ρu r ) 2πr h = −ρu r + ∆(ρu r ) 2πr h m

(2.2)

El signo negativo es considera porque el flujo es en dirección opuesta a la dirección r. La masa en el volumen de control es, £ ¤ Vr = 2π (r + ∆r )2 − r 2 h

(2.3)

Vr = 2πr ∆r h

(2.4)

Desarrollando,

Entonces la la masa de fluido en el volumen de control es,

m = 2πr ∆r hφρ

(2.5)

Sustituyendo las ecs. 2.2 a 2.5 en la ecuación de continuidad, se tiene: £ ¤ (2πr ∆r hφρ)t +∆t − (2πr ∆r hφρ)t +∆t −ρu r 2πh(r + ∆r ) − −ρu r + ∆(ρu r ) 2πr h = £ ∆t ¤ £ ¤ 2πr ∆r h (φρ)t +∆t − (φρ)t +∆t −ρu r 2πh(r + ∆r ) − −ρu r + ∆(ρu r ) 2πr h = ∆t

(2.6)

Dividiendo la ec. 2.6 por 2πr ∆r h,

−

£ ¤ ρu r 2πh(r + ∆r ) − −ρu r + ∆(ρu r ) 2πr h

2πr ∆r h

=

£ ¤ (φρ)t +∆t − (φρ)t +∆t

∆t

(2.7)

simplificando £ ¤ (φρ)t +∆t − (φρ)t +∆t ρu r ∆(ρu r ) − − = r ∆r ∆t

(2.8)

La ec. 2.8, también puede ser escrita como, · ¸ 1 ∆(ρu r ) ∆(φρ) ρu r + r =− r ∆r ∆t

(2.9)

6

Tomando límites cuando ∆r y ∆t tienden a 0, · ¸ 1 ∂(ρu r ) ∂(φρ) ρu r + r =− r ∂r ∂t

(2.10)

Si se deriva el siguiente producto de funciones, · ¸ ∂r ∂(r ρu r ) ∂(ρu r ) = ρu r +r ∂r ∂r ∂r

(2.11)

Por lo tanto, considerando la ec. 2.11 en la ec. 2.9, se tiene, 1 ∂(r ρu r ) ∂(φρ) =− r ∂r ∂t

(2.12)

La ec. 2.12 es la ecuación de continuidad en coordenadas radiales y representa la conservación de masa. Sustituyendo la ley de Darcy, el cual establece que la velocidad es proporcional al gradientes de potencial (o presión) con la velocidad. Incluyendo la ec 1.3 para flujo horizontal en la ec. 2.12, se llega a , · µ ¶¸ 1 ∂ k ∂p ∂(φρ) rρ − =− r ∂r µ ∂r ∂t

(2.13)

Considerando la viscosidad, µ y la permeabilidad, k, constantes, eliminando signos negativos en ambos lados y aplicando la regla de la cadena en el lado derecho de la ec. 2.13, se tiene, · µ ¶¸ µ ¶ 1 ∂ ∂p µ ∂(ρ) ∂(φ) rρ = φ +ρ r ∂r ∂r k ∂t ∂t

(2.14)

El lado izquierdo de la ec. 2.14 es expandido y la regla de la cadena es usada en ambos lados, · ¸ µ ¶ ρ ∂ ∂p ∂p ∂ρ ∂p µ 1 ∂ρ 1 ∂φ ∂p r + = φρ + r ∂r ∂r ∂r ∂p ∂r k ρ ∂p φ ∂p ∂t

(2.15)

Definiendo la compresibilidad de la roca c r ,

cr =

1 ∂φ φ ∂p

(2.16)

Sustituyendo la ecuación de estado 1.1 y la ec. 2.16, en la ec. 2.15, · ¸ µ ¶2 ¡ ¢ ∂p ρ ∂ ∂p ∂p µ r +cf ρ = φρ c f + c r r ∂r ∂r ∂r k ∂t

(2.17)

7

Dividiendo entre la densidad, · ¸ µ ¶2 ¢ ∂p 1 ∂ ∂p ∂p µ ¡ r +cf = φ c f + cr r ∂r ∂r ∂r k ∂t

Considerando que la compresibilidad total es c t = c f + c r , y que el termino c f queño considerado con el resto de terminos, se obtiene lo siguiente, · ¸ ∂p φµc t ∂p 1 ∂ r = r ∂r ∂r k ∂t

(2.18) ³

´ ∂p 2 ∂r

es pe-

(2.19)

La ec. 2.19 representa la ecuación de difusión en coordenadas radiales para un fluido con compresibilidad pequeña y constante, donde los cambios del gradiente de presión deben de ser pequeños Como es recomendado en la refrencias técnicas, se tendrá que considerar las premisas para las cuales, esta, fue obtenida. Fundamentalmente es la base de la teoría de k pruebas de presión. la constante φµc , es conocido como término de difusividad hidraúlica. t Con el propósito de validar las soluciones numéicas que se desarrollen, se presentará a continuación la solución numérica para la Ecuación de Difusividad en Coordenadas Radiales.

8

3 S OLUCIÓN A NALÍTICA DE LA E CUACIÓN DE D IFUSIÓN PARA YACIMIENTOS N O F RACTURADOS Una solución analítica es, por mucho, mejor que una solución numérica. Aunque las soluciones analíticas son obtendidas para geometrías regulares y algunas condiciones de frontera, la ec. 2.19 puede ser resuelta analíticamente porque es lineal y el término de difusividad hidraúlica es constante e independiente de la presión. Siendo lineal, esta, puede ser resuelta analíticamente para diferentes condiciones iniciales y de frontera. La Figura 3.1 muestra un diagrama de la geometria del problema a resolver y sus diferentes condiciones de frontera.

Figure 3.1: Condiciones de Fronter Coordenadas Radiales En el área de pruebas de presión, la ec. 2.19 ha sido resuelta para las siguientes condiciones: • Gasto constante • Yacimiento infinito • Yacimiento finito, frontera cerrada

9

• Yacimiento finito, con frontera a presión constante

3.1 S OLUCIÓN F UENTE L INEAL En ausencia de almacenamiento y daño, el transiente de presión ( respuesta de presión que se podría medir en la cara del pozo) en un yacimiento, en flujo radial, presentando comportamiento infinito, y considerando al pozo como una línea fuente (source line, radio del pozo despreciable) bajo la condición de gasto constante, es (Matthews and Russel, 1967, pp. 130) # " r D2 1 PD = − Ei − 2 4t D

(3.1)

donde las variables adimensional es estan definidas por: ¤ 2πkh £ p i − p r,t qB µ kt

PD = tD = rD =

2 φµc t r w

(3.2)

r rw

En la obtención de la 3.1 se han realizado ciertas suposiciones y es recomendable mantenerlas en mente para conocer sus limitaciones. Las suposiciones son las siguientes: • Difusividad constante, esto es, k, φ, c t , µ constantes • Geometría radial • Espesor de la formación constante • Gasto, qB , constante • Yacimiento finito, con frontera a presión constante Sustituyendo las ecs. definidas en 3.2 en 2.18,

p r,t

· ¸ φµc t r 2 qB µ Ei − = pi + 4πkh 4kt

(3.3)

La ec. 3.3 determina la presión en cualquie radio, r, y a cualquier tiempo, t., incluido el radio del pozo, r w . Precaución se deberá tener en el cálculo de la integral exponencial, E i (−x).

10

3.1.1 I NTEGRAL E XPONENCIAL La integral exponencial, E i (−x), definida por ∞

Z

E i (−x) =

x

e −u du u

(3.4)

La ec. 3.4 pueder ser aproximada por la siguiente serie, (Craft and Hawkin, 1959, pp. 311) ∞

Z

−E i (−x) =

x

¸ · e −u x2 x3 x4 xn ∼ d u = ln x + 0.5772 − x + − + − .... + u 2×2 3×3 4×4 n ×n

(3.5)

Considerando, los dos primeros términos de la serie, ∞

Z

−E i (−x) =

x

e −u du ∼ = [ln x + 0.5772] u

(3.6)

La ec. 3.6 puede ser usada para valores de x 10.9, la E i (−x) puede ser considera 0 para la mayoria de los casos prácticos en ingeniería de yacimientos.

11

Table 3.1: Información pi

k 2

h

µ

ct −1

(lb/pg )

(mD)

(Pie)

Cp

(lb/pg)

4000

60

15

1.5

12x10−6

B

φ

(bbl/STB)

(Fracción)

1.251

0.15

s

0

3.1.2 E JEMPLO DE A PLICACIÓN Un pozo de aceite produce un gasto constante de 200 STB/D. Las propiedades de la formación productora se muestran en la tabla 3.1 Calcular la presión para los radios de 0.25, 5, 10, 50, 100, 500, 1000, 1500, 2000, y 2500 pie. Para los tiempos de 1, 24 y 48 hrs. Grafique los resultados de la siguiente forma: • Gráfica normal de radio (horizontal) contra presión (vertical) • Gráfica semilog del radio (horizontal) contra presión (vertical) • Comente los resultados Solución Paso 1. Sustituir los valores mostrados en la tabla 3.1 en la ec. 3.3. Precaución se deberá tener en el manejo de los signos. La 3.3 presenta un signo negativo en el argumento en la E i . Este signo es útil cuando se esta calculando la E i (−x) graficamente, como se muestra en la Figura 3.2 En nuestro caso, donde se esta utilizando alguna aproximación, el argumento de la integral exponencial debe introducirse con signo positivo, E i (x) y el resultado será considerado negativo. De esta manera, la ec. 3.3 es reescrita como:

p r,t = p i −

¸ · qB µ φµc t r 2 Ei 4πkh 4kt

(3.10)

Paso 2. Se realizan los cálculos con ayuda de un código escrito en MatLab. Se recomienda usar el SI para evitar el manejo de constantes de unidades o cualquier otro sistema de unidades consistente. En este problema se una el SI, aunque los resultados sepresentan en las unides que se proporcionan en los datos de entrada. Paso 3. La Tabla 3.2 muestra los resultados para los tiempos y radios requeridos. Paso 4. La Figura 3.3 muestra los resultados en forma gráfica (normal y semilog) para los tiempos y radios requeridos.

12

Table 3.2: Resultados ejemplo de aplicación, tiempos 1, 12, y 24 hr

1 2 3 4

Tiempo

Radio

Argumento x

Ei(x)

p(r,t)

hr

pie

1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

0.25 5.00 10.00 50.00 100.00 500.00 1000.00 1500.00 2000.00 2500.00

2.667e-006 1.067e-003 4.266e-003 1.067e-001 4.266e-001 1.067e+001 4.266e+001 9.599e+001 1.707e+002 2.667e+002

12.26 6.27 4.88 1.76 0.66 0.00 0.00 0.00 0.00 0.00

3459.13 3723.46 3784.49 3922.13 3970.89 4000.00 4000.00 4000.00 4000.00 4000.00

12.0 12.0 12.0 12.0 12.0 12.0 12.0 12.0 12.0 12.0

0.25 5.00 10.00 50.00 100.00 500.00 1000.00 1500.00 2000.00 2500.00

2.222e-007 8.888e-005 3.555e-004 8.888e-003 3.555e-002 8.888e-001 3.555e+000 8.000e+000 1.422e+001 2.222e+001

14.74 8.75 7.37 4.15 2.79 0.27 0.01 0.00 0.00 0.00

3349.49 3613.86 3675.02 3816.67 3876.68 3988.29 3999.71 4000.00 4000.00 4000.00

24.0 24.0 24.0 24.0 24.0 24.0 24.0 24.0 24.0 24.0

0.25 5.00 10.00 50.00 100.00 500.00 1000.00 1500.00 2000.00 2500.00

1.111e-007 4.444e-005 1.778e-004 4.444e-003 1.778e-002 4.444e-001 1.778e+000 4.000e+000 7.111e+000 1.111e+001

15.44 9.44 8.06 4.84 3.47 0.63 0.07 0.00 0.00 0.00

3318.90 3583.27 3644.44 3786.28 3846.87 3972.05 3997.05 3999.83 4000.00 4000.00

psi

radio=[0.25 5,10,50,100,500,1000,1500,2000,2500]; %hr radio=radio*u.ft; %m time=[1 12 24]; %hrs time=time*u.hr; %seg

5 6

Nt=length(time);

13

Figure 3.2: Integral Exponencial

7

Nr=length(radio);

14

8 9 10 11 12 13 14 15 16 17

for i=1:Nt for j=1:Nr term1=Qo*Bo*vis/(4*pi*perm*h); x=por*vis*ct*radio(j)*radio(j)/(4*perm*time(i)); %llamado a funcion de calculo de Ei [Ei]=−IntExp(x); p(i,j)=pini+term1*(Ei); end end

La Figura 3.3 muestra el comportamiento de la presión con respecto al radio. Se observa en la grafica normal que las caídas de presión son mayores a medidad que se acerca al pozo. La caída de presión, con respecto a la presión incial, es mayor a medidad que el tiempo avanza. Se observa, aunque hay pocos puntos, que el transiente de presión aun no alcanza las frontera. Aqui surge una pregunta, si hacemos el cálculo de la presión a tiempos grandes, se sentirán los efectos de frontera? En la gráfica semolog se observa que el comportamiento de la presión con respecto al radio muestra un comportamiento lineal. En otras palabras, la caída de presión es proporcional al log de r.

3.2 C ÁLCULO DE PRESIÓN EN LA CARA DE LA FORMACIÓN , p w f Al evaluar la ec. 3.3 en la cara del pozo, esto es, en r = r w , la ec. 3.6 puede ser usada, en la cara del pozo,

r rw = =1 rw rw rD 1 x= = 4t D 4t D

rD =

(3.11)

Entonces, la ec. 3.6 se convierta a, 1 + 0.5572 −E i (−x) ∼ = [ln x + 0.5772] = ln 4t D

(3.12)

Arreglando el termino en el (ln)

1 + 0.5571 ∼ −E i (−x) ∼ = ln t D − 0.8091 = ln 4t D

(3.13)

Sustituyendo la ec. 3.13 en la ec. 3.1,

PD =

1 [ln t D − 0.8091] 2

(3.14)

Nótese la eliminación del signo negativo.

15

Tiempo (hr): 1, 12 y 24 4000 3950

Presión (Psia)

3900 3850 3800 3750 3700 3650 3600

100

200

300

400

500 600 Radio (pie)

700

800

900

1000

Tiempo (hr):1, 12, y 24 4000 3900 3800

Presión (Psia)

3700 3600 3500 3400 3300 3200 3100 3000 -1 10

0

10

1

2

10

10

3

10

4

10

Radio (pie)

Figure 3.3: Presión vs radio de drene

16

3.2.1 FACTOR DE D AÑO En el caso donde la formación se encuentre alterada en su permeabilidad, ya sea por por fluidos de perforación ( disminuación de la permeabilidad ) o por una estimulación (mejoramiento de la permeabilidad), la ec. 3.1 no puede ser aplicada. Debido a que se supuso una permeabilida de la formación constante en todo el dominio. Por lo tanto, para que la ec. 3.1 se usada en dichas formaciones es neceario considerar un término extra donde se considera la alteración de la formación.

Figure 3.4: Zona alterada en la vecindad del pozo (Hawkins, 1959, pp.319) menciona que si la zona alterada es considerada como una region de radio r s con diferente permeability,k s , esto es, Figura 3.4, entonces la caída de presión se puede modelar usando la ecuación de flujo en régimen estacionario, (Dake, 1998, pp.139)

q=

2πkh B µ ln( rrwe )

∆p

(3.15)

Arreglando para ∆p,

17

∆p =

re qB µ ln( ) 2πkh rw

(3.16)

Por lo tanto, en relación a la Figura 3.4, la caída de presión adicional resultante de una formación alterada (cerca del pozo) con permeabilidad, k s , es · ¸ · ¸ qB µ qB µ rs rs ∆p s = − ln ln 2πkh rw 2πk s h rw

(3.17) · ¸ · ¸ rs qB µ k − 1 ln = 2πkh k s rw

Al no exister zona alterada, k s = k, entonces la caída de presión adicional, ∆p s = 0. La caída de presión adimensional debida al daño es definida como,

∆p sD =

2πkh∆p s qB µ

(3.18)

Escribiendo la ec. 3.17 en forma adimensional, ¸ · ¸ rs k − 1 ln ∆p sD = ks rw ·

(3.19)

3.2.2 C AÍDA DE P RESIÓN T OTAL EN LA CARA DEL POZO Para determinar la caída de presión total se combinan la ec. 3.19 y la ec. 3.3. En forma adimensional, la ecuación para determinar las caídas de presión a cualquier radio y tiempo es, " # r D2 1 + ∆p sD (3.20) PD = − Ei − 2 4t D Se deberá recordar que la alteración de la formación (zona dañada) sucede en las vecindades del pozo, cerca de r w . Si la ec. 3.20 es usada para calcular las caídas de presión en el interior del yacimiento, a un r >> r w , el segundo término de la ecuación es nulo. Debido a que las mediciones de presión se realizan en la cara del pozo, esto es, en r = r w . Es importante obtener una expresión analítica para su cálculo. En radios de drene pequeños, incluidos el radio del pozo, r w , la aproximación logarítmica de la integral exponencial es válida, esto es, al sustituir la ec.3.14, obtenida para r = r w , en la ec. 3.20, se obtiene lo siguiente,

18

PD =

1 [ln t D + 0.80907] + ∆p sD 2

(3.21)

Sustituyendo, también, el término que representa la caída de presión por el daño a la formación, ∆p sD , esto es, la ec. 3.19, se tiene,

PD =

· ¸ · ¸ 1 k rs − 1 ln [ln t D + 0.80907] + 2 ks rw

Como es complicado conocer el radio de daño, r s y su permeabilidad, k s , el término, es considerado el factor de daño, s. De esta manera, la ec. 3.22, se convierte en,

PD =

1 [ln t D − 0.8091] + s 2

(3.22) h

k ks

i h i − 1 ln rrws ,

(3.23)

Así, el factor de daño, s, es considerado como una caída de presión (adimensional) extra. Introduciendo el factor de daño en el interior del corchete y sustituyendo la definición de cada variable adimensional, esto es, ecs.3.2 en la ec. 3.23, se llega a, · ¸ ¤ 1 2πkh £ kt ln pi − p w f = − 0.8091 + 2s 2 qB µ 2 φµc t r w ¸ · £ ¤ qB µ 1 kt − 0.8091 + 2s ln pi − p w f = 2 2πkh 2 φµc t r w

(3.24)

Nótese que se usa p w f , debido a que estamos en r = r w . Aislando el ln t y considerando la definición, ln x = 2.3025 log x. La ec. 3.24 se convierte en, · ¸ k qB µ 1 −0.8091 + 2s pi − p w f = 2.3025 log t + log + 2 2πkh 2 2.3025 φµc t r w · ¸ £ ¤ 2.3025qB µ 1 k −0.8091 + 2s pi − p w f = log t + log + 2 2πkh 2 2.3025 φµc t r w · ¸ £ ¤ 1.151qB µ 1 k pi − p w f = log t + log − 0.3514 + 0.8686s 2 πkh 2 φµc t r w £

¤

(3.25)

La ec. 3.25 esta escrita en unidades consistente, por ejemplo, para el Sisteme Internacional de unidades, las variables tiene las siguientes unidades, p w f , Presión de fondo fluyendo, (Pa) p i , Presión inicial, (Pa) k , Permeabilidad, (m2 ) h, Espesor,(m) µ, Viscosidad,(Pa seg)

19

c t , Compresibilidad,(Pa−1 q, Gasto, (m3 /seg) B , Factor de Volumen, (m3 /Sm3 ) φ, Porosidad,(Fracción) t , Tiempo, (Seg) s, Daño a la formación, (adim) La ec. 3.25 en unidades del sistema ingles es,

pw f

· ¸ 162.6qB µ k = pi − log t + log − 3.2274 + 0.8686s 2 kh φµc t r w

(3.26)

Donde: p w f , Presión de fondo fluyendo, (lb/pg2 ) p i , Presión inicial, (lb/pg2 ) k , Permeabilidad, (mD) h, Espesor,(Pie) µ, Viscosidad,Cp c t , Compresibilidad,(lb/pg)−1 q, Gasto, STB/Día B , Factor de Volumen, (bbl/STB) φ, Porosidad,(Fracción) t , Tiempo, Días s, Daño a la formación, (adim)

3.2.3 E JEMPLO DE A PLICACIÓN

20

4 S OLUCIÓN N UMÉRICA DE LA E CUACIÓN DE D IFUSIÓN PARA YACIMIENTOS NO F RACTURADOS En esta sección se aproximará, mediante diferencias finitas, la solución de la ecuación de difusión para formaciones no fracturadas, es decir, yacimientos homogeneos. Posteriormente se planteará el caso de yacimientos naturalmente fracturados. Se ha discutido ampliamente que al considerar las propiedades del fluido dependientes de la presión, le ecuación diferencial es no lineal y, en general, para su solución, se debe de recurrir a métodos numéricos. El caracter continuo de la solución de la Ecuación en Derivadas Parciales es cambiado a un caracter discreto tanto en espacio como en tiempo. Reescribiendo la ec. 2.13, en lugar de la ec. 2.19, e introduciendo el factor de volumen de aceite, B , se tiene lo siguiente, 1 ∂ r ∂r

·

¸ µ ¶ k ∂p ∂ φ r = µB ∂r ∂t B

(4.1)

Nótese que la ec. 2.13 ha sido dividida por la densidad @ c.e., donde se ha aplicado la definición del fector de volumen y se han eliminado ambos signos negativos. ρ 1 = ρ c.e B

(4.2)

La ec. 4.1, junto con sus condiciones de frontera, describe el flujo monofásico de un fluido a través de un medio poroso, fluyendo en dirección radial y horizontal. Esta ecuación tiene aplicación en analizar el comportamiento del flujo radial de un solo pozo, y el diseôsar pruebas de decremento e incremento de presión. La formulación propuesta para solucionar numéricamente el flujo de fluidos en medios porosos considera almacenamiento en el pozo y daño a la formación en la vecindad del pozo. La Figura 4.1 muestra un diagrama del sistema superficie-pozo-yacimiento. En ella, se establecer una convención para enumerar las celdas numéricas, las fronteras e identificar las contribuciones a el gasto total que se tiene en la superficie. Es importante establecer que la celda 1 esta asignada al pozo (espacio libre). La asignación de los nodos al medio poroso será, entonces, del nodo 2 al nodo Imax. Para plantear el problema superficie-pozo-yacimiento, primero se iniciará, en términos de gastos y presiones, el análisis, para la interacción entre la superficie, el pozo y el primer celda numérica del yacimiento.

21

qsup

qsup = qsf + qws

qws pozo

Nodos

1 2

3

4

Imax

qsf 1

2

3

4

Imax

Fronteras

Figure 4.1: Discretización del sistema superficie-pozo-yacimiento

4.1 C ONTRIBUCIONES AL G ASTO T OTAL EN S UPERFICIE Se observa que el gasto total en la superficie,q sup , es la contribución del gasto de la formación, q s f , y el gasto que aportará, por un tiempo corto, la expansión de los fluidos almacenados en el pozo, q w s . El gasto total expresado a condiciones de yacimiento,

q sup B = q w s + q s f

(4.3)

4.1.1 C ONTRIBUCIÓN AL G ASTO POR ALMACENAMIENTO EN EL P OZO, q w s El almacenamiento del pozo puede ser causado, principalmente, por dos efectos. Una es almacenamiento por expansión del fluido almavenado en el pozo, y la otra por cambio en el nivel del fluido dentro del pozo. Expansión del Fluido en el Pozo Al abrir un pozo, empacado, a producción, Figura 4.2 (a) , el pozo sufre una caída de presión, esta viaja a través de la columna del fluido que contiene el pozo. Antes de que el disturbio de

22

presión alcance la cara del pozo, la caída de presión causa una expansión en los fluidos del pozo. El tiempo en que el yacimientos siente el disturbio de presión es relativamente corto. La producción total en superficie, en el tiempo que dura el efecto de almacenamiento, es, por lo tanto, la contribución de la expansión del fluido, primero, y la producción del yacimiento después.

qsup = qws + qsf

qsup = qws + qsf

qws qws

qsf

(a)

qsf

(b)

Figure 4.2: Efecto de almacenamiento. (a) Expansión de Fluido. (b) Cambio de Nivel. A medida que el fluido se expande, el volumen del pozo es drenado y la formación, como respuesta, inicia a aportar progresivamente. Este es el almacenamiento por expansión del fluido. Puesto en otras palabras, el gasto total que se mide en la superficie, q sup es la suma del gasto proveniente del pozo (volumen almacenado) q w s y el gasto de la formación q s f (wb significa wellbore y s f sandface). La capacidad del pozo para almacenar o descargar fluidos es conocido como almacenamiento del pozo. Cambio en el Nivel del Pozo El segundo tipo de efecto de almacenamiento, más común, es debido al cambio de nivel en

23

la columna de líquido. Esto sucedo en pozos donde no se ha colocado un empacador como lo muestra la Figura 4.2 (b). Cuando un pozo es abierto, la reducción en presión causa una caída del nivel de líquido en el espacio anular. El líquido del espacio anular se une con el gasto proveniente de la formación. La suma representa el gasto total que se mide en la superficie. La caída del nivel del fluido generalmente abastece mucho mas fluido que la formación creando un efecto de almacenamiento mucho mayor que el almacenamiento por la expansión del fluido. El coeficiente de almacenamiento C, es un parÃametro ˛ usado para cuantificar el efecto de almacenamiento de un pozo. Como se establecio anteriormente, el efecto de almacenamiento se presenta a tiempos cortos y tiene su relevancia en pruebas de presión. El coeficiente de almacenamiento, C, se define como el volumen del fluido que el pozo aportará debido a una caÃŋda de presiøsn, esto es,

C=

∆V Vi ni − V (t ) = ∆p p i ni − p w f (t )

(4.4)

Asumiendo un coeficiente constante de almacenamiento, C, y derivando con respecto al tiempo, d [p i ni − p w f (t )] d [Vi ni − V (t )] =C dt dt d p w f (t ) dV (t ) − = −C dt dt

(4.5)

La ec. 4.5 puede ser representada de esta manera,

qw s = C

d p w f (t ) dt

(4.6)

Nótese que en la ec. 4.6, la contribución al gasto total debido al efecto de almacenamiento, q w s , decrece a medidad que la p w f se estabiliza. Asimismo, la ec. 4.6 resultará en un signo negativo, debido a que las p w f , disminuye con respecto al tiempo a medidad que el pozo produce un gasto q. Este efectto se visualizará al momento de aproximar la ecuación en diferencias finitas. 4.1.2 C ONTRIBUCIÓN AL G ASTO POR A PORTE DEL YACIMIENTO, q s f Hagamos usos de la Figura 4.3 para representar la contribución del gasto de la formación, q s f , al gasto total q sup . Para este anaálisis, se raelizarán dos suposiciones: • El regimen que se presenta del nodo 2 (yacimiento) al nodo 1 (pozo) es un regimen pseudoestacionario (conocido, también, como psudopermanante, semi estacionario; Psudo Steady State, PSS, en ingles).

24

1

2

Fronteras

qsf

Figure 4.3: Convención usada para malla radial • La celda 2 (yacimiento) es la que eventualmente presentará daño a la formación. Consideremos la solución de la Ecuación de difusión para regimen pseudoestacionario, considerando la caída extra de presión causada por el factor de daño,

pe − p w f

· ¸ qµ re 1 =− ln − +s 2πkh rw 2

(4.7)

En la ec. 4.7 se ha considerado, por convención, que el gasto que se produce es negativo (-) y el gasto que se inyecta es positivo (+). Considernado el análisis entre el nodo 1 (pozo) y nodo 2 (primera celda del yacimiento), se tiene los siguiente, pe ⇒ p2 p w f ⇒ p1 r e ⇒ r 2+ 1 2

(4.8)

r w ⇒ r1 q ⇒ qs f

25

Por lo tanto, la ecs. 4.6 se convierte en,

qw s = C

d p w f (t ) dt

⇒C

d p 1 (t ) dt

(4.9)

y la ec. 4.7, se convierte en, · r 1 ¸ 2+ 2 1 ln − +s p2 − p1 = − 2πkh r1 2

qs f µ

(4.10)

En términos del gasto de la formación, q s f ,

qs f = −

2πkh p2 − p1 h r 1 i 2+ µ ln 2 − 1 + s r1

(4.11)

2

Con el propósito de facilitar la integración de la ec. 4.11 con el yacimiento, se considerará la siguiente aproximación, µ

ln

r 2+ 1 2

r1

¶ µ ¶ µ ¶ µ ¶ r2 − r1 ∆r 1 ∼ r2 − r w − = = = 2 rw r1 r 1+ 21

(4.12)

Entonces, haciendo la consideración de la ec. 4.12, la ec. 4.11, se reescribe como,

qs f

2πkh p 2n+1 − p 1n+1 i h =− B µ ¡ ∆r ¢ 1 + s r

(4.13)

1+ 2

Sustituyendo la ecs. 4.6 y 4.11 en la ec. 4.3, se tiene que el gasto total en superficie a condiciones de yacimiento esta definido como,

qs f B = C

d p 1 (t ) 2πkh p2 − p1 h¡ ¢ i − ∆r dt µ 1 +s r

(4.14)

1+ 2

Expresada a condiciones de superficie,

qs f =

C d p 1 (t ) 2πkh p2 − p1 h¡ ¢ i − ∆r B dt Bµ 1 +s r

(4.15)

1+ 2

Aproximando la derivada en el primer término de la ec. 4.14 mediante diferencias finitas y considerando que, en el segundo término, la p 2 y p 1 serán evaluadas al tiempo n + 1, se

26

obtiene la ecuación que acopla el gasto en superficie con el almacenamiento, el factor de daño y la primer celda del yacimiento:

q sup =

C p 1n+1 − p 1n 2πkh p 2n+1 − p 1n+1 i h − B ∆t B µ ¡ ∆r ¢ 1 + s r

(4.16)

1+ 2

Nótese que en la ec. 4.16, la cual representa la ecuación para el nodo 1 y su interacción con el nodo 2, únicamente hay variables (incógnitas), como se esperaba, del nodo 1 y nodo 2. Regresando a la Figura 4.1, hasta ahora se tiene una expresión que relaciona el gasto en la superficie, el gasto debido al almacenamiento y el gasto de la propia formación. Esta relación, como se establecio anteriormente, será designada al nodo 1 (pozo). A continuación se obtendrán las ecuaciones del resto de los nodos, los cuales pertenecen al yacimiento (medio poroso) y posteriormente se acoplaran para resolverse simultaneamente.

27

4.2 D ISCRETIZACIÓN DE LA E CUACIÓN DE D IFUSIÓN Para migrara del medio continuo al medio discreto, la ec. 4.1 se aproxima mediante diferencias finitas. La Figura 4.4 muestra las caracteríticas del dominio físico a considerar para realizar la discretización. Nótese que la celda numérica que representa el pozo en la celda 1 y la distribución logarítmica de los nodos. Las celdas son mas refinadas en las cercanias del pozo.

Imax

4

2 1 22 3

3

3

Fronteras Imax

4

2

4

Imax

3

1 2

4

Imax

Figure 4.4: Discretización Dominio Radial A continuación se presenta un par de algoritmos para determinar la posición de los nodos, las fronteras y el volumen de cada celda numérica en coordenadas radiales.

4.3 M ALLA R ADIAL La convención que se utilizará para generar la malla radial es similar a la utilizada en la generación de la malla cartesiana. La Figura 4.5 muestra la simbología a utilizar. Supongase que los radios de los nodos son relativos a un origen, esto es r = 0.

28

∆𝒓

∆𝒓

𝟏 𝒊−𝟐

𝟏 𝒊+𝟐

i

i-1

∆𝒓𝒊−𝟏

i+1

∆𝒓𝒊

𝒓

𝒓

𝟏 𝒊−𝟐

𝒊+

𝟏 𝟐

i

i-1 𝒓𝒊−𝟏

∆𝒓𝒊+𝟏

i+1 𝒓𝒊+𝟏

𝒓𝒊

Figure 4.5: Convención usada para malla radial

4.3.1 B LOQUES D ISTRIBUIDOS Pendiente... 4.3.2 B LOQUES C ENTRADOS En la distribución de nodos centrados para malla radial, los nodos son espaciados logarítmicamente a partir del pozo y hacia la frontera externa. En el apéndice A se demuestra el porqué de este espaciamiento.

re ψ= rw ·

¸

1 I max

r1 = r w (4.17)

ψ ln ψ r2 = r1 (ψ − 1) r i = ψr i −1

para i = 3,4,...Imax

29

Las fronteras de las celdas son definidas mediante un promedio logarítmico,

ri + 1 = 2

r i +1 − r i ³ ´ ln r ri +1 i

para i = 2,3,...Imax

(4.18)

El volumen de la roca esta dado por

¡ ¢ Vi = π r i2+1 − r i2 ∆z

para i = 2,3,...Imax

(4.19)

4.3.3 E JERCICIO : G ENERACIÓN DE M ALLAS R ADIALES Colocar Tarea 18

30

4.4 A PROXIMACIÓN DE LA E CUACIÓN DE D IFUSIÓN MEDIANTE D IFERENCIAS F INITAS La ec. 4.1 , 1 ∂ r ∂r

·

¸ µ ¶ k ∂p ∂ φ r = µB ∂r ∂t B

y sus condiciones iniciales y de frontera: p(r, t = 0) = p 0 ¸ ∂p qB µ r =− ∂r r w 2πkh · ¸ qB µ ∂p =− =0 r ∂r r e 2πkh ·

(4.20)

son discretizada mediante el método de diferencias finitas. Llamaremos al término de lado izquierdo de la ec. 4.1 , Término de flujo y al término del lado derecho, Término de acumulación. Aproximando la ec. 4.1 en el nodo i al nivel de tiempo n + 1, se tiene lo siguiente, 1 ∂ r ∂r

·

k ∂p r µB ∂r

¸n+1 i

µ ¶ ∂ φ n+1 = ∂t B i

para i = 2,3,...Imax

(4.21)

n =1,2,... con sus respectivas condiciones iniciales y de frontera, p(r, t = 0) = p 0 ∂p r ∂r ·

·

r

∂p ∂r

¸n+1 1+ 21

=−

¸n+1 I max+ 12

=−

qB µ 2πkh

(4.22)

qB µ =0 2πkh

Convencionalmente, se considera con signo negativo a un gasto que se produce y positivo a un gasto que se inyecta 4.4.1 D ISCRETIZACÍON DEL T ÉRMINO DE F LUJO Realizando la siguiente sustitución: k ∂p r U= µB ∂r ·

¸

(4.23)

31

Con la ayuda de la Figura 4.5, aplicando diferencias finitas centrales, la discretización de la ec. 4.23 en el nodo i al nivel del tiempo n + 1 es, · ¸ U n+1 −U n+1 i − 12 1 ∂U n+1 ∼ 1 i + 12 = r ∂r i ri ∆r i

(4.24)

En la ec. 4.24 , se esta considerando un punto adelante del nodo i , esto es,(i + 21 ) y un punto atras del nodoi , esto es,(i − 12 ) . y la distancia que los separa, ∆r i . Sustituyendo U por su definición original, esto es, ec. 4.23, se tiene,

1 ri

U n+1 −U n+1 i + 12 i − 12 ∆r i

h

=

1 ri

in+1 h in+1 k ∂p k ∂p r − r µB ∂r i + 1 µB ∂r i − 1 2 2

(4.25)

∆r i

Factorizando r i + 1 , r i − 1 y ∆r i , 2

2

1 ri

n+1 U n+1 1 −U 1 i+ 2

i− 2

∆r i

2

Separando movilidad, esto es,

1 ri

n+1 U n+1 1 −U 1 i− 2

i+ 2

∆r i

r i + 1 µ k ∂p ¶n+1 r i − 1 µ k ∂p ¶n+1 2 2 − = r i ∆r i µB ∂r i + 1 r i ∆r i µB ∂r i − 1

k µB

y la

(4.26)

2

∂p ∂r , se tiene

r i + 1 µ k ¶n+1 µ ∂p ¶n+1 r i − 1 µ k ¶n+1 µ ∂p ¶n+1 2 2 = − r i ∆r i µB i + 1 ∂r i + 12 r i ∆r i µB i − 1 ∂r i − 12 2

2

∂p 1 ∂r , evaluada en los puntos i + 2

En la ec. 4.27, aparece se obtiene lo siguiente,

µ

µ

∂p ∂r

¶n+1

∂p ∂r

¶n+1

i + 21

i − 21

∼ =

(4.27)

e i − 12 . Aplicando, derivadas centrales,

p in+1 − p in+1 +1 ∆r i + 1

2

∼ =

(4.28)

p in+1 − p in+1 −1 ∆r i − 1

2

Sustituyendo la ec.4.28 en la ec. 4.27, la aproximación mediante diferenicas finitas centrales, del lado izquierdo de la ec. 4.1 es, 1 ∂ r ∂r

·

k ∂p r µB ∂r

¸n+1 i

∼ =

à ! à ! r i + 1 µ k ¶n+1 p n+1 − p n+1 r i − 1 µ k ¶n+1 p n+1 − p n+1 i +1 i i i −1 2 2 − r i ∆r i µB i + 1 ∆r i + 1 r i ∆r i µB i − 1 ∆r i − 1 2

2

2

2

(4.29) para i = 2,3,...Imax; n =1,2,...

32

4.4.2 D ISCRETIZACÍON T ÉRMINO DE A CUMULACIÓN El término del lado derecho de la ec. 4.21, es el Término de acumulación. Aplicando diferencias regresivas en tiempo, µ ¶ n+1 n ∂ φ n+1 ∼ (φb)i − (φb)i = ∂t B i ∆t

(4.30)

donde, b = (1/B ) También, el término de acumulación se puede expander como, µ ¶ µ ¶ ∂(φb) 1 ∂b 1 ∂φ ∂p ∂ φ = = φb + ∂t B ∂t b ∂p φ ∂p ∂t

(4.31)

Sustituyendo las definiciones de compresibilidad del fluido y de la roca, se tiene, µ ¶ ¡ ¢ ∂p ∂ φ ∂(φb) = φb c f + c r = ∂t B ∂t ∂t µ ¶ ∂(φb) ∂ φ ∂p = = φbc t ∂t B ∂t ∂t

(4.32)

Aproximando mediante diferencias regresivas la ec. 4.32, (p in+1 − p in ) ∂p ∼ n+1 φbc t = (φbc t ) ∂t ∆t

(4.33)

Las ecuaciones 4.30 y 4.33 son aproximaciones equivalente. La primera es usada en la formulación totalmente implícita y la segunda para formulaciones solo implicita en presión. Sustituyendo las aproximaciones en diferencias finitas de los Términos de Flujo, ec. 4.29, y de acumulación, ec.4.30, en la ecuación original 4.21, se tiene lo siguiente,

à ! à ! r i + 1 µ k ¶n+1 p n+1 − p n+1 r i − 1 µ k ¶n+1 p n+1 − p n+1 (φb)in+1 − (φb)ni i +1 i i i −1 2 2 − = r i ∆r i µB i + 1 ∆r i + 1 r i ∆r i µB i − 1 ∆r i − 1 ∆t 2

2

2

(4.34)

2

para i = 2,3,...Imax; n =1,2,... La ec. 4.21 será usada en la formulación Totalmente Implícita, TI. Sustituyendo las aproximaciones en diferencias finitas de los Términos de Flujo, ec. 4.29, y de acumulación, ec.4.33, en la ec. 4.21, se tiene lo siguiente,

33

à ! à ! r i + 1 µ k ¶n+1 p n+1 − p n+1 r i − 1 µ k ¶n+1 p n+1 − p n+1 (p n+1 − p in ) i +1 i i i −1 2 2 − = (φbc t )n+1 i r i ∆r i µB i + 1 ∆r i + 1 r i ∆r i µB i − 1 ∆r i − 1 ∆t 2

2

2

2

(4.35) para i = 2,3,...Imax; n =1,2,... La ec. 4.35 será la elguida para la formulación implícita en presión, la cual se desarrollará a continuación. Multiplicando la ec. 4.35 por el volumen de roca de la celda, en general, para la cela i , se tiene que, Vr,i = r ∆r ∆θ∆z, para uno ángulo de 360o , es decir, 2π, en radianes, el volumen es,Vr,i = 2πr ∆r ∆z, por lo tanto, se tiene, ! ! Ã Ã (2πr ∆r ∆z)i r i + 1 µ k ¶n+1 p n+1 − p n+1 (2πr ∆r ∆z)i r i − 1 µ k ¶n+1 p n+1 − p n+1 i +1 i i i −1 2 2 − = r i ∆r i µB i + 1 ∆r i + 1 r i ∆r i µB i − 1 ∆r i − 1 2

2

2

2

(p n+1 − p in ) (2πr ∆r ∆z)i (φbc t )in+1 i ∆t

(4.36)

La Figura 4.6 muestra, en general, el volumen de cada celda considerada. La ec. 4.36 puede ser simplificada,

µ

2πr ∆z ∆r

¶ i + 21

µ

k µB

¶n+1 i + 12

¡

µ ¶ µ ¶ (p in+1 − p in ) ¢ ¢ 2πr ∆z k n+1 ¡ n+1 n+1 n+1 n+1 p in+1 − p − p − p = (V bc ) p t i +1 i i i −1 ∆r ∆t i − 21 µB i − 1 2

(4.37) Definiendo el término de Transmisibilidad, µ

Ti ± 1 = 2

A ∆r

¶ i ± 12

(λ)i ± 1

2

(4.38)

Donde el Área y la movilidad estan definidas como, A i ± 1 = 2πr i ± 1 ∆z 2 µ ¶2 k λi ± 1 = 2 µB i ± 1

(4.39)

2

Considerando la ec. 4.38 en la ec. 4.37 ,

34

D𝒓

𝒓 − D𝒓/𝟐

𝒓

𝒓 + D𝒓/𝟐

Dz

𝑽𝒓 = 𝝅 (𝒓 + D𝒓/𝟐)𝟐 − (𝒓 − D𝒓/𝟐)𝟐 D𝒛

𝑽𝒓 = 𝟐𝝅𝒓 D𝒓D𝒛

Figure 4.6: Volumen de roca para nodo i

(p in+1 − p in ) ¡ n+1 ¢ ¡ n+1 ¢ n+1 n+1 n+1 n+1 Tin+1 p − p − T p − p = (V bc ) p t i i +1 i i i −1 + 12 i − 12 ∆t

(4.40)

para i = 2,3,...Imax; n =1,2,... Tratando el termino de transmisibilidad y el término, (Vp bc t ) al nivel de tiempo n, es decir, evaluados explícitamente a la p n , se tiene lo siguiente, (p in+1 − p in ) ¡ ¢ ¡ n+1 ¢ n+1 n n+1 n Tin+ 1 p in+1 − p p − p − T = (V bc ) 1 p t +1 i i i −1 i i− 2 ∆t 2

(4.41)

para i = 2,3,...Imax; n =1,2,... De esta manera estamos linealizando, de forma pedestre, la ec.4.40. Esta ecuación es la que se estará utilizando para construir el modelo numérico en coordenas radiales. Aunque podemos utilizar la formulacióm Totalmente implícita, esto es, ec. 4.34, eso se realizará como

35

proyecto final. Incluyendo la ecuación para el nodo 1, esto es, ec. 4.16, y desarrollando la ec. 4.41 para cada nodo i en el medio poros, esto es, para 2,3,...I , se tiene lo siguiente,

n+1 n C p 1 −p 1 B ∆t

i=1

p n+1 −p 1n+1 i +s r 1+ 1

h 2 − 2πkh B µ ¡ ∆r ¢

=

q sup

2

i=2

Tn

¡

i=3

Tn

¡

2+ 21

¢ p 3n+1 − p 2n+1 − T n

¡

p 2n+1 − p 1n+1

¢

=

(Vp bc t )n2

(p 2n+1 −p 2n ) ∆t

¢ p 4n+1 − p 3n+1 − T n

¡

p 3n+1 − p 2n+1

¢

=

(Vp bc t )n3

(p 3n+1 −p 3n ) ∆t

2− 21

3+ 21

3− 21

...

...

=

...

...

...

=

...

...

...

=

...

i = Imax T n

¡

I max+ 12

¢ n+1 n p n+1 I max+1 − p I max − T

¡

I max− 12

n+1 p n+1 I max − p I max−1

¢

= (Vp bc t )nImax

n (p n+1 I max −p I max ) ∆t

(4.42)

Este sistema de ecuaciones se deben resolverse simultáneamente. Nótese que las presiones al nivel de tiempo n + 1 son las icógnitas a resolver, esto es, p in+1 , para i = 1, 2, 3...I max. Los términos evaluados al tiempo n son conocidos. 4.4.3 A COPLAMIENTO DE C ONDICIONES DE F RONTERA El acoplamiento de las condiicones de frontera, ecs. 4.22, se realizará en las fronteras donde se espcifica la condición. La Figura 4.7 muestra la celda donde se acoplaran ´ dichas condiciones. Condición de Frontera Externa. La condición de frontera exterior en r = r e fue identificada en la ec. 4.22c. Aproximando la derivada mediante diferencias finitas centrales: ∂p r ∂r ·

¸n+1 I max+ 21

· ¸ p I max+1 − p I max n+1 ∼ = r I max+ 1 2 r I max+1 − r I max

(4.43)

Por lo tanto, sustituyendo la ec. 4.43 en la ec. 4.22c, se llega a, ·

r I max+ 1 2

p I max+1 − p I max r I max+1 − r I max

¸n+1

qB µ =− 2πkh ·

¸ I max+ 21

=0

(4.44)

36

Imax Nodos

1

2

𝐼𝑚𝑎𝑥 +

1 2

2 1/2 Imax

1 1/2 2

Fronteras

1

Figure 4.7: Condiciones de Frontera

37

Resolviendo para el gasto q, se tiene lo siguiente, Ã

q =−

2πr I max+ 1 h 2

!µ

k µB

¶

¡

p I max+1 − p I max

¢n+1

=0 r I max+1 − r I max I max+ 21 µ ¶ µ ¶ ¡ ¢n+1 A k q =− p I max+1 − p I max =0 ∆r I max+ 12 µB I max+ 1

(4.45)

2

Recordando la definición de transmisibilidad, ecs. 4.39; la ec. 4.45, se puede expresar como, ¡ ¢ n+1 q = −T Inmax+ 1 p n+1 (4.46) I max+1 − p I max = 0 2

Nótese que en la ec. 4.46, se esta evaluando la Transmisibilidad al nivel de tiempo n, es decir, explícitamente. La ec. 4.46 debe ser acoplada al sistema de ecuaciones 4.42. Onservamos que en el nodo i = I max se tiene un término identico a la ec. 4.46. Por lo tanto, ese término, de acuerdo a la condición de frontera impuesto deber ser eliminado. Condición de Frontera Interna. Similarmente para la frontera interna, esto es, r = r w , se realiza el mismo procedimiento. Llegando al resultados siguiente, ·

r 1+ 1

2

p2 − p1 r2 − r1

¸n+1

·

=−

qB µ 2πkh

¸ 1+ 12

(4.47)

Resolviendo para el gasto q, que en realidad para este frontera es el gasto que aporta el yacimiento, esto es, q s f , se tiene lo siguiente, Ã ! 2πr 1+ 1 h µ k ¶ ¡ ¢n+1 2 qs f = − p2 − p1 r2 − r1 µB 1+ 1 (4.48) µ ¶ 2 µ ¶ ¡ ¢n+1 k A qs f = − p2 − p1 ∆r 1+ 21 µB 1+ 1 2

En términos de la transmisibilidad, ¡ n+1 ¢ n q s f = −T1+ − p 1n+1 1 p2 2

(4.49)

La ec. 4.49 debe ser acoplada al sistema de ecuaciones 4.42. Observamos que en nodo i = 1 se tiene un término identico a la ec. 4.49. Por lo tanto, este término, de acuerdo a la condición de frontera impuesto deber ser considerado como q s f para satisfacer la condición de frontera interna. Sustituyendo las dos condiciones de frontera, esto es, ecs. 4.46 y 4.49 en el sistema de ecuaciones 4.42, se tiene lo siguiente:

38

n+1 n C p 1 −p 1 B ∆t

i=1

p n+1 −p 1n+1 i +s r 1+ 1

h 2 − 2πkh B µ ¡ ∆r ¢

=

q sup

2

Tn

i=2 i=3

2+ 12

Tn

¡

3+ 12

¡

¢ p 3n+1 − p 2n+1 + q s f

¢ p 4n+1 − p 3n+1 − T n

¡

3− 12

p 3n+1 − p 2n+1

¢

=

(Vp bc t )n2

(p 2n+1 −p 2n ) ∆t

=

(Vp bc t )n3

(p 3n+1 −p 3n ) ∆t

...

...

=

...

...

...

=

...

...

...

=

...

0−Tn

i = Imax

¡

I max− 21

n+1 p n+1 I max − p I max−1

= (Vp bc t )nImax

¢

(4.50)

n (p n+1 I max −p I max ) ∆t

Para considerar el factor de daño en la vecindad del pozo, esto es, nodo 2. Se considerará la expresión que se obtuvo en la ec. 4.16 como expresión para q s f en el sistema 4.50. De esta manera, habiendo acoplado las condiciones de frontera, el sistema de ecuaciones a solucionar, considerando almacenamiento y daño a la formación es,

n+1 n C p 1 −p 1 B ∆t

i=1 i=2 i=3

Tn

¡

2+ 12

Tn

3+ 12

¡

p n+1 −p 1n+1 i +s r 1+ 1

h 2 − 2πkh B µ ¡ ∆r ¢

¢ p 3n+1 − p 2n+1 − 2πkh Bµ

¢ p 4n+1 − p 3n+1 − T n

3− 12

¡

=

2 p n+1 −p n+1 h¡ 2 ¢ 1 i ∆r +s r 1+ 1 2

p 3n+1 − p 2n+1

¢

q sup

=

(Vp bc t )n2

(p 2n+1 −p 2n ) ∆t

=

(Vp bc t )n3

(p 3n+1 −p 3n ) ∆t

...

...

=

...

...

...

=

...

...

...

=

...

i = Imax

−T n

I max− 21

¡

n+1 p n+1 I max − p I max−1

¢

= (Vp bc t )nImax

(4.51)

n (p n+1 I max −p I max ) ∆t

39

4.4.4 C ONSTRUCCIÓN DEL S ISTEMA DE E CUACIONES Ax = b El desarrollo del sistema de ecuaciones 4.51 para cada nodo de la malla numérica se muestra en la siguiente página. Se puede observar lo siguiente: • Cada uno de los nodos genera 3 incógnitas, excepto el primer nodo y el último, los cuales generan dos. • Cada uno de los nodos genera 3 coeficientes, excepto el primer nodo y el último, los cuales generan dos. • El nodo en cuestión tiene relación únicamente con sus nodos vecinos Definiendo los coeficientes, uno para cada incógnita de presión, esto es, a i , b i , y c i , para p in+1 , p in+1 , p in+1 , respectívamente; y otro extra para el término del lado derecho, el cual con−1 +1 tiene todos los términos conocidos, esto es, d i , se tiene los siguiente, Para el nodo i = 1, a 1 =0 



2πkh  C  h¡ ¢ i b1 =  + B ∆t B µ ∆r r 1+ 1 + s 2

  c1 = − 



Bµ

d 1 =q sup +

2πkh h¡ ¢ ∆r r

1+ 21

+s

(4.52)

 i

C pn B ∆t 1

Expresando la Transmisibilidad, T1+ 1 , de la siguiente manera, 2

Bµ

2πkh h¡ ¢ ∆r r

1+ 12

µ ¶ 2πhr 1+ 1 µ k ¶ ∆r 2 i= ∆r 1+ 1 B µ 1+ 1 ∆r + r s 1+ 21 +s 2 2 µ ¶ ∆r =T1+ 1 2 ∆r + r s 1+ 21

(4.53)

40

41

... ...

...;

...;

i = Imax;

!

T n 1 p in+1 −1 i− 2

2

2

!

p 1n+1

p 1n+1

h¡ 2πkh i ¢ B µ ∆r +s r 1+ 1

h¡ 2πkh i ¢ B µ ∆r +s r 1+ 1

+

...

Ã

C B ∆t

...;

i >= 3;

i = 2;

i = 1;

Ã

− + 2

h¡ 2πkh i ¢ B µ ∆r +s r 1+ 1

2

+

h¡ 2πkh i ¢ B µ ∆r +s r 1+ 1

i

´ Vp bc t n ∆t 2

p 2n+1 ³

!

!

p 2n+1

i− 2

...

...

...

p n+1 Tn I max− 12 I max−1

i+ 2

µ ³ ´ ¶ Vp bc t n n n p in+1 − T 1 + T 1 + ∆t

Tn 1 2+ 2

Ã

−

Ã

i max− 21

µ − Tn

+

³

´ Vp bc t n ∆t i max

...

...

...

i+ 2

+T n 1 p in+1 +1

2+ 2

+T n 1 p 3n+1

¶

p n+1 I max

= −

=

=

=

=

=

=

³

³

³

...

...

...

(4.54)

p nImax

´ Vp bc t n n p ∆t i i

´ Vp bc t n n p ∆t 2 2

´ Vp bc t n ∆t I max

−

−

q sup + BC∆t p 1n

Además, multiplicando solo este renglón, (i = 1), por (-1). Los coeficientes se puedes expresar de la siguiente manera,

a 1 =0 Ã

! µ ¶ C ∆r b1 = − + T1+ 1 2 B ∆t ∆r + r s 1+ 12 µ

c 1 =T1+ 1 2

∆r ∆r + r s

(4.55)

¶ 1+ 12

C d 1 = − q sup + pn B ∆t 1 µ

¶

Entonces, para el nodo i = 2, a 2 =c 1 n c 2 =T2+ 1 2

µ

d2 = −

Vp bc t

¶n

∆t

2

(4.56)

p 2n

µ ¶ d2 b2 = − c2 + a2 − n p2

Para el resto de los nodos i = 3, 4, 5, ..., I max − 1, a i =c i −1 c i =Tin+ 1 2

di = −

(Vp bc t )ni ∆t "

(4.57)

p in

di bi = − ci + ai − n pi

#

Para el nodo i = I max,

42

a I max =c I max−1 c I max =0 d I max = −

(Vp bc t )nImax ∆t

p nImax

"

d I max b I max = − a I max − n p I max

(4.58)

#

Considerando los coeficientes 4.52, 4.56, 4.57 y 4.58, en la ec. 4.54, se tiene, n+1 a i p in+1 + c i p in+1 −1 + b i p i +1 = d i

para i = 1,2,3,...Imax

(4.59)

n =1,2,...

Nótese que el sistema 4.59 es similar al obtenido en la solución del problema en 1D, horizontal en coordenadas cartesianas. Escribiendola como un sistema de ecuaciones, se tiene lo siguiente:



b1  a  2               

c1 b2 a3

n 

c2 b3 a4

c3 b4 a5

c4 b5 a6

c5 b6 a7

c6 b7 a8

c7 b8 a9

c8 b9 a 10

                c9  b 10

                

p1 p2 p3 p4 p5 p6 p7 p8 p9 p 10

n+1                 

         =        

d1 d2 d3 d4 d5 d6 d7 d8 d9 d 10

n                 

(4.60)

Los coeficientes a 0 y c 10 en este se usan en algunas problemas para acoplar algunas condiciones de frontera. En este caso se consideran 0. El sistema de ecuaciones 4.60, es un sistema tridiagonal (disperso) y se puede solucionar con el algoritmo de Tomas para matrices tridiagonales.

4.5 A LGORITMO DE S OLUCIÓN la siguiente secuencia de instrucciones es una propuesta de algoritmo de solución para el sistema de ecuaciones 4.60,

43

1. Lectura de información 2. Definición de Condiciones Iniciales 3. Definición de Condiciones de Frontera 4. Verificar consistencia de Unidades 5. Definir un ∆t , por ejemplo, 1 sec. 6. Obtención de los coeficientes a i , b i , c i , d i 7. Resolver el Sistems tridiagonal con el Algoritmo de Thomas 8. Reasignar o actualizar las nuevas presiones p in ⇐ p in+1 , para i = 1, 2, 3, ...I max. 9. Elegir otro ∆t , por ejemplo, podemos incrementar el ∆t anterior, esto es, ∆t = 1.2 ∗ ∆t , para inciar otro calculo de p in+1 . 10. Regresar al paso 6. 11. Parar hasta que se cumpla el tiempo de simulación total.

El paso 6 se puede realizar procesando nodo por nodo, es decir, construir la matriz de coeficientes, A, ec. 4.60 renglón por renglón. Para evitar el uso de arreglos, se suguiere renombrar las transmisibilidades como se muestra en la Figura 4.8. Para cada nodo se tendrán dos transmisibilidades, esto es, T1 y T2. De esta manera, con el uso de instrucciones de condicion (i f ...el se..end ) se podrá, facilmente, definir cuando una transmisibilidad, especialmente en los nodos de las fronteras, es cero. Por ejemplo, en el nodo 1 se determinará únicamente la transmisibilidad en T2, ya que la transmisibilidad T1 es nula. En el ultimo nodo, I max, la transmisibilidad T2 será igual a 0 (no existe flujo en la frontera externa!). El resto de los nodos tendrán las dos transmisibilidades. De esta manera no hay necesidad de crear vectores para almacenar variables. Se analiza nodo por nodo y se obtienen los coeficientes. El siguiente código en MatLab ilustra un poco más el algoritmo de solución,

1 2 3 4 5

clc; clear all; close all; %Structure containing the SI conversion unit factors u=units;

6 7

% Lectura de Datos

8 9

% Cambio de Unidades al SI

10

44

r  rw

i

i 𝑇𝑥1 ← 𝑇

1 𝑖−2

𝑇

1 𝑖+2

→ 𝑇𝑥2

Figure 4.8: Tratamiento de Transmisibilidad

11 12 13 14 15

% Domensionamiento de arreglos rCell=zeros(Imax,1); bCell=zeros(Imax,1); vCell=zeros(Imax,1); % etc..

16 17 18 19

%Initial Condition PrsOld(:,1)=PrsIni; PrsNew(:,1)=PrsIni;

20 21 22

%Radial Grid [ rCell,bCell,vCell ] = RGrid(Re,Rw,Imax,dz);

23 24 25

Tacum=0.0; %Cumulative time while Tacum < Tsim

26 27 28 29 30 31 32

% Para designar Transmisibilidades T1 y T2 for i=1:Imax if(i==1) Tx1=0.0; else Area=2.0*pi*bCell(i−1)*dz;

45

Dr=rCell(i)−rCell(i−1); Tx1=(Area/Dr)*(Perm*fvfo/Vis);

33 34

end

35 36

if(i==Imax) Tx2=0.0; else Area=2.0*pi*bCell(i)*dz; Dr=rCell(i+1)−rCell(i); Tx2=(Area/Dr)*(Perm*fvfo/Vis); end

37 38 39 40 41 42 43 44 45

a(i)= b(i)= c(i)= d(i)=

46 47 48 49 50

end

51 52 53

end

54 55 56

%Thomas Algorithm [PrsNew] = thomas(Imax,a,b,c,d );

57 58 59

%Cumulative timestep Tacum=Tacum+Dt;

60 61

% Impresiones de tiempos y presiones

62 63 64 65

% Actualizaci\'on de presiones para otro Dt PrsOld(:,1)=PrsNew(:,1);

66 67 68

%Increse of timestep Dt=Dt*1.2;

69 70

end %while

46

5 S OLUCIÓN N UMÉRICA PARA YACIMIENTOS N ATURALMENTE F RACTURADOS El siguiente planteamiento para la solución numérica de flujo de una fase, en dirección radial, yacimientos naturalmente fracturados, YNF, es una extensión del flujo monofásico en Yacimientos no Fracturados. Para conceptualizar a un YNF y estar en posición de modelarlo veamos la Figura 5.1 , la cual muestra una afloramiento real de un YNF. En general, un YNF es conceptualizado como dos sistemas que interactuan en espacio y tiempo. Estos dos sistemas que estan presentes e interactuan son separados, abstractamente, para poder modelarlos. Como se muestra en la parte inferior de la Figura 5.1, un sistema corresponde a la matriz de la roca (porosidad primaria) donde se almacena la mayor parte del fluido y presenta una permeabilidad baja. El otro sistema es la red de fracturas (porosidad secundaria), con alta permeabilidad pero donde se tiene un almacenamiento pequeño de fluidos. Los bloques de matriz se encuentran desconectados, rodeados por una red de fracturas ortogonales continuas . A esta descripción, en la literatura técnica, se le conoce como modelo de doble porosidad-una permeabilidad o modelo de Warren & Root.

5.1 C ONCEPTUALIZACIÓN DE UN YACIMIENTOS N ATURALMENTE F RACTURADO En la idealización de un YNF, con el propósito de modelar numéricamente el flujo de fluidos, han surgido varias conceptualizaciones. La Figura 5.2 muestra algunos modelos que se han sido utilizados para modelar el flujo de fluidos en YNF. Para entender mejor, se hace uso de un esquema con tres celdas, como las usadas en la solución numérica para representar, esquemáticamente, el flujo de fluidos en y entre las celdas, y luego, hacia el pozo productor. En la parte superior de la figura se muestra, esquematicamente, la conceptualización de un yacimiento no fracturado, homogeneo. Para ser consistentes con la nomenclatura, a este modelo se pudiera identificar como 1φ − 1k. Solo un sistema esta presente, esto es, matriz homogenea, no fracturada, con cierta porosidad y permeabilidad presentes. Este es la conceptualizacin ´ que se ha venido trabjando hasta ahora. En la parte intermedia de la Fig. 5.2, se presenta la conceptualización del modelo de doble porosidad - doble permeabilida, 2φ − 2k. En este modelo, los dos medios que han sido identificados se consideran continuos, es decir, ambos sistemas transportan fluido del y a través del yacimiento hacia los pozos productores. Por último, se presenta la conceptualización de un modelo de doble porosidad-una permeabilida, esto es, 2φ − 1k. En este modelo, el sistema 1, generalmente esta representado como la red de fracturas y es considerado comoel medio continuo. A través de este sistema el fluido fluye hacia los pozos. El Sistema 2, los bloques de matriz, alimentan, localmente, con fluido a las fracturas. Este es el medio descontinuo, esta es la razón de por que no hay conexión entre

47

Figure 5.1: Conceptualización de una Formación Naturalmente Fracturada

48

celdas numéricas. Existen otas conceptializaciones, por ejemplo, triple porosidad, fractura discreta, etc. Dependiendo del mecanismo de producción en un YNF, esto es, drene gravitacional, imbibición, difusión etc. y de las características estructurales y geológicas del yacimiento, la conceptualización es eleguida. Para el propósito de este curso y por simplicidad, el modelo de doble porosidad, una permeabilidad, 2φ − 1k, sera desarrollado para flujo en una sola fase, en forma horizontal. Es importante hacer la anotación para explicar cuando se menciona 1φ o 1k, esto no quiere indicar que únicamente existe un valor de porosidad o un valor de permeabiliad. Generalmente se presenta una distribución de estos parámetros. Lo que se quiere indicar es que solo se ha identificado un solo sistema, "un tipo de roca" con ciertas propiedades y que se pueden caracterizar mediante un solo sistema. En otras palabras, no hay sistemas contrastantes que modifiquen por si mismos el flujo (k) y el almacenamiento del fluido (φ).

5.2 F LUJO DE M ONOFÁSICO DE ACEITE EN YACIMIENTOS N ATURALMENTE F RACTURADOS Para comparar las Ecuaciones en Derivadas Parciales, EDP, que describen el flujo radial, horizontal y monofásico de un fluido en un YNF, se presentan, a continuación, las tres conceptualizaciones mostrados en la Figura 5.2.

5.2.1 M ODELO DE YACIMIENTO N O F RACTURADO. 1φ − 1k La EDP que describe el flujo de una fase es la que se ha desarrollado hasta ahora, esto es, ec. 4.1: 1 ∂ r ∂r

·

¸ µ ¶ k ∂p ∂ φ r = µB ∂r ∂t B

(5.1)

Obviamente con sus respectivas CI y CF.

5.2.2 M ODELO DE YACIMIENTO F RACTURADO. D OBLE P OROSIDAD - D OBLE P ERMEABILIDAD. 2φ − 2k Como se establecio anteriormente. Este modelo considera, ambos sitemas identificados, continuos en todo el dominio del yacimiento. Los dos medios identificados transportan fluidos hacia los pozos productores. Como se discutirá posteriormente, existen rangos de contribución de cada medio al transporte de fluidos. Debido a esto se han publicado una clasificación de tipos de YNF., Figura

49

(a)

i-1

i

i+1

pozo

i-1

i

i+1

pozo

(b) pozo

i-1

i

i+1

pozo

(c)

Sistema 1 (fractura) Flujo S1-S2 (misma celda)

Sistema 2 (matriz) Flujo entre Celdas

Figure 5.2: Modelos de Yacimientos. (a) Homogeneo. (b) Doble Porosidad-Doble permeabilidad. (c) Doble Porosidad - Una Permeabilidad

50

Figure 5.3: Gráfica de Nelson. Clasificación de Yacimientos Naturalmente Fracturados.

Por simplicidad, considerando las propiedades PVT constantes, la EDP que describe el flujo de un mode de 2φ − 2k son las siguientes: Medio Continuo - Fracturas (f)

µ ¶ ∂ φf = ∂t B

(5.2)

¸ µ ¶ k m ∂p m ∂ φm ˜ r + qm f = µB ∂r ∂t B

(5.3)

1 ∂ r ∂r

·

1 ∂ r ∂r

·

kf µB

r

∂p f ∂r

¸

− q˜m f

Medio Continuo - Matriz (m)

Nótese el parámetro incluyendo la transferencia de fluidos entre los dos medios, esto es, q˜m f . Este representa el intercambio , por causa de un gradiente de presiones, de fluidos entre los dos sistemas. En la ec. 5.2 es negativo y la ec. 5.3 es positivo. Es simplemente un

51

balance de masa. El fluido que gana un medio, lo pierde el otro. También, es necesario distinguir las propiedades petrofísicas para cada medio, esto es lo que causa el adjetivo de heterogeneidad y la necesidad de recurrir a modelos relativamente más sofisticadas, obligando, por lo tanto, a la distición de medios o sistemas de flujo. Nótese, también, las distinción de presiones. Las propiedades PVT, como se mencionó, se mantendran constantes. Observando las ecuaciones 5.2 y 5.3, son similares a la ec. 5.1, la cual, ya ha sido resuelta en forma analítica y numérica. Esto simplificará el trabajo y la curva de aprendizaje para comprender el flujo de fluidos en YNF. 5.2.3 M ODELO DE YACIMIENTO F RACTURADO. M ODELO DE D OBLE P OROSIDAD - U NA P ERMEABILIDAD. 2φ − 1k Como se establecio anteriormente. Este modelo considera, de los dos sitemas identificados, que solo un sistema es continuo en todo el dominio del yacimiento. El segundo sistema es discontinuo y en la mayoria de los casos, este sirve como almacenador y provedor de fluidos hacia el sistema continuo. Generalmente se considera el sistema continuo como la red de fracturas y el medio descontinuo como los bloques de matriz. Por simplicidad, considerando las propiedades PVT constantes, la EDP que describe el flujo de un mode de 2φ − 1k son las siguientes: Medio Continuo - Fracturas (f)

1 ∂ r ∂r

·

kf µB

r

∂p f

¸

∂r

− q˜m f =

µ ¶ ∂ φf ∂t B

(5.4)

Medio Discontinuo - Matriz (m)

+q˜m f =

µ ¶ ∂ φm ∂t B

(5.5)

Observando la ec. 5.5, esta faltando el primer término, el cual se establecio anteriormente que es el término de flujo, y tiene sentido, debido a que el sistema de bloque de matriz es considerado como sistema discontinuo, de ahi el nombre de doble porosidad - una peremabilidad, 2φ − 1k. En otras palabras, el sistema de bloques de matriz, no transporta fluidos a grandes distancias, estos solo intercambian, localmente, fluidos con el sistema de fracturas, las cuales lo transportan hacia los pozos productores.

52

Habiendo discutido las dos conceptualizaciones más comunes para modelar, numéricamente, el flujo de fluidos a través de medios porosos fracturados, se procederá plantear numéricamente el modelo de doble porosidad - una permeabilidad, 2φ − 1k.

5.3 D ISCRETIZACIÓN DEL M ODELO DE D OBLE P OROSIDAD - U NA P ERMEABILIDAD Como ya se ha discretizado, desarrollado y solucionado el problema de 1φ − 1k, procederemos a incluir lo siguiente:

• Incluir el el término de intercambio de fluidos entre la matriz y la fractura, discretizarlo y platear su solución. • Discretizar la EDP en el sistema de matriz para posteriormente acoplarla a la formulación numérica de la fractura • Plantear la solución la formulación matemática del problea a resolver, es el siguiente, 1 ∂ r ∂r

·

kf µB

r

∂p f ∂r

¸

− q˜m f +q˜m f

µ ¶ ∂ φf = ∂t B µ ¶ ∂ φm = ∂t B

(5.6)

y sus condiciones iniciales y de frontera:

p m = p f (r, t = 0) = p 0 ¸ · ∂p f qB µ =− r ∂r r w 2πk f h · ¸ ∂p f qB µ r =− =0 ∂r r e 2πk f h

(5.7)

Nótese que la ecuación de la matriz, ec. 5.6b, necesita solo una CI. Eliminando en la ec. 5.6a, el subindice ( f ) y el término de intercambio de fluidos matriz-fractura, q˜m f , esta, ya ha sido planteda y resuelta numéricamente. Por lo tanto, únicamente nos esforzaremos en trabajar con el término de intercambio de fluidos matrz-fractura, q˜m f , lo cual será a continuación.

5.3.1 T ÉRMINO DE T RANSFERENCIA DE F LUIDOS M ATRIZ -F RACTURA Existen varias funciones de transferencia de fluidos matriz-fractura publicadas, la mas simple es la propuesta por Warren & Root. Esta es,

53

k =σ µB ·

q˜m f

¸

(p f − p m )

(5.8)

m

donde; σ = Factor de forma, (m −2 ) k m = Permeabilidad de matriz, (m 2 ) p f = Presión de fractura, (P a) µ = Viscosidad, (P a − Sec) B = factor de Volumen , (m 3 /Sm 3 ) La ec. ec. 5.8 muestra que el intercambio es función de la diferencia de presiones entre el sistema de fracturas y el sistema de bloques de matriz. En la ec. 5.8 la tilde en, q˜m f , significa que este intercambio de fluidos esta expresado por unidad de volumen de roca. Por lo tanto, en la etapa de discretización de la ec. 5.2, al momento de multiplicarla por el volumen de roca, Vr se tendrá lo siguiente:

q m f = q˜m f Vr = σVr

·

k µB

¸

(p f − p m )

(5.9)

m

Analizando las unidades de la ec. 5.9, estas, están en términos de unidades de gasto a c.e.:

Sm 3 seg ,

"

q m f = σ(m

−2

3

)Vr (m )

#

k(m 2 )

(p f (P a) − p m (P a)) ⇒

m3

µ(P a ∗ sec)B Sm 3

m

Sm 3 seg

(5.10)

La ec. 5.9 será utilizada en la discretización del problema.

5.3.2 D ISCRETIZACIÓN DE MEDIO C ONTINUO - F RACTURAS Para simplificar el trabajo, haciendo uso de la discretización del problema homogeneo, no fracturado, esto es, ec. 4.41, e incluiremos, únicamente , el término de transferencia matrizfractura. La ec. 4.41 es

(p in+1 − p in ) ¡ ¢ ¡ n+1 ¢ n+1 n n+1 n Tin+ 1 p in+1 − p − T p − p = (V bc ) p t i +1 i i i −1 i − 12 ∆t 2 para i = 2,3,...Imax; n =1,2,...

54

Nótese que solo trabajaremos en la parte del medio poroso. La ecuación para el pozo será la misma. Entonces, considerando la ec. 4.41 como base, modificando los subíndices para distinguir entre propiedades de fractura y matriz, y colocando el tiempo al cual se evaluara el término de tranferencia matriz-fractura, se tiene lo siguiente,

T fn,i + 1 2

³

n+1 p n+1 f ,i +1 − p f ,i

´

− T fn,i − 1 2

³

n+1 p n+1 f ,i − p f ,i −1

´

− σVr

·

k µB

¸n m,i

n+1 (p n+1 f ,i − p m,i )

= (Vp bc t )nf,i

(p n+1 − p nf,i ) f ,i

(5.11)

∆t

para i = 2,3,...Imax; n =1,2,... Para el medio discontinuo, es decir, el sistema de bloques de matriz, la discretización esta hecha. Tomemos la ec. 5.11. Eliminemos los términos de flujo, cambiemos el signo al término de tranferencia matriz-fractura y cansideremos el subíndice de la matriz donde aplique, se tiene lo siguiente: ·

+σVr

k µB

¸n m,i

n+1 n (p n+1 f ,i − p m,i ) = (V p bc t )m,i

n+1 n (p m,i − p m,i )

∆t

(5.12)

para i = 2,3,...Imax; n =1,2,... El acoplamiento de las condiciones de frontera se obviarán. Es un procedimiento similar al hecho anteriormente. Ahora, estaria aplicandose al medio continuo, esto es, al sistema de red de fracturas. 5.3.3 A COPLAMIENTO DE LA E CUACIÓN DE LA M ATRIZ EN LA E CUACIÓN DE LA F RACTURA Observando las ecs. 5.11 y la ec. 5.12, podemos apreciar que únicamente existen incógnitas de presiones, todo el resto de variables y términos son conocidos al tiempo (n). Por lo tanto, n+1 la p n+1 , en la ec. 5.12 puede ser puesta en términos de la p m,i y ser sustituida en la ec. 5.11. f ,i De esta manera, creando las siguientes variables,

¸ km n µB i · ¸ Vp,m bc t ,m n

Ψi = σVr Ωm,i =

Ω f ,i =

·

∆t i · ¸n Vp, f bc t , f ∆t

(5.13)

i

55

La ec. 5.12 pasa a ser la siguiente, n+1 n+1 n +Ψi (p n+1 f ,i − p m,i ) = Ωm,i (p m,i − p m,i )

(5.14)

n+1 Solucionando para p m,i ,

n+1 p m,i

=

n Ψi p n+1 + Ωm,i p m,i f ,i

(5.15)

Ψi + Ωm,i

Sustituyendo la ec. 5.15 en la ec. 5.11

T fn,i + 1 2

³

n+1 p n+1 f ,i +1 − p f ,i

´

− T fn,i − 1 2

³

n+1 p n+1 f ,i − p f ,i −1

´

Ã

− Ψi

p n+1 f ,i −

n Ψi p n+1 + Ωm,i p m,i f ,i

!

Ψi + Ωm,i

(5.16)

n = Ω f ,i (p n+1 f ,i − p f ,i )

Nótese que la ec. 5.16 presenta, ahora, incógnitas únicamente del sistema de fracturas, esto es, p n+1 , para i = 2, 3, ..., I max f ,i Considerando el nodo 1, donde se localiza el pozo, y acoplando, similarmente, las condiciones de frontera, como se hizo en el caso de yacimientos homogeneos, esto es, no fracturados, se llegará un sistema de ecuaciones similar a la ec. 4.51, esto es,

56

57

³

2

p n+1 −p n+1 f ,2 f ,1 i +s 1 r 1+

h − 2πkh B µ ¡ ∆r ¢

I max− 12

³

µ

´ n+1 p n+1 − p − ΨI max p n+1 − f ,I max f ,I max−1 f ,I max

...

...

i = Imax −T n

...

...

3− 2

...

3+ 2

µ ΨI max +Ωm,I max

n ΨI max p n+1 +Ωm,I max p m,I max f ,I max

µ µ ¶¶ ´ n n+1 n+1 Ψ2 p n+1 +Ωm,2 p m,2 f ,2 2πkh hp f ,2 −p f ,1 i n+1 n+1 n+1 p − p − p − − Ψ ¢ ¡ 2 Bµ Ψ2 +Ωm,2 ∆r f ,2 f ,3 f ,2 2+ 12 +s r 1+ 1 2 µ µ ¶¶ ³ ´ ³ ´ n Ψ3 p n+1 +Ωm,3 p m,3 f ,3 n+1 n n+1 n+1 n+1 T n 1 p n+1 − p − T p − p − Ψ p − 3 1 Ψ3 +Ωm,3 f ,4 f ,3 f ,3 f ,2 f ,3

Tn

n+1 n C p f ,1 −p f ,1 B ∆t

...

i=3

i=2

i=1

¶¶

...

...

...

Ω f ,3 (p n+1 − p nf,3 ) f ,3

Ω f ,2 (p n+1 − p nf,2 ) f ,2

q sup

(5.17)

= Ω f ,I max (p n+1 − p nf,I max ) f ,I max

=

=

=

=

=

=

Si eliminamos el subindice ( f ) y el término de transferencia de fluidos matriz-fractura, este sistema de ecuaciones, es similar al obtenido en el caso de yacimiento homogeneo. Por lo tanto, el modelo de 2φ − 2k, puede ser visto como una extensión al modelo de 2φ − 1k. Continuando con el desarrollo, se agruparán términos semejantes y para estar en posición de escribirlo como un sistema de ecuaciones de la forma Ax = d , similar a los sistemas anteriormente descritos.

58

5.3.4 C ONSTRUCCIÓN DEL S ISTEMA DE E CUACIONES Ax = b El desarrollo, para el caso de 2φ − 1k, es similar al caso de 1φ − 1k, excepto que se deben de identificar las presiones para cada medio, esto es, p m y p f e incluir el término de transferencia de fluidos matriz-fractura. Procediendo de forma similar al caso anterior, el desarrollo del sistema de ecuaciones 4.51 para el modelo 2φ − 21k se muestra en la siguiente página. Se puede observar lo siguiente: • El nodo en cuestión tiene relación únicamente con sus nodos vecinos en el medio continuo, esto es, para el sistema de fracturas. • También se debe notar la distinción entre la permeabilidad ,k, porosidad, φ, y compresibilidad, c t para ambos medios. • Las propiedades del fluido serán consideradas constantes. • Cada uno de los nodos genera 3 incógnitas, excepto el primer nodo y el último, los cuales generan dos. • Cada uno de los nodos genera 3 coeficientes, excepto el primer nodo y el último, los cuales generan dos. Definiendo los coeficientes, uno para cada incógnita de presión, esto es, a i , b i , y c i , para p n+1 , p n+1 , p n+1 , respectívamente; y otro extra, para el término del lado derecho, el cual f ,i −1 f ,i f ,i +1 contiene todos los términos conocidos, esto es, d i , se tiene los siguiente, Procediendo de forma similar, al modelo 1φ − 1k, donde las ec. para el nodo 1 no cambia, a 1 =0 Ã

! ¶ µ C ∆r b1 = − + T1+ 1 2 B ∆t ∆r + r s 1+ 12 µ ¶ ∆r c 1 =T1+ 1 2 ∆r + r s 1+ 12 µ ¶ C p nf,1 d 1 = − q sup + B ∆t

(5.18)

Solo este renglón, (i = 1), se ha multiplicado por (-1) y en la ec. 5.18, se considerado la Transmisibilidad, T1+ 1 , de la siguiente manera, 2

59

Bµ

2πk f h h¡ ¢ ∆r r

1+ 12

+s

i=

2πhr 1+ 1 µ k f ¶ 2 ∆r 1+ 1 2

µ

=T1+ 1

2

Bµ

∆r ∆r + r s

µ

1+ 21

∆r ∆r + r s

¶ 1+ 21

(5.19)

¶ 1+ 21

Definiendo la siguiente variable, αi , para i = 2, 3, 4, ...I max,

αi =

Ψi Ωm,i Ψi + Ωm,i

(5.20)

Entonces, para el nodo i = 2, se tiene, a 2 =c 1 n c 2 =T2+ 1 2

(5.21)

¡ ¢ b 2 = − c 2 + a 2 + Ω f ,2 + α2 n d 2 = − (Ω f ,2 p nf,2 + α2 p m,2 )

Para el resto de los nodos i = 3, 4, 5, ..., I max − 1, a i =c i −1 c i =Tin+ 1 2

(5.22)

¡ ¢ b i = − c i + a i + Ω f ,i + αi n d i = − (Ω f ,i p nf,i + αi p m,i )

60

61

... ...

...;

i− 2

...;

i = I;

!

!

p n+1 f ,1

p n+1 f ,1

T n 1 p n+1 f ,i −1

2

2

2πk h h¡ ¢ f i B µ ∆r +s r 1+ 1

+

2πk h h¡ ¢ f i +s B µ ∆r r 1+ 1

...

Ã

C B ∆t

...;

i >= 3;

i = 2;

i = 1;

Ã

−

+ 2

2πk h h¡ ¢ f i B µ ∆r +s r 1+ 1

2

!

p n+1 f ,2 !

+ Ωm,2 + α2 p n+1 f ,2

2πk h h¡ ¢ f i +s B µ ∆r r 1+ 1

i+ 2

...

...

...

T n 1 p n+1 f ,I −1 I−2

i− 2

µ ¶ n n − T 1 + T 1 + Ωm,i + αi p n+1 f ,i

Tn 1 2+ 2

Ã

−

Ã

I − 12

µ − Tn

¶ + Ωm,I + αI p n+1 f ,I

...

...

...

i+ 2

+T n 1 p n+1 f ,i +1

2+ 2

+T n 1 p n+1 f ,3

q sup + BC∆t p nf,1

=

=

=

=

=

(5.23)

n −Ω f ,I p nf,I − αI ∗ p m,I

...

...

...

n −Ω f ,i p nf,i − αi ∗ p m,i

n = −Ω f ,2 p nf,2 − α2 ∗ p m,2

=

Para el nodo i = I max, a I max =c I max−1 c I max =0 (5.24) ¡ ¢ b I max = − a I max + Ω f ,I max + αI max n d I max = − (Ω f ,I max p nf,I max + αI max p m,I max )

Considerando los coeficientes 5.18, 5.21, 5.22 y 5.24, en la ec. 5.23, se tiene, n+1 n+1 a i p n+1 f ,i −1 + b i p f ,i + c i p f ,i +1 = d i

para i = 1,2,3,...Imax

(5.25)

n =1,2,...

Nótese que el sistema 5.25 es similar al obtenido en la solución del modelo 1φ − 1k en 1D, horizontal en coordenadas cartesianas. Escribiendola como un sistema de ecuaciones, se tiene lo siguiente:



b1  a  2               

c1 b2 a3

n 

c2 b3 a4

c3 b4 a5

c4 b5 a6

c5 b6 a7

c6 b7 a8

c7 b8 a9

c8 b9 a 10

                c9  b 10

                

p f ,1 p f ,2 p f ,3 p f ,4 p f ,5 p f ,6 p f ,7 p f ,8 p f ,9 p f ,10

n+1                 

         =        

d1 d2 d3 d4 d5 d6 d7 d8 d9 d 10

n                 

(5.26)

Los coeficientes a 0 y c 10 se usan en algunas problemas para acoplar algunas condiciones de frontera. Por el momento, y por simplicidad, no se considerarán. El sistema de ecuaciones 5.26, es un sistema tridiagonal (disperso) y se puede solucionar con el algoritmo de Tomas para matrices tridiagonales.

5.4 A LGORITMO DE S OLUCIÓN PARA EL M ODELO DE 2φ − 1k la siguiente secuencia de instrucciones es una propuesta de algoritmo de solución para el sistema de ecuaciones 5.26, 1. Lectura de información

62

2. Definición de Condiciones Iniciales 3. Definición de Condiciones de Frontera 4. Verificar consistencia de Unidades 5. Definir un ∆t , por ejemplo, 1 sec. 6. Obtención de los coeficientes a i , b i , c i , d i 7. Resolver el Sistems tridiagonal con el Algoritmo de Thomas para las incógnitas de las para i = 1, 2, 3, ..., I max. fracturas, esto es, p n+1 f ,i 8. Habiendo obtenidos las incógnitas de las fracturas, se soluciona para las incógnitas de n+1 la matriz, esto es, p m,i para i = 2, 3, ..., I max, usando la ec. 5.15. (para i = 1 no se tiene medio poroso) 9. Reasignar o actualizar las nuevas presiones paa ambos medios, esto es, p nf,i ⇐ p n+1 , f ,i n+1 n , para i = 2, 3, ...I max ⇐ p m,i para i = 1, 2, 3, ...I max y p m,i

10. Elegir otro ∆t , por ejemplo, podemos incrementar el ∆t anterior, esto es, ∆t = 1.2 ∗ ∆t , n+1 para inciar otro calculo de p n+1 y de p m,i f ,i 11. Regresar al paso 6. 12. El procedimiento finaliza hasta que se cumpla el tiempo de simulación total.

El paso 6 se puede realizar procesando nodo por nodo, es decir, construir la matriz de coeficientes, A, ec. 4.60 renglón por renglón como se seguirio en el modelo 1φ − 1k El siguiente código en MatLab ilustra un poco más el algoritmo de solución,

1 2 3 4 5

clc; clear all; close all; %Structure containing the SI conversion unit factors u=units;

6 7

% Lectura de Datos

8 9

% Cambio de Unidades al SI

10 11 12 13 14 15

% Domensionamiento de arreglos rCell=zeros(Imax,1); bCell=zeros(Imax,1); vCell=zeros(Imax,1); % etc..

16 17 18

%Initial Condition ambos medios, matriz y fractura PrsOldf(:,1)=PrsIni;

63

19

PrsOldm(:,1)=PrsIni;

20 21 22

%Radial Grid [ rCell,bCell,vCell ] = RGrid(Re,Rw,Imax,dz);

23 24 25

Tacum=0.0; %Cumulative time while Tacum < Tsim

26

% Para designar Transmisibilidades T1 y T2 (fractura, medio continuo) for i=1:Imax if(i==1) Tx1=0.0; else Area=2.0*pi*bCell(i−1)*dz; Dr=rCell(i)−rCell(i−1); Tx1=(Area/Dr)*(Permf*fvfo/Vis); end

27 28 29 30 31 32 33 34 35 36

if(i==Imax) Tx2=0.0; else Area=2.0*pi*bCell(i)*dz; Dr=rCell(i+1)−rCell(i); Tx2=(Area/Dr)*(Permf*fvfo/Vis); end

37 38 39 40 41 42 43 44

a(i)= b(i)= c(i)= d(i)= end

45 46 47 48 49 50 51 52

end %Thomas Algorithm [PrsNewf] = thomas(Imax,a,b,c,d );

53 54 55

%Cumulative timestep Tacum=Tacum+Dt;

56 57

% Impresiones de tiempos y presiones

58 59

% Obtencion de presion de matriz

60 61 62 63 64 65

for i=2:Imax PSI=sigma*vCell(i)*Permm/(Vis*FVFo); OMEm=vCell(i)*Porom*fvfo*Ctm/Dt; PrsNewm(i)=(PSI*PrsNewf(i)+OMEm*PrsOldm(i))/(PSI+OMEm); end

66 67 68 69

% Actualizaci\'on de presiones para otro Dt PrsOldf(:,1)=PrsNewf(:,1); PrsOldm(:,1)=PrsNewm(:,1);

70 71 72

%Increse of timestep Dt=Dt*1.2;

64

73 74

end %while

65

6 S OLUCIÓN A NALÍTICA PARA EL M ODELO DE D OBLE P OROSIDAD U NA PERMEABILIDAD Considerando que un yacimiento naturalmente fracturado es conceptualizado con el modelo de doble porosidad - una permeabilidad, 2φ−1k, para un flujo monofásico, en un yacimiento infinito, y gasto de aceite constante, Warren and Root desarrollaron una solución analítica al modelo de doble porosidad, una permeabilidad, esto es, ecs. 5.6, 1 ∂ r ∂r

·

kf µB

r

∂p f ∂r

¸

µ ¶ ∂ φf ∂t B µ ¶ ∂ φm = ∂t B

− q˜m f = +q˜m f

(6.1)

y sus condiciones iniciales y de frontera: p f (r, t = 0) = p m (r, t = 0) = p 0 · ¸ ∂p f qB µ r =− ∂r r w 2πk f h

(6.2)

lim p f (r, t ) = lim p m (r, t ) = p i ni

r →∞

r →∞

Utilizando las siguientes variables adimensionales, 2πk f h £

PD = tD = rD =

qB µ

¤

kf t 2 (φm c t ,m + φ f c t , f )µr w

r rw

λ=σ ω=

p i − p r,t

(6.3)

km 2 r kf w φ f ct , f

φ f c t , f + φm c t ,m

Las ec. 6.1 es, · ¸ ∂p f D ∂P f D 1 ∂ ∂P mD rD =ω + (1 − ω) r D ∂r D ∂r D ∂t D ∂t D ∂P mD λ(P f D − P mD ) = −(1 − ω) ∂t D

(6.4)

66

Las condiciones iniciales y de frontera 6.2 se convierten en, P f D (r D , 0) = P mD (r D , 0) = 0 · ¸ ∂P f D = −1 ∂r D r D =1

(6.5)

lim p f D (r D , t D ) = lim p mD (r D , t D ) = 0

r D →∞

r D →∞

Usando la Transformada de Laplcae, Warren and Root solucionaron las ec. 6.4 junto con sus C.I. y C.F. ecs. 6.5, esto es, P f D (r D = 1, t D ) = =

µ · µ ¶ ¶¸ λt D 1 λt D ln t D + 0.80907 + E i − − − Ei − − +s 2 ω(1 − ω) (1 − ω)

(6.6)

Se observa en la ec. 6.6 que la respuesta del transiente de presión de un pozo en un yacimiento naturalmente fracturado esta caracterizado, básicamente, por dos parámetros, esto es, λ y ω, los cuales realcionan las propiedades de la matriz y la fractura. El parámetro ω, es la relación de almacenamiento de las fracturas conrespecto al total del sistema. El parámetro λ refleja el contraste de permabilidades entre matriz y fracturas, se le conoce como parametro de fluj interposo. Analizando la ec. 6.6 para tiempos cortos, las E i pueden ser sustituidas por su aproximación logarítmica, esto es, ec.3.13, reduciéndose a, · ¸ tD 1 ln + 0.80907 + s P f D (r D = 1, t D ) = 2 ω

(6.7)

Con la definición de t D dada por la ec. 6.3. Para tiempo largos, en la ec. 6.6 las E i llegan a ser cero, se conviertiendose en,

P f D (r D = 1, t D ) =

1 [ln t D + 0.80907] + s 2

(6.8)

Las ecs. 6.7 y 6.8, representan dos lines rectas en una gráfica semilog de t vs p w f , con una separación dada por,

∆Sep =

· ¸ 1 1 2 ω

(6.9)

Cuando el medio continuo, la fractura, primer recta semilog, ha dejado de manifestarse, entonces, entra en acción el medio discontinuo, esto es, la matriz, alimentando a las fracturas con fluido y soportando la presión, o sea, la segunda recta semolig. A este lapso de tiempo se

67

le conoce como periódo de transición del modelo de doble porosidad. El punto d einflexión de esta zona de transición esta dado por, 1 1−ω ln 1.26 2 λ

Pf D =

(6.10)

Un Yacimiento Naturalmente Fracturado puede ser caracterizado mediante el uso de la función derivada, esto es, tomando la derivada de la ec. 6.6, se tiene los siguiente

tD

dPf D d tD

· µ ¶ µ ¶¸ 1 λt D λt D = 1 + exp − − − exp − − 2 ω(1 − ω) (1 − ω)

(6.11)

A tiempos cortos y tiempos largos, la funcón derivada, para un comportamiento de doble porosidad, toma el valor de 0.5, esto es,

tD

dPf D d tD

=

1 2

(6.12)

Es posible, a partir de la ec. 6.11 obtener el punto de inflexión de la función dada por la ec. 6.6, de la siguiente manera,

tD tD

dPf D d tD dPf D d tD

=0 · µ ¶ µ ¶¸ 1 λt D λt D = 1 + exp − − − exp − − =0 2 ω(1 − ω) (1 − ω)

(6.13)

Esto arroga el siguiente t D mínimo,

t D,mi n =

ω 1 ln λ ω

(6.14)

Con el siguiente punto en la función derivada, ·

tD

dPf D d tD

¸

= mi n

i 1 ω 1h 1 + ω 1−ω − ω 1−ω 2

(6.15)

Por lo tanto, graficando estos dos puntos, ec. 6.14 y ec.6.15, mostraran el punto mínimo en la función derivada. La Figura 6.1 muestra el transiente de presión para tiempos cortos y tiempos largos en un yacimientos naturalmente fracturado, conceptualizado con el modelo de doble porosidad; asi como también el periódo de transición entre las dos rectas semilog y la función derivada. Observando la Figura 6.1 podemos hacer las siguientes afirmaciones:

68

Modelo de Doble Porosidad 10 Tiempor Cortos W-R Tiempor Largos

8

6

tD

4

Inflexión 2

0

-2

-4 -4 10

-2

0

10

2

10

10

4

6

10

10

PD Función Derivada-Modelo de Doble Porosidad 1 Derivada 0.9 0.8

tD*d[PD/tD]

0.7 0.6 0.5 0.4 0.3 0.2

Mínimo 0.1 0 -4 10

-2

10

0

2

10

10

4

10

6

10

tD

Figure 6.1: (a) Modelo de Warren & Root. (b) Función Derivada. ω = 0.01, λ = 1e-3.

69

• Durante tiempos cortos, cuando la ec. 6.7 es válida, esta, presentará una linea recta en una gráfica semilog, mostrando el efecto de las fracturas exclusivamente. • Durante tiempos largos, cuando la ec. 6.8 es válida, esta, presentará una linea recta en una gráfica semilog, mostrando le influencia del sistema completo, esto es, fractura mas matriz. • Durante tiempos intermedios, esto es, durante el periódo de transición, la capacidad de aporte de la fractura disminuirá y causará que la matriz aporte fluido y recarge, alimente, al medio fracturado, hasta que el sistema fractura-matriz se refleje totalmente como un sistema total. • Durante tiempos cortos de producción, en general, la respuesta es eclipsada por el alamacenamiento del pozo. A continuación se presenta el algoritmo de Warren and Root. Este, esta en variables adimiensionales, para usar variables reales, solo hay que usar las definiciones y despegar las variables reales, es decir, p w f y t 1 2

%% Radial Flow Single Phase. Warren and Root Model %

3 4 5 6

clc; clear all; close all;

7 8 9 10 11 12

% Some definition to be used with real variable % %Omega = Porof*Ctf /(Porof*Ctf + Porom*Ctm); %Lambda = sigma*km*(rw^2)/kf; %Aux1 = kf/( (Porof*Ctf+Porom*Ctm)*Vis*(rw^2)); %tD

13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

Omega = 0.01; Lambda=1e−3; Cont=0; for i=−3:5 for j=1:9 Cont=Cont+1; % tD generation tD = j*10.0D0^i; TT(Cont)=tD; Aux2 = Lambda*tD/(Omega*(1.0D0 − Omega)); Auxa = −Aux2; if(Aux2 < 100) [Expo1]= IntExp(Aux2); Expo1 = −Expo1; else Expo1=0.0D0; end

32

70

33 34 35 36 37 38 39 40

Aux3 = Lambda*tD/(1.0D0−Omega); Auxb = −Aux3; if(Aux3 < 100) [Expo2]= IntExp(Aux3); Expo2 = −Expo2; else Expo2 = 0.0; end

41 42 43

%Early Times pwD1 (Cont)=0.5*(log(tD)+0.80907 +log(1/Omega)) ;

44 45 46

%Presion Adimensional pwD2 (Cont)=0.5*(log(tD)+0.80907 +Expo1−Expo2) ;

47 48 49

%Late Times pwD3 (Cont)=0.5*(log(tD)+0.80907) ;

50 51 52 53 54

%Derivative Function arg2=−(Lambda*tD)/(1−Omega); arg1=arg2/Omega; FDer(Cont)=0.5*(1+exp(arg1)−exp(arg2));

55 56 57 58 59 60 61

%For Real variable %Presion Real (q in negative (for simulation). So, Pini+ %Pwf(Cont)= Pini+(q*B*Vis*pwD(Cont))/(2*pi*kf*h); %A1=kf*h/(2*pi*q*B*Vis); %DPwD=0.5*(1.0 +exp(Auxa)−exp(Auxb)); %Dpwd=DPwD/A1;

62 63 64

end end

65 66 67

Dist=0.5*log(1/Omega); PfDD=0.5*log(1.26*(1−Omega)/Lambda)

68 69 70

%Minimo en PwD vs tD PfD= 0.5*log(1.26*(1−Omega)/Lambda);

71 72 73 74

%Minimos en Funciø sn Derivada tdMin=(Omega/Lambda)*log(1/Omega) FDerMin=0.5*(1+Omega^(1/(1−Omega))−Omega^(Omega/(1−Omega)));

75 76 77 78 79 80 81 82 83

figure(1); semilogx(TT,pwD1,'−k', TT,pwD2,'bo',TT,pwD3,'−r','LineWidth',2) grid on title('Modelo de Doble Porosidad') ylabel('tD') xlabel('PD') legend('Tiempor Cortos','W−R','Tiempor Largos' ) hold on

84 85 86

figure(1); semilogx(tdMin,PfD,'ro','LineWidth',4)

71

87 88 89 90 91 92 93 94 95 96 97 98

figure(2); semilogx(TT,FDer,'−b','LineWidth',2) grid on title('Funciø sn Derivada−Modelo de Doble Porosidad') ylabel('tD*d[PD/tD]') xlabel('tD ') legend('Derivada' ) hold on

99 100 101

figure(2); semilogx(tdMin,FDerMin,'o','LineWidth',4)

72

A CKNOWLEDGMENTS Para todos aquellos que no le han vendido su alma al diablo... References [?].[?].

73