Simulación de Yacimientos: Modelamiento de Flujo Monofásico (Gas) en 1 Dimensión Di mensión – Solución Solución Método Psor. Andrés Felipe Villota Rosada – Rosada –
[email protected] Jorge Alejandro González Firigua – Firigua –
[email protected] [email protected] Aron Kennedy Pantoja -
[email protected]
Resumen En este documento se abordarán el modelamiento del flujo de un yacimiento de gas en una dimensión. En la simulación de un yacimiento son varias las variables que influyen en el comportamiento de dicho flujo, propiedades de la roca de yacimiento, propiedades PVT del o los fluidos presentes, así, como la geometría del yacimiento; el modelamiento del flujo se inicia planteando un modelo matemático en el que se encuentran las expresiones para el flujo en una dimensión, así como las condiciones limite y de frontera para dar solución a las expresiones correspondientes a cada fase. Seguidamente se elabora un modelo numérico numérico que permita elaborar un código usando el software Matlab y de esta manera conocer el comportamiento de las presiones en el yacimiento, este último último se divide en e n múltiples bloques de los cuales se tiene conocimiento de las propiedades ya mencionadas; seguidamente se hace uso de los volúmenes finitos finitos para encontrar a partir del modelo matemático, expresiones con las que se pueda describir el comportamiento del flujo en cada uno de los bloques de la malla elaborada; Lo anterior permite encontrar un sistema de ecuaciones para el comportamiento del flujo en el yacimiento al que se da solución mediante métodos numéricos. Finalmente, como como se mencionó mencionó anteriormente, anteriormente, dada la la complejidad de las ecuaciones, se escribe un código de computadora que permite conocer el valor de las propiedades a diferentes tiempos, generar el sistema de ecuaciones y a partir su solución, encontrar el comportamiento de las presiones en el yacimiento.
Introducción La simulación de yacimientos es un área muy importante en la industria del petróleo. Juega un papel fundamental en el e l análisis del comportamiento de un yacimiento y de las variaciones que se presentan debido a cambios en las diferentes variables que influyen sobre el sistema, tales como la ubicación de los pozos y su papel en el desarrollo desarrollo del de l yacimiento, también los cambios que se puedan presentar en las propiedades de los fluidos presentes y de las propiedades de la roca. El objetivo en la simulación simulación de un yacimiento yacimiento es mediante mediante las ecuacio ecuaciones nes de difusividad difusividad y haciendo uso de software computacional describir de manera realista los procesos que ocurren al interior del yacimiento, de manera que se permitan obtener buenas aproximaciones para aspectos de interés interés como las reservas reservas que se pueden recuperar y de la manera óptima óptima de hacerlo. A continuación, se presentan las propiedades del yacimiento, del fluido que en este caso es gas, el modelo matemático, numérico y computacional para describir el flujo en una
dimensión. También se presenta el diagrama del procedimiento a seguir para la elaboración del simulador. Finalmente se presentan los resultados y su respect ivo análisis.
Modelo físico está conformado por la información de la malla necesaria para modelar el flujo, en este caso monofásico (gas) en una sola dimensión. En la tabla 1 se listan los parámetros constantes tanto del yacimiento (ancho, espesor, presión inicial, presión de fondo, compresibilidad de poro) como del fluido (densidad a condiciones estándar) y parámetros de entrada.
Parámetro Ancho (ft) Espesor (ft) Compresibilidad (1/psia) Presión Inicial (psia) Densidad (lbm/ft3) Delta de tiempo (días) rea de flujo (ft2) Pwf (psia)
Magnitud 100 90 2,1x10-6 2100 0,059186
9000 2000 Tabla 1. Parámetros de simulación
Bloque 1 2 3 4 5 6 7 8 9 10 11
Profundidad (ft)
Espesor (ft)
Permeabilidad (mD)
3930.000 3887.400 3851.670 3841.150 3806.350 3800.060 3806.350 3824.150 3851.670 3887.400 3930.000
309.510 15.475 325.800 16.290 342.250 17.148 361.000 18.050 380.000 19.000 400.000 20.000 380.000 19.000 361.000 18.050 342.950 17.148 325.800 16.290 309.510 15.476 Tabla 2. Información de la malla
Porosidad 0.153 0.146 0.139 0.132 0.126 0.120 0.126 0.132 0.139 0.146 0.153
Figura 1. Diagrama de la malla
Los pozos inyectores se encuentran localizados en los bloques 1 y 11 mientras que el pozo productor se encuentra localizado en el bloque 6. Para las propiedades PVT se requiere una función que modele su comportamiento, en el caso de la viscosidad, la función se encontró por medio de regresión lineal y la función que modela el comportamiento del factor volumétrico se encontró por medio de regresió n potencial.
Factor volumétrico
Viscosidad
0,003
0,04
0,0025
] p c 0,03 [
] f c 0,002 s / B 0,0015 R [ g 0,001 B
d a d i 0,02 s o c s i 0,01 V
0,0005 0
0 0
1000
2000
3000
4000
5000
0
1000
Presión [psi]
20.052−.
2000
3000
4000
5000
Presión [psi]
4060.0106
Figura 2. Propiedades Pvt del fluido
Modelo matemático: Para el modelo matemático se parte de la ecuación general de difusividad y se restringe para flujo monofásico (gas) en una sola dimensión (x).
Donde
∝
∝ ∝ ∝ ( ) ∝̃ ( ) ̃ ( ) () (̃)
1
=Area de flujo
2 3 4
Definición de condiciones iniciales En
0
la presión es todos los bloques es igual a la presión inicial
0 2100 ∀
Definición de condiciones de frontera
5
Se asume un yacimiento que tiene una condición de frontera tipo Von Newman (barrera de no flujo):
0
6
0 ∗ ()
Definición de tasa terminal constante, sea de inyección o producción:
Donde
es el factor de forma:
ℎ l2 og(/) 0,14 √∆ ∆ ∅ ()∅°1∅° ∅ ∅
Para un bloque rectangular el radio equivalente
8
es
Definición de la porosidad como una función de la presión
Modelo numérico:
7
9 10 11
Realizando la discretización por volúmenes finitos, se toma cada término de la ecuación y luego se reúne cada uno de los términos para armar la ecuación ya discretizada.
Primer término (acumulación)
Integramos
Segundo término (transporte)
∅
12
∫∆ ∫ ∅ ∫∆ ∫∆ ∅ + ∅ ∅ ∫∆ + ∅ ∅
13
∫∆ ∫∆ ∅
14
Tercer término (Fuente)
∫∆ ∆ ∅ ∅ + − ∆
15
∫∆ ∫∆ ∫∆ ∆ ∆ + +∆+ ∆−
16 17
Luego evaluamos el potencial en dirección x
−
Definiendo la transmisibilidad como:
18 19
−
+ ∆ + − ∆ − ∆∅ + ∅ ∅ ∆ ∗∆∅ + (
Reemplazando:
Además
Entonces
20 21 )
(22)
+ ∅ ++ − − ∆ ( )
23
24 25
+(++ +)+ −(+ −+) − + ∅ ∆ ( ) 26
Utilizando factor común
∅ ∅ + + + − -+ − ∆ ++ = -∆ + ∗ + − ∗ − 27 −
Se determina
∅ ∅+ ∅′ ° ∅ ′∅ ∅+ + ∅ ∅ 1 1 + 1 + ++ + −+ ++ + −+
28 29
=
30
Como lo que nos interesa saber es el cambio de la presión con el tiempo, partimos de (27), debemos agrupar los términos que contengan , , , de manera que usando la notación de stencils la ecuación queda planteada de la siguiente manera:
,
Donde y representan los srencil west, central y east respectivamente, termino independiente.
31 representa el
−∗ ⁄ +∗ ⁄ ′ ∗ ∗ +⁄ −⁄ ∆ ′ ∗ ∗ +⁄+ −⁄ − ∗ ∆ Condiciones iniciales
Para
32 33 34 35 36
Condiciones límite
− = 0 + = 0 Para
37
Para
38
luego el sistema de ecuaciones queda planteado de la siguiente manera:
Método de Solución
++ + −+ ∀ 2:10 + −+ 1 ++ + 11
39 40 41
Para dar solución al sistema de ecuaciones y encontrar las presiones en cada bloque para un tiempo determinado se usó el método numérico PSOR, la ecuación general para un elemento i del vector P n+1 es:
− + + 1 ∗ , ∗,∗ ,∗ = =+
42
w es la constante de relajación que acelera la convergencia del método de Gauss-Seidel. Para la convergencia se usó el criterio de convergencia de error absoluto:
=:max| |≤
1
43
Cuando se cumple el anterior criterio las presiones se actualizan, es decir y es necesario volver a calcular las propiedades de esta manera llegando a un nuevo sistema de ecuaciones. Para tener la solución al tiempo n se debe cumplir:
=:max|+ |≤
44
Lo anterior se debe validar con el balance de materiales EBM:
= + ∑ ∑= , 1 ≤
45
Diagrama de flujo
Resultados Presión en los bloques con el avance del tiempo Bloques
dia15
dia30
dia45
dia60
dia90
dia120
(Inyector) 1
2222,80
2233,50
2234,80
2235,00
2235,00
2235,00
2
2195,60
2206,00
2207,20
2207,40
2207,40
2207,40
3
2168,50
2178,30
2179,50
2179,70
2179,70
2179,70
4
2141,70
2150,50
2151,60
2151,70
2151,70
2151,70
5
2115,00
2122,50
2123,40
2123,50
2123,60
2123,60
(Productor) 6
2088,50
2094,30
2095,10
2095,20
2095,20
2095,20
7
2115,00
2122,50
2123,40
2123,50
2123,60
2123,60
8
2141,70
2150,50
2151,60
2151,70
2151,70
2151,70
9
2168,50
2178,30
2179,50
2179,70
2179,70
2179,70
10
2195,60
2206,00
2207,20
2207,40
2207,40
2207,40
(Inyector) 11
2222,80
2233,50
2234,80
2235,00
2235,00
2235,00
Tabla 1. Resultados del modelamiento del flujo en el yacimiento durante 120 días
Análisis de resultados La tabla 3 nos arroja un resultado de presiones para un periodo de 15 a 120 días para los 11 bloques que conforman la malla. Se puede notar que las presiones son simétricas para cada tiempo en los pozos inyectores (P (1)=P (11), P (2)=P (10)…) con excepción al pozo productor que se encuentra en el bloque 6. La información de la tabla 3 nos arroja la posición de los pozos inyectores y el pozo productor, ya que en los bloques 1 y 11 se presenta un aumento continuo de la presión debido a condición de tasa terminal constante, mientras que en el bloque 6 la presión disminuye con el tiempo. Los pozos inyectores, presurizan su bloque (1 y 11) y se genera una propagación de la onda de presión a través de los bloques cercanos, debido a la simetría, los bloques 1 y 11 propagan el aumento de presión por los bloques 2 y 10 y así respectivamente para los demás. Este comportamiento se mantiene en el tiempo, en cuanto al comportamiento en el tiempo del pozo productor, se observa que la presión es menor que en los bloques vecinos, dando como resultado un comportamiento esperado, ya que la presión siempre se mantiene por debajo de la presión inicial, lo que nos confirma que el comportamiento es el esperado para un pozo productor.
Figura 3. Comportamiento y Avance de la Presión en el tiempo, pozos 1:11
CONCLUSIONES La simulación realizada para las presiones tanto del pozo productor como para los pozos inyectores nos arrojó un resultado coherente, esto debido a que en los pozos inyectores presentan un aumento continuo de la presión debido a condición de tasa terminal constante, generando una propagación de la onda de presión que va de los bloques 1 y 11 y se propagan a los bloques vecinos, y en el pozo productor la presión es menor que la de los bloques vecinos y siempre se encuentra por debajo de la presión inicial.
Anexos INFORMACI�N DE LA MALLA espesorX = [309.51 325.80 342.95 361.00 380.00 400.00 380.00 361.00 342.95 325.80 309.51]; anchoY = ones(1,length(espesorX))*100; espesorZ = ones(1, length(espesorX))*90; Prof = [3930.00 3887.40 3851.67 3824.15 3806.35 3800.00 3806.35 3824.15 3851.67 3887.40 3930.00]; K = [15.476 16.290 17.148 18.050 19.000 20.000 19.000 18.050 17.148 16.290 15.476]; Poroini = [0.153 0.146 0.139 0.132 0.126 0.120 0.126 0.132 0.139 0.146 0.153]; Cp = 2.1E-6; Pi = ones(1, length(espesorX))*2100; INFORMACI�N DE LOS POZOS Pwf=2000;
%------------------------------------------------------------------------% PARAMETROS DE LA SIMULACI�N dt=0.5; Tol=10e-6; %------------------------------------------------------------------------% INFORMACI�N DEL FLUIDO Densid = 0.059186; %-----------------------------------------------------------------------------%parametros invariantes en el tiempo %volumen de los bloques V=espesorX.*anchoY.*espesorZ/5.615; %area de las caras Area=anchoY.*espesorZ; %DeltaZ dZright(length(Prof))=0; dZleft(1)=0; for i=1:length(Prof)-1 dZright(i)=Prof(i+1)-Prof(i); dZleft(i+1)=Prof(i+1)-Prof(i); end %factor de forma del pozo Bc=0.001127; req=0.14*sqrt(espesorX(6)^2+anchoY(6)^2); Gw=(2*pi*Bc*K(6)*espesorZ(6))/(log(req/0.33)); tol=10e-5; itmax=1000; errormax=20; t=0; Pn=Pi; j=1; Pcal(1,:)=Pi; errorEVB(j)=1; EBM=1; while (t<120) t=t+dt %primera iteracion de tiempo %suponemos Pk presion a la sigueinte nivel de tiempo Pk=Pn+10; gamma = 0.21584e-3*Densid*32.174; d=derivadaPhiBg(Cp,Poroini,Pi,Pk); %Qv=Qv(Pn,Gw); Qv=zeros(1,length(Pk)); Qv(1)=5000000; Qv(11)=5000000; Qv(6)= -Gw*(Pk(6)-Pwf)/(Bg(Pk(6))*u(Pk(6))); Qv=5.614583*Densid*Qv; Qvv(j)=Qv(6); errormax=10; while errormax>1 T=Transmisibilidad_cara(K,Area,espesorX,Pk); Tleft=T(1,:); Tright=T(2,:);
%Matriz A A=zeros(length(K)); for i=1:length(K) A(i,i)=-Tright(i)-Tleft(i)-V(i)*d(i)/(dt*5.614583); end for i=1:length(K)-1 A(i,i+1)=Tright(i); A(i+1,i)=Tleft(i+1); end %vector b b=(-5.64583*Densid*Qv - V.*d.*Pn/(dt*5.614583) + Tright.*gamma.*dZright - Tleft.*gamma.*dZleft)'; Pk=Pk'; %[P,error,it]= sor(A,b,1.2,Pk,tol,itmax); P=inv(A)*b; errorvector=P-Pk; errormax=max(abs(errorvector)); Pk=P'; end Pn=Pk; j=j+1; Pcal(j,:)=Pn; %criterio de balance de materiales %Z=Prof+min(Prof); Z=0; num=V.*poro(Cp,Poroini,Pi,(Pcal(j,:)gamma*Z))./(5.615*dt*Bg(Pcal(j,:)-gamma*Z)) V.*poro(Cp,Poroini,Pi,(Pcal(j-1,:)-gamma*Z))./(5.615*dt*Bg(Pcal(j-1,:)gamma*Z)); den=Qv; numm=0; denn=0; for l=1:length(K) numm=numm+num(i); denn=denn+den(i); end EBM=abs(numm-denn);
end x=[1 2 3 4 5 6 7 8 9 10 11]; plot(x,Pcal(1,:)) axis([1 11 min(min(Pcal)) max(max(Pcal))]) pause(1) for i=1:15 plot(x,Pcal(i,:)) axis([1 11 min(min(Pcal)) max(max(Pcal))]) pause(4/241) end for i=15:241 plot(x,Pcal(i,:)) axis([1 11 min(min(Pcal)) max(max(Pcal))])
pause(4/241) end for i=1:15 plot(x,Pcal(i,:)) axis([1 11 min(min(Pcal)) max(max(Pcal))]) hold on pause(4/241) end for i=1:7 j=30*i plot(x,Pcal(j,:)) axis([1 11 min(min(Pcal)) max(max(Pcal))]) pause(4/241) end EEEE=max(errorEVB)
Calaculo de transmisibilidades function Transmi = Transmisibilidad_cara(K,Area,espesorX,P) %trasnmicibilidad centro bloque T=(0.001127*K.*Area)./(u(P).*Bg(P).*espesorX); %trasnsmisibilidad cara por porm,edio ponderado por distancia Transmi(1,1)=0; Transmi(2,length(Area))=0; for i=1:length(Area)-1 %caras a la izquierda (left) Transmi(1,i+1)=(T(i)*espesorX(i)+T(i+1)*espesorX(i+1))/(espesorX(i)+espes orX(i+1)); %caras a la derecha (right) Transmi(2,i)=(T(i)*espesorX(i)+T(i+1)*espesorX(i+1))/(espesorX(i)+espesor X(i+1)); end
Calculo de la solución SOR function [P,error,it]= sor(A,b,w,P0,tol,itmax) %rutina que implementa un metodo sor basico it=0; error=1000; D=diag(diag(A)); E=-tril(A,-1); F=-triu(A,+1); res=norm(A*P0-b); error=res/norm(b); while ((it
tol)) it=it+1; P=(D-w*E)\(w*F*P0+(1-w)*D*P0+w*b); res = norm((A*P)-b); error = res/norm(b); P0=P; end
