Correlación de Datos Experimentales Aplicado al Flujo Multifásico

INTEVEP, S.A. Centro de Investigaci´ on y Apoyo Tecnol´ ogico. Filial de PDVSA. Los Teques, Edo. Miranda. VENEZUELA.

CORRELACION DE DATOS EXPERIMENTALES APLICADO AL FLUJO MULTIFASICO ANDRES L. GRANADOS M. (EPPR/322).

Departamento de Producci´ on (EPPR). Secci´ on Manejo de Crudos (EPPR/3). Unidad de transporte de Crudos (EPPR/32).

PROYECTO (STE) No. 6555. REPORTE TECNICO No. INT-EPPR/3-NT-92-004.

MARZO, 1992.

CONFERENCIA SOBRE ESTADO DEL ARTE EN MECANICA DE FLUIDOS COMPUTACIONAL. Auditorium de INTEVEP. Los Teques, 27-29 Mayo. 1991. Edo. Miranda, Venezuela.

NUEVAS CORRELACIONES PARA FLUJO MULTIFASICO

ANDRES L. GRANADOS M. INTEVEP, S.A. Departamento de Producci´ on (EPPR/32). Centro de Investigaci´ on y Apoyo Tecnol´ ogico. Filial de PDVSA. Los Teques, Estado Miranda. Venezuela.

RESUMEN.-

Se han desarrollado correlaciones y expresiones anal´ıticas de las mismas con la finalidad de sistematizar de una forma m´ as eficiente los algoritmos de c´ alculos implementados en los paquetes de simulaci´ on. Muchos de los m´ odulos de los paquetes de simulaci´ on del flujo multif´ asico emplean correlaciones basadas en lectura directa de tablas o gr´ aficos, y en la actualidad lo que hace el algoritmo de c´ alculo es que tiene en memoria, asignando a variables, los datos num´ericos de dichas lecturas. Luego realiza una b´ usqueda binaria para identificar el o los intervalos donde se aplicar´ a posteriormente una interpolaci´ on polin´ omica. Este proceso es sumamente ineficiente por la gran cantidad de operaciones l´ ogicas involucradas dentro de la programaci´ on computacional del algoritmo. Esto se acent´ ua todav´ıa m´ as cuando se trabajan con tablas o gr´ aficos que representan funciones de dos o m´ as variables. En este trabajo se est´ an presentando tres correlaciones desarrolladas a partir de datos tabulados o gr´ aficos. Estas correlaciones son: Grado API y Gravedad Espec´ifica, Coeficiente de Arrastre para Esferas, Correlaci´ on de Eaton para la Fracci´ on de l´ıquido y Viscosidad de Gases Naturales. La primera correlaci´ on esta basada en los datos tabulados para la correcci´ on del grado API medido a la temperatura ambiente para llevarlo a

60o F. La segunda, la tercera y la cuarta correlaciones se hallaron con datos obtenidos de lectura directa sobre gr´aficos que representaban las curvas de las caracter´ısticas se¯ naladas. Todos los ajustes de curvas (Correlaci´ on 2 y 3) y de superficies (Correlaci´ on 1 y 4) fueron realizados empleando un programa de regresi´ on m´ ultiple lineal y no lineal, denominado ADJUST. Este programa realiza los ajustes lineales de forma autom´ atica, es decir, no hay que definir a priori la funci´ on a ajustar. El programa lo hace autom´ aticamente. Los resultados de todos estos ajustes son bastante buenos num´ericamente hablando y los errores en todos los casos son muy inferiores a los errores experimentales que se generaron durante la obtenci´ on de los datos originales.

INTRODUCCION

Se han desarrollado correlaciones y expresiones anal´ıticas de las mismas con la finalidad de sistematizar de una forma m´as eficiente los algoritmos de c´alculos implementados en los paquetes de simulaci´ on. Muchos de los m´odulos de los paquetes de simulaci´ on del flujo multif´ asico emplean correlaciones basadas en lectura directa de tablas o gr´ aficos, y en la actualidad lo que hace el algoritmo de c´alculo es que tiene en memoria, asignando a variables, los datos num´ericos de dichas lecturas. Luego realiza una b´ usqueda binaria para identificar el o los intervalos donde se aplicar´ a posteriormente una interpolaci´ on polin´ omica. Este proceso es sumamente ineficiente por la gran cantidad de operaciones l´ ogicas involucradas dentro de la programaci´ on computacional del algoritmo. Esto se acent´ ua todav´ıa m´as cuando se trabajan con tablas o gr´aficos que representan funciones de dos o m´as variables. En este trabajo se est´ an presentando tres correlaciones desarrolladas a partir de datos tabulados o gr´aficos. Estas correlaciones son: • Grado API y Gravedad Espec´ifica. • Coeficiente de Arrastre para Esferas. • Correlaci´on de Eaton para la Fracci´ on de l´ıquido. • Viscosidad de los Gases Naturales. La primera correlaci´ on esta basada en los datos tabulados para la correcci´ on del grado API medido a la temperatura ambiente para llevarlo a 60o F. La segunda, la tercera y la cuarta correlaciones se hallaron con datos obtenidos de lectura directa sobre gr´ aficos que representaban las curvas de las caracter´ısticas se¯ naladas. Todos los ajustes de curvas (Correlaci´on 2 y 3) y de superficies (Correlaci´ on 1 y 4) fueron realizados empleando un programa de regresi´ on m´ ultiple lineal y no lineal, denominado ADJUST. Este programa realiza los ajustes lineales de forma autom´ atica, es decir, no hay que definir a priori la funci´ on a ajustar. El programa lo hace autom´ aticamente. Los resultados de todos estos ajustes son bastante buenos num´ericamente hablando y los errores en todos los casos son muy inferiores a los errores experimentales que se generaron durante la obtenci´on de los datos originales.

2

METODOS NUMERICOS

POLINOMIOS DE INTERPOLACION

Sean (x0 , f (x0 )), (x1 , f (x1 )), (x2 , f (x2 )) hasta (xn , f (xn )) n + 1 puntos discretos que representan a una funci´ on y = f (x). Como se sabe, existe un u ´ nico polinomio y = Pn (x) de grado n que pasa por los n + 1 puntos mencionados. Estos polinomios son adecuados para realizar estimaciones de la funci´ on y = f (x) para un valor x cualquiera perteneciente al intervalo [x0 , x1 , x2 , x3 , . . . , xn ] que contiene a todos los puntos, estando los valores xi no necesariamente ordenados. A este proceso se le denomina “Interpolaci´ on”. Si el valor x est´a fuera del intervalo de los puntos entonces el proceso se denomina “Extrapolaci´on”.

DIFERENCIAS DIVIDIDAS Las diferencias divididas [1] simbolizadas por f [ · ] se definen de forma recurrente de la siguiente forma f [x0 ] = f (x0 )

f [x1 , x0 ] =

f [x2 , x1 , x0 ] =

f [x3 , x2 , x1 , x0 ] =

f [xn , xn−1 , . . . , x1 , x0 ] =

(1.a)

f [x1 ] − f [x0 ] x1 − x0

(1.b)

f [x2 , x1 ] − f [x1 , x0 ] x2 − x0

(1.c)

f [x3 , x2 , x1 ] − f [x2 , x1 , x0 ] x3 − x0

f [xn , xn−1 , . . . , x2 , x1 ] − f [xn−1 , xn−2 , . . . , x1 , x0 ] xn − x0

(1.d)

(1.e)

Las diferencias divididas cumplen con la propiedad f [xn , xn−1 , . . . , x1 , x0 ] = f [x0 , x1 , . . . , xn−1 , xn ]

∀n ∈ IN

(2)

Esta propiedad, expresada de otra forma, lo que significa es que, sin importar el orden en que est´ an los valores xi dentro de una diferencia dividida, el resultado es siempre el mismo. Una forma de expresar todas las diferencias divididas posibles de generar mediante, por ejemplo, un conjunto de puntos x0 , x1 , x2 , x3 y x4 , no necesariamente ordenados, es lo que se denomina el Diagrama 3

Romboidal de diferencias divididas. Para el ejemplo propuesto se tiene que el diagrama romboidal se representa como

x0

f [x0 ]

x1

f [x1 ]

x2

f [x2 ]

x3

f [x3 ]

x4

f [x4 ]

f [x0 , x1 ] f [x1 , x2 ] f [x2 , x3 ] f [x3 , x4 ]

f [x0 , x1 , x2 ]

f [x0 , x1 , x2 , x3 ]

f [x1 , x2 , x3 ]

f [x0 , x1 , x2 , x3 , x4 ]

f [x1 , x2 , x3 , x4 ]

f [x2 , x3 , x4 ]

(3)

Manipulaci´ on algebra´ıca de la diferencias de o´rdenes crecientes conlleva, mediante inducci´ on, a una forma sim´etrica similar para la n-´esima diferencia dividida, en t´ermino de los argumentos xi y de los valores funcionales f (xi ). Esta forma sim´etrica puede ser escrita de manera compacta como

f [x0 , x1 , x2 , . . . , xk−1 , xk ] =

k
i=0

k

f (xi )

j=0 (xi j
=i

− xj )

(4)

POLINOMIOS DE NEWTON Los polinomios de Newton de grado n en diferencias divididas [1], como se dijo antes, permiten hacer estimaciones de la funcion y = f (x) en la forma f (x) = Pn (x) + Rn (x)

(5)

donde Pn (x) es el polinomio de grado n Pn (x) =f [x0 ] + (x − x0 ) f [x0 , x1 ] + (x − x0 )(x − x1 ) f [x0 , x1 , x2 ] + · · · + (x − x0 )(x − x1 )(x − x2 ) · · · (x − xn−1 ) f [x0 , x1 , x2 , . . . , xn−1 , xn ]

=

n k−1


(x − xj ) f [x0 , x1 , x2 , . . . , xk−1 , xk ]

(6)

k=0 j=0

y Rn (x) es el error cometido en la interpolaci´ on

Rn (x) =

n  j=0

(x − xj )

f (n+1) (ξ) (n + 1)! 4

ξ ∈ [x0 , x1 , . . . , xn−1 , xn ]

(7)

Naturalmente Rn (xi ) = 0 para i = 1, 2, 3, . . . , n, ya que el polinomio pasa por cada uno de los puntos (xi , f (xi )). Cuando el l´ımite superior de una productoria es menor que el l´ımite inferior, el resultado de dicha productoria es la unidad.

POLINOMIOS DE LAGRANGE Los polinomios de Lagrange son otra forma de expresar los mismos polinomios Pn (x) de la expresi´on (5), pero se tiene que [2]

Pn (x) =

n


Li (x) f (xi )

(8)

i=0

donde

Li (x) =

n  (x − xj ) (xi − xj ) j=0

(9)

j
=i

El error Rn (x) contin´ ua siendo el mismo que la expresi´on (7).

METODO DE MINIMOS CUADRADOS

Sea un conjunto de p valores x1 , x2 , x3 , . . ., xp , donde cada xi representa un a (m + 1)-upla de la forma xi = (x1 , x2 , x3 , . . . , xm , xm+1 )i

(10)

con todas sus componentes independientes entre s´ı. Sea y = f (x) una funci´ on escalar que expresa una de las componentes de la (m + 1)-upla en funci´ on de las restantes. Es decir, por ejemplo, que xm+1 = f (x1 , x2 , x3 , . . . , xm )

(11)

Esto se puede hacer sin p´erdida de generalidad, puesto que siempre se puede hacer una transformaci´on del tipo ˜ = H(x) x ˜ i tengan todos sus componentes independientes entre s´ı, al igual que los xi en (10). tal que los x 5

(12)

Dados los valores xi y definida la funci´ on y = f (x), se puede ahora tratar de encontrar una funci´ on de aproximaci´ on Y = F (x, c) dependiente, no s´ olo de los x, sino tambi´en de n par´ ametros cj , los cuales son expresados en la forma de una n-upla como c = (c1 , c2 , c3 , . . . , cn )

(13)

Estos par´ ametros c se escogen tales que la funci´on Yi = F (xi , c) se aproxime lo mejor posible a la funci´ on yi = f (xi ), para todos los valores xi , con i = 1, 2, 3, . . . , p. El m´etodos de los m´ınimos cuadrados en particular lo que trata es de encontrar los valores de los p´ arametros c de la funci´ on Y = F (x, c), tales que el valor

S=

p


(Yi − yi )2

(14)

i=1

sea el m´ınimo posible. El valor S representa la sumatoria de todas las desviaciones, entre la funci´on definida para los puntos y la funci´ on de aproximaci´ on encontrada, al cuadrado. El valor S se puede interpretar de dos formas posibles. Se puede interpretar como un funcional de la funci´ on F , es decir, S = S(F (x, c), donde a su vez la funci´ on F depende de unos par´ ametros c que forman parte de la misma. En este caso el m´etodo de m´ınimos cuadrados se convierte en un problema variacional. El valor S tambi´en se puede interpretar como una funci´ on de los par´ ametros, es decir, S = S(c), asumiendo una funci´ on de aproximaci´ on ya encontrada. Esta u ´ltima forma es la que vamos a interpretar aqu´ı. Con la aclaratoria precedente, entonce la definici´on (14) se puede expresar como

p
S(c) = [F (xi , c) − f (xi )]2

(15)

i=1

En algunos casos se alteran la desviaciones con una funci´ on de peso W (x) para hacer que el ajuste de los par´ ametros tienda hacer la aproximaci´on mejor para unos valores de xi que para otros. Esto es

S(c) =

p


W (x)[F (xi , c) − f (xi )]2

(16)

i=1

Sin embargo, todas las deducciones se har´ an para el m´etodo de los m´ınimos cuadrados expresado como est´a en (15). Extender estos resultados a como est´a expresado el m´etodo en (16) es muy sencillo. 6

El m´etodo de m´ınimos cuadrados es en s´ı un procedimiento para encontrar el valor m´ınimo de la funci´ on (15) de S = S(c). Para ello debe encontrarse los valores de cj , tales que hagan las derivadas de S(c) todas nulas. En otras palabras,

p  ∂F 
∂S =2 [F (xi , c) − f (xi )] =0 ∂cj ∂cj xi i=1

(17)

Las expresiones (17) se denominan “Ecuaciones Normales” y deben cumplirse simult´ aneamente para j = 1, 2, 3, . . . , n. Su nombre se debe a que la derivadas son calculadas para una hipersuperficie donde las direcciones cj son ortogonales entre s´ı y est´an evaluadas en un punto c donde todas son nulas. Las direcciones son ortogonales debido a que los par´ ametros cj son todos independientes entre s´ı. Es bueno hacer notar que lo que se halla mediante este procedimiento es un m´ınimo de la funci´ on escalar S(c) y no un m´ aximo, puesto que la funci´ on Yi = F (xi , c) puede estar tan alejada de los valores xi como se quiera, variando los valores de los par´ ametros cj .

AJUSTE LINEAL. • Series de Funciones Bases. El ajuste lineal es el m´as empleado y el m´ as reportado en la literatura. Su nombre se debe a que la funci´ on de aproximaci´ on posee una expresi´on lineal de la forma [2]

F (x, c) =

n


ck gk (x)

(18)

k=1

lo que no es m´as que una serie de funciones gj (x) todas diferentes entre s´ı, por lo que son consideradas que forman parte de una base de un espacio de funciones. Para el caso particular de la funci´ on de aproximaci´ on definida por (18) se tiene que

∂F = gj (x) ∂cj 7

(19)

Substituyendo este resultado en la expresi´on (17), junto con la definici´ on (18) e intercambiando las sumatorias de k con la sumatoria de i, se obtiene  p 
n
∂S =2 ck gk (xi ) − f (xi ) gj (xi ) = 0 ∂cj i=1

(20.a)

k=1

p
n


ck gk (xi )gj (xi ) =

i=1 k=1 p n 



p


f (xi )gj (xi )

(20.b)

i=1

 p
gj (xi )gk (xi ) ck = gj (xi )f (xi )

k=1 i=1

(20.c)

i=1

Al final queda un sistema de ecuaciones lineales de la forma n


Ajk ck = bj

[A]c = b

(21)

k=1

donde los elementos de la matriz del sistema y el vector independiente se expresan como

Ajk =

p


gj (xi )gk (xi )

(22.a)

gj (xi )f (xi )

(22.b)

i=1

bj =

p
i=1

• Series de Funciones Ortogonales. Como ejemplos [3] de funciones de aproximaci´ on m´ as utilizados se tienen las series de funciones polin´ omicas F (x, c) =

n


ck xk−1

(23)

ck cos[(k − 1)x]

(24)

k=1

y la serie de funciones trigonom´etricas

F (x, c) =

n
k=1

Tambi´en existen series de funciones racionales, hiperb´ olicas, polinomios de Chebyshev, polinomios de Legendre, etc. Tambi´en se pueden tener combinaciones de estas funciones. 8

AJUSTE NO LINEAL. En el ajuste no lineal de los par´ ametros cj la funci´ on de aproximaci´ on F (x, c) tiene una expresi´on distinta a la expresi´ on (18), por consiguiente, lo que se obtiene es un sistema de ecuaciones no lineales en las variable cj que puede ser resuelto con cualquier m´etodo para tales tipo de sistemas, como, por ejemplo, el m´etodo de Newton-Raphson. Sin embargo esto trae como consecuencia que el procedimiento de m´ınimos cuadrados se vuelva m´ as complicado, ya que hay que calcular la matriz jacobiana del sistema de funciones no lineales. Para evitar el inconveniente mencionado se han desarrollado varios m´etodos, dentro los cuales est´an: - M´etodo del m´ aximo descenso. - M´etodo de Gauss-Newton. - M´etodo de Marquardt. Todos estos m´etodos se derivan del siguiente an´ alisis. Sea la expansi´ on en series de Taylor hasta el t´ermino de primer orden de la funci´ on de aproximaci´on F (xi , c) alrededor de un valor estimado c∗ de los par´ ametros. Esto es,

F (xi , c) = F (xi , c∗ ) +

n 
∂F ∗ k=1

∂ck

xi

∆c∗k + O(∆c∗ 2 )

(25.a)

donde ∆c∗k = ck − c∗k

(25.b)

Substituyendo este resultado en la expresi´on (17) e intercambiando las sumatorias de k con la sumatoria de i, se obtiene un sistema de ecuaciones de la forma

p  p n 
 ∂F 

∂F   ∂F ∗ ∗ [F (xi , c∗ ) − f (xi )] ∆ck = − ∂cj xi ∂ck xi ∂cj xi i=1 i=1

(26)

k=1

on del punto valor c∗k al valor ck , entonces Si se supone que las derivadas ∂F/∂cj no sufren gran variaci´ la expresi´on (26) se podr´ıa reescribir aproximadamente como

n


A∗jk ∆c∗k = b∗j

k=1

9

(27.a)

donde p 
∂F ∗  ∂F ∗ ∂cj xi ∂ck xi i=1

(27.b)

 ∂F ∗ [F (xi , c∗ ) − f (xi )] ∂cj xi i=1

(27.c)

A∗jk =

b∗j = −

p


Con base en este an´alisis, entonces se pueden aplicar los diferentes m´etodos que se explican a continuaci´on.

• M´etodo de Gauss-Newton. El m´etodo de Gauss-Newton consiste en un procedimiento iterativo que se origina a partir de las expresiones (27) junto con la definici´on (25.b). De esta forma resulta el siguiente algoritmo iterativo con s como indicador del n´ umero de la iteraci´on n


Asjk ∆csk = bsj

(28)

k=1

donde p 
∂F s  ∂F s ∂cj xi ∂ck xi i=1

(29.a)

 ∂F s [F (xi , cs ) − f (xi )] ∂cj xi i=1

(29.b)

Asjk =

bsj = −

p


y luego se obtiene cs+1 = cs + ∆cs

(30)

La expresi´on (28) representa un sistema de ecuaciones lineales que se resuelve en cada iteraci´on s, conociendo los par´ ametros cs . Despu´es se substituye este resultado en la expresi´on (30) para obtener los valores cs+1 de los par´ ametros en la siguiente iteraci´on. El procedimiento se contin´ ua hasta obtener convergencia hacia la soluci´on c de las ecuaciones normales (17), aplicando un criterio de parada de la forma ∆cs  < εmax en el error local de las variables cs y donde εmax es la tolerancia permitida para dicho error local. 10

(31)

Frecuentemente es recomendable alterar el algoritmo relaj´andolo en la forma cs+1 = cs + ω∆cs

(30 )

para asegurar la convergencia del proceso iterativo. Aqu´ı ω es el factor de relajaci´ on y en cada iteraci´on se altera ω  = ρω

ρ1

(32 )

Normalmente se emplean los valores de ρ = 0.5 y τ = 2, para producir el efecto de una b´ usqueda del ω ´ optimo mediante la bisecci´on consecutiva de los intervalos [csk , csk + ∆csk ], comenzando con un ω = 1. Cuando las derivadas de las expresiones (29) se hacen complicadas de calcular, estas pueden ser obtenidas num´ericamente de la siguiente forma  ∂F s ∂cj

xi

∼ =

F (xi , cs ) − F (xi , cs−1 (j) ) csj − cs−1 j

(34.a)

donde s−1 s s s F (xi , cs−1 , . . . , csn ) (j) ) = F (xi , c1 , c2 , c3 , . . . , cj

(34.b)

• M´etodo del M´ aximo Descenso. El m´etodo del m´ aximo descenso est´a basado en el hecho de que ∂S

s

= −bsj ∂cj

(35)

on Es decir, que S(cs ) se incrementa en la direcci´on indicada por el gradiente (35). Si se escoge una direcci´ ∆cs opuesta a este gradiente tal que ∆csj = ωbsj 11

(36)

se obtendr´ a el m´aximo descenso de la funci´ on S(c). La expresi´on (36) se puede reescribir como

n


s Djk ∆csk = bsj

(37)

s Djk = δjk

(38)

cs+1 = cs + ω∆cs+1

(39)

k=1

donde

y se tiene que

s Sin embargo, el m´etodo puede ser modificado de manera tal que la matriz Djk tenga dimensiones

acorde con la funci´ on S(c), y, por consiguiente, se puede hacer s Djk = As δjk

(40)

El valor de ω se modifica de igual forma que el m´etodo de Gauss-Newton para asegurar la convergencia, pero por el contrario el m´etodo del m´ aximo descenso converge muy lentamente y por lo tanto no es recomendable su uso.

• M´etodo de Marquardt. La formula algor´ıtmica del m´etodo de Marquardt [4] es la siguiente

n


s (Asjk + λDjk )∆csk = bsj

(41.a)

k=1

cs+1 = cs + ∆cs

(41.b)

donde el factor λ funciona similar a un factor de relajaci´on y le da al m´etodo de Marquardt un car´ acter hibrido donde existe un compromiso entre el m´etodo del m´ aximo descenso y el m´etodo de Gauss-Newton. Cuando λ → 0, la direcci´ on del m´etodo se dirije hacia el m´etodo de Gauss-Newton. Cuando λ → ∞, la direcci´on del m´etodo se dirije hacia el m´etodo del m´ aximo descenso. Los estudios de Marquardt indican que el m´etodo posee un ´angulo promedio entre los m´etodos de Gauss-Newton y M´ aximo Descenso de 90o . La selecci´on de un λ entre 0 e ∞ produce una direcci´on intermedia. 12

Para efectos de garantizar la convergencia en cada iteraci´ on se altera el factor λ de la forma λ = λ/ρ

ρ1

(42 )

N´otese que incrementar λ en el m´etodo de Marquardt es equivalente a disminuir ω en el m´etodo de Gauss-Newton. Normalmente, se toman los valores de λinicial = 10−3 , ρ = 0.1 y τ = 10. Cuando en varias iteraciones consecutivas el m´etodo mejora su convergencia, es decir se cumple la relaci´on (33), entonces λ → 0, y esencialmente se estar´a empleando el m´etodo de Gauss-Newton. Si la convergencia no mejora, por el contrario, λ se incrementa y se estar´a usando pr´ acticamente el m´etodo del M´ aximo Descenso.

EVALUACION DEL AJUSTE. La evaluaci´on del ajuste viene dada mediante el an´ alisis de de ciertos factores que permiten, por un lado comparar cuan bueno es un ajuste en relaci´ on a otro, y por otro lado comparar cuando un ajuste reproduce bien el conjunto de puntos de datos. Estos coeficientes son los siguientes: Suma de los cuadrados de las desviaciones con respecto a la funci´on de aproximaci´ on.

S(c) =

p
[F (xi , c) − f (xi )]2

(44)

i=1

Media de la variable dependiente.

1
f (xi ) p i=1 p

fm =

13

(45)

Suma de los cuadrados de las desviaciones con respecto a la media de la variable dependiente.

Sm =

p


[f (xi ) − fm ]2

(46)

i=1

Desviaci´on est´ andar (σ) o´ varianza (σ 2 ) con respecto a la funci´ on de aproximaci´ on.  S p−n

σ=

(47)

2 Desviaci´on est´ andar (σm ) o´ varianza (σm ) con respecto a la media fm .

 Sm p−1

σm =

(48)

Coeficiente de determinaci´o (r2 ) o´ coeficiente de correlaci´on (r). r2 indica el porcentaje de la incertidumbre inicial que ha sido disminuido usando la funci´ on de aproximaci´ on.

r2 =

Sm − S Sm

(49)

En algunas literaturas definen el coeficiente de determinaci´ on (R2 ) o´ coeficiente de correlaci´on (R) de la siguiente forma. R2 =

σm − σ σm

(50)

σm fm

(51)

Coeficiente de variaci´on. Cv =

Desviaci´on RMS (Root of the Mean Square).  δrms = 14

S p

(52)

Desviaci´on m´ axima. δmax = max |F (xi , c) − f (xi )| 1≤i≤p

(53)

En la desviaci´on est´ andar σ, la cantidad S est´a dividida por (p − n), debido a que n par´ ametros (c1 , c2 , c3 , . . . , cn ), derivados de los datos originales (x1 , x2 , x3 , . . . , xp ), fueron usados para computar S(c). De aqu´ı que se hallan perdido n grados de libertad en la probabilidad. En la desviaci´on est´ andar σm , la cantidad Sm est´a dividida por (p − 1), debido a que la media de o de los datos originales (x1 , x2 , x3 , . . . , xp ), fu´e usada para la variable dependiente, fm , la cual se deriv´ computar Sm . De aqu´ı que se halla perdido un grado de libertad [5]. La desviaci´on est´ andar σm debe ser mayor que σ, de otra forma no se justifica el uso de la funci´ on de aproximaci´on, y la media fm da una mejor aproximaci´ on a los datos que la funci´ on de aproximaci´ on propuesta. La funci´ on de aproximaci´ on que mejor se ajusta a los datos originales (x1 , x2 , x3 , . . . , xp ), no es aquella que ofrece un menor valor de S, sino aquella que brinda una menor desviaci´on est´ andar σ, con respecto a la funci´on de aproximaci´ on. Esto implica que el coeficiente de determinaci´on R2 sea m´as cercano a la unidad. El coeficiente de variaci´on Cv nos brinda una medida normalizada de cual es la dispersi´ on de los datos originales y normalmente se da en forma porcentual [5]. Cuando la dispersi´ on de los datos es muy grande significa que los puntos est´an muy dispersos y si se grafican formar´ an una nube ancha alrededor de cualquier correlaci´on que se trate de hallar. En este caso, la mejor correlaci´on la ofrecer´ıa la media fm . La desviaci´on RMS y la desviaci´on m´ axima dependen del ajuste particular que se est´a realizando. La desviaci´on RMS se puede interpretar como una desviaci´on promedio del ajuste. La desviaci´ on m´ axima acota cuanto va a hacer el mayor error cometido con el ajuste. Entre mayor sea la diferencia entre estas dos desviaciones, mejor ser´a el ajuste por si mismo.

15

GRADO API Y GRAVEDAD ESPECIFICA

Se ha desarrollado una correlaci´ on para la correcci´ on del grado API obtenido de una muestra de on de laboratorio a una temperatura diferente a 60o F y mediante un procedimiento indirecto de medici´ la gravedad espec´ıfica. La base de datos consisti´o en una tabla de correcci´on del grado API obtenida de un ap´endice del manual TREATING OIL EMULSIONS [6], el cual es la reproducci´ on de una porci´ on de “Table 5, Reduction of Observed API Gravity to API Gravity at 60o F ”, PETROLEUM MEASUREMENT TABLES, American Edition (Philadelphia, Pa.: American Society for Testing and Materials, 1952). La correlaci´on se obtuvo empleando la t´ecnica num´erica de ajuste de par´ ametros con el m´etodo de m´ınimos cuadrados aplicado a una serie de funciones bases, o regresi´on lineal como tambi´en se le conoce. Para esto se utiliz´o el programa de computadora denominado ADJUST, el cual realiza ajustes de superficies (es decir, ajustes de funciones de varias variables) de forma autom´atica. El programa prueba una infinidad de funciones bases combinadas o compuestas de diferentes formas y de todos los ajustes efectuados escoge el mejor o los mejores. El usuario del programa s´olamente define las funciones elementales que intervienen (no es necesario definir sobre que, se aplica la funci´on elemental) y los rangos de n´ umero de t´erminos y factores a probar. Inicialmente se intent´ o hacer el ajuste m´as simple de la forma S ∗ = S ∗ (T, S) donde: S ∗ = Gravedad espec´ıfica corregida a 60o F . S = Gravedad espec´ıfica a la temperatura T . T = Temperatura (o F ). Con

S = 141.5/(131.5 +o AP I)

y

S ∗ = 141.5/(131.5 +o AP I ∗ )

Como no se lograron resultados satisfactorios, luego se intent´o ajustar la desviaci´ on S − S ∗ con respecto a S y la desviaci´ on de la temperatura a los 60o F , normalizada en 100oF . Es decir, se hizo

T − 60 S − S ∗ = f S, 100 16

y el ajuste obtenido fue el siguiente √ 1 + S (S ∗ − S) = A + B



T − 60 100



 −C

T − 60 100

2

donde el par´ ametro A di´ o un valor despreciable, por lo que luego se elimin´ o, y el coeficiente de determinaci´on di´ o aproximadamente 98.855347 %. La eliminaci´on de A tambi´en se justifica debido a que hace la expresi´on encontrada consistente con el comportamiento esperado, es decir, para T = 60o F , S ∗ = S. Finalmente, una vez eliminado el t´ermino independiente A, se procedi´ o a hacer el ajuste definitivo de la forma     2  √ T − 60 T − 60 S =S+ B / 1+S −C 100 100 ∗

En este caso el coeficiente de determinaci´on di´ o aproximadamente 98.855748 %, levemente superior al ajuste anterior. Los par´ ametros dieron B = 4.936603951 × 10−2 C = 1.480129729 × 10−3 La expresi´on anal´ıtica desarrollada aqu´ıpara la gravedad espec´ıfica sirve tambi´en para los grados API si se substituyen en dicha expresi´ on las definiciones de S y S ∗ en funci´ on de AP I y AP I ∗ , respectivamente. La validez de la expresi´on hallada se prob´ o haciendo conversiones de grado API a gravedad espec´ıfica y viceversa, y comprobando con la tabla de correcci´on. De esta forma no se observaron errores mayores que 0.05oAP I en los rangos de temperaturas y grados API 10o ≤ AP I ≤ 42o 45o F ≤ T ≤ 180oF Por u ´ ltimo, es bueno mencionar que se prob´ o el ajuste con el t´ermino c´ ubico en T , pero los resultados fueron levemente inferiores al ajuste finalmente reportado, lo que se puede ver comparando el coeficiente de determinaci´on que en este caso di´o 98.854875 %.

17

COEFICIENTE DE ARRASTRE PARA ESFERAS

Se desarroll´ o una correlaci´on para el coeficiente de arrastre en el flujo sobre esferas que cubre pr´ acticamente todos los reg´ımenes: Laminar, turbulento y desprendimiento de la capa l´ımite. En el estudio del flujo multif´ asico es importante esta correlaci´on para el modelaje de los mecanismos de separaci´on de dos fases inmiscibles, por ejemplo, en separadores bif´ asicos y trif´ asicos, o para el modelaje de los mecanismos de partici´on y coalescencia en la mezcla de dos fases inmiscibles en r´egimen turbulento, por ejemplo, en mezcladores est´aticos o din´ amicos. Los datos empleados en el ajuste de la correlaci´ on para el coeficiente de arrastre, CD , en funci´ on del n´ umero de Reynolds, Re, se obtuvieron de la gr´ afica correspondiente publicada en el libro de texto de Vennard y Street [7], mediante una lectura directa sobre la misma.El coeficiente de arrastre CD se define como CD =

F/A 1 2 2 ρU

donde F =Fuerza de arrastre actuando sobre la esfera. A=Area transversal de la esfera perpendicular al flujo A = πD. D=Di´ ametro de la esfera. U =Velocidad del flujo libre no perturbado. y el n´ umero de Reynolds se define como Re =

UD ν

El problema de ajustar una sola expresi´on anal´ıtica al coeficiente de arrastre para esferas se atac´o en dos etapas. En la primera etapa se ajustaron los datos para el n´ umero de Reynolds de hasta 2 × 105 en una expresi´on anal´ıtica denominada CD1 . Esta expresi´ on cubre los reg´ımenes laminar y turbulento. Mediante el ajuste automatizado con el programa ADJUST se realiz´ o el siguiente ajuste ln(CD1 Re/24) = f (X) = f [ln(1 + Re)] donde X = ln(1 + Re). As´ı se obtuvo una expresi´on de la forma

CD1 =

A1 exp(B1 X 3/2 + C1 X 3 ) Re

donde 18

A1 =24 B1 =0.1491674466 C1 =0.001110998948 Aqu´ı puede observarse que la expresi´ on es consistente en la medida que Re → 0, ya que CD1 → 24/Re. En la segunda etapa se ajustaron los datos para el n´ umero de Reynolds desde 2 × 105 hasta 2 × 106 en una expresi´on anal´ıtica denominada CD2 . Esta segunda expresi´on cubre la zona donde existe el desprendimiento de la capa l´ımite. Realizando un ajuste directo con el programa ADJUST para una par´ abola se obtuvo la expresi´ on CD2 = A2 − B2 X + C2 X 2 donde A2 =13.70070120 B2 =2.108301309 C2 =0.08162948441 Para enlazar estos dos resultados parciales se emple´o la funci´ on escal´on continua de la forma

CD = CD1 + (CD2 − CD1 )

A 1+A

donde A = exp(B(X − C)/D) B=5.3 C=12.35 D=0.47 La figura 1 muestra los resultados para CD1 y CD2 en l´ıneas punteadas, y muestra tambi´en el resultado final de la correlaci´on para CD en l´ınea continua. Como puede observarse el ajuste es bastante bueno.

19

Figura 1. Correlaci´ on para el coeficiente de arrastre sobre una esfera.

20

CORRELACION DE EATON PARA LA FRACCION DE LIQUIDO

La correlaci´on de Eaton [8] es muy empleada en la predicci´ on de la fracci´ on de l´ıquido para el flujo bif´ asico. Sin embargo, es bueno mencionar que normalmente se emplea esta correlaci´on haciendo la lectura directa de puntos sobre una gr´afica y luego realizando una interpolaci´ on polin´ omica, lo cual ocupa mucho espacio (memoria) y tiempo (c´ alculo) por parte del usuario de la correlaci´ on o por parte de la computadora. Por la raz´ on antes mencionada se ha desarrollado una expresi´on anal´ıtica que permita r´ apidamente realizar los c´alculos de la fracci´on de l´ıquido. Debido a la forma muy particular de la curva en la gr´ afica, se decidi´ o realizar con el programa ADJUST un ajuste no lineal empleando la funci´ on y = tanh(x). El estudio de los resultados de varios de estos ajustes di´o como la mejor alternativa la expresi´ on

Y =

   exp(A ln X + B) 1 A ln X + B = 1 + tanh 1 + exp(A ln X + B) 2 2

donde A=0.9377978363 B=1.034705148 Aqu´ı Y representa a la fracci´on de l´ıquido y X representa al n´ umero adimensional de Eaton. El coeficiente de determinaci´on fue en este caso del 99.942265 % y la desviaci´on m´ axima fue de 0.0175 en valor absoluto. El figura 2 muestra la correlaci´on de Eaton representando con puntos los valores tomados por lectura directa de la gr´ afica original y con l´ınea continua la expresi´on anal´ıtica desarrollada aqu´ı.

21

Figura 2. Correlaci´ on de Eaton para la fracci´ on de l´ıquido.

22

VISCOSIDAD DE LOS GASES NATURALES

En esta parte se ha desarrollado una correlaci´ on para determinar la viscosidad µg (cP ) a partir on de la presi´on P (psia) y la temperatura T (o F ). Aqu´ı se ha empleado polinomios de interpolaci´ dependientes de una variable y de dos variables hasta de s´eptimo y tercer ´ordenes, respectivamente. Las propiedades cr´ıticas se calculan de polinomios de primer grado ajustados a las curvas de la figura 3, las cuales fueron originadas por Brown y otros [9]. As´ı Pc = 708.75 − 57.5ˆ γg Tc = 169 + 314ˆ γg

(psia) (o R)

donde γg =

ρg

(aire → γg = 1)

ρa(C.N.)

Las propiedades reducidas se calculan de acuerdo a su definici´ on como

Pr =

P Pc

Tr =

T Tc

El peso molecular del gas se obtiene aplicando la ley de Avogrado e igualando la relaciones de las densidades y de los pesos moleculares entre el gas y el aire, por consiguiente, se tiene que Mg = 29γg

(Ma = 28.97Kg/Kmol)

La viscosidad del gas a presi´on atmosf´erica µog se calcula empleando polinomios de interpolaci´ on de s´eptimo orden en funci´ on del peso molecular y polinomios de primer orden en funci´ on de la temperatura obtenidos de la figura 4 superior, presentada originalmente en el trabajo de Carr y otros [10]. As´ı resulta que µog =

T˜ − 40 (B − D) + D 360

donde B=

7
i=0

23

bi Mgi

b0 = 0.56029480984621560080587295514067244∗10−7 b1 = 0.27312398997391739841697617712590227∗10−2 b2 = −0.18132032693001915414663570202538451∗10−3 b3 = 0.59390863711497936749989944058703826∗10−5 b4 = −0.10923922746785512786170687790770486∗10−6 b5 = 0.11435283316789526225894393568977734∗10−8 b6 = −0.63606853259926439127595300600699914∗10−11 b7 = 0.1459517617526926305515487822092447∗10−13

y

D=

7


di Mgi

i=0

d0 = −0.23583553216856789920903745788993099∗10−6 d1 = 0.19066214788694697028554710122475509∗10−2 d2 = −0.13697529931119144030375366373765708∗10−3 d3 = 0.48148244554067375318995242580238245∗10−5 d4 = −0.94948365870792018998160113031591755∗10−7 d5 = 0.10660996339278927171939987508140218∗10−8 d6 = −0.63600504977699266594757441336664721∗10−11 d7 = 0.15633092783799045959673459937429638∗10−13

La viscosidad del gas a presi´on atmosf´erica µog puede sufrir una correcci´on como lo indica la figura ogeno 4 superior debido a impuresas como nitr´ogeno (N2 ), di´oxido de carbono (CO2 ) y sulfuro de hidr´ (H2 S). El contenido de estas impurezas en el gas se expresa como una fracci´on molar y(·) . Asi se obtiene on atmosf´erica una viscosidad corregida µ ˜og a presi´

µ ˜og = µog + yN2 (α1 log γˆg + β1 ) + yCO2 (α2 log γˆg + β2 ) + yH2 S (α3 log γˆg + β3 ) 24

con α1 = 8.48∗ 10−3

β1 = 9.59∗ 10−3

α2 = 9.08∗ 10−3

β2 = 6.29∗ 10−3

α3 = 8.49∗ 10−3

β3 = 3.73∗ 10−3

y donde γˆg es la gavedad espec´ıfica del gas expresada a condiciones est´andar. La viscosidad asi obtenida a su vez se vuelve a corregir en funci´on de la presi´on y la temperatura empleando la figura 4 inferior. De aqu´ı se deriva la siguiente relaci´ on que contiene polinomios de interpolaci´ on en dos direcciones: la presi´ on reducida y la temperatura reducida. Esto es,

µg exp A = µ ˜og Tr

donde A=

3


aij Pri Tri

1,j=0

a00 = −2.46211820 a10 = 2.97054714 a20 = −2.86264054∗10−1 a30 = 8.05420522∗10−3 a01 = 2.80860949 a11 = −3.49803305 a21 = 3.60373020∗10−1 a31 = −1.04432413∗10−2 a02 = −7.93385684∗10−1 a12 = 1.39643306 a22 = −1.49144925∗10−1 a32 = 4.41015512∗10−3 a03 = 8.39387178∗10−2 a13 = −1.86408848∗10−1 a23 = 2.03367881∗10−2 a33 = −6.09579263∗10−4 25

Finalmente, la viscosidad del gas se obtiene como

˜og µg = µ

µg exp A =µ ˜og µ ˜og Tr

La correlaci´on para µ ˜og es valida s´ olamente para el rango de temperatura 40o F < T < 400o F . La correlaci´on para µg /˜ µog tiene la restricci´on de que fu´e obtenida apartir de datos de presi´ on reducida en el rango 1 < Pr < 20.

26

Figura 3. Propiedades pseudocr´ıticas de los gases naturales.

27

Figura 4. Correlaci´ on de las viscosidad de los gases naturales.

28

REFERENCIAS

[1] Carnahan, B.; Luther, H. A.; Wilkes, J. O. Applied Numerical Methods. John Wiley & Sons, 1969. [2] Gerald, C. F. Applied Numerical Analysis. 2nd Edition. Addison-Wesley, 1970. [3] Burden R. L.; Faires, J. D. Numerical Analysis. 3rd Edition. PWS. Boston, 1985. [4] Marquardt, D. “An Algorithm for Least Squares Estimation of Non-Linear Parameters”. SIAM J. Appl. Math.. Vol. 11, 1963, pp. 431-441. [5] Chapra, S. C.; Canale, R. P. Numerical Methods for Engineers with Personal Computer Applications. McGraw-Hill Book Company, 1985. [6] Petroleum Extension Service and American Petroleum Institute. “Treating Oil Field Emulsions”. 3rd Edition. Dallas, Texas, 1974. [7] Vennard, J.K.;Street, R.L. “Elementary Fluid Mechanics”. 6th Edition, John Wiley and Sons. New York, 1982. [8] Eaton, B. A. “The Prediction of Flow Patterns, Liquid Holdup and Pressure Losses Occurring During Continuous Two-phase Flow in Horizontal Pipelines”. Trans. AIME, pp.815. 1967. [9] Brown, G. G.; et al. “Natural Gasoline and Volatile Hydrocarbons”. N.G.A.A. 1948. [10] Carr, N. L.; et al. “Viscosity of Hydrocarbon Gases under Pressure”. Trans. AIME, pp.264. 1954.

29