INTRODUCCIÓN

El transporte de crudo en el Ecuador constituye una actividad imperiosa, sobre todo si se considera que casi la totalidad del crudo proviene de la región amazónica. El país cuenta con dos oleoductos: el Oleoducto de Crudos Pesados (OCP), que es privado, destinado a la transportación de crudo pesado; y, el SOTE, que es un oleoducto de propiedad del estado ecuatoriano, diseñado en los años 70 para transportar crudos livianos. En la actualidad, surge la necesidad de buscar alternativas que permitan la transportación de crudos pesados por la infraestructura del SOTE, para lo cual se conocen dos opciones: la primera que consiste en calentar el crudo para disminuir su viscosidad, y la segunda que es la de diluir el crudo pesado con otro de menor viscosidad o con diluyente. De otro modo, recurrir a la experimentación física de este método de transporte en el SOTE no es posible, debido a los altos riesgos de ruptura del oleoducto, lo que implica serios problemas de contaminación ambiental y su posterior remediación; lo cual afecta la economía y la imagen del país. Con estos antecedentes, disponer de un simulador matemático que permita experimentar en el computador el transporte de crudos pesados por el SOTE, se vuelve imprescindible. Las caídas de presión, los cambios de temperatura y los esfuerzos que se producen en la tubería, se simulan mediante la ejecución de programas computacionales elaborados en base a algoritmos desarrollados de complicados métodos numéricos de aproximación de soluciones de ecuaciones en derivadas parciales. Esto es factible gracias a la modelización matemática fundamentada en las leyes de la mecánica de fluidos, la hidráulica de tuberías y bombas centrífugas; así también, a los avances permanentes del análisis numérico de ecuaciones en derivadas parciales y al procesamiento digital de imágenes.

METODOLOGÍA

La construcción del modelo matemático se ha fundamentado teóricamente en las ecuaciones de movimiento y de continuidad; la primera determina la condición de un objeto en función de la velocidad, aceleración, masa y otras variables-función que afecta su movimiento; es decir define la evolución temporal de un sistema físico en un espacio; y la segunda expresa la conservación de la masa de manera numérica y en forma diferencial (aplicada al caso de investigación).

Un modelo matemático de flujo inestable e isotérmico en una tubería, se obtiene mediante la aplicación de la segunda ley de Newton, la conservación de masa, leyes de comportamiento y de información experimental. La dinámica del fluido en la tubería se expresa en términos de presión p (x, t) y velocidad v (x, t) que son funciones de la posición y del tiempo. Los modelos que se obtienen, se encuentran basados en ecuaciones diferenciales unidimensionales, que son simplificaciones de modelos tridimensionales, estas se aproximan, usando un esquema numérico, con el cual se elabora un algoritmo computacional que simula el transporte del crudo por el SOTE.

Los objetivos obedecen al problema de investigación y son los siguientes:

1. 1. Diseñar un modelo matemático basado en las ecuaciones de movimiento y de continuidad.
2. 2. Aproximar la solución del modelo matemático, mediante el esquema numérico fundamentado en diferencias finitas centrales.
3. 3. Elaborar un programa computacional, que bajo determinadas condiciones como las de viscosidad del crudo, caudal, condiciones iniciales de bombeo y estructurales del Oleoducto, permita simular el transporte de crudo a través del SOTE.

El proceso constructivo del modelo matemático se encuentra descrito en una metodología sistemática, lógica y deductiva, que refleja su integralidad y alcance.

ECUACION DE MOVIMIENTO

Se considera un volumen de control Ω (t) constituido por una porción de fluido en la tubería de sección transversal de área A y longitud δx. En este dominio se toman dos puntos referenciales notados con 1 y 2 de Ω(t), como se muestra en la figura 1.

Las fuerzas que obran en este elemento de fluido son: F₁ fuerza ejercida sobre la sección 1 debida a la presión, esto es, F₁=ρA donde ρ es la presión y A es el área de la sección transversal; F₂ es la fuerza debida al peso del fluido, la misma que se expresa como F₂=mg=ρAδxg, donde m es la masa del fluido, δ su densidad y g=9,8m/s² el valor de la gravedad. Interesa su componente paralela al eje central de la tubería, que se le nota F_2,x y que está definida como F_2,x =ρAδxgsen(θ), donde θ es el ángulo de inclinación del tubo respecto de la horizontal. A continuación, la fuerza F₃ ejercida sobre la sección 2 debida al cambio de presión respecto de la posición, definida como F₃=pA+ ∂/∂x(pA) δx. La fuerza debida a la fricción F₄ = Γ₀ πDdx, donde Γ₀ es el esfuerzo cortante, D el diámetro.

Por la fórmula de Darcy-Weisbach , con f el coeficiente de fricción. La fuerza ejercida por variación de la sección transversal F₃ está definida como . Para mayor información referirse a [2], [6], [7].

Por la segunda Ley de Newton, se tiene:

Remplazando los resultados precedentes, se obtiene la ecuación siguiente:

Luego de algunas simplificaciones, la ecuación precedente se escribe como

Como

Reemplazando este resultado en la igualdad precedente y dividiendo para Aδx, resulta

Para δx suficientemente pequeño, el término es suficientemente pequeño por lo que es despreciable, y la ecuación precedente se reduce a:

Además y dividiendo para ρ, se obtiene la siguiente ecuación de Euler:

En la práctica la ecuación (1) es reemplazada por la ecuación 2 ([2], [6], [7], [8]):

donde 1<λ≤2

En el caso de flujo estable, ∂v/∂t= 0, con lo que la ecuación reduce a (Ec. 3):

Adicionalmente, si se considera la densidad ρ y el área A constantes, lo que implica ∂v/∂x= 0, y de esta resulta que la velocidad v es constante. La ecuación (3) se transforma en

de donde: con Δx = x₂ – x₁

Esta ecuación es equivalente a la ecuación de Darcy-Weisbach ([2], [6], [7], [8]): en la que se considera el perfil del oleoducto. La ecuación (4) es la ecuación de flujo de estado estable generalmente utilizada como la condición inicial en los modelos de flujo inestable e isotérmico.

En la figura 2 se muestra la cabeza piezométrica H o línea de nivel hidráulico sobre un dato arbitrario, que puede reemplazarse por la presión p=ρg(H-z), donde z es función de X y representa la elevación del eje central del tubo. Resulta

Si se considera que v ∂v/∂x≈ 0 que es válido para flujos inestables de baja velocidad, la ecuación (1) se escribe como

Remplazando Ec. (5) en la ecuación precedente, se obtiene

De esta ecuación, se obtiene la ecuación 6

La ecuación (6) es válida para fluidos poco compresibles que fluyen a baja velocidad.

ECUACIÓN DE CONTINUIDAD

Se considera nuevamente el volumen de control Ω(t) y denotamos con Γ(t) su frontera. En Ω(t) fluye un fluido de densidad ρ. Se tiene

Luego, por el principio de conservación de la masa dm/dt=0. Nótese que el volumen de la porción de fluido contenida en Ω(t) es A(x,t) dx, siendo A(x,t) el área de la sección transversal que es función de X y de t. Resulta

y por el teorema de transporte de Reynolds, tenemos

donde n(s,t) es el vector normal exterior a Γ(t) y v(s,t) es la velocidad del fluido. Por el teorema de la divergencia de Gauss (ver [2], [15]),

Luego

Como se privilegia el movimiento del fluido en la dirección del eje del tubo, los componentes en las direcciones ortogonales a esta se consideran nulas, o sea v(x,t)=(v(x,t),0,0). La ecuación precedente se escribe como (Ec. 7):

que es la ecuación de continuidad unidimensional, la misma que se la escribe como sigue

dividiendo la ecuación precedente para ρA, se deduce la siguiente:

El primer término de la ecuación (8) conduce a la elasticidad de la pared del tubo y su tasa de deformación con la presión, luego

El término dp/pdt toma en consideración la compresibilidad del fluido. La presión se introduce a través de la definición del módulo de elasticidad del fluido. Así:

entonces k dp/p=dp

de donde dp/dt=kdp/pdt

Reemplazando estos dos últimos resultados en (8) se tiene

realizando las simplificaciones necesarias se reduce a (Ec. 9):

donde .Tomando en consideración que se sigue que:

La ecuación (10) se conoce como ecuación unidimensional de conservación de masa para fluidos ligeramente compresibles. En el caso de flujo estable ∂p/t =0 con lo que (10) se escribe como (Ec. 11):

Para flujos con velocidad pequeñas v∂p/∂x=0 con lo que (10) se reduce a la siguiente:

y tomando en consideración que, se sigue

Esta última ecuación conjuntamente con la de Euler, forman el modelo matemático conocido como de la columna elástica que se utiliza para simular el conocido problema de golpe de ariete (para más detalles referirse a [2], [3], [4], [6], [7], [8]). En el caso de estudio, este problema no será tratado debido al interés de modelos de flujo completamente desarrollados, incompresibles en los que la elasticidad de la tubería no es considerada o dicho de otro modo es de efecto despreciable. Para el SOTE, este problema en principio, no es prioritario. Para una segunda fase este modelo tiene mucho interés fundamentalmente por vincularse con el sistema SCADA.

EXISTENCIA Y UNICIDAD DE LA SOLUCIÓN DEL MODELO ESTACIONARIO E ISOTÉRMICO

El modelo estacionario de flujo incompresible y completamente desarrollado en una tubería está gobernado por el siguiente sistema de ecuaciones diferenciales:

donde la velocidad v y la presión p son las variables a analizar, que en este modelo intervienen y son funciones definidas en [L₁,L₂], donde xЄ[L₁,L₂] diámetro interno D(x) y consecuentemente el área de la sección transversal es densidad del fluido p(x), perfil del oleoducto dz(x)/dx y , viscosidad cinemática v(x), factor de fricción f, donde la función se define por:

Nótese que se define a f₁, f₂ son los factores de fricción a 2000 y 4000 Reynolds respectivamente. En este intervalo se tiene la transición de flujo laminar a turbulento. El número de Reynolds está definido como

La fricción

El diámetro interno del oleoducto D, la densidad del fluido δ, la viscosidad cinemática v son funciones definidas en [L₁, L₂]. estrictamente positivas, constantes a trozos y acotadas, es decir que D, p, v ∈ H¹ (L₁, L₂). Luego existen D_min, D_max, p_min, p_max, v_min, v _max ∈ R⁺ tales que para todo x∈ (L₁, L₂),

Además v₀>0. De la ecuación

se sigue que:

que muestra que y en consecuencia

De la ecuación

Analicemos cada término:

1)

2)

3) La función factor de fricción es una función positiva y acotada,

Por lo tanto,

En conclusión, las funciones v, p ∈ H₁ (L₁, L₂) arriba definidas son solución del sistema de ecuaciones diferenciales. Es muy conocido en los cursos de análisis que estas funciones son la única solución del sistema de ecuaciones diferenciales propuesto.

APROXIMACIÓN NUMÉRICA (POSICIÓN DEL PROBLEMA)

Sean L₁>0, L₂>0, tales que L₁<L₂. La ecuación de continuidad se escribe como

y la ecuación de Euler

En el caso estacionario, ∂/∂t (ρA)=0 y ∂v/∂t=0. Luego, las dos ecuaciones precedentes se escriben como el sistema no lineal de ecuaciones diferenciales

donde v y p son las funciones incógnitas definidas en [L₁; L₂]. Suponemos que estas funciones son de clase C¹.

DISCRETIZACIÓN

El modelo matemático se discretiza con diferencias finitas centrales.

Sea n∈Z⁺ y τ(n)={L₁=x₀<⋯<x_n=L₂} una partición de [L₁;L₂], donde x_j-1<x_j , j=1,2,⋯,n. Ponemos ∆x_j=x_j-x_j-1, ∆x=(L₂-L₁)/n. En el caso de mallado uniforme, tenemos ∆x=∆x_j; j=1,2,⋯,n. Se designa con ρ_j, A_j, v_j, p_j aproximaciones de ρ(x_j), A(x_j), v(x_j), p(x_j), j=1,2,⋯,n.

El esquema numérico resultante se escribe como sigue:

donde

que se denomina sistema de ecuaciones discreto, donde v_j y p_j son las incógnitas, j=1,2,...,n. Con el propósito de facilitar el cálculo, se define. Entonces

De la primera ecuación del sistema discreto obtenemos

Luego:

En esta última ecuación se calculap_j,j=1,2,…

que permite calcular la presión en cada punto x_j. Como el área de la sección transversal de la tubería es entonces con D_i el diámetro interno que se escribirá D_j en x_j, el esquema numérico precedente se escribe como:

CÁLCULO DEL COEFICIENTE DE FRICCIÓN. MÉTODO DE NEWTON

La ecuación de Colebrook ([1], [2], [9]) relaciona el coeficiente de fricción f con el número de Reynolds la rugosidad relativa ε/D mediante la ecuación:

Donde D denota el diámetro interno de la tubería, y ε es la rugosidad que para oleoductos cumple con normas internacionales (0.045 mm para acero comercial). Se observa que el factor de fricción f no puede calcularse directamente, por lo tanto, es necesario recurrir a los métodos numéricos, particularmente el método de Newton. La ecuación de Colebrook se expresa en forma equivalente como la siguiente:

Se define la función θ de ]0, 1[ en R como:

Entonces

El método de Newton genera una sucesión f_n definida como:

Una aproximación inicial de f₀ se obtiene de la ecuación de S. E. Haaland:

La aproximación inicial f₀ calculada con la ecuación de Haaland difiere a lo mucho en aproximadamente 5% con el valor real, por lo que el método de Newton converge en pocas iteraciones para la precisión fijada ε . En la práctica ε=10^-10. Este algoritmo es aplicado cuando el número de Reynolds Re≥4000. Cuando 0<Re≤2000, el factor de fricción se calcula como f=64/Re. Tomando en consideración que Re=ρv/v y que son funciones dependientes de la posición x en el oleoducto, en cada punto x del oleoducto debe calcularse f en la forma propuesta. Para el caso térmico, la viscosidad es también función de la temperatura θ. Escribiremos al número de Reynolds Re como Re=(x,ρ(x),v(x),v(x,θ)).

MODELOS DE TRANSPORTE DE PETRÓLEO POR TUBERÍAS PARA MEZCLAS

La parte principal en el problema de transporte de crudos pesados en un oleoducto es el control sobre la viscosidad de este fluido. Uno de los métodos de transporte de crudos pesados es mezclar con otro crudo más ligero llamado diluyente, de modo que dé como resultado un crudo intermedio que tenga una densidad y viscosidad razonables. Con este propósito se determina la fracción de dilución, la densidad de la mezcla, así como la viscosidad a la temperatura de operación.

Sean m_d la masa del diluyente, m_c la masa de crudo (pesado) y se designe con m_m a la masa de la mezcla. En lo que sigue los subíndices o superíndices d, c, m se referirán al diluyente, al crudo y a la mezcla. Tenemos entonces que m_m=m_c+m_d=ρ_cV_c+ρ_dV_d, donde ρ_c, ρ_d son las densidades del crudo y del diluyente, y V_c, V_d son los respectivos volúmenes. Las fracciones de crudo y de diluyente se designan con χ_d y χ_c , definidas como:

con V_m=V_d+V_c el volumen de la mezcla. Resulta χ_d+χ_c=1 de donde

Con 0<χ_d<1, 0<χ_c<1. Se define ρ_m=ρ_d χ_d+ρ_c χ_c. Luego

Esta relación será utilizada sistemáticamente en el modelo de transporte de crudos con diluyente. En términos de masas, se tiene:

La gravedad específica de un fluido se designa con S_pgr (specific gravity) y se define como S_pgr=ρ/ρa , con ρa la densidad del agua a 60ºF. Entonces

y de esta relación obtenemos

Puesto que

la ecuación precedente se escribe como

de donde

Esta es la ecuación de la densidad de la mezcla en unidades de campo a utilizarse en el modelo. Esta expresión es corregida con la temperatura si se toma en consideración que

donde T_oper es la temperatura de operación.

Para el cálculo de la viscosidad de la mezcla se aplica la ecuación ASTM ([3], [4]) para la viscosidad cinemática a la temperatura de operación medida en cst. Para el efecto se asume que la viscosidad del crudo es V_c a una temperatura T_c, la viscosidad del diluyente es V_d a la temperatura T_d. Resultados experimentales muestran que la viscosidad de la mezcla V_m es:

Esta viscosidad es la que se incorpora en el modelo de transporte de crudo por tuberías. Para el caso específico de petróleo del oriente ecuatoriano se dispone de la siguiente función para el cálculo de la viscosidad de una mezcla de crudos de diferentes densidades y viscosidades:

La investigación desarrollada, constituye un invaluable aporte a la ciencia y una contribución de alto impacto para el estado ecuatoriano, la comunidad y el medio ambiente, por efecto de las siguientes razones:

• Posibilita simular el transporte de crudo pesado a través del SOTE, de lo cual se deriva, un importante ahorro económico para el Estado Ecuatoriano.

• El impacto medio ambiental se reduce notablemente, dado que las pruebas se las efectiviza en el simulador y no en la estructura física del SOTE.

• El estudio coadyuva a la generación de aprendizajes significativos, consolida el conocimiento y fortalece los principios de la Comunidad Científica.

RESULTADOS

Se desea transportar un caudal de 200000 BPD de petróleo de 13 ºAPI, cuya viscosidad a 80 ºF es de 80 cst y a 120 ºF es de 55 cst, usando como diluyente un crudo de 29 ºAPI, cuya viscosidad a 80 ºF es de 25 cst y a 120 ºF es de 13 cst a condiciones de bombeo iniciales de temperatura 110 ºF, rugosidad absoluta de 0.0045 cm y una presión de succión de 120 psi.

Si se desea que el API de la mezcla sea 16.5 ºAPI a 60 ºF, el programa calcula la fracción de diluyente que se tiene que emplear, en este caso es de 0.237, lo que significa un caudal de diluyente de 62200.7 BPD, lo que da como resultado un caudal de la mezcla de 262200.7 BPD.

La Figura 3, muestra las curvas de viscosidad de los fluidos en función de la temperatura, donde se puede apreciar, que a mayor temperatura, la viscosidad disminuye. El simulador permite configurar las estaciones de bombeo, como por ejemplo la estación de Lumbaqui, que cuenta con seis bombas, cada una con seis etapas, funcionando a 3400 RPM.

La Figura 4, muestra la curva que corresponde a una sola etapa de una bomba para agua. El simulador corrige automáticamente la curva para la mezcla a ser transportada mediante un modelo matemático.

La Figura 5, muestra el cálculo del punto óptimo de operación de la bomba, donde se puede visualizar que la bomba opera en condiciones óptimas.

En la figura 6, se presenta el perfil del oleoducto en color verde, en color rojo se tiene la máxima presión admisible (M.P.A.), mientras que en color negro, se representa el gradiente hidráulico a lo largo del oleoducto. Es claro que el gradiente hidráulico tiene que estar entre la M.P.A. y el perfil del oleoducto. Si la curva del gradiente hidráulico está por encima de la M.P.A. la tubería se rompería, mientras que si se encuentra por debajo del perfil del oleoducto, significaría que la presión es insuficiente para que el fluido se mueva. En este caso se puede observar que es factible transportar un caudal de 200000 BPD de 13 ºAPI usando 62200 BPD de diluyente de 29 ºAPI.

CONCLUSIONES Y RECOMENDACIONES

Mediante la modelización matemática y la simulación numérica, se generan escenarios virtuales para ver si es factible transportar un determinado caudal de crudo pesado usando un diluyente, desde la estación de Lago Agrio, hasta Balao, sin correr el riesgo de tener un colapso de la infraestructura.

El diseño del modelo matemático isotérmico, se ajusta a las condiciones de viscosidad del crudo, caudal, condiciones iniciales de bombeo y estructurales del SOTE.

El esquema numérico de diferencias finitas, ha permitido hallar una solución aproximada del sistema de ecuaciones diferenciales, que permite identificar la presión y la velocidad en cualquier punto del Oleoducto.

Los resultados obtenidos mediante el simulador, han sido comparados con los resultados de presión y velocidad, arrojados por el Sistema SCADA (controla y supervisa el proceso de transporte de crudo a lo largo del oleoducto), presentando un margen de error relativo menor al 5%.

Se recomienda investigar y desarrollar modelos matemáticos térmicos para el transporte de crudo pesado, que permita simular el transporte de crudo por el oleoducto, añadiendo estaciones de calentamiento virtuales, que en el caso de determinarse su factibilidad en la simulación, se pueda sugerir su construcción.

BIBLIOGRAFÍA

[1] Bergada, J. M. (2015). Mecánica de Fluidos. Iniciativa Digital Politécnica, Universidad Politécnica de Cataluña, España.

[2] Streeter, V., Wylie, B. & Bedford, K. (2000). Mecánica de Fluidos, Novena Edición, Editorial McGraw-Hill, Bogotá.

[3] Mohitpour, M., Golshan, H., & Murray, A. (2003). Pipeline Design & Construction. A Practical Approach, Second Edition, Editorial American Society of Mechanical Engineers ASME Press, New York.

[4] Mohitpour, M., Szabo, J. & Van Hardeveld, T. (2005). Pipeline Operation & Maintenance. A Practical Approach, Editorial American Society of Mechanical Engineers (ASME Press), New York.

[5] Pozrikidis, C. (2017). Fluid Dinamics, Editorial Springer, Berlín.

[6] Huilgol, R. (2015). Fluid Mechanics of Viscoplasticity, Editorial Springer, Berlín.

[7] Rodriguez, L., Vega, L. (2013). Cálculos Térmicos e Hidrodinámicos de un oleoducto, Editorial Académica Española.

[8] Larock, B., Jeppson, R. & Watters, G. (2000). Hydraulics of Pipeline Systems, Editorial CRC Press, Boca Raton.

[9] Benalcázar Gómez, H.(2008). Análisis Numérico, Preprinter, Quito.

[10] Wesseling, P. (2001). Principles of Computational Fluid Dynamics, Editorial Springer-Verlag, Berlín.

[11] Tannehill, J., Anderson, D. & Pletcher, R. (1997). Computational Fluid Mechanics and Heat Transfer, Editorial Taylor & Francis, Philadelphia.

[12] Hirsch, C. (2007). Numerical Computation of Internal & External Flows, Editorial Elsevier, Amsterdam.

[13] Chung, T. (2002). Computational Fluid Dynamics, Editorial Cambridge University Press, Cambridge.

[14] Teman, R. & Miranville, A. (2000). Mathematical Modeling in Continuum Mechanics, Editorial Cambridge University Press, Cambridge, UK.

Notas de autor

hbenalcazarg@uce.edu.ec

Información adicional

Cómo citar: Benalcázar Gómez, H. G., Albuja Proaño, G. A., Naula Reina, I. C., Carrillo Flores, R. A., & Izurieta Cabrera, C. F. (2017). Modelación y simulación numérica del transporte de petróleo por el Sistema de Oleoducto Transecuatoriano. FIGEMPA: Investigación y Desarrollo, 3(1), 22–30. https://doi.org/10.29166/revfig.v1i1.52