Resolución del modelo matemático para transporte de gas, considerando un sistema isotérmico.

Resolution of the mathematical model for gas transport, considering an isothermal system.

El proyecto de gasificar el país y ofertar gas natural a otros países (presentado por el gobierno nacional como plan económico y energético para los próximos años) exige el desarrollo de una red de gasoductos que resuelva el problema de desconexión gasífera que existe tanto en el país como en algunos países vecinos.Es un hecho bien establecido que el flujo en los gasoductos es inestable, (Thorley y col., 1987). Las condiciones son siempre cambiantes con el tiempo, no importa cual pequeño sea el cambio, (Zhouy col., 2000, Villarreal, 2009, Benjamin y col., 1978, Martins y col., 2017,Thorley y col., 1987, Jibrin Helma y col.,2018). Por tal motivo el diseño de tuberías para el transporte de gas natural bajo condiciones de operatividad eficientes y seguras, requiere del conocimiento de la caída de presión, cambios de temperatura (caso no isotérmico) y flujo de masa distribuida en todo el sistema.

Debido a la complejidad del modelo matemático, resulta conveniente desarrollar métodos numéricos para su resolución, y con esta información puede ser utilizada durante el diseño del gasoducto, definición de los manuales de operación y capacitación de los operadores del sistema. Bajo condiciones isotérmicas, las ecuaciones de conservación de masas, conservación de la cantidad de movimiento y la ecuación de estado del gas natural, (Zhou y col., 2000, Villarreal 2009, Martins y col., 2017, Vásquez, 2003, Benjamin y col.,1978)constituyen las ecuaciones que describe el flujo transitorio del gas natural en gasoductos.

Generalmente, se supone considerar flujo isotérmico, aplicado a fricciones de estado estable, y despreciando expansión o contracción de la pared de la tubería por efectos de la presión. Cabe destacar que, para la simulación del flujo transitorio del gas natural en gasoductos, se está considerando el término de inercia en la ecuación de cantidad de movimiento.

Sobre cada paso del tiempo el modelo matemático se comportará hiperbólico,(LeVequey col.,1990) y se resolverá usando el método de separación (Splitting) (Dyachenko y col., 2016); es decir, dividiéndolo en un modelo no lineal,homogéneohiperbólico y un sistema de ecuaciones diferenciales ordinarias, donde las soluciones exactas de este último se usaran como un operador, y así obtener la solución del sistema completo no lineal y no homogéneo.

El concepto del método de Disminución de la Variación Total (TVD),(Zhouy col., 2000, Harten1997), ha conducido a una gran gama de métodos numéricos de primer y segundo orden de exactitud, la cual es de gran utilidad para ámbitos teóricos y aplicaciones ingenieriles.En el presente estudio, se implementará un método (TVD) explícito de segundo orden de exactitud, en 5 nodos, y es formulado para resolver el sistema no lineal, homogéneo hiperbólico de las leyes de conservación, que describe el transporte transitorio del gas natural por gasoductos inclinados. Por otro lado, se debe prestar atención en las condiciones numéricas de frontera en la entrada y salida del gasoducto; para esto, usaremos el método de las características,(Zhou y col.,2000).

Además, necesitamos un método numérico adicional para el tratamiento de las variables conservativas de interés en el segundo y penúltimo nodo. Para este propósito usaremos el esquema explícito de Godunov, (Zhouy col., 2000,Godunov 1959), de primer orden de exactitud, en tres nodos. A continuación se presenta el modelo, conformado por las ecuaciones de conservación de masas, cantidad de movimiento y la ecuación de estado del gas en coordenadas cartesianas, en una dimensión y bajo la mecánica clásica las cuales representarán el transporte de gas en gasoductos lineales inclinados, considerando condiciones de borde de densidad al comienzo y masa de flujo al final del gasoducto.

Si a la figura1 le aplicamos estos conceptos físicos, se tiene el sistema de ecuaciones de Navier – Stokes, (Vásquez 2003.Benjamin y col.,1978):

Donde es el caudal másico, la densidad del gas, el área transversal del gasoducto, el esfuerzo cortante, el diámetro del gasoducto, la aceleración de gravedad, es la inclinación del gasoducto, la aceleración del fluido y el intervalo espacial. La ecuación representa la ley de conservación de masas, y la ecuación es la aplicación de la sumatoria de fuerzas en el volumen de control representado por la figura ,(Villarreal2009). De la ecuación se puede obtener

Por otra parte, si se considera la ecuación de Darcy- Weisbach,(Villarreal, 2009, White 2011),y se multiplica por el inverso del área y del intervalo de espacio, de la ecuación se puede deducir

Reescribiendo la derivada total de la velocidad en función del tiempo y acomodando convenientemente los factores se tiene

Considerando que las variables a estudiar serán y, y por la ley de gas semi-ideal, se puede reescribir el sistema de ecuaciones de la siguiente forma:

Donde es la densidad del gas, es el flujo de masa, es la presión, es la compresibilidad del gas, es la constante del gas transportado, es la temperatura del gas, es la distancia a lo largo del gasoducto, es el tiempo, es el factor de fricción, es el diámetro interno del gasoducto, es la aceleración de la gravedad, es la velocidad del sonido en el fluido y es el ángulo de inclinación del gasoducto. Observemos que en la ecuación, el valor absoluto fue agregado para considerar el sentido de la velocidad del fluido,( Villarreal 2009,. Martins y col., 2017, Benjamin y col., 1978, Wanga y col., 2018). El problema debe resolverse para el caso transitorio, utilizando los métodos numéricos expuestos anteriormente.

Los valores iniciales para cada punto del gasoducto serán obtenidos del caso estacionario, entonces del sistema de ecuaciones (6) tenemos

para todo x $\in$ [0,L]. Y además

Integrando sobre (0,$∆$x) a la ecuación diferencial de variables separables (19) se obtiene la solución general implícita, dada por:

Donde ;derivado la función (20) respecto a la variable $\rho (∆\chi )$ se obtiene la ecuación:

Se obtendrá el valor de ${\rho }_{∆\chi }$

A través de la ecuación (22), se pueden obtener los valores iniciales de la densidad a lo largo del gasoducto. A continuación plantearemos el uso del método Splitting al sistema de ecuaciónes (4), para separar en un sistema hiperbólico no lineal homogéneo y un sistema de ecuaciones diferenciales ordinarias. En la forma matricial, en sistema (4) puede reescribirse como sigue:

Donde

El cual puede separarse en el sistema de ecuaciones (24) Villarreal 2017]:

Para la aplicación del método, es necesario obtener la solución numérica a las variables conservativas con la ecuación (24. A), y posteriormente utilizar estos en la solución exacta de la ecuación (24. B) para obtener la solución final. Vemos que el sistema de ecuaciones (24.b) tiene la siguiente forma:

Agrupando términos semejantes y operando convenientemente, se obtiene

Luego, haciendo , entonces k = 1 si m $\ge$ 0 ó k = -1 si m <0. Además y $\omega$ =$\rho gsin\alpha$

Ahora se considerarán dos casos (Y siempre un gasoducto con inclinación positiva),el primero sera el caso en que m$\ge$ 0. La ecuación diferencial quedaría

Integrando la ecuación diferencial sobre el intervalo (0,$∆t$), tenemos que

Por otro lado, si m<0 y mediante pasos análogos al caso anterior, se obtiene la solución

Donde $\psi =e^{2\sqrt{\gamma \omega ∆t}}$

Se empezará formulando el método de las características para el primer y el últmo nodo. A continuación se presentan las ecuaciones (8) y (10) discretizadas, usando diferencias finitas progresivas sobre el prier nodo y regresivas en el ultimo, respectivamente:

Para los puntos intermedios, el modelo matemático no lineal homogéneo (24. a), será resuelto mediante el esquema de Harten explicito TVD de segundo orden de exactitud sobre cinco puntos. A continuación se presenta la matriz Jacobiana para las ecuaciones del gas dinámico en caso isotérmico y en una dimensión

Donde la velocidad promedio $u$ esta definida por:

Esta ecuación a menudo se llama rho-promedio o (Roepromedio) de la velocidad, (Zhou y col., 2000). Observemos que la matriz $\stackrel{\_}{A}$ es la matriz Jacobiana $\frac{af\left(u\right)}{au}$ $\overline{u_1}$=$\overline{u}$ -c, $\overline{u_2}$ $\overline{u}$ + c y $\overline{R_1}$= (1,$\overline{u}$ - c); $\overline{R_2}$ (1,$\overline{u}$ + c). Por otro lado, de la ecuación (12), se tiene que $∆U_{j+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ = (U₁, U₂) = (${\rho }_2$ $-{\rho }_1$, m₂-m₁)

Sustituyendo los valores de $\overline{R^k}$

$({\rho }_2$ $-{\rho }_1$,m₂-m₁) = ( ${\alpha }^1$+${\alpha }^2$,${\alpha }^1$($\overline{u}$ $-c$)+$a^2$($\overline{u}$+c)

resolviendo este sistema, obtenemos:

Sustituyendo la ecuación (33), los auto valores y autovectores; además las ecuaciones (35) y (36) sobre el esquema TVD de Harten, para el sistema hiperbólico general homogéneo de leyes de conservación, se obtiene el esquema explicito TVD de segundo orden de exactitud sobre cinco puntos para el modelo hiperbólico sin fuente, dado por

Observemos que $R^k$(2) es la segundo componente de los auto vectores $R^1$ y $R^2$. Por otro lado, los parámetros de estas ecuaciones están dadas por:

Donde l= j-2,1j,j+1

Además $u_n^t$ $=\frac{m_n^t}{{\rho }_n^t}$ y $\lambda =\frac{∆t}{∆\chi }$. Por últimose presenta el esquema de Godunov, explícito de primer orden de exactitud para el segundo y el penúltimo nodo, sobre tras puntos:

Sustituyendo en las ecuaciones anteriores, las componentesdel vector de variables conservativas $\overrightarrow{U}$

Donde j=1,2,n-1, y $\lambda =\frac{∆t}{∆x}$. sustituyendo estos valores en la ecuación (49), obtenemos los valores de densidad del gas y la masa de flujo, como sigue:

Como parte de la experimentación numérica, se usarán los métodos explícitos de Goduno v de primer de orden de exactitud y el esquema TVD de Harten de segundo orden de exactitud, para simular el problema de cierre instantáneo de una válvula al final del gasoducto. Cabe destacar que esta condición severa se usa para probar la estabilidad y consistencia de dichos métodos numéricos, para calcular los valores de interés en el problema de transporte de gas en gasoductos inclinados.

A continuación se presentan los datos de entrada del gasoducto, parámetros presentes en el proceso de transporte de gas en gasoductos inclinados, así como las condiciones de borde en el gasoducto, esto es: Longitud del gasoducto: L=1000m, diámetro interno del gasoducto D=0,5m. fricción nterna del gasoducto $f_g$=3,1x10^-3, velocidad del sonido en el fluido c=360 m/s, longitud entre nodos espaciales consecutivos $∆\chi =10m$, pasos temporales $∆t=$0,01s, densidad del gas al comienzo del gasoducto $\rho \left(0\right)=50\raisebox{1ex}{$kg$}\!\left/ \!\raisebox{-1ex}{$m^3$}\right.$ y el ángulo de inclinación $\alpha$

Primero presentamos los resultados gráficos obtenidos al resolver el sistema de ecuaciones en derivadas parciales hiperbólico no lineal, no homogéneo en caso estacionario; cabe destacar que estos resultados serán usados como condiciones iniciales del problema transitorio de transporte de gas en gasoductos inclinados.

Ahora bien, para efectos de comprobar que el método numérico presentado en este trabajo es consistente y estable para este problema de transporte de gas, se dejara correr el problea transitorio sin imponer en la condición de borde al final del gasoducto $m\left(x\frac{t}{n}\right)$=0, por un tiempo de 0$\le t\le$

Como podemos observar la densidad del gas continua manteniéndose de la misma forma; es decir, la densidad del gas a lo largo del gasoducto inclinado está convergiendo al resultado esperado. Por otro lado, vemos que el vector flujo de masa presenta una pequeña baja temporal ya que, debido a los ajustes en la caída de densidad, el flujo de masa del gas vuelve a sus valores de la condición inicial. Ahora tenemos la gráfica de densidad del gas y flujo de masa a lo largo del gasoducto,para el tiempo t=100 segundos.

Ahora se puede observar como las variables conservativas de densidad y flujo de masa del gas ya han convergido a la solución, obteniendo resultados gráficos consistentes y estables. Por otro lado, consideremos la prueba mayor de estabilidad y consistencia de los métodos presentados en este trabajo, mediante la simulación de un cierre instantáneo de válvula al final del gasoducto; es decir, haciendo $m\left(x_n^t\right)=0$, para el tiempo 120<t$\le$180, segundos.

Como se puede observar en las gráficas anteriores, se ha obtenido el valor flujo de masa del gas al final del gasoducto igual a cero; es decir, se ha cerrado instantáneamente la válvula al final del gasoducto, teniendo como resultado el aumento de la densidad al final del gasoducto, por efectos del confinamiento de los volúmenes fijos de gas al final. A continuación se presenta los resultados obtenidos de densidad y flujo de masa del gas a lo largo del gasoducto para el tiempo t = 141 y t =171 segundos.

En las gráficas anteriores se observa como las variables conservativas de interés son calculadas; es decir, el método numérico se ha comportado aceptablemente bien, o en otras palabras, se sigue calculando las variables sin presentar oscilaciones en sus gráficas. Además, en la figura se observa como las ondas de propagación de masa de flujo se devuelven aguas arriba, dando así aumento de densidad a lo largo del gasoducto.

Cabe destacar que estas gráficas convergerán a rectas horizontales después de cierto tiempo, ya que cuando el flujo de masa sea cero en todo el gasoducto, entonces no habrá caída de densidad; es decir, la presión es una recta horizontal. Por otro lado, como otra prueba al método numérico presentado en este trabajo, haremos apertura instantánea de la válvula al final del gasoducto; es decir, impondremos nuevamente $m\left(x_n^t\right)=110$, a partir del tiempo t>180 segundos, obtenido los siguientes resultados gráficos para la densidad del gas y la mas de flujo a lo largo del gasoducto.

a se puede observar nuevamente, como el flujo de masa del gas vuelve a ser como la condición de frontera al inicio del problema transitorio del transporte de gas, dando como resultado la caída instantánea de la densidad al final del gasoducto; tal como se esperaba. A continuación se presentan los resultados gráficos de las variables conservativas en este trabajo para el tiempo t = 181 y t = 250 segundos.

En la figura 11 se puede observar como las gráficas de las variables conservativas van convergiendo al caso estacionario obtenido transitoriamente. Cabe destacar, que si seguimos obteniendo los resultados gráficos para las variables conservativas, se obtendrán como en el caso estacionario del problema transitorio de transporte de gas en gasoductos inclinados; y los errores para densidad del gas y flujo de masa del gas, convergerán a los siguientes valores prácticamente ceros, Error_p = 0.001 y Error_m = 0.01.

Por último en la figura 12 se obtuvieron resultados gráficos semejantes al caso estacionario tal como se esperaba sin presentar problemas de oscilaciones en dichas gráficas; en otras palabras, el método numérico presentado en este trabajo se comportó de manera estable y consistente a la hora del cálculo de las variables de interés, para este problema durante todo el tiempo de ejecución. Más aun para consolidar lo dicho anteriormente, a continuación se presentan los resultados gráficos de las variables conservativas en este trabajo para el tiempo t = 400 segundos.

En los resultados obtenidos de la experimentación numéricae pudo observar que, cuando se comienza a calcular el problema transitorio las gráficas presentaron modificaciones, respecto a los resultados gráficos del caso estable o estacionario. Esto se deba una pequeña perturbación inicial que conduce a una onda amortiguada para los primeros nodo espaciales, la cual después de cierto tiempo convergieron al estado estable proporcionalmente como las condiciones iniciales. Se asume esto al método Splitting implementado en este trabajo, el cual será sometido a pruebas y experimentos para su mejoramiento.

Por otro lado, estos resultados gráficos de las variables conservativas de interés en este trabajo, tuvieron un comportamiento aceptable durante todo el proceso transitorio del transporte de gas en gasoductos inclinados; es decir, se mantuvieron estables y consistentes durante todo el proceso de cálculo. Luego del inicio de las perturbaciones programadas, de cierre y apertura instantánea de válvula al final, los resultados obtenidos de forma numérica y gráfica no presentaron ninguna forma de inconsistencia o inestabilidad numérica, y están acordes a lo que muestra la literatura.

Cabe destacar que la longitud de un gasoducto inclinado en la vida real son de longitudes cortas en comparación con la longitud del gasoducto considerado en la experimentación numérica de este trabajo. Además, tenemos que el cierre de válvula al final, en la experimentación numérica fue de forma instantánea, haciendo el flujo de masa al final del gasoducto igual a cero, y en la vida real no es así; es decir, en un cierre de válvula, el flujo de masa al final del gasoducto se va haciendo cero, linealmente entre 10 ª 20 segundos. Considerando estas dos observaciones, de la longitud del gasoducto inclinado y el cierre progresivo de la válvula al final del gasoducto, se obtendrán resultados gráficos excelentes, sin modificaciones durante el proceso transitorio. Además, que dichos resultados serán más acertados a la realidad de un proceso de transporte de gas natural.

Además, se pueden resaltar varios puntos importantes. Es necesaria una mejor experimentación empleando el método Splitting para tener conclusiones más acertadas sobre su uso y precisión. Una validación del método sería esencial para emplear este modelo en aplicaciones de la vida real, a fin de que se verifique la veracidad de los resultados obtenidos.

Por último, el método numérico aplicado a la resolución del modelo de transporte de gas en gasoductos inclinados presentado en este trabajo demostró ser efectivo, estable, consistente y estar acorde a los resultados numéricos y gráficos que se encuentran en la literatura, por lo cual se seguirá trabajando con el modelo en futuras investigaciones y se empleará para resolver distintos problemas de transporte de gas.

El presente trabajo ha sido posible gracias al apoyo brindado por el Departamento de Cálculo, de la Escuela Básica de la Facultad de Ingeniería, de la Universidad de Los Andes. Mérida – Venezuela.

El presente trabajo ha sido posible gracias al apoyo brindado por el Departamento de Cálculo, de la Escuela Básica de la Facultad de Ingeniería, de la Universidad de Los Andes. Mérida – Venezuela.