INTEVEP, S.A.

Centro de Investigación y Apoyo Tecnológico.

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ón (EPPR).

Sección 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ón (EPPR/32).

Centro de Investigación y Apoyo Tecnológico. 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ás eficiente los algoritmos de cálculos implementados en los paquetes de simulación. Muchos de los módulos de

los paquetes de simulación del flujo multifásico emplean correlaciones basadas en lectura directa de tablas o gráficos, y en

la actualidad lo que hace el algoritmo de cálculo es que tiene en memoria, asignando a variables, los datos numéricos de

dichas lecturas. Luego realiza una búsqueda binaria para identificar el o los intervalos donde se aplicará posteriormente una

interpolación polinómica. Este proceso es sumamente ineficiente por la gran cantidad de operaciones lógicas involucradas

dentro de la programación computacional del algoritmo. Esto se acentúa todavı́a más cuando se trabajan con tablas o gráficos

que representan funciones de dos o más variables. En este trabajo se están presentando tres correlaciones desarrolladas a

partir de datos tabulados o gráficos. Estas correlaciones son: Grado API y Gravedad Específica, Coeficiente de Arrastre

para Esferas, Correlación de Eaton para la Fracción de lı́quido y Viscosidad de Gases Naturales. La primera correlación

esta basada en los datos tabulados para la corrección 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áficos que
representaban las curvas de las caracterı́sticas sen̄aladas. Todos los ajustes de curvas (Correlación 2 y 3) y de superficies

(Correlación 1 y 4) fueron realizados empleando un programa de regresión múltiple lineal y no lineal, denominado ADJUST.

Este programa realiza los ajustes lineales de forma automática, es decir, no hay que definir a priori la función a ajustar. El

programa lo hace automáticamente. Los resultados de todos estos ajustes son bastante buenos numéricamente hablando y

los errores en todos los casos son muy inferiores a los errores experimentales que se generaron durante la obtención de los

datos originales.
 INTRODUCCION

Se han desarrollado correlaciones y expresiones analı́ticas de las mismas con la finalidad de sistem-
atizar de una forma más eficiente los algoritmos de cálculos implementados en los paquetes de simulación.
Muchos de los módulos de los paquetes de simulación del flujo multifásico emplean correlaciones basadas
en lectura directa de tablas o gráficos, y en la actualidad lo que hace el algoritmo de cálculo es que tiene
en memoria, asignando a variables, los datos numéricos de dichas lecturas. Luego realiza una búsqueda
binaria para identificar el o los intervalos donde se aplicará posteriormente una interpolación polinómica.
Este proceso es sumamente ineficiente por la gran cantidad de operaciones lógicas involucradas dentro
de la programación computacional del algoritmo. Esto se acentúa todavı́a más cuando se trabajan con
tablas o gráficos que representan funciones de dos o más variables. En este trabajo se están presentando
tres correlaciones desarrolladas a partir de datos tabulados o gráficos. Estas correlaciones son:

• Grado API y Gravedad Específica.

• Coeficiente de Arrastre para Esferas.

• Correlación de Eaton para la Fracción de lı́quido.

• Viscosidad de los Gases Naturales.

La primera correlación esta basada en los datos tabulados para la corrección 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áficos que representaban las curvas de las caracterı́sticas
sen̄aladas. Todos los ajustes de curvas (Correlación 2 y 3) y de superficies (Correlación 1 y 4) fueron
realizados empleando un programa de regresión múltiple lineal y no lineal, denominado ADJUST. Este
programa realiza los ajustes lineales de forma automática, es decir, no hay que definir a priori la función a
ajustar. El programa lo hace automáticamente. Los resultados de todos estos ajustes son bastante buenos
numéricamente hablando y los errores en todos los casos son muy inferiores a los errores experimentales
que se generaron durante la obtención 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ón y = f (x). Como se sabe, existe un ú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ón
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ón”. Si el valor x está fuera del intervalo de los puntos entonces el proceso se denomina
“Extrapolación”.

DIFERENCIAS DIVIDIDAS

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

f [x1 ] − f [x0 ]
f [x1 , x0 ] = (1.b)
x1 − x0

f [x2 , x1 ] − f [x1 , x0 ]
f [x2 , x1 , x0 ] = (1.c)
x2 − x0

f [x3 , x2 , x1 ] − f [x2 , x1 , x0 ]
f [x3 , x2 , x1 , x0 ] = (1.d)
x3 − x0

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

f [xn , xn−1 , . . . , x1 , x0 ] = (1.e)
xn − x0

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án 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 ]
f [x0 , x1 ]
x1 f [x1 ] f [x0 , x1 , x2 ]
f [x1 , x2 ] f [x0 , x1 , x2 , x3 ]
x2 f [x2 ] f [x1 , x2 , x3 ] f [x0 , x1 , x2 , x3 , x4 ] (3)
f [x2 , x3 ] f [x1 , x2 , x3 , x4 ]
x3 f [x3 ] f [x2 , x3 , x4 ]
f [x3 , x4 ]
x4 f [x4 ]

Manipulación algebraı́ca de la diferencias de órdenes crecientes conlleva, mediante inducción, a

una forma simétrica similar para la n-ésima diferencia dividida, en término de los argumentos xi y de los
valores funcionales f (xi ). Esta forma simétrica puede ser escrita de manera compacta como

k
f (xi )
f [x0 , x1 , x2 , . . . , xk−1 , xk ] = k (4)
i=0 j=0 (xi − xj )
j
=i

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ón


n
f (n+1) (ξ)
Rn (x) = (x − xj ) ξ ∈ [x0 , x1 , . . . , xn−1 , xn ] (7)
j=0
(n + 1)!

4
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ón (5), pero se tiene que [2]

n
Pn (x) = Li (x) f (xi ) (8)
i=0

donde

n
(x − xj )
Li (x) = (9)
j=0
(xi − xj )
j
=i

El error Rn (x) continúa siendo el mismo que la expresión (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ón escalar que expresa una de las componentes de la (m + 1)-upla en
función de las restantes. Es decir, por ejemplo, que

xm+1 = f (x1 , x2 , x3 , . . . , xm ) (11)

Esto se puede hacer sin pérdida de generalidad, puesto que siempre se puede hacer una transformación
del tipo
x̃ = H(x) (12)

tal que los x̃i tengan todos sus componentes independientes entre sı́, al igual que los xi en (10).

5
 Dados los valores xi y definida la función y = f (x), se puede ahora tratar de encontrar una
función de aproximación Y = F (x, c) dependiente, no sólo de los x, sino también de n parámetros cj , los
cuales son expresados en la forma de una n-upla como

c = (c1 , c2 , c3 , . . . , cn ) (13)

Estos parámetros c se escogen tales que la función Yi = F (xi , c) se aproxime lo mejor posible a la función
yi = f (xi ), para todos los valores xi , con i = 1, 2, 3, . . . , p.

El métodos de los mı́nimos cuadrados en particular lo que trata es de encontrar los valores de
los párametros c de la función Y = F (x, c), tales que el valor

p
S= (Yi − yi )2 (14)
i=1

sea el mı́nimo posible. El valor S representa la sumatoria de todas las desviaciones, entre la función
definida para los puntos y la función de aproximación encontrada, al cuadrado.

El valor S se puede interpretar de dos formas posibles. Se puede interpretar como un funcional
de la función F , es decir, S = S(F (x, c), donde a su vez la función F depende de unos parámetros
c que forman parte de la misma. En este caso el método de mı́nimos cuadrados se convierte en un
problema variacional. El valor S también se puede interpretar como una función de los parámetros, es
decir, S = S(c), asumiendo una función de aproximación ya encontrada. Esta última forma es la que
vamos a interpretar aquı́.

Con la aclaratoria precedente, entonce la definición (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ón de peso W (x) para hacer que el ajuste de los
parámetros tienda hacer la aproximación mejor para unos valores de xi que para otros. Esto es

p
S(c) = W (x)[F (xi , c) − f (xi )]2 (16)
i=1

Sin embargo, todas las deducciones se harán para el método de los mı́nimos cuadrados expresado como
está en (15). Extender estos resultados a como está expresado el método en (16) es muy sencillo.

6
 El método de mı́nimos cuadrados es en sı́ un procedimiento para encontrar el valor mı́nimo de
la función (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,

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

Las expresiones (17) se denominan “Ecuaciones Normales” y deben cumplirse simultáneamente

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án evaluadas en un punto c donde todas son nulas.
Las direcciones son ortogonales debido a que los parámetros 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ón
escalar S(c) y no un máximo, puesto que la función Yi = F (xi , c) puede estar tan alejada de los valores
xi como se quiera, variando los valores de los parámetros cj .

AJUSTE LINEAL.

• Series de Funciones Bases.

El ajuste lineal es el más empleado y el más reportado en la literatura. Su nombre se debe a que
la función de aproximación posee una expresión lineal de la forma [2]

n
F (x, c) = ck gk (x) (18)
k=1

lo que no es más 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ón de aproximación definida por (18) se tiene que

∂F
= gj (x) (19)
∂cj

7
Substituyendo este resultado en la expresión (17), junto con la definición (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 (20.a)
∂cj i=1 k=1

p

n

p
ck gk (xi )gj (xi ) = f (xi )gj (xi ) (20.b)
i=1 k=1 i=1

n 


p 

p
gj (xi )gk (xi ) ck = gj (xi )f (xi ) (20.c)
k=1 i=1 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

p
Ajk = gj (xi )gk (xi ) (22.a)
i=1

p
bj = gj (xi )f (xi ) (22.b)
i=1

• Series de Funciones Ortogonales.

Como ejemplos [3] de funciones de aproximación más utilizados se tienen las series de funciones
polinómicas

n
F (x, c) = ck xk−1 (23)
k=1

y la serie de funciones trigonométricas

n
F (x, c) = ck cos[(k − 1)x] (24)
k=1

También existen series de funciones racionales, hiperbólicas, polinomios de Chebyshev, polinomios de

Legendre, etc. También se pueden tener combinaciones de estas funciones.

8
AJUSTE NO LINEAL.

En el ajuste no lineal de los parámetros cj la función de aproximación F (x, c) tiene una expresión
distinta a la expresión (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étodo para tales tipo de sistemas, como, por
ejemplo, el método de Newton-Raphson. Sin embargo esto trae como consecuencia que el procedimiento
de mı́nimos cuadrados se vuelva más 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étodos, dentro los cuales
están:

- Método del máximo descenso.

- Método de Gauss-Newton.

- Método de Marquardt.

Todos estos métodos se derivan del siguiente análisis.

Sea la expansión en series de Taylor hasta el término de primer orden de la función de aproxi-
mación F (xi , c) alrededor de un valor estimado c∗ de los parámetros. Esto es,

n 

∂F ∗
F (xi , c) = F (xi , c∗ ) + ∆c∗k + O(∆c∗ 2 ) (25.a)
∂ck xi
k=1

donde
∆c∗k = ck − c∗k (25.b)

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

n 
p 

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

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

n
A∗jk ∆c∗k = b∗j (27.a)
k=1

9
donde
p 

∂F ∗  ∂F ∗
A∗jk = (27.b)
i=1
∂cj xi ∂ck xi

p  ∂F ∗
b∗j = − [F (xi , c∗ ) − f (xi )] (27.c)
i=1
∂cj xi

Con base en este análisis, entonces se pueden aplicar los diferentes métodos que se explican a
continuación.

• Método de Gauss-Newton.

El método de Gauss-Newton consiste en un procedimiento iterativo que se origina a partir de las

expresiones (27) junto con la definición (25.b). De esta forma resulta el siguiente algoritmo iterativo con
s como indicador del número de la iteración

n
Asjk ∆csk = bsj (28)
k=1

donde
p 

∂F s  ∂F s
Asjk = (29.a)
i=1
∂cj xi ∂ck xi

p  ∂F s
bsj = − [F (xi , cs ) − f (xi )] (29.b)
i=1
∂cj xi

y luego se obtiene
cs+1 = cs + ∆cs (30)

La expresión (28) representa un sistema de ecuaciones lineales que se resuelve en cada iteración
s, conociendo los parámetros cs . Después se substituye este resultado en la expresión (30) para obtener
los valores cs+1 de los parámetros en la siguiente iteración. El procedimiento se continúa hasta obtener
convergencia hacia la solución c de las ecuaciones normales (17), aplicando un criterio de parada de la
forma
∆cs  < εmax (31)

en el error local de las variables cs y donde εmax es la tolerancia permitida para dicho error local.

10
 Frecuentemente es recomendable alterar el algoritmo relajándolo en la forma

cs+1 = cs + ω∆cs (30 )

para asegurar la convergencia del proceso iterativo. Aquı́ ω es el factor de relajación y en cada iteración
se altera
ω  = ρω ρ<1 (32)

de manera de garantizar que dentro de esa iteración se cumpla que

S(cs+1 ) < S(cs ) (33)

Recuérdese que lo que se está buscando es el valor de c para que la función S(c) se haga mı́nima.

Si la relación (33) se cumple en una iteración, entonces en la siguiente iteración se permite un

incremento de ω en la forma
ω = τ ω τ >1 (32 )

Normalmente se emplean los valores de ρ = 0.5 y τ = 2, para producir el efecto de una búsqueda
del ω óptimo mediante la bisección 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éricamente de la siguiente forma

 ∂F s F (xi , cs ) − F (xi , cs−1

∼ (j) )
= (34.a)
∂cj xi csj − cs−1
j

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

• Método del Máximo Descenso.

El método del máximo descenso está basado en el hecho de que

∂S s
= −bsj (35)
∂cj

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

11
se obtendrá el máximo descenso de la función S(c).

La expresión (36) se puede reescribir como

n
s
Djk ∆csk = bsj (37)
k=1

donde
s
Djk = δjk (38)

y se tiene que
cs+1 = cs + ω∆cs+1 (39)

s
Sin embargo, el método puede ser modificado de manera tal que la matriz Djk tenga dimensiones
acorde con la función S(c), y, por consiguiente, se puede hacer

s
Djk = As δjk (40)

El valor de ω se modifica de igual forma que el método de Gauss-Newton para asegurar la

convergencia, pero por el contrario el método del máximo descenso converge muy lentamente y por lo
tanto no es recomendable su uso.

• Método de Marquardt.

La formula algorı́tmica del método de Marquardt [4] es la siguiente

n
(Asjk + λDjk
s
)∆csk = bsj (41.a)
k=1

cs+1 = cs + ∆cs (41.b)

donde el factor λ funciona similar a un factor de relajación y le da al método de Marquardt un carácter

hibrido donde existe un compromiso entre el método del máximo descenso y el método de Gauss-Newton.
Cuando λ → 0, la dirección del método se dirije hacia el método de Gauss-Newton. Cuando λ → ∞, la
dirección del método se dirije hacia el método del máximo descenso.

Los estudios de Marquardt indican que el método posee un ángulo promedio entre los métodos
de Gauss-Newton y Máximo Descenso de 90o . La selección de un λ entre 0 e ∞ produce una dirección
intermedia.

12
 Para efectos de garantizar la convergencia en cada iteración se altera el factor λ de la forma

λ = λ/ρ ρ<1 (42)

hasta que se cumpla dentro de la misma iteración que

S(cs+1 ) < S(cs ) (33)

Una vez satisfecha la relación anterior en una iteración se puede disminuir λ en la siguiente iteración de
manera que
λ = λ/τ τ >1 (42 )

Nótese que incrementar λ en el método de Marquardt es equivalente a disminuir ω en el método de

Gauss-Newton.

Normalmente, se toman los valores de λinicial = 10−3 , ρ = 0.1 y τ = 10.

Cuando en varias iteraciones consecutivas el método mejora su convergencia, es decir se cumple

la relación (33), entonces λ → 0, y esencialmente se estará empleando el método de Gauss-Newton. Si
la convergencia no mejora, por el contrario, λ se incrementa y se estará usando prácticamente el método
del Máximo Descenso.

EVALUACION DEL AJUSTE.

La evaluación del ajuste viene dada mediante el análisis de de ciertos factores que permiten, por
un lado comparar cuan bueno es un ajuste en relación 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ón de aproximación.

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

Media de la variable dependiente.

1

p
fm = f (xi ) (45)
p i=1

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

p
Sm = [f (xi ) − fm ]2 (46)
i=1

Desviación estándar (σ) ó varianza (σ 2 ) con respecto a la función de aproximación.

S
σ= (47)
p−n

2
Desviación estándar (σm ) ó varianza (σm ) con respecto a la media fm .

Sm
σm = (48)
p−1

Coeficiente de determinació (r2 ) ó coeficiente de correlación (r). r2 indica el porcentaje de la

incertidumbre inicial que ha sido disminuido usando la función de aproximación.

Sm − S
r2 = (49)
Sm

En algunas literaturas definen el coeficiente de determinación (R2 ) ó coeficiente de correlación

(R) de la siguiente forma.

σm − σ
R2 = (50)
σm

Coeficiente de variación.
σm
Cv = (51)
fm

Desviación RMS (Root of the Mean Square).

S
δrms = (52)
p

14
 Desviación máxima.
δmax = max |F (xi , c) − f (xi )| (53)
1≤i≤p

En la desviación estándar σ, la cantidad S está dividida por (p − n), debido a que n parámetros
(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ón estándar σm , la cantidad Sm está dividida por (p − 1), debido a que la media de
la variable dependiente, fm , la cual se derivó de los datos originales (x1 , x2 , x3 , . . . , xp ), fué usada para
computar Sm . De aquı́ que se halla perdido un grado de libertad [5].

La desviación estándar σm debe ser mayor que σ, de otra forma no se justifica el uso de la función
de aproximación, y la media fm da una mejor aproximación a los datos que la función de aproximación
propuesta.

La función de aproximación 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ón estándar σ,
con respecto a la función de aproximación. Esto implica que el coeficiente de determinación R2 sea más
cercano a la unidad.

El coeficiente de variación Cv nos brinda una medida normalizada de cual es la dispersión de los
datos originales y normalmente se da en forma porcentual [5]. Cuando la dispersión de los datos es muy
grande significa que los puntos están muy dispersos y si se grafican formarán una nube ancha alrededor
de cualquier correlación que se trate de hallar. En este caso, la mejor correlación la ofrecerı́a la media
fm .

La desviación RMS y la desviación máxima dependen del ajuste particular que se está realizando.
La desviación RMS se puede interpretar como una desviación promedio del ajuste. La desviación máxima
acota cuanto va a hacer el mayor error cometido con el ajuste. Entre mayor sea la diferencia entre estas
dos desviaciones, mejor será el ajuste por si mismo.

15
 GRADO API Y GRAVEDAD ESPECIFICA

Se ha desarrollado una correlación para la corrección del grado API obtenido de una muestra de
laboratorio a una temperatura diferente a 60o F y mediante un procedimiento indirecto de medición de
la gravedad especı́fica.

La base de datos consistió en una tabla de corrección del grado API obtenida de un apéndice
del manual TREATING OIL EMULSIONS [6], el cual es la reproducción de una porción 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ón se obtuvo empleando la técnica numérica de ajuste de parámetros con el método

de mı́nimos cuadrados aplicado a una serie de funciones bases, o regresión lineal como también se le
conoce. Para esto se utilizó el programa de computadora denominado ADJUST, el cual realiza ajustes
de superficies (es decir, ajustes de funciones de varias variables) de forma automática. 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ólamente define las funciones
elementales que intervienen (no es necesario definir sobre que, se aplica la función elemental) y los rangos
de número de términos y factores a probar.

Inicialmente se intentó hacer el ajuste más 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ó ajustar la desviación S − S ∗ con respecto
a S y la desviación 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

   2
√ T − 60 T − 60
1 + S (S ∗ − S) = A + B −C
100 100

donde el parámetro A dió un valor despreciable, por lo que luego se eliminó, y el coeficiente de determi-
nación dió aproximadamente 98.855347 %. La eliminación de A también se justifica debido a que hace
la expresión encontrada consistente con el comportamiento esperado, es decir, para T = 60o F , S ∗ = S.

Finalmente, una vez eliminado el término independiente A, se procedió a hacer el ajuste definitivo
de la forma
    2 
∗ T − 60 T − 60 √
S =S+ B −C / 1+S
100 100

En este caso el coeficiente de determinación dió aproximadamente 98.855748 %, levemente superior al

ajuste anterior. Los parámetros dieron

B = 4.936603951 × 10−2

C = 1.480129729 × 10−3

La expresión analı́tica desarrollada aquı́para la gravedad especı́fica sirve también para los gra-
dos API si se substituyen en dicha expresión las definiciones de S y S ∗ en función de AP I y AP I ∗ ,
respectivamente.

La validez de la expresión hallada se probó haciendo conversiones de grado API a gravedad

especı́fica y viceversa, y comprobando con la tabla de corrección. 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 último, es bueno mencionar que se probó el ajuste con el término cúbico en T , pero los
resultados fueron levemente inferiores al ajuste finalmente reportado, lo que se puede ver comparando el
coeficiente de determinación que en este caso dió 98.854875 %.

17
 COEFICIENTE DE ARRASTRE PARA ESFERAS

Se desarrolló una correlación para el coeficiente de arrastre en el flujo sobre esferas que cubre
prácticamente todos los regı́menes: Laminar, turbulento y desprendimiento de la capa lı́mite. En el estu-
dio del flujo multifásico es importante esta correlación para el modelaje de los mecanismos de separación
de dos fases inmiscibles, por ejemplo, en separadores bifásicos y trifásicos, o para el modelaje de los
mecanismos de partición y coalescencia en la mezcla de dos fases inmiscibles en régimen turbulento, por
ejemplo, en mezcladores estáticos o dinámicos.

Los datos empleados en el ajuste de la correlación para el coeficiente de arrastre, CD , en función

del número de Reynolds, Re, se obtuvieron de la gráfica correspondiente publicada en el libro de texto de
Vennard y Street [7], mediante una lectura directa sobre la [Link] coeficiente de arrastre CD se define
como
F/A
CD = 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ámetro de la esfera.

U =Velocidad del flujo libre no perturbado.

y el número de Reynolds se define como

UD
Re =
ν

El problema de ajustar una sola expresión analı́tica al coeficiente de arrastre para esferas se atacó
en dos etapas. En la primera etapa se ajustaron los datos para el número de Reynolds de hasta 2 × 105
en una expresión analı́tica denominada CD1 . Esta expresión cubre los regı́menes laminar y turbulento.
Mediante el ajuste automatizado con el programa ADJUST se realizó el siguiente ajuste

ln(CD1 Re/24) = f (X) = f [ln(1 + Re)]

donde X = ln(1 + Re). Ası́ se obtuvo una expresión de la forma

A1
CD1 = 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ón es consistente en la medida que Re → 0, ya que CD1 → 24/Re.

En la segunda etapa se ajustaron los datos para el número de Reynolds desde 2 × 105 hasta
2 × 106 en una expresión analı́tica denominada CD2 . Esta segunda expresión cubre la zona donde existe
el desprendimiento de la capa lı́mite. Realizando un ajuste directo con el programa ADJUST para una
parábola se obtuvo la expresión
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ó la función escalón continua de la forma

A
CD = CD1 + (CD2 − CD1 )
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én el
resultado final de la correlación para CD en lı́nea continua. Como puede observarse el ajuste es bastante
bueno.

19
Figura 1. Correlación para el coeficiente de arrastre sobre una esfera.

20
 CORRELACION DE EATON PARA LA FRACCION DE LIQUIDO

La correlación de Eaton [8] es muy empleada en la predicción de la fracción de lı́quido para el

flujo bifásico. Sin embargo, es bueno mencionar que normalmente se emplea esta correlación haciendo
la lectura directa de puntos sobre una gráfica y luego realizando una interpolación polinómica, lo cual
ocupa mucho espacio (memoria) y tiempo (cálculo) por parte del usuario de la correlación o por parte de
la computadora.

Por la razón antes mencionada se ha desarrollado una expresión analı́tica que permita rápidamen-
te realizar los cálculos de la fracción de lı́quido. Debido a la forma muy particular de la curva en la gráfica,
se decidió realizar con el programa ADJUST un ajuste no lineal empleando la función y = tanh(x). El
estudio de los resultados de varios de estos ajustes dió como la mejor alternativa la expresión

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

donde

A=0.9377978363

B=1.034705148

Aquı́ Y representa a la fracción de lı́quido y X representa al número adimensional de Eaton. El coeficiente

de determinación fue en este caso del 99.942265 % y la desviación máxima fue de 0.0175 en valor absoluto.

El figura 2 muestra la correlación de Eaton representando con puntos los valores tomados por
lectura directa de la gráfica original y con lı́nea continua la expresión analı́tica desarrollada aquı́.

21
Figura 2. Correlación de Eaton para la fracción de lı́quido.

22
 VISCOSIDAD DE LOS GASES NATURALES

En esta parte se ha desarrollado una correlación para determinar la viscosidad µg (cP ) a partir
de la presión P (psia) y la temperatura T (o F ). Aquı́ se ha empleado polinomios de interpolación
dependientes de una variable y de dos variables hasta de séptimo y tercer órdenes, 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 (psia)

Tc = 169 + 314γ̂g (o R)

donde
ρg
γg = (aire → γg = 1)
ρa(C.N.)

Las propiedades reducidas se calculan de acuerdo a su definición como

P
Pr =
Pc

T
Tr =
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ón atmosférica µog se calcula empleando polinomios de interpolación
de séptimo orden en función del peso molecular y polinomios de primer orden en función de la temperatura
obtenidos de la figura 4 superior, presentada originalmente en el trabajo de Carr y otros [10]. Ası́ resulta
que

T̃ − 40
µog = (B − D) + D
360
donde
7


B= bi Mgi
i=0

23
 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

7


D= 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ón atmosférica µog puede sufrir una corrección como lo indica la figura
4 superior debido a impuresas como nitrógeno (N2 ), dióxido de carbono (CO2 ) y sulfuro de hidrógeno
(H2 S). El contenido de estas impurezas en el gas se expresa como una fracción molar y(·) . Asi se obtiene
una viscosidad corregida µ̃og a presión atmosférica

µ̃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ándar.

La viscosidad asi obtenida a su vez se vuelve a corregir en función de la presión y la temper-

atura empleando la figura 4 inferior. De aquı́ se deriva la siguiente relación que contiene polinomios de
interpolación en dos direcciones: la presión reducida y la temperatura reducida. Esto es,

µg exp A
=
µ̃og Tr

donde

3


A= 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

µg exp A
µg = µ̃og = µ̃og
µ̃og Tr

La correlación para µ̃og es valida sólamente para el rango de temperatura 40o F < T < 400o F .
La correlación para µg /µ̃og tiene la restricción de que fué obtenida apartir de datos de presión reducida
en el rango 1 < Pr < 20.

26
Figura 3. Propiedades pseudocrı́ticas de los gases naturales.

27
Figura 4. Correlación 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 Emul-
sions”. 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