Instituto Politécnico Nacional
Escuela Superior De ingeniería y
Arquitectura.
Unidad Ticoman ciencias de la tierra.
ingeniería Petrolera.
Simulación Numérica de Yacimientos.
Dr. Octavio Cazares Candia
Trabajo final análisis de simulación en
FORTRAN
De Vicente Pérez Giovanni
Boleta: 2014041802
Introducción
A diario se ve cómo crece la tecnología, un ejemplo de ello es la creciente
capacidad y actualización de las computadoras, la inmensa investigación en el
campo de la Ciencia de la Computación que otorgan nuevas herramientas, para
apoyar el proceso de la toma de decisiones en diversas disciplinas y áreas de
diseño en la industria. La Simulación es una de las herramientas que puede
aproximar cualquier comportamiento dinámico de un yacimiento, de esta manera se
pueden prevenir eventos indeseables y se corrigen con tiempo de manera que se
alcancen con éxito los proyecto.
La demanda de energéticos es cada vez más exigente, por lo que se debe hacer
hincapié y redoblar esfuerzos para generar y hacer de la industria petrolera un
apoyo energético y económico, a través de la SNY, en la cual se pueda realizar un
análisis más profundo en el estudio de los yacimientos y se obtengan factores de
recuperación altos.
En la industria petrolera lo mas importante es la producción de hidrocarburos ya
que esto significa mayores ganancias y para poder llegar a la meta de mayor
producción se buscaron herramientas donde se puedan analizar todas las variables
que se pueden presentar en un yacimiento implementando modelos matemáticos
que cubran las diferentes características particulares de cada yacimiento, en estos
programas de simulación podemos hacer distintas pruebas para encontrar la mejor
manera de explotar este yacimiento.
Planteamiento del problema
Se tiene un dominio unidimensional Ω (ver Figura 1), el cual representa a un medio
poroso totalmente saturado por una fase fluida monocomponente. Por tal, para
analizar la hidrodinámica de flujo en el dominio se requiere plantear un modelo
matemático y su respectivo modelo numérico, el cual posteriormente pueda ser
implementado (modelo computacional) y resuelto numéricamente utilizando
lenguaje de programación (FORTRAN). El balance de masa para flujo monofásico
mono componente en un medio poroso está dad o de la siguiente manera:
𝝏(𝝓𝝆𝒐)
1
𝝏𝒕 + 𝛁 ∙ (𝝆𝒐𝒖𝒐) = 𝒒𝒐
Donde 𝜙 es la porosidad, 𝜌"y 𝑢"son la densidad y la velocidad de la fase oleica,
respectivamente, y 𝑞" es el término fuente/sumidero. En la ecuación (1) la velocidad
de la fase oleica está dada por la ley de Darcy:
𝑲
2
𝒖𝒐 = − 𝝁𝒐 (𝛁𝒑𝒐 − 𝝆𝒐𝐠𝛁𝒛)
aquí 𝑲 es la permeabilidad absoluta, 𝜇" y 𝑝" son la viscosidad y la presión de la
fase oleica, respectivamente; y 𝐠 es la aceleración de la gravedad.
Fig. I -
Los parámetros de simulación de la matriz porosa y del fluido que la satura son las
siguientes:
Porosidad (∅) = 0.2
Densidad del fluido (ρ) = 1000 𝑘𝑔/𝑚3
Viscosidad cinemática (𝜈) = 1𝐸 − 5 𝑚2/𝑠
Compresibilidad total ( 𝐶𝑇 ) = 20𝐸 − 6 𝑝𝑠𝑖-1
Longitud del dominio (𝑙) = 400 𝑓𝑡
La condición inicial y las condiciones de frontera que se aplican al dominio Ω son
las siguientes:
Condición inicial
donde 𝒑𝟎 = 2000 psi, y representa la presión inicial en el dominio Ω
Condiciones de frontera
En 𝑥 = 0
En 𝑥 = 𝑙
Aquí representan el valor de la variable en cada frontera izquierda (en 𝑥 = 0) del
dominio Ω.
Objetivos del ejercicio
Los objetivos del ejercicio son de carácter representativo para desarrollar un
modelo unidimensional, para entender los fundamentos que se dieron en el
presente trabajo así como hacer uso de las bases que se plantearon en el trabajo;
se ocuparon los principios básicos tales como son las características del
Yacimiento, el significado mismo de un Yacimiento. Una vez conocidas las
características del yacimiento, se iniciará con la selección más acorde con el
Yacimiento es decir seleccionamos el modelo que se simulara con la Ruta crítica de
la SNY. Conocido el modelo a simular se desarrolla la EDP del ejemplo en su forma
de diferencias finitas para dar paso a su solución con este método, conociendo que
toda EDP necesita condiciones iniciales y de frontera
Desarrollo
Puntos a desarrollar:
Plantear un modelo matemático y numérico para el modelado de flujo monofásico
monocomponente ligeramente compresible en un medio poroso indeformable, con
porosidad, viscosidad y permeabilidad constante.
a. El modelo numérico debe plantearse mediante el uso del método de
diferencias finitas, utilizando un esquema implícito.
b. La malla numérica debe ser en no contacto con la frontera.
Modelo matematico
Para plantear el modelo matemático es necesario primero plantear las siguientes
suposiciones para un modelo tener un ánalisis más sencillo:
Medio poroso indeformable
No hay generación de masa
Flujo darciano
Flujo en condiciones isotérmicos
Equilibrio térmico y termodinámico
El efecto de la gravedad no se considera(g=0)
Flujo unidimensional
Porosidad, viscosidad y permeabilidad constante
Ligeramente compresible
Flujo monofásico y monocomponente
Podemos partir de la ecuación de difusión:
Trabajando solo en el eje x ya que solo se trata de una dimensión podemos
determinar que :
Sabiendo que la velocidad de darcy para un flujo monofásico se define por:
Se sustituye la ecuación de darcy en la ecuación, y despreciamos la gravedad, por
lo tanto:
Tomando en cuenta que el diferencial de presión es sobre el posicionamiento, esta
se puede escribir de la siguiente manera
Factorizando / porque k=cte y 𝜇 =cte , por lo que tenemos que:
Se deriva la ecuación en función a la posición que están, aplicando la siguiente
formula (UV):
Siendo U la densidad y V la derivada parcial de la presión con respecto a su
dirección, queda de la siguiente manera las derivadas:
Quedando de la siguiente manera la ecuación:
Simplificando y factorizando la densidad:
Aplicando la “Regla de la cadena”, se tiene que:
Considerando que la compresibilidad está dada por
Y que la expansión térmica está dada por
Se multiplica nuestra ecuación por
Por lo tanto, queda de la siguiente manera:
Se deriva la ecuación :
Y luego utilizamos la regla de la cadena, por lo tanto, la ecuación nos queda:
Conociendo las ecuaciones de Compresibilidad y Expansión térmica, multiplicamos
la ecuación por y hacemos la correspondiente sustitución, así que la ecuación se
rescribe de la siguiente manera.
Se procede a sustituir en la ecuación y nos queda de la siguiente manera
Como porosidad es constante y tenemos poca compresibilidad y como el fluido es
isotérmico y considerando que no hay transferencia de masa
Se multiplica la ecuación por nos quema:
Sin embargo, se deben usar unidades de campo, las cuales son:
Por lo tanto:
Sin embargo:
Por lo tanto, se debe convertir a [𝑝𝑠𝑖/𝑓𝑡8].
Sabiendo que:
Entonces:
Por lo tanto, para usar la ecuación en unidades de campo, se debe usar la
constante, teniendo que:
Modelo numérico
Partiendo del modelo matemático para unidades de campo:
Considerando
Sabiendo que las ecuaciones de diferencias finitas nos dicen que, la derivada de
tiempo es igual a:
La derivada espacial está dada por:
Por lo tanto:
Entonces:
Reorganizando los términos, se tiene que:
Factorizando
Considerando
Nos queda:
Multiplicamos por -1 y reordenamos:
2. Utilizando lenguaje de programación FORTRAN, desarrollar un código
numérico que permita obtener la solución del problema descrito previamente.
Se conocen los datos que se deben introducir a FORTRAN ya que son los que
se mencionan anteriormente, con excepción de los siguientes:
• Convertimos la viscosidad cinemática a dinámica:
• Obtenemos el valor de C2 resolviendo la derivada parcial, que es 0 porque
el t es constante.
Para este caso de estudio se usará un total de 9 nodos (NX= 9), un DT= 2 días
y un criterio de convergencia de 1E-6.
A criterio propio, opté por utilizar:
• NX= 20
• DT=3
• Criterio de convergencia de: 1E-5
Ahora bien, para la elaboración del código se debe hacer considerando el
siguiente diagrama de flujo (Figura 2)
3. Obtener la solución transitoria del problema hasta alcanzar el estado
permanente.
4. Elaborar reporte.
Diagrama de flujo
INICIO
Parámetros de entrada
Generación de malla computacional
Generación de malla
Condición
Condición inicial
𝑦
Adivinar el valor de Px0 y Pxold
Cálculo de coeficientes y condiciones de
• frontera
•
Solver (Px)
𝑃𝑥 )
=𝑃OLD
Pxold=𝑃𝑥
Px
Criterio en paso del tiempo
No
12
𝜀 = 10
𝟎 𝟏
𝒑
𝑃0 =
= 𝒑𝑃1 SI
NO Criterio en
estado
permanente
𝜀 = 10’-12
SIsi
Imprimir resultados ados FIN
Fig. II - Diagrama de flujo.
PROGRAM
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Programa para resolver la Ecuación de Difusividad unidimensional.
DIMENSION a(105),b(105),c(105),d(105),p(105),pn(105) ! SE CREA UN ARCHIVO PARA GUARDAR LOS
RESULTADOS OBTENIDOS open(12,file='[Link]')
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!características del yacimiento
pini = 1000. xlong= 1500. c1=1000. c2=1000. area=100 N=16 delx = xlong /(N-1) 128 poro=0.15 perm=15
visc=20 bo=1.118 comp=35.E-6 tiempo=90 delt=15 qp= 3500 qi=2800
alfa=(poro*visc*comp*delx**2)/(0.00633*(perm*delt)) q1=qp*visc*bo/(perm*delx*area)
q2=qi*visc*bo/(perm*delx*area)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! CONDICIONES INICIALES DO 10 i=1,N p(i)= pini pn(i)= pini 10
CONTINUE
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! PUNTOS INTERNOS DO 20i=2,N-1 a(i)=-1 b(i)=2+alfa c(i)=-1 20 CONTINUE
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !
CONDICIONES DE FRONTERA
b(1)=1. c(1)=0. d(1)=c1 a(N)=0. b(N)=1. d(N)=c2
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !
COMIENZAN LOS PASOS DE TIEMPO
t=0
30 CONTINUE
t=t+delt
DO 40 i=2,N-1
d(i)=pn(i)*alfa
! SE INTRODUCEN LOS GASTOS DE LOS POZOS EN EL ALGORITMO PARA UN POZO PRODUCTOR Y UNO
INYECTOR EN
&X=500, X=1000
d(6)=d(6)-q1 d(11)=d(11)+q2 40 2
El algoritmo descrito es para un arreglo de 16 nodos y para un arreglo con más
celdas es necesario modificar el número de nodos, se modularon las operaciones
para poder modificar cada una de las condiciones si así se desea. Los resultados
son guardados en un archivo de datos WordPad que se muestra a continuación.
Resultados
Para verificar que el código es correcto, realizar algunas pruebas, por lo que la validación se demuestra
según lo siguiente:
Con el código FORTRAN proporcionado por el profesor, se procede a correr la simulación con los mismos
valores que se tienen en el libro y los graficamos, al ser los mismos datos y mismas escalas, las gráficas
tienen que ser las mismas, con lo cual se comprobará el correcto funcionamiento de este.
–
interpretación y solución
Se realiza el gráfico de la figura 5 con los datos obtenidos a partir de las iteraciones
de 20 nodos con un paso de tiempo de 3 días, al graficar cada perfil de presión a
distintos días se puede observar que a un mayor número de días los perfiles de
presión toman un comportamiento lineal, siendo a partir del día 60 que la diferencia
es prácticamente nula, por lo tanto se puede decir que la presión en el yacimiento
es independiente de la posición, con lo cual, de esta manera, el sistema llega a un
estado permanente. Sin embargo, se opta por considerar que llega al estado
permanente a los 60 días y no hasta el último día de la simulación ya que la presión
cambia muy poco, por lo que al tener un resultado prácticamente igual, puede
resultar ser una mejor opción escoger un tiempo de 60 días para ahorrarle trabajo
al programa, es decir, para optimizar.
.
En el gráfico de la figura 6 se puede observar que los perfiles con un paso
temporal (DT) de 0.375, 0.1875 y 0.09375 se sobreponen uno sobre el otro, lo cuál
se puede interpretar que estos son independientes del paso temporal, por lo que
una mejor opción es cambiar el paso temporal (DT) a un valor de 0.09375, ya que
de esta manera lso resultados serán más precisos, sin embargo, al ser
prácticamente igual a DT= 0.375, se decide por DT= 0.375 debido a que requiere
menor trabajo computacional y el valor de la presión es muy similar, es decir, para
optimizar el proceso..
la figura se puede observar que las aproximaciones convergen en un paso
temporal (DT) de 0.375 para todas las configuraciones del número de nodos (NX),
por lo que se opta por cualquiera de ellas, sin embargo, lo conveniente, esto según
optimización, es usar un dx= 21.05 (NX= 20) y así observar que la presión para un
paso temporal de DT=0.375 tiene una variación mínima, pues esta apenas varía 10
psi a comparación de la presión para un DT= 0.09375, así que es más conveniente
optar por el DT= 0.375, para así de esta manera obtener que la solución sea
independiente del tamaño de malla y a su vez independiente del paso temporal. En
conclusión, se puede optar por un DT= 0.375 y un dx= 21.05 (NX=20) con el fin de
que sea la configuración más adecuada para realizar la simulación y llegar a la
solución exacta o al menos a la más aproximada.
la grafica de la figura se observa que, mientras el paso temporal (DT) sea mayor,
la presión continuará variando, esto de acuerdo a la posición, sin embargo, de igual
manera se puede observar que los últimos tres perfiles, es decir los que utilizan un
paso temporal de 0.375, 0.1875 y 0.09375, los cuáles son DT más pequeños, son
quienes tienen una pendiente prácticamente nula, es decir, de 0, lo que significa
que la presión no cambia o lo hacen en menor medida sin importar la posición, , es
decir, es independiente de la posición; de igual manera, como en el gráfico que
muestra la figura 6, se puede confirmar que utilizar un paso temporal de (DT) 0.375
es lo más adecuado.
el gráfico de la figura , se compraran los perfiles de presión para la configuración
del paso temporal (DT) y DX inicial, así como el paso temporal (DT) y DX que
resulta ser el optimo, y así poder decir que de manera general, la segunda opción
del perfil de presión es mejor debido a la excatitud que esta tiene, o mejor, dicho la
aproximación es mayor.
IX -
días.
Finalmente con la gráfica que muestra la figura 10, se observa el
comportamiento que tiene la presión, es decir, la variación que tiene la presión para
los distintos perfiles variando a su vez los días hasta llegar a t=60, justo cuando el
sistema alcanza el estado permanente.
Fig. X - Gráfica que muestra la comparación de la presión respecto al
tiempo a un x =21.05.
Conclusión
Para hacer el modelado de un Yacimiento Petrolero es necesario realizar
un análisis previo, el cual permita determinar la confiabilidad de la
información con la que se hará el modelado, ya que esta información
brinda el tipo de modelado que se construirá, por lo que los resultados
dependen de la validez de la información. Cuando se dispone de poca
información, la construcción de un modelo simple resulta ser la opción
más viable, pues al carecer de información, está puede ser inferida o
correlacionada con algún campo vecino, por lo que no es una opción
construir un modelo detallado que no sea del todo representativo del
Yacimiento que se desea estudiar.
Este ejercicio tiene como fin ser una herramienta académica para poder
familiarizarse con el simulador y por esta razón el nivel de complejidad es
inferior a lo que realmente se presenta en el campo laboral sin embargo
en el ejercicio pudimos notar que entre más es el número de nodos mas
cercanos estamos de la realidad del yacimiento, pero llega un punto en
donde aun que aumentemos la cantidad de nodos estos ya no afectaran
a nuestro sistema por lo que lo más recomendable es identificar el
numero optimo de nodos que podemos emplear esto dependiendo de las
variables que utilicemos en nuestra simulación como lo son las
dimensiones de las mallas del caso que estamos estudiando y un paso
temporal independiente.
.
Bibliografía
Roberto C. M .C. (2014). T E S I S QUE PARA OPTAR POR EL GRADO DE: MAESTRO EN CIENCIAS
DE LA TIERRA “Simulación Numérica de Inyección de Agua en Yacimientos Petroleros
empleando el Método de Líneas de Corriente” UNAM. 2014.
HOSUÉ ALBERTO ROSALES RANGEL, (2009); T E S I S, QUE PARA OBTENER EL GRADO DE
INGENIERO PETROLERO “MODELOS SIMPLES DE SIMULACIÓN NUMÉRICA PARA SOLUCIONES
PRÁCTICAS EN LA INGENIERÍA DE YACIMIENTOS” UNAM 2009
Victor Hugo Arana Ortiz, David Trujillo Escalona, Juventino Sánchez Vela, (2003); “APUNTES DE
SIMULACIÓN NUMÉRICA DE YACIMIENTOS” UNAM 2003
