Wag.pdf

Proyecto: "Inyección alternada de agua y gas (WAG) como sistema de recuperación mejorada" CIMAT Octubre 2013

1

Introducción

Este proyecto fue dividido en cuatro etapas: ˆ En la primera etapa se hizo una revisión de diferentes ecuaciones que han sido propuestas para modelar un flujo monofásico que se difunde a través de un medio poroso que, por sus propiedades, producen una difusión anómala. ˆ En la segunda etapa se eligó la ecuación de Tarasov para modelar la difusión anómala en un dominio que puede tener uno o varios pozos. Se aplicó el método de los elementos finitos basados en volúmenes de control (MEFVC) para resolver numéricamente la ecuación diferencial asociada. ˆ En la tercera etapa se planteó el modelo transitorio de difusión usando el modelo de Tarasov. MEFVC fue aplicado para resolver numéricamente el problema. ˆ En la última etapa se formuló un problema para poder estimar los parámetros fractales de la ecuación transitoria, a partir de un conjunto de mediciones de la caída de presión registrada en un punto particular del dominio del problema, tomando varias mediciones en diferentes instantes de tiempo.

2

Difusión para pozos con geometría radial

Existen varias propuestas para modelar el comportamiento subdifusivo de un fluido en un medio poroso. En algunos enfoques se usan derivadas fracciona1

rias. Para el caso de pozos con geometría radial se estudiaron varias ecuaciones. Una de ellas fue la ecuación de Metzler-Glökle-Nonnenmacher (MGN) [Camacho-Velázquez 08]: ‚ Œ ∂p ∂α pD 1 ∂ D β (1) r = r β+θ ∂rD D ∂rD ∂t α D

D

donde θ es el índice de conductividad y β = dmƒ −θ−1, siendo dmƒ la dimensión fractal de masa, y p la presión en un punto del dominio. La solución de esta ecuación se puede expresar en forma analítica, en términos de la transformada de Laplace de la presión. Por ejemplo, para el problema pt =

1

∂ €

r D−1 ∂r p(r, 0) = 0, ∂p = −1, ∂r r=r ∂p = 0, ∂r r=re

Š r d pr ,

r ≤ r ≤ re , t > 0.

(2)

r ≤ r ≤ re .

(3)

t > 0.

(4)

t > 0.

(5)

Se puede ver que la solución de esta ecuación se puede expresar en términos de la transformada de Laplace P¯ de p. Para el caso de una reserva acotada con producción bajo el régimen de tasa constante (condiciones de frontera (4) y (5)), la solución está dada por  α/ 2 Θ   α/ 2 Θ  s r s r eD eD α/ 2 Θ K (s r / Θ) + K ν (sα/ 2 r Θ / Θ)  νΘ ν ν−1 ν−1 Θ Θ r (6) P¯ = α/ 2  α/ 2 Θ   α/ 2 Θ  s r s r s eD eD Kν−1 (sα/ 2 / Θ)ν−1 − ν−1 (sα/ 2 / Θ)Kν−1 Θ Θ donde Kν , ν son las funciones modificadas de Bessel de orden ν [Watson 96], y Θ = (D − d + 1)/ 2, 1−d ν = . D−d+1 El cálculo de la transformada inversa requiere el uso de un método numérico [Davies 79, de Hoog 82]. En la Fig. 1 se muestran el comportamiento de la caída de presión, del flujo y del flujo acumulativo medidos en la pared de pozo r = 1 que tiene un radio exterior re = 500, al variar el parámetro α del modelo, que tiene que ver con el orden de la derivada temporal fraccionaria. 2

104

10−1

q(1, t)

Pressure drop

103

10

100

α = 1.00 α = 0.90 α = 0.80 α = 0.70

2

10−2 10−3

101

10−4

100

10−5

101

102

103

104

105 106 Time

107

108

109

α = 1.00 α = 0.90 α = 0.80 α = 0.70 101

102

103

104

105 106 Time

107

108

109

108 107

Q(t)

106 105 104 103

α = 1.00 α = 0.90 α = 0.80 α = 0.70

102 101 101

102

103

104

105 106 Time

107

108

109

Fig. 1: Cantidades calculadas a partir del modelo MGN, cuando se varía el parámetro α del modelo. En el intervalo de tiempo (0, 109 ) se muestra el comportamiento de la caída de presión, el flujo q(1, t) en la pared del pozo y el flujo acumulativo Q(t).

3

Modelo de difusión transitorio

El modelo de Tarasov introduce integrales fraccionarias sobre un medio continuo que permiten incorporar las propiedades fractales de un medio poroso por medio de ciertos parámetros de la ecuación que se deriva de las leyes de conservación. Estos parámetros corresponden a la dimensión fractal de las regiones bidimensionales y la dimensión fractal de sus fronteras, y al variarlos es posible producir un efecto subdifusivo.

3.1

El modelo de Tarasov bidimensional

En el modelo de Tarasov [Tarasov 05] se reemplaza el medio fractal por un medio continuo con masa de dimensión fractal, la cual es descrita por integrales fraccionarias que toman en cuenta la fractalidad del medio. La masa del 3

medio obedece una relación de potencia, M(R) = kRD

(D < 2),

donde M es la masa del medio, R es el radio de un círculo que contiene a la masa, y D es la dimensión fractal de la masa. Sea W una región en R2 cuyo punto medio es . Denotamos por A(W ) al área de W y MD (W ) a su masa. Se define la integral fraccional en R2 como Z ( D ρ)(r 0 ) =

ρ(r) dAD .

(7)

W

donde dAD = c3 (D, r, r0 )d1 d2 y c3 (D, r, r0 ) =

22−D (D/ 2)

kr − r 0 k

D−2

2

kr − r 0 k =

,

2 X

(k − k0 )2 .

k=1

El punto r 0 ∈ W es el punto inicial de la integral fractional. Entonces MD (W) = ( D ρ)(r 0 ) =

22−D

Z

(D/ 2)

ρ(r)kr − r 0 kD−2 d1 d2 .

(8)

W

Para obtener la ecuación de continuidad para el medio fractal, se considera una región W cuya masa tiene dimensión D y con frontera ∂W de dimensión d. Entonces Z Z Z d dMD (W) = q dAD i.e. ρ(R, t) dAD = q dAD (9) dt dt 



W

donde R = r − r 0 . Suponemos que la región W sePmueve con el medio. El campo de velocidad es denotado por  = (R, t) = k k ek . Para una función α dada, la derivada total temporal se define como Z Z Z d ∂α α dAD = dAD + α · n dd . (10) dt ∂t W

donde n =

P

k



∂W

nk ek es el vector normal a ∂W, y dd = c2 (d, R) d1 ,

c2 (d, R) = 4

21−d (d)

kRkd−1 .

Entonces, por el teorema de la divergencia Z

α · n dd =

∂W

Z

αc2 (d, R) · n d1 =

∂W

=

Z

Z

∇ · [αc2 (d, R)] dA2

W

c−1 (D, R)∇ · [αc2 (d, R)] dAD 3

W

Sustituyendo en (9) se tiene Œ Z Z ‚ d ∂α ∇ · [αc2 (d, R)] + α dAD = dAD dt ∂t c3 (D, R)

(11)



W

Por (9), si tomamos α = ρ(R, t), (11) debe ser igual a la integral de q sobre W, para todo W centrado en r 0 , por lo que ∂ρ ∂t

+

∇ · [ρc2 (d, R)] c3 (D, R)

=q

Esto es, ∂ρ ∂t

+

( D ) 2D−d−1 2 (d)

kRk2−D ∇ · (kRkd−1 ρ ) = q,

(12)

donde 1 < d < D < 2. Por otra parte, se define la compresibilidad del fluido, cƒ , como 1 ∂ρ , cƒ = ρ ∂p T en una temperatura fija T. Integrando la expresión anterior con respecto a p, se obtiene ρ = ρ0 ecƒ (p−p ) , 0

donde ρ0 es la densidad en la presión de referencia p0 . Con los primeros dos términos de la expansión de Taylor de la función anterior, tenemos la siguiente aproximación ρ ≈ ρ0 [1 + cƒ (p − p0 )]. Entonces

5

∂ρ

≈ cƒ ρ0

∂p

. ∂t ∂t Adicionalmente, si suponemos que el flujo cumple la ley de Darcy,

(13)

 = −κ∇p, al sustituir esta expresión y (13) en (12), obtenemos la ecuación transitoria de Tarasov cƒ ρ

0

∂p ∂t

−

) 2D−d−1 ( D 2 (d)

k −  0 k2−D ∇ · (k −  0 kd−1 ρ κ∇p) = q,

(14)

donde  0 ∈ Ω está dado, ƒ (, t) es un término fuente y d < D ≤ 2. En principio, al variar los parámetros D y d de la ecuación anterior podemos controlar la forma en que la difusión ocurre.

3.2

Solución numérica del modelo de Tarasov

Para resolver la ecuación del modelo de Tarasov se utilizó el método de los elementos finitos basados en volúmenes de control (MEFVC) [Chen 07]. Para obtener una idea del error en los cálculos, se resolvío el modelo de Tarasov para el caso de un pozo acotado con geometría radial, para poder comparar la solución numérica con la que se obtiene a partir de (6). En este caso, en la ecuación 2D (14), tenemos que q = 0 y  0 = 0. por lo que la podemos reescribir como g()

∂p ∂t

− ∇ · ( ()∇p ) = ƒ (, t)

(15)

con g() =

cƒ ρ0 (d) ( D ) 2D−d−1 2

k −  0 kD−2 ,

() = k −  0 kd−1 ρ κ, (d) ƒ (, t) = k −  0 kD−2 q(, t). D D−d−1 ( 2 ) 2

El dominio Ω se discretiza en triángulos, como se muestra en la Fig. 2(a). A cada nodo de la discretización se le puede asociar un volumen de control, como se muestra en la Fig. 2(b).

6

4

4

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

−2

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

● ● ●

●

●

●

●

● ●

●

(a) Discretización

●

●

●

●

●

●

● ●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

−4

●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

●

●

●

●

−2

●

●

● ●

−4

●

●

●

●

●

●

● ●

●

●

●

● ●

●

● ●

●

● ●

●

●

●

●

● ●

●

●

●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

●

●

●

● ●

● ●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

● ●

●

●

0

0

●

−4

nodes[, 2]

●

● ●

●

●

●

●

●

●

● ●

2

● ●

●

●

●

●

−2

●

● ●

●

●

●

● ●

●

nodes[, 2]

2

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

● ●

●

●

● ●

●

●

●

●

●

●

●

●

●

(b) Volúmenes de control

0

2

4

−4

−2

0

2

4

Fig. 2: Ejemplo de la triangulación del dominio y los volúmenes de control asociados nodes[, 1] a la discretización. nodes[, 1]

La ecuación diferencial se integra en cada volumen de control. Para cada volumen se obtiene una ecuación en la que aparecen como variables los valores de la presión en los nodos de la discretización. A partir de estas ecuaciones y usando el esquema de Crank-Nicholson para resolver el problema, obtenemos el sistema de ecuaciones de la forma Œ ‚ Œ ‚ Δt Δt Δt A pk+1 = (Fk+1 + Fk + 2b) + G − A pk , (16) G+ 2 2 2 donde pk es la solución en el instante de tiempo k-ésimo, el término b tiene que ver con las condiciones de frontera del problema, G = diag(G1 , ..., Gn ), Fk = (F1k , ..., Fnk )> , con Z

G ≈

g d,

Fk

V

Z

≈

ƒ (, tk ) d,

V

y las entradas de la matriz A = [Aj ] tienen que ver con las trasmisibilidades Tj asociadas a los nodos  y j, calculadas sobre las aristas de los volúmenes de control: ¨ −Tj ,  6= j, Aj = P k∈S Tk ,  = j, S índices de los nodos vecinos al nodo . 7

Las consideraciones que tomamos para construir y resolver (16) son las siguientes: ˆ Usar cómputo paralelo para calcular las entradas de las matrices A, G y de los vectores Fk+1 y Fk . ˆ Resolver el sistema lineal de ecuaciones (16) usando el algoritmo de gradiente biconjugado [Press 07], que es un algoritmo iterativo y que la convergencia a la solución depende de los valores que se dan al vector p al inicio del algoritmo. ˆ Usar la solución pk del instante k para inicializar el algoritmo de solución, porque si el incremento en el tiempo Δt es pequeño, el algoritmo hará pocas iteraciones para converger a pk+1 .

Para el caso del problema transitorio 2D en un dominio anular, modelado por la ecuación (??), en la Fig. 3 se muestran los perfiles de la caída de presión registrada en un punto de la frontera interior, para diferentes valores de los parámetros D y d, así como los valores de la caída de presión en todo el dominio para dos instantes de tiempo. Por ejemplo, si fijamos r = 1, re = 10 y usamos una discretización de Ω que tiene las siguientes características: 29800 59004 88804

nodos elementos triangulares aristas

Así, en cada paso de tiempo hay que resolver un sistema de ecuaciones de 29800 variables. En la Fig. 4 se muestran las solución numérica y la analítica para D = 2.00 y d = 1.00, que es el caso en que no hay difusión anómala. Para la resolución de la gráfica, las dos curvas quedan superpuestas. Para resaltar las diferencias, calculamos el error relativo entre las dos soluciones, |pk − p(  , tk )| |p(  , tk )|

,

(17)

donde pk es el valor la solución numérica en el nodo   = (1, 0) en el instante de tiempo tk , y p(  , tk ) el valor de la solución analítica correspondiente. Los errores son del orden de 10−4 en todo el intervalo de tiempo. De este modo, la solución numérica reproduce adecuadamente los valores de la solución analítica. 8

1.0

1.2

D=2.00, D=2.00, D=2.00, D=2.00,

0.8

d=1.00 d=1.00 d=1.00 d=1.00

d=1.00 d=1.15 d=1.30 d=1.45

0.6 0.0

0.0

0.2

0.2

0.4

Presion

0.4

0.6

Presion

0.8

1.0

D=2.00, D=1.85, D=1.70, D=1.55,

0.0

0.5

1.0

1.5

2.0

0.0

0.5

Tiempo

1.0

1.5

2.0

Tiempo

t = 2.0

t = 20.0

Fig. 3: Ejemplos de la solución del problema transitorio 2D.

0.00030 0.00025 0.00020 0.00015

5

Error relativo

10

Solucion analitica Solucion numerica

2

Caida de presion

20

Otro ejemplo se muestra la Fig. 5, en donde los parámetros seleccionados fueron D = 1.85 y d = 0.75. En este caso se resuelve el problema en el

10

20

50

100

200

500

1000

0

Tiempo

200

400

600

800

1000

Tiempo

Fig. 4: Soluciones en (1, 0) para D = 2.0 y d = 1.0, y gráfica del error relativo.

9

10 2

0.00010

0.00020

20

Error relativo

50

Solucion analitica Solucion numerica

5

Caida de presion

100 200

intervalo de tiempo [0, 104 ] usando un tamaño de paso Δt = 0.5. Aunque incrementamos el intervalo de tiempo y el tamaño de paso, el error relativo se mantiene dentro del rango anterior. Los tiempos registrados para obtener los resultados usando cuatro procesadores se muestran en la Tabla 1.

20

50 100

500

2000 5000

0

Tiempo

2000

4000

6000

8000

10000

Tiempo

Fig. 5: Soluciones en (1, 0) para D = 1.85 y d = 0.75, y gráfica del error relativo entre ellas.

Tiempo para construir la matriz del sistema Tiempo para construir los vectores independientes Tiempo para resolver los sistemas de ecuaciones Tiempo total registrado

0.008 313.673 4051.214 4466.225

segundos segundos segundos segundos

Tabla 1: Tiempos registrados usando 4 procesadores.

Además de medir el error que se comete en un punto de la frontera interior, podemos el error en todos los nodos   de la discretización del dominio para un instante tiempo dado. Para ello fijamos el valor del tiempo tk = 10000 y para cada r = k  k calculamos el valor de p(r , tk ) de la solución analítica a partir de la ecuación de la transformada inversa (6), y lo comparamos con el valor pk de la solución numérica en t = tk . En la Fig. 6 se muestra la solución en tk = 10000 y el error relativo en todo el dominio. Se puede ver que los errores más grándes quedan distribuidos en anillos y siguen siendo del mismo orden que el error medido en el punto (0, 1), por lo que podemos esperar que el error medido en ese punto es representativo de todo el dominio. En otra prueba, cambiamos el dominio de modo que el radio exterior re = 150. Generamos una discretización con 60000 elementos y 30300 nodos.

10

Solución

Error relativo

Fig. 6: Solución en el instante t = 10000 para D = 1.85 y d = 0.75, y error relativo medido en todo el dominio, para re = 10.

Queremos ver si esto afecta en el cálculo de la solución del problema transitorio. Lo resolvemos en el intervalo de tiempo [0, 10000] tomando un tamaño de paso Δt = 0.5. La comparación entre las soluciones numérica y la analítica en el punto (1, 0), así como el error relativo entre ellas se muestran en la Fig. 7. El error es mayor al del caso anterior, pero aun es pequeño. Sin embargo, si medimos el error relativo en todo el dominio en el instante t = 10000 obtenemos la g´rafica que se muestra en la Fig. 8. Vemos que los errores más grandes ocurren cerca de la frontera exterior, donde los triángulos son más grandes. Por tanto, hay que refinar más la discretización para mejorar la aproximación en esa zona, aunque cerca de la frontera interior el error es pequeño.

4

Determinación de los parámetros fractales del modelo

Si queremos estimar los parámetros d y D del modelo a partir de un conjunto de mediciones {p1 , p2 , ..., pm } de la caída de presión en un punto dado del dominio, registradas en los instantes de tiempo t1 < t2 < ... < tm , tenemos que proponer valores para d y D, y para una pareja de éstos, hay que resolver la ecuación diferencial para obtener los valores de la solución en los instantes t y compararlos con los valores p . Hay dos inconvenientes: ˆ Si se tienen que hacer varias pruebas de los valores de d y D, hay que resolver muchas veces la ecuación diferencial, y esto puede ser computacionalmente costoso.

11

1e−03 Error relativo

8e−04

7 6 5 4

4e−04

3

Caida de presion

6e−04

8 9

Solucion analitica Solucion numerica

20

50 100

500

2000 5000

0

Tiempo

2000

4000

6000

8000

10000

Tiempo

Fig. 7: Soluciones en (1, 0) para D = 1.85 y d = 0.75, y gráfica del error relativo entre ellas, para re = 150.

Solución

Error relativo

Fig. 8: Solución en el instante t = 10000 para D = 1.85 y d = 0.75, y error relativo medido en todo el dominio, para re = 150.

ˆ Si tm es grande, tenemos que discretizar el intervalo [0, tm ] finamente, para que cuando resolvamos numéricamente la ecuación diferencial, el error acumulado en cada paso de tiempo sea pequeño. Esto también hace que se eleve el costo computacional para resolver la ecuación diferencial para cada pareja (d, D).

Puesto que usualmente la presión que se puede medir en un pozo es la pared del mismo, sólo deberíamos considerar que esa es la única información disponible para hacer la estimación de los parámetros. 12

Para simplificar el problema aún más, podemos considerar la ecuación 1D, que es más sencilla y podemos tener la expresión de la solución analítica en términos de la transformada de Laplace de la presión, como se explica en la siguiente sección. La ventaja de esto es que podemos simular estos valores más rápido y podemos enfocarnos en el método de estimación de parámetros, el cual no esta asociado a la forma particular en que se generaron los datos.

4.1

Solución del problema para geometría radial

Consideremos el caso particular en que dominio tiene forma de anillo de radio interior r y radio exterior re . El problema transitorio 1D a resolver está dado por las ecuaciones (2), (4) y (5), cuya solución está dada por (6).

4.2

Planteamiento del problema

Supongamos que p representa el valor de la caída de presión medida en la pared de pozo r en el instante de tiempo t . Queremos determinar el valor de los parámetros fractales de la ecuación diferencia a partir de un conjunto formado m registros de la caída de presión, { ( t , p ) }m , =1

t1 < t2 < ... < tm .

50

x x

20

x

10

x x x x x x x x x x x x x x x

x

x

2

5

Caída de presión

100

200

En la Fig. 9 se muestra un ejemplo de este conjunto. A partir de esté queremos estimar los valores D y d, de modo que el valor p(t ; d, D) de la

1e+01

1e+03

1e+05 Tiempo

Fig. 9: Las marcas indican las muestras ( t , p ) que tenemos para resolver el problema.

13

solución de la ecuación en el instante t se aproxime a p . Una forma de medir las diferencias en todo el conjunto de valores es mediante el error cuadrático de las diferencias, m X

[p(t ; d, D) − p ]2 ,

=1

100

sin embargo, como vemos en la Fig. 9, los valores de la caída de la presión aumentan conforme aumenta el tiempo, de modo que hay un sesgo y puede ocurrir que se prefieran valores D y d que reduzcan las diferencias [p(t ; d, D) − p ]2 para  grandes, que son los errores más signitificativos, pero no hagan un buen ajuste para valores pequeños de . En algunos casos los últimos valores p de la secuencia pueden ser los que más pueden ayudar en la determinación de los valores D y d. Por ejemplo, supongamos que los verdaderos valores son D = 1.65 y d = 0.88. En la gráfica de la izquierda de la Fig. 10 se ve que si tenemos correcto el valor d y sólo tenemos muestras para t < 104 hay varias soluciones de la ecuación que pasan cerca de la curva que queremos reconstruir, por lo que es más dificil determinar el valor de D. En cambio, si tenemos muestras para t > 105 , las curvas se diferencian mejor y podemos dar una mejor estimación de D. Por otro lado, en la gráfica de la derecha de la Fig. 10 se ve que si tuvieramos correcto el valor de D y sólo tenemos que deteminar d, es mejor tener muestras para t < 105 , ya que para para ese intervalo las curvas se diferencian mejor.

D=1.65, D=1.65, D=1.65, D=1.65, D=1.65,

20

d=0.80 d=0.84 d=0.88 d=0.92 d=0.96

5

5

10

20

Caída de presión

50

d=0.88 d=0.88 d=0.88 d=0.88 d=0.88

10

Caída de presión

50

D=1.59, D=1.62, D=1.65, D=1.68, D=1.71,

1e+02

1e+03

1e+04

1e+05

1e+06

Tiempo

1e+02

1e+03

1e+04

1e+05

1e+06

Tiempo

Fig. 10: Curvas de la caída de presión obtenidas al variar uno de los parámetros, D o d.

14

10

20

d=0.85 d=0.91 d=0.88 d=0.85 d=0.91

5

Caída de presión

50

D=1.62, D=1.62, D=1.65, D=1.68, D=1.68,

1e+02

1e+03

1e+04

1e+05

1e+06

Tiempo

Fig. 11: Curvas de la caída de presión obtenidas al variar los parámetros D y d.

En general cuando variamos los dos parámetros D y d se obtienen curvas como las que muestran en la Fig. 11. En este caso no es claro cuales muestras son más importantes. Así, debemos normalizar las diferencias para que todas contribuyan por igual en la determinación de D y d. Dados D y d, calculamos la solución p(t ; d, D) de la ecuación en el instante t y medimos el error global mediante la suma de los cuadrados de los errores relativos, E(d, D) =

m – p(t ; d, D) − p ™2 X  

p

=1

.

(18)

De este modo, el problema que queremos resolver es (D∗ , d∗ ) sujeto a

= rg min E(d, D) = d,D

m – p(t ; d, D) − p ™2 X   =1

p(t; d,D) solución de (??) , 0 < d < D ≤ 2, 1 < D.

p

(19)

Un ejemplo de la gráfica de la función de error E(d, D) se muestra en la Fig. 12. Para generarla se tomaron m = 100 distribuidas en el intervalo [50, 106 ]. Hay un óptimo global, que se identificar con más facilidad en la Fig. 13 cerca

15

del punto d = 0.88 y D = 1.65, que son los verdaderos valores de los parámetros fractales que se usaron para generar las observaciones.

12 10

12

8

10

6

8

4

6

2

4

02 1.9 1.8 D

1.7 1.6 1.50.6

0.8

1

1.2 d

1.4

1.5 1.6 1.7 D 1.8 1.9

2

1.6

0 1.6

1.4

1.2 d

1

0.8

0.6 2

Fig. 12: Gráfica de la función de error.

2

10

8

1.8

6

1.7

4

1.6

2

D

1.9

1.5

0 0.6

0.8

1

1.2

1.4

1.6

d

Fig. 13: Curvas de nivel de la función de error. El óptimo se encuentra en D = 1.65 y d = 0.88.

En este ejemplo se puede ver que cerca del óptimo la superficie cae rápidamente a cero. Esto se aprecia mejor en la Fig. 14, con los perfiles de la gráfica de error, E(d, 1.65) y E(0.88, D). Cerca del óptimo la pendiente de la 16

curva cambia abrutamente. Es por esto que para resolver el problema (19) que están básados en derivadas pueden tener problemas de convergencia. E(d, 1.65)

E(0.88, D)

7

5

6 4 5 3 Error

Error

4

3

2

2 1 1

0

0 0.6

0.8

1

1.2

1.4

1.6

1.5

1.6

1.7

d

1.8

1.9

2

D

Fig. 14: Perfiles de la función de error.

Debido a que sólo son dos parámetros que hay que calcular, y que para no tener problemas con métodos de optimización que requieren la estimación de derivadas, aplicamos el método de Nelder-Mead [Press 07]. Este método tiene las siguientes características: ˆ Sólo requiere evaluaciones de la función a optimizar, sin que ésta tenga que cumplir alguna condición de suavidad. ˆ Puede aplicarse para optimizar funciones con cualquier número de variables, pero trabaja mejor en bajas dimensiones. ˆ No hay garantía de que converga al óptimo global de la función, de modo que depende de la inicialización. ˆ Hay estrategias para reinicializarlo de forma automática de modo que pueda hacer una mejor exploración del espacio de soluciones.

Para el caso en que estamos interesados, la función a optimizar no tiene varios mínimos locales y el espacio de soluciones está acotado, por lo que no es difícil inicializarlo. Para poder aplicarlo al problema (19) la versión que implementamos maneja las restricciones 0 < d < D ≤ 2.

17

4.3

Experimentos realizados

En esta sección desarrollamos varias pruebas para saber la dificultad que tiene encontrar una buena estimación de los parámetros del modelo. Para ello fijamos los valores de d y D, y generamos el conjunto de datos del problema { ( t , p ) }m , =1

t1 < t2 < ... < tm ,

Queremos ver como se comporta el algoritmo cuando ˆ se tienen varias combinaciones de los parámetros d y D, ˆ se cambia el número m de muestras, y ˆ se cambia el intervalo [t1 , tm ] en donde tenemos las muestras.

Las pruebas las dividimos en dos casos: en uno de ellos tenemos los datos sin ruido y en el otro agregamos una cierta cantidad de ruido.

4.4

Pruebas con los datos sin ruido

Tomamos un pozo con radio exterior re = 500 y las muestras las generamos dentro del intervalo de tiempo [50, 106 ]. Los vértices del triángulo inicialmente están en los puntos (1.2, 1.95),

(1.2, 1.75),

(1.45, 1.95).

En la Tabla 2 se muestra los verdaderos valores d∗ y D∗ con que se generaron los datos, que son los óptimos de la función en cada caso. También se muestra los valores dn y Dn a los que el algoritmo converge, el número n de iteraciones requeridas y el valor de la función de error E en esos puntos. En todos los casos, con menos de 40 iteraciones, el método converge muy cerca del óptimo. Los errores relativos están por debajos de 0.21%. Se puede decir no hay una región en el rectángulo 0.6 ≤ d∗ ≤ 1.4, 1.6 ≤ D∗ ≤ 1.9 que presente alguna dificultad particular que afecte la convergencia del método. Por esto, para la siguientes pruebas sólo tomamos una combinación de valores de d∗ y D∗ . Fijamos d∗ = 0.6, D∗ = 1.6 y el radio exterior re = 500. En la Tabla 3 se muestran los resultados de variar la cantidad de muestras y el intervalo de tiempo en que se encuentran distribuidas. Vemos que aunque las muestras estén distribuidas en [50, 106 ], o al principio de este intervalo o al final del mismo, el método siempre converge al mismo resultado, sin importar demasiado el número de muestras que se tomen. De este modo, se puede decir que el método es robusto a estas variaciones. 18

Tabla 2: Prueba con datos sin ruido y varios valores d∗ y D∗ . d∗ 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6

D∗ 1.9 1.9 1.9 1.9 1.9 1.9 1.9 1.9 1.9 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8 1.8

n 23 23 22 23 25 27 28 33 36 28 24 25 25 25 27 28 27 39

dn 1.4000 1.2999 1.2000 1.1000 0.9983 0.8999 0.7999 0.6999 0.5999 1.3999 1.3001 1.2000 1.0999 0.9984 0.9000 0.8000 0.7000 0.6000

Dn E(dn , Dn ) 1.9000 0.0009 1.8999 0.0007 1.8999 0.0009 1.9000 0.0013 1.9019 0.0374 1.9001 0.0012 1.9000 0.0006 1.9000 0.0004 1.9000 0.0005 1.7999 0.0014 1.7999 0.0011 1.8000 0.0008 1.8000 0.0012 1.8014 0.0376 1.7999 0.0010 1.7999 0.0011 1.7999 0.0009 1.7999 0.0004

d∗ 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6

D∗ 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.7 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6

n 27 26 28 29 33 28 29 32 27 32 29 28 26 32 24 31 27 26

Tabla 3: Prueba con datos sin ruido, variando el número en donde se encuentran. Intervalo: [50, 106 ] d∗ D∗ m n dn Dn 0.6 1.6 25 26 0.59994 1.60003 0.6 1.6 50 26 0.60004 1.59995 0.6 1.6 100 26 0.60004 1.59995 0.6 1.6 200 26 0.60004 1.59995 Intervalo: [50, 103 ] ∗ ∗ d D m n dn Dn 0.6 1.6 25 32 0.59978 1.60049 0.6 1.6 50 34 0.60009 1.59983 0.6 1.6 100 32 0.60010 1.59979 0.6 1.6 200 30 0.60008 1.59985 Intervalo: [105 , 106 ] ∗ ∗ d D m n dn Dn 0.6 1.6 25 31 0.60007 1.59995 0.6 1.6 50 31 0.60007 1.59995 0.6 1.6 100 32 0.60000 1.59999 0.6 1.6 200 32 0.60000 1.59999

19

dn 1.3999 1.2999 1.1999 1.0999 0.9984 0.8999 0.8000 0.7000 0.5999 1.4000 1.3000 1.2000 1.0999 1.0021 0.8999 0.8000 0.6999 0.6000

Dn E(dn , Dn ) 1.7000 0.0007 1.7000 0.0011 1.6999 0.0008 1.7000 0.0010 1.7011 0.0376 1.7000 0.0011 1.7000 0.0016 1.6999 0.0005 1.7001 0.0011 1.6000 0.0010 1.5999 0.0011 1.5999 0.0016 1.6000 0.0010 1.5996 0.0376 1.6000 0.0011 1.5999 0.0011 1.6000 0.0008 1.5999 0.0007

de muestras y el intervalo

E(dn , Dn ) 0.0004 0.0005 0.0007 0.0010 E(dn , Dn ) 0.0002 0.0002 0.0004 0.0005 E(dn , Dn ) 0.0004 0.0006 0.0005 0.0008

4.5

Pruebas con los datos con ruido

Repetimos algunas de las pruebas anteriores agregando ruido a los datos. Si p es el valor de la presión en el instante  sin ruido, perturbamos estos valor del problema de modo que el valor perturbado pε cumpla con p − pε  <ε p

para un ε dado. De este modo, la cantidad en que se perturban los datos se va escalando de acuerdo a la magnitud de los datos sin ruido. Por ejemplo, en la Fig. 15 se muestran el conjunto de datos {p } y {pε } cuando ε = 0.4. Reemplazamos en la función de error los valores p por pε , de modo que E(d, D) =

m – p(t ; d, D) − pε ™2 X  

pε

=1

.

Al resolver el problema (19) usando los datos perturbados es de esperar que entre más grande ser ε peor será la estimación de los valores d∗ y D∗ .

Presion

102

101

105

106 Tiempo

Fig. 15: Ejemplo del conjunto de datos del problema. Los que no tienen ruido están marcado con ’+’ y los que datos perturbados con 40% de ruido se indican con ’o’.

20

Tabla 4: Prueba con datos con ruido, variando el número de muestras y el intervalo en donde se encuentran con ε = 0.2. Intervalo: [50, 106 ] ∗ ∗ d D m n dn Dn E(dn , Dn ) 0.6 1.6 25 28 0.58502 1.63380 0.5230 0.6 1.6 50 26 0.61256 1.59286 0.7962 0.6 1.6 100 23 0.61941 1.57982 1.1143 0.6 1.6 200 27 0.61503 1.59456 1.6846 3 Intervalo: [50, 10 ] d∗ D∗ m n dn Dn E(dn , Dn ) 0.6 1.6 25 36 0.88483 1.01314 0.5251 0.6 1.6 50 35 0.54019 1.74054 0.8013 0.6 1.6 100 37 0.47938 1.85346 1.1375 0.6 1.6 200 28 0.60716 1.61777 1.6943 Intervalo: [105 , 106 ] d∗ D∗ m n dn Dn E(dn , Dn ) 0.6 1.6 25 29 0.56482 1.62585 0.5232 0.6 1.6 50 29 0.60164 1.60259 0.8016 0.6 1.6 100 31 0.61518 1.59460 1.1376 0.6 1.6 200 32 0.60187 1.60496 1.6941

En la Tabla 4 se tienen los resultados para ε = 0.2. Lo que se observa es la estimación de los parámetros es mejor cuando se tienen muestras de la caída de presión para tiempos grandes. Los errores más grandes ocurren cuando se tienen pocas muestras y están se toman en un intervalo corto al inicio de la simulación. Si este es el caso, solo aumentando el número de muestras se puede mejorar la estimación. En la Tabla 5 se repite la prueba anterior usando ε = 0.4, y debido a que el ruido es mayor, se observa que los resultados de la estimación de los parámetros empeora. El mejor de los casos ocurre cuando se tiene la mayor cantidad de muestras en instantes de tiempo grandes. Para entender lo que implica equivocarse en la determinación de parámetros, generamos las curvas de la caída de presión cuando se usan los verdaderos valores de los parámetros (d∗ = 0.6 y D∗ = 1.6) y uno de las estimaciones que obtuvimos en la prueba anterior, por ejemplo, cuando d = 0.6153 y D = 1.6252. En la Fig. 16 se comparan estas curvas. Hay diferencias entre ellas y puede ocurrir que para tiempos más largos aumenten éstas. Hay que tomar en cuenta que dependiendo de otros parámetros, como el radio exterior re , el mismo nivel de ruido en las mediciones puede ser más destructivo para cierto rando de valores de re . Por ejemplo, las pruebas anteriores se hicieron para re = 500 y vemos que podemos obtener buenas aproximaciones para ε = 0.4 con m = 200 muestras. Si cambiamos el radio exterior, 21

Tabla 5: Prueba con datos con ruido, variando el número de muestras y el intervalo en donde se encuentran con ε = 0.4. Intervalo: [50, 106 ] ∗ ∗ d D m n dn Dn E(dn , Dn ) 0.6 1.6 25 29 0.59326 1.67818 1.0575 0.6 1.6 50 24 0.64473 1.59837 1.6265 0.6 1.6 100 27 0.67166 1.56008 2.2742 0.6 1.6 200 28 0.65294 1.60464 3.4523 3 Intervalo: [50, 10 ] d∗ D∗ m n dn Dn E(dn , Dn ) 0.6 1.6 25 37 0.92408 1.02894 1.0882 0.6 1.6 50 36 0.54603 1.85965 1.6336 0.6 1.6 100 69 0.52115 1.89378 2.3619 0.6 1.6 200 32 0.68302 1.60893 3.4665 Intervalo: [105 , 106 ] d∗ D∗ m n dn Dn E(dn , Dn ) 0.6 1.6 25 37 0.53804 1.67291 1.0525 0.6 1.6 50 34 0.61520 1.61777 1.6339 0.6 1.6 100 38 0.65018 1.59861 2.3431 0.6 1.6 200 32 0.61529 1.62516 3.4651

de modo que re = 2000, la Tabla 6 se puede ver como se afecta a la estimación de los parámetros cuando el nivel de ruido va en aumento. Sin ruido, la estimación es precisa, pero para ε = 0.15 ya hay diferencias importantes. Se requiere aumentar el número de muestras de 100 a 400 para mejorar la estimación de los parámetros cuando ε = 0.2. De este modo, para re = 500 un ruido con ε = 0.2 es moderado, mientras que para re = 2000, ε = 0.2 es un ruido destructivo si no se tienen suficientes muestras. Tabla 6: Prueba con datos con ruido, variando el número de muestras y el intervalo en donde se encuentran con ε = 0.4. Intervalo: [120, 106 ] ∗ ∗ d D m ε n dn Dn E(dn , Dn ) 0.6 1.6 100 0.00 42 0.60002 1.59995 0.0003 0.6 1.6 100 0.05 42 0.59541 1.60781 0.2823 0.6 1.6 100 0.10 30 0.63696 1.50195 0.5691 0.6 1.6 100 0.15 27 0.65150 1.47050 0.8514 0.6 1.6 100 0.20 25 0.66302 1.45052 1.1358 0.6 1.6 200 0.20 31 0.66611 1.45642 1.6651 0.6 1.6 400 0.20 33 0.60321 1.61471 2.3123

22

100 20 10 5

Caída de presión

50

d=0.600, D=1.600 d=0.615, D=1.625

1e+02

1e+03

1e+04

1e+05

1e+06

Tiempo

Fig. 16: Comparación entre la curva de la caída de presión generada con los verdaderos valores de los parámetros (en rojo) y la curva obtenida con los parámetros estimados (en azul).

Referencias [Camacho-Velázquez 08] Rodolfo Camacho-Velázquez, Gorgonio Fuentes-Cruz & Mario Vásquez-Cruz. Decline curve analysis of fractured reservoirs with fractal geometry. SPE Res. Eval. & Eng., vol. 11, no. 3, pages 606–619, 2008. SPE104009-PA. [Chen 07]

Zhangxin Chen. Reservoir simulation: mathematical techniques in oil recovery. CBMS-NSF regional conference series in applied mathematics; 77. SIAM, 2007.

[Davies 79]

B. Davies & B. L. Martin. Numerical inversion of Laplace transform: A survey and comparison of methods. Journal of Computational Physics, vol. 33, no. 1, pages 1–32, 1979.

[de Hoog 82]

F. R. de Hoog, J. H. Knight & A. N. Stokes. An improved method for numerical inversion of the Laplace trans23

forms. SIAM J. Stat. Comput., vol. 3, no. 3, pages 357–366–, 1982. [Press 07]

W. H. Press, S. A. Teukolsky, W. T. Vetterling & B. P. Flannery. Numerical recipes: The art of scientific computing. Cambridge University Press, 3rd edition edition, 2007.

[Tarasov 05]

Vasily E. Tarasov. Fractional hidrodinamics equations for fractal media. Annals of Physics, vol. 318, pages 286–307, 2005.

[Watson 96]

G. N Watson. A treatise on the theory of bessel functions. Cambridge University Press, 1996.

24