UNAM-FI
Tarea 3 Simulación
de modelos de
yacimientos con
fluidos
compresibles:
Simulación Matemática de Yacimientos
Villegas Tapia Edgar
19-5-2018
Introducción
Los simuladores de gas se utilizan para llevar a cabo los ajustes de la historia de
datos de presión y producción, así como la predicción del comportamiento de un
pozo y/o yacimiento de gas. Los estudios para este tipo de yacimientos y pozos tal
vez no siempre son los más fáciles de resolver, dependiendo de la simulación o
problemas reales que se traten de resolver como son la complejidad geológica del
yacimiento, geometría del pozo, etc. Los principales parámetros que se pueden
obtener de este tipo de simulador son entre otros: el volumen de gas inicial, G, el
gasto de producción de gas, qg, el comportamiento la distribución de las presiones,
p, y/o pseudopresiones de los gases reales, m(p).
Un pozo vertical es perforado es un yacimiento de gas con un radio de 600 pies. El
yacimiento es discretizado en 9 celdas en la dirección r. El yacimiento es horizontal,
homogéneo e isotrópico con un espesor de 30 pies, porosidad es de 0.20 y la
permeabilidad de 0.1 md. La presión inicial del yacimiento es de 4500 [psia], no hay
flujo en la frontera externa. El pozo tiene un rw de 0.25 pies, produce 50,000 SCF/D.
El gas que se produce tiene una densidad relativa de 0.7 a una temperatura de
yacimiento de 150°F. La ecuación que representa el flujo de fluidos a través del
medio poroso es:
1𝜕 𝑘𝑟 𝜕𝑝 𝑉𝑏 𝜕 𝜙
[𝑟𝛽𝑐 ]= ( )
𝑟 𝜕𝑟 µ𝑔 𝐵𝑔 𝜕𝑟 𝛼𝑐 𝜕𝑡 𝐵𝑔
Calcular la distribución de presión en el yacimiento para un tiempo de 1, 10, 50, 100,
500 y 1000 días, usando el programa de computo desarrollado (pueden utilizar
cualquier método de programación para llegar a la solución); imprimir una gráfica
de presión vs log(distancia) donde se grafiquen los resultados anteriores y entregar
una tabla de las presiones de cada una de las celdas a cada tiempo especificado.
2. Imprimir una gráfica de (Pi – Pwf) vs tiempo (desde 0.1 días hasta los 1000 días).
3. Entregar el código impreso de su programa desarrollado.
La presión de fondo fluyendo la pueden calcular del siguiente modelo de pozo:
2𝜋𝛽𝑐 𝑘𝑟1 ℎ1
𝑞𝑠𝑐 𝑟 (𝑝1 − 𝑝𝑤𝑓 )
𝐵1 µ1 𝑙𝑛 (𝑟1 )
𝑤
Desarrollo
Se desarrolla el programa de cómputo para la resolución del problema en
diferencias finitas para un yacimiento de fluidos compresibles.
! información de entrada empleada:
4500,P(x=L,t>0) !psi (presion en frontera)
0.0 ,P(x=0,t>0) !psi (pwf)
4500 ,P(x,t=0) !psi (pws)
600 , !distancia en pies
0.1 , !perme en mD
0.20 ,! porosidad
0.001127 ,!Bc
5.615, !ac
30 , !espesor [ft]
0.7, !gravdead específica
150, !temperatura °F
0.25, !rw
50000, !producción
!Código de programación. !En el primer paso declaramos las variables
!a utilizar
program solucion
real,dimension (6):: t
real,dimension (6):: t,dist,pcl,ddd
real, dimension(102,11) ::p
real, dimension(102)::posicionx,posnew,a,c,d,pe,w,g,be
integer imax, b, new, n, j, ka, i,h
double precision po, pl, pr, L, k, Porosidad, visc, deltax, sum,bc,ac,densr,temp,rw,prod,bg,xmg,ymg,k1,k2,k3,z
double precision alpha, nhe, beta,ppc,tpc,tcr,pcr,ddz,cz,bz,az,az1,az2,az3,bg
real::pws=4500
t=(/1.0,10.0,50.0,100.0,500.0,1000.0/)
!Con esta función abrimos un block de notas, que es un archivo .txt
!desde el programa, se manda a llamar con el comando open.
open(unit=10,file="Datos.txt",status="unknown",action="read")
!Con read, leemos la información de entrada.
read (10,*) po
read (10,*) pl
read (10,*) pr
read (10,*) L
read (10,*) k
read (10,*) Porosidad
read (10,*) bc
read (10,*) ac
read (10,*) h
read (10,*) densr
read (10,*) temp
read (10,*) rw
read (10,*) prod
!Se cierra el archivo abierto en el paso anterior archivo
close(10)
call system("color 04")
!Se despliega un menu, donde se pregunta cuantos nodos
!o celdas trabajara el programa
call system("color 02")
print*,"el analisis estipulado sera de"
print*,''
print*,' 9 nodos'
print*,''
imax=9
new=1
do b=2,10,1
deltax=L/(imax-1)
posicionx(b)=(new-1)*deltax
new=new+1
end do
!Aquí se finaliza la construcción del mallado con un ciclo do.
call system("pause")
call system("color 03")
call system("cls")
print*,'La solucion que se aplicara es la sigguiente'
print*,''
print*,'Implicito'
call system("pause")
call system ("cls")
!Después se emplea correlacion del gas
!de Lee- Gonzalez para el cáculo de viscosidad
!para cada celda
xmg=3.448+(986.4/(temp+460))+(0.2897*densr)
ymg=2.447-(0.2224*xmg)
k1=(9.4+(0.5794*densr))*((temp+460)**1.5)
k2=209+(550.4*densr)+(temp+460)
k3=k1/k2
visc=(k3/10000)*(exp(xmg*(((densr*28.9)/62.428)**ymg)))
!Para el cálculo de z por celda se utiliza la correlacion
! de Brill and beggs que queda de la siguiente forma.
ppc=756.8-(131*densr)-(3.6*densr*densr)
tpc=169.2+(349.5*densr)-(74*densr*densr)
tcr=(temp+460)/tpc
pcr=pr/ppc
ddz=10**(0.3106-(0.49*tcr)+(0.1824*tcr*tcr))
cz=0.132-(0.32*log10(tcr))
bz1=(0.62-(0.23*tcr))*pcr
bz2=((0.066/(tcr-0.86))-0.037)*pcr*pcr
bz3=(0.32/(10**(9*(tcr-1))))*(pcr**6)
bz=bz1+bz2+bz3
az1=(tcr-0.92)**0.5
az2=(0.36*tcr)
az3=0.101
az=(1.39*az1)-az2-az3
z=az+((1-az)/exp(bz))+(cz*(pcr**ddz))
!el cálculo de bg (factor de volumen de gas) viene
!dado por la ecuación de estado de los gases reales
!puede caluclarse un Bg por cada celda
bg=(.02827*temp*z)/pr
!Solución numérica a la ecuación de difusión por metodo implícito
!En yacimientos de fluidos compresibles (yacimientos de gas)
!por medio de diferencias finitas.
print*,'Esquema de solución implicito'
p=0.0D0
nhe=(L*visc*L*h*Porosidad*bg)/(k*bc*ac*(.02817*temp*z))
do n=1,6,1
do j=1,imax,1
posnew(j)=posicionx(j+1)
end do
do j=1,imax,1
p(j,1)=po
101 end do
do i=1,10,1
p(1,i)=pl
p(imax,i)=pr
end do
alpha= (t(n)*(nhe)*(L+(deltax/2)))/(deltax**2)
beta = (t(n)*(nhe)*(L-(deltax/2)))/(deltax**2)
do i=1,imax
a(i)=alpha
be(i)=alpha+beta
c(i)=beta
end do
new=1
do j=1,imax,1
if (j==1) then
d(1)=-p(j,n)-a(1)*pl
else
if (j==imax) then
d(imax)=-p(j,n)-c(imax)*pr
else
d(j)=-p(j,n)
end if
end if
end do
!Con este metodo obtenemos los coeficientes de una matriz
!tridiagonal, esta se resuelven por medio de el algoritmo
!de Thomas para la resolución de matrices tridiagonales.
!Algoritmo de Thomas para la solución de matrices tridiagonales
w(1)=c(1)/be(1)
g(1)=d(1)/be(1)
do i=2,imax
w(i)=c(i)/(be(i)-a(i)*w(i-1))
g(i)=(d(i)-a(i)*g(i-1))/(be(i)-a(i)*w(i-1))
end do
pe(imax)=g(imax)
do i=imax-1,1,-1
pe(i)=g(i)-w(i)*pe(i+1)
end do
!Fin thomas
do ka=2,imax-1
if (n.ne.10) p(ka,n+1)=pe(ka)
end do
!Ahora se generan los datos para el ejercicio 1
print*,''
print*,"para tiempo",t(n)
print*,''
print*,' NODO distanca Presion[psi]'
print*,''
if (n==1)then
open(unit=11,file="EsquemaImplicito_Resultados1.txt",status="unknown",action="write")
do b=2,imax+1,1
150 if (b==4) then
write(11,*) posicionx(b),(1*(p(b-1,n)))
print*, new,posicionx(b),(1*(p(b-1,n)))
new=new+1
else if (b==6)then
write(11,*) posicionx(b),(1*(p(b-1,n)))
print*, new,posicionx(b),(1*(p(b-1,n)))
else if (b==8)then
write(11,*)posicionx(b),(1*(p(b-1,n)))
print*, new,posicionx(b),(1*(p(b-1,n)))
else
write(11,*) posicionx(b),(p(b-1,n))
print*, new,posicionx(b),(p(b-1,n))
end if
new=new+1
end do
close(11)
else if (n==2)then
open(unit=12,file="EsquemaImplicito_Resultados2.txt",status="unknown",action="write")
do b=2,imax+1,1
if (b==4) then
write(12,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
new=new+1
else if (b==6)then
write(12,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else if (b==8)then
write(12,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else
write(12,*) posicionx(b),(1*(p(b-1,n)))
print*, new,posicionx(b),(p(b-1,n))
end if
new=new+1
end do
close(12)
else if (n==3)then
open(unit=13,file="EsquemaImplicito_Resultados3.txt",status="unknown",action="write")
do b=2,imax+1,1
if (b==4) then
write(13,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
new=new+1
else if (b==6)then
write(13,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else if (b==8)then
200 write(13,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else
write(13,*) posicionx(b),(1*(p(b-1,n)))
print*, new,posicionx(b),(p(b-1,n))
end if
new=new+1
end do
close(13)
else if (n==4)then
open(unit=14,file="EsquemaImplicito_Resultados4.txt",status="unknown",action="write")
do b=2,imax+1,1
if (b==4) then
write(14,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
new=new+1
else if (b==6)then
write(14,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else if (b==8)then
write(14,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else
write(14,*) posicionx(b),(1*(p(b-1,n)))
print*, new,posicionx(b),(p(b-1,n))
end if
new=new+1
end do
close(14)
else if (n==5)then
open(unit=15,file="EsquemaImplicito_Resultados5.txt",status="unknown",action="write")
do b=2,imax+1,1
if (b==4) then
write(15,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
new=new+1
else if (b==6)then
write(15,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else if (b==8)then
write(15,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else
write(15,*) posicionx(b),(1*(p(b-1,n)))
print*, new,posicionx(b),(p(b-1,n))
end if
new=new+1
end do
251 close(15)
else if (n==6)then
open(unit=16,file="EsquemaImplicito_Resultados6.txt",status="unknown",action="write")
do b=2,imax+1,1
if (b==4) then
write(16,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
new=new+1
else if (b==6)then
write(16,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else if (b==8)then
write(16,*) posicionx(b),(-1*(p(b-1,n)))
print*, new,posicionx(b),(-1*(p(b-1,n)))
else
write(16,*) posicionx(b),(1*(p(b-1,n)))
print*, new,posicionx(b),(p(b-1,n))
end if
new=new+1
end do
close(16)
end if
end do
call system ("gnuplot grafico3.txt")
!Para el ejercicio 2, se grafica (PI-PWF) vs tiempo
open(unit=21,file="EsquemaImplicito_Resultados1.txt",status="unknown",action="read")
open(unit=22,file="EsquemaImplicito_Resultados2.txt",status="unknown",action="read")
open(unit=23,file="EsquemaImplicito_Resultados3.txt",status="unknown",action="read")
open(unit=24,file="EsquemaImplicito_Resultados4.txt",status="unknown",action="read")
open(unit=25,file="EsquemaImplicito_Resultados5.txt",status="unknown",action="read")
open(unit=26,file="EsquemaImplicito_Resultados6.txt",status="unknown",action="read")
read (21,*) dist(1),pcl(1)
read (22,*) ddd(1)
read (22,*) dist(2),pcl(2)
read (23,*)ddd(2)
read (23,*) dist(3),pcl(3)
read (24,*)ddd(3)
read (24,*) dist(4),pcl(4)
read (25,*)ddd(4)
290 read (25,*) dist(5),pcl(5)
read (26,*)ddd(4)
read (26,*) dist(6),pcl(6)
call system("pause")
call system("color 09")
call system("cls")
print*,''
print*,'Ejercicio 2 PWS-Pwf vs tiempo'
print*,''
print*,' PWS-PWF Tiempo '
print*,''
open(unit=8,file="ejercicio2.txt",status="unknown",action="write")
do nn=1,6,1
pwf=((prod*bg*visc*log((L/9)/rw))/(2*3.141516*bc*k*h))+pcl(nn)
Print*,(pws-pwf),t(nn)
write(8,*) t(nn),(pws-pwf)
end do
call system ("gnuplot grafico4.txt")
end program valores
Fig. 1 – Solución a 9 celdas del problema.
Figura 1, 2 y 3. – Distribución de presiones por distancia
Fig. 4. - Gráfica de presión vs distancia. Fig. 4 – Perfiles de
presiones vs
distancia.
Fig. 5 – Resolución de ejercicio 3, Gráfico de (Pws – Pwf) vs tiempo.
Conclusiones
La resolución de este problema muestra que las fronteras se alcanzan con prontitud,
pero esto es normal, ya que debido a que se trata de un yacimiento de fluidos
compresibles (yacimientos de gas), se está acostumbrado a presenciar este
comportamiento en este tipo de yacimientos, que alcanzan su frontera rápidamente
comparados con los yacimientos de fluidos ligeramente compresibles.
