1
Desempeño en CFD del modelo euleriano calculando el coeficiente de arrastre a
partir del modelo Drift-Flux para flujos liquido-gas en tuberías verticales
Esteban Guerrero Andrioli
Proyecto de Grado 2014-20
Departamento de Ingeniería Química
Universidad de los Andes
Bogotá, Colombia
Resumen
El flujo bifásico tiene una gran aplicación a nivel industrial. Además la caracterización de éste resulta
de gran utilidad para aumentar la eficiencia y rentabilidad de los procesos. Por ejemplo, las industrias
petroleras necesitan conocer la cantidad de gas natural y crudo provenientes de un reservorio para
calcular la rentabilidad de la producción. Sin embargo, describir correctamente las propiedades que
definen este tipo de flujo requiere de mecanismos eficientes como son las herramientas
computacionales. Actualmente, la mecánica de fluidos computacional (CFD) es la técnica más
utilizada para estos análisis. Para modelar el fluido bifásico comúnmente se utiliza en esta herramienta
el modelo de volumen de fluido (VOF) y el modelo euleriano. La formulación matemática de estos
modelos genera una diferencia en su convergencia, rapidez y exactitud. El modelo VOF no discrimina
las velocidades de las fases y el modelo euleriano permite acoplar diferentes modelos como el Drift-
Flux para mejorar la precisión. Este trabajo demuestra el desempeño de acoplar el modelo euleriano
y el modelo Drift-Flux. Con el fin de llevar a cabo este objetivo se debe seleccionar el mejor modelo
de Drift-Flux, encontrar las diferencias entre el modelo euleriano y el modelo VOF, y validar estos
modelos con resultados experimentales. En este proyecto se realizó la modelación de seis
experimentos obteniendo que el modelo euleriano presenta un error cuadrático medio (13.86%)
inferior al modelo VOF (19.04%) para flujos con baja fracción de vacío. Por otro lado, se evidenció
que el modelo euleriano es independiente de la malla permitiendo un tiempo de simulación menor al
del modelo VOF. Adicionalmente, se encontró que el mejor modelo de Drift-Flux reportado es el de
Maier et al. (1997). Finalmente, se obtuvo que el modelo euleriano presenta divergencia al calcular
el coeficiente de arrastre con solo un término del modelo de Drift-flux.
Palabras clave: CFD, Flujo bifásico, Drift-flux, Coeficiente de Arrastre, Fracción de vacío.
Nomenclatura
𝐶 Número de Courant [-]
𝐶0 Parámetro de distribución [-]
𝐶𝐷 Coeficiente de arrastre [-]
𝑑 Diámetro de burbuja [m]
𝐷 Diámetro de la tubería [m]
𝑔 Gravedad [= 9.81 m/s2]
𝑖, 𝑗 Índices de fases [-]
𝑃 Presión [Pa]
𝑅𝑒 Número de Reynolds [-]
𝑡 Tiempo [s]
𝑢𝐷 Velocidad de deriva [m/s]
𝑢𝑀 Velocidad de mezcla [m/s]
𝑢𝑔 Velocidad superficial del gas [m/s]
𝑢𝑙 Velocidad superficial del líquido [m/s]
𝑢𝐺 Velocidad total del gas [m/s]
𝑊 Flujo másico [kg/s]
𝛼 Fracción de vacío [-]
𝜌 Densidad [kg/m3]
𝜎 Tensión superficial [N/m]
𝜏 Esfuerzo molecular [Pa]
𝜏𝑡 Esfuerzo turbulento [Pa]
2
1. INTRODUCCIÓN
El flujo multifásico es una condición de operación frecuente en la industria. Este tipo de flujo se
puede encontrar en sistemas de generación de energía, transporte de material, transferencia de calor,
equipos de separación o reacción y equipos de control ambiental (Ishii & Hibiki, 2011).
Adicionalmente, el flujo multifásico se presenta en sistemas externos a la industria como la
circulación de sangre y oxígeno en el cuerpo humano.
1.1 Planteamiento del problema La industria petrolera ha encontrado una gran importancia en caracterizar el flujo multifásico debido
a la alta presencia en sus tres sectores: upstream, midstream y downstream. Estos sectores se encargan
respectivamente de la producción, el transporte y la refinación del petróleo. En la primera parte de
esta cadena de operación se obtiene una mezcla de agua, crudo y gas del reservorio, haciendo uso de
tuberías verticales, horizontales o inclinadas. En pozos maduros es común realizar un levantamiento
artificial mediante la inyección de vapor, agua o gas para incrementar la producción del pozo (Zhang
et al., 2003). En este orden de ideas, el sector de upstream trabaja con flujo bifásico de gas-líquido.
Las composiciones de la mezcla que definirán las propiedades del flujo bifásico dependerán de la
naturaleza del reservorio. Además, se debe tener en cuenta que a medida que aumenta el tiempo de
producción del pozo las composiciones del flujo cambian. Frecuentemente la composición del agua
aumenta con el paso del tiempo. Sin embargo, la variación en la composición lleva a severos
problemas operacionales, desde una producción incontinua hasta el apagado del proceso (Abdulkadir,
2011). Adicionalmente, esta variación tiene efectos económicos en la producción, esto se debe a que
según los valores comerciales la rentabilidad del pozo depende principalmente de la cantidad de
barriles de crudo producidos más no de la cantidad de gas natural extraído. Por esta razón, la industria
petrolera tiene como necesidad caracterizar con precisión las velocidades y fracciones de gas
presentes en el flujo bifásico.
1.2 Patrones de flujo
Las condiciones de operación del flujo bifásico en una tubería permiten que éste obtenga una
hidrodinámica particular. Este comportamiento se ve afectado principalmente por las diferentes
velocidades que presentan la fase líquida y gaseosa. A este tipo de hidrodinámica se le conoce como
patrones de flujo. De esta manera, el diseño de herramientas para caracterizar el flujo bifásico deben
considerar las diferentes configuraciones que presentan las fases para conocer las propiedades de
interés.
Los patrones de flujo difieren dependiendo de la orientación de la tubería. En este proyecto se
estudiará el flujo bifásico en tuberías verticales para la producción de petróleo. Los únicos patrones
de flujo en esta orientación son: flujo de burbuja dispersa, flujo tapón, flujo agitado, flujo anular y
flujo neblina. La configuración de las fases en estos patrones se muestra en la Figura 1. El flujo
neblina se omitió en la figura debido a que presenta la misma apariencia que el flujo de burbuja
dispersa a diferencia que las fases están invertidas (Abdulkadir, 2011). Sin embargo, se debe tener en
cuenta que estos patrones de flujo se vuelven inestables cuando el diámetro de la tubería supera los
50 mm de diámetro (Sharaf & Luna-Ortiz, 2014).
3
Figura 1. Patrones de flujo en tuberías verticales. a) Flujo burbuja dispersa. b) Flujo tapón. c) Flujo agitado. d) Flujo
anular (Bratland, 2010).
Un bajo flujo de gas en una tubería donde se presenta un flujo bifásico gas-líquido causa que la fase
gaseosa se comporte como burbujas aisladas en la fase continua. Estas burbujas se caracterizan por
ser esféricas y significativamente menores al diámetro del tubo (Thome, 2004). La configuración
mencionada anteriormente se conoce como flujo de burbuja dispersa (Figura 1.a). Una vez el flujo de
gas empieza a incrementar, la cantidad de burbujas aumenta llevándolas a coalescer hasta alcanzar
un tamaño cercano al área transversal de la tubería. Estas burbujas presentan una forma ondulada y
plana en sus extremos superior e inferior respectivamente. El incremento del tamaño en las burbujas
lleva a que se presente en contraflujo una película de líquido entre la burbuja de gas (conocidas como
burbujas de Taylor o en forma de bala) y la pared de la tubería. Este patrón (Figura 1.b) se conoce
como flujo tapón (Shang, 2014). Incrementando aún más el flujo de gas se desestabiliza el sistema
generando un flujo agitado (Figura 1.c). Este patrón consiste en que la película de líquido del flujo
tapón se desestabiliza cayendo a la bala líquida entre burbujas, donde vuelve a tomar una dirección
ascendente hasta volver a desestabilizarse (Abdulkadir, 2011). El patrón de flujo anular (Figura 1.d)
se obtiene cuando el incremento del flujo de gas es tal que el estrés interfacial es mayor a los efectos
de gravedad, expulsando el líquido del centro de la tubería. El aspecto de este patrón consiste en una
película delgada de líquido pegada a la pared mientras en el centro fluye una corriente de gas, la
interfaz líquido-gas de este patrón se mantiene ondulada (Thome, 2004). Por último, el flujo neblina
se genera cuando las ondulaciones de flujo anular crecen de manera significativa debido al incremento
del flujo de gas. Estas ondulaciones producen que gotas de la fase líquida entren de manera dispersa
en la fase gaseosa hasta terminar con la película de líquido (Shang, 2014).
Las condiciones de operación de un flujo bifásico pueden ubicarse en mapas de patrones de flujo para
predecir en cuál de los patrones descritos anteriormente se encuentra el proceso. En la literatura se
expone una gran variedad de mapas empíricos. Comúnmente, estos consisten de ejes dependientes de
parámetros como la densidad, viscosidad, tensión superficial, velocidad superficial, entre otros. En la
Figura 2 se muestra el mapa propuesto por Hewitt et al. (1986) para tuberías verticales.
1.3 Herramientas experimentales
La caracterización del flujo bifásico se puede realizar mediante el uso de métodos experimentales.
Principalmente se utiliza tomografías eléctricas, sondas de agujas, sondas ópticas, sondas ultrasónicas
y sensor de malla de alambre por conductancia o conductividad. Estos métodos tienen la capacidad
de medir en una posición específica la fracción de vacío y las velocidades de cada una de las fases.
De acuerdo a estos resultados es posible reportar el patrón de flujo en el que se encuentra el sistema.
Sin embargo, estos métodos presentan inconvenientes en la experimentación. Las tomografías
4
eléctricas tienen baja resolución, las sondas de agujas y el sensor de malla de alambre deben
encontrarse en contacto con el fluido causando su desgaste, y las sondas ópticas deben tener una
fracción de gas pequeña y un acceso óptico (Sharaf, 2014).
Otro problema relacionado al uso de estos sistemas es que sus metodologías están basadas en datos
experimentales tomados en laboratorios que comúnmente utilizan líquidos con viscosidades menores
a 0.02 Pa.s y bajas presiones de 1 a 3 atm. Por el contrario, la producción de petróleo maneja en sus
procesos viscosidades mayores a 10 Pa.s y presiones altas alrededor de 31 atm (Jeyachandra et al.,
2012). Por otro lado, el gas del flujo bifásico a nivel de laboratorio presenta un comportamiento ideal
debido a la baja presión. No obstante, la producción de crudo se opera a altas presiones y caídas de
presión que alteran este comportamiento. Por este motivo, el uso de estas herramientas en la
caracterización del flujo debe validarse en operaciones industriales o ser sustituidas por herramientas
más eficaces.
1.4 Herramientas computacionales
Las limitaciones experimentales están siendo superadas por el manejo de la Mecánica de Fluidos
Computacional (CFD, por sus siglas en inglés – Computational Fluid Dynamics). Esta herramienta
usa procesadores para resolver a una gran velocidad los modelos matemáticos que describen un
sistema. De esta manera es posible caracterizar un flujo bifásico utilizando ecuaciones de mecánica
de fluidos. La mayoría de los modelos usados en CFD son un conjunto de ecuaciones que presentan
derivadas parciales en su estructura. Por este motivo, esta herramienta discretiza estos términos
utilizando el método de volúmenes finitos. La principal ventaja de usar CFD es la posibilidad de
simular diferentes ambientes y condiciones de operación sin tener en cuenta los requisitos rigurosos
y costos que presentan las herramientas experimentales (Abdulkadir, 2011). De esta manera, es
posible plantear un modelo para describir cualquier sistema una vez se haya validado con resultados
experimentales a unas condiciones de operación específicas.
El modelamiento y simulación en CFD de un flujo bifásico se realiza principalmente con el modelo
euleriano o el modelo de volumen de fluido (VOF, por sus siglas en inglés – Volume Of Fluid). Estos
modelos se diferencian por sus sistemas de ecuaciones, generando ventajas y desventajas el uno sobre
el otro. El modelo euleriano cuenta con un set de ecuaciones de conservación de momento y de masa
para cada una de las fases presentes. Sin embargo, es necesario definir la velocidad con la cual el
coeficiente de fuerza de arrastre va a ser evaluado. En este proyecto se propone utilizar el modelo de
flujo de deriva (Drift-Flux de ahora en adelante) para extraer la correlación que describe la velocidad
de deriva. Este término hace referencia a la velocidad del gas cuando la fase líquida se mantiene en
reposo. El modelo de Drift-Flux sigue la estructura mostrada en la Ecuación (1).
𝑢𝐺 = 𝐶0𝑢𝑀 + 𝑢𝐷 (1)
Donde 𝑢𝐺 es la velocidad total del gas, 𝑢𝑀 es la velocidad de mezcla que resulta siendo la suma de
las velocidades superficiales del gas y líquido, 𝐶0 es el parámetro de distribución y 𝑢𝐷 es la velocidad
de deriva. Los dos últimos parámetros dependen de la fracción volumétrica de la fase gaseosa en la
posición evaluada, termino conocido como fracción de vacío. Por esta razón, existen diversas
correlaciones según el patrón de flujo en el que se encuentre el sistema. No obstante, la mayoría de
los modelos de Drift-Flux planteados en la literatura son basados en datos experimentales (Choi et
al., 2012). Debido a esto la robustez del modelo euleriano puede reducirse. Por otro lado, el modelo
VOF únicamente tiene un set de ecuaciones de conservación de momento y de masa para las dos
fases. De esta manera las propiedades del fluido deben ser calculadas usando la teoría de propiedades
parciales y debe agregarse una ecuación de conservación de fracción de volumen. Sin embargo, este
modelo matemático presenta problemas en el cálculo de la velocidad de deslizamiento entre fases
debido a que comparte las ecuaciones de momento en la interface y no discrimina entre fases. Por
5
ejemplo, en la literatura se evidencia que este método presenta inconvenientes en capturar la película
de líquido presente en el flujo tapón (Horgue et al., 2012). Por esta razón, usar el modelo VOF
requiere de una mallada más fina (mayor tiempo computacional) que el modelo euleriano con el fin
de identificar correctamente las interfaces.
1.5 Objetivos del proyecto
1.5.1 Objetivo general
El trabajo presentado demuestra el desempeño en CFD del modelo euleriano calculando el coeficiente
de arrastre a partir de correlaciones de Drift-Flux para flujos bifásicos en tuberías verticales. La
evaluación del desempeño se discutirá a partir de dos criterios. En primer lugar, se realizaron una
comparación con el tiempo computacional gastado por el modelo VOF, comúnmente usado en las
investigaciones actuales. En segundo lugar, se evaluó la exactitud del modelo al validar los resultados
computacionales con resultados experimentales.
1.5.2 Objetivos específicos
I. Seleccionar de la base de datos de Pineda et al. (2014) y Coddington & Macian (2002) el
modelo Drift-Flux más adecuado para las condiciones de operación estudiadas.
II. Establecer las diferencias en la modelación de flujo bifásico utilizando el modelo euleriano
y el modelo VOF.
III. Validar los modelos físicos utilizando los trabajos de Krepper et al. (2005), Sun et al. (2004),
y Westende (2008). Con el fin de realizar esta comparación, las simulaciones se realizaron
para un flujo bifásico de agua-aire.
2. MATERIALES Y MÉTODOS
La simulación de flujo bifásico en CFD para tuberías verticales debe considerar ciertos aspectos para
desarrollarse adecuadamente. En este capítulo se menciona las dimensiones de las tuberías y
condiciones de operación a modelar. Adicionalmente, se determina cómo generar el mejor mallado y
especificar un paso de tiempo congruente a éste. Por último, se explica la formulación matemática de
los modelos analizados en este proyecto.
2.1 Geometrías y condiciones de operación
El desempeño de los modelos en CFD para la caracterización de flujos bifásicos se validó utilizando
resultados experimentales llevados a cabo con diferentes técnicas de medición. Sun et al. (2004)
utilizaron un anemómetro láser Doppler para medir la fracción de vacío y los flujos de las dos fases.
Krepper et al. (2005) calcularon las fracciones de vacío mediante señales de corriente que fueron
medidas a través de un sensor de malla de alambre por conductividad. Por último, Westende (2008)
midió las velocidades utilizando un anemómetro de fase Doppler. Con el fin de replicar cada uno de
los experimentos en CFD se utilizó el software computacional STAR-CCM+ v9.02 de la compañía
CD-adapco. Las dimensiones de los tubos y condiciones de operación se muestran en la Tabla 1,
donde 𝑢𝑖 representa la velocidad superficial de la fase 𝑖, 𝑧 es la altura del tubo y la relación z/D
describe la ubicación de la herramienta de medición en la tubería.
6
Tabla 1. Geometrías y condiciones de operación.
Autor Diámetro (m) Longitud (m) Caso ul (m/s) ug (m/s) Sensor (z/D)
Sun et al.
(2004) 0.0508 3.81
A 0.6150 0.049 51
B 3.4580 0.318 51
Krepper et al.
(2005) 0.0512 4.00
C 1.0000 0.220 60
D 1.0000 0.340 60
Westende
(2008) 0.0500 8.00
E 0.0394 12.200 39
F 0.0411 21.200 100
Los patrones de flujo trabajados en este proyecto pueden ser predichos a partir de un mapa de patrones
de flujo. En este caso, se utilizó la gráfica reportada por Hewitt et al. (1986) debido a que maneja las
mismas sustancias. Además, ésta se utiliza para diámetros similares a los trabajados en los casos de
estudio (≈0.05 m). Los experimentos que se simularon, descritos en la Tabla 1, se ubicaron en este
mapa de patrones de flujo como se muestra en la Figura 2. De esta manera, es posible predecir que
en este proyecto se manejan flujos de burbuja dispersa y flujos anulares. En otras palabras, se modeló
flujos bifásicos con altas y bajas fracciones de vacío.
Figura 2. Ubicación de casos de estudio en el mapa de patrones de flujo (Hewitt et al., 1986).
2.2 Generación del mallado
La solución del sistema de ecuaciones que describe un flujo bifásico requiere de la discretización de
la geometría para resolver las derivadas parciales. Esta discretización se conoce como mallado y se
conforma de múltiples unidades llamadas celdas. Las dimensiones y el arreglo espacial de estas
permite establecer una gran variedad de mallados para una misma geometría. No obstante, la
convergencia, rapidez y exactitud de la solución dependen de la calidad de la malla realizada.
7
En este proyecto se utilizaron geometrías simples de forma cilíndricas. Sin embargo, los modelos
físicos tratados en este sistema aumentan la complejidad del problema, causando que la calidad del
mallado sea de gran importancia. El trabajo de Hernandez et al. (2010) determinó que en una tubería
la mejor distribución de mallado en la sección transversal es de forma mariposa. Haciendo uso de este
tipo de mallado la simulación presenta una mayor eficiencia que el uso de distribuciones rectangulares
(Mallado-H), circulares (Mallado-O) o amorfas. En la Figura 3 se ilustra la distribución de cada uno
de los tipos de malla mencionados anteriormente.
Figura 3. Tipos de distribución de malla. a) Circular (Mallado-O), b) Rectangular (Mallado-H), c) Mariposa, d) Amorfa.
La selección del mallado se realizó a partir de un test de independencia de mallado con el fin de
eliminar la influencia en la solución y disminuir el tiempo computacional requerido. De esta manera,
se simuló el Caso D utilizando cuatro mallados con los siguientes números de celdas: 43400, 228780,
312800 y 415140. El criterio de selección fue la exactitud con respecto a los datos experimentales y
el menor tiempo computacional gastado. Este test se realizó para el modelo euleriano y el modelo
VOF.
2.3 Criterio de estabilidad
El comportamiento de los patrones de flujo no permite alcanzar un sistema en estado estacionario.
Por esta razón, la modelación de flujos bifásicos se realiza en estado transitorio, en donde se observa
la evolución y fluctuaciones del flujo con respecto al tiempo. La convergencia de este modelo depende
crucialmente del paso de tiempo establecido. La inestabilidad del sistema ocurre cuando el paso de
tiempo presenta una magnitud elevada en relación a la velocidad, tal que el flujo recorre grandes
espacios de celdas sin resolver puntos intermedios. Esta situación produce que en las celdas sin
solución se supongan valores para la siguiente iteración, llevando el sistema a divergir (Abdulkadir,
2011). Estos problemas son evitados al calcular un paso de tiempo apropiado mediante la definición
del número de Courant. Este número representa la relación entre el paso de tiempo de simulación y
el tiempo de residencia del fluido en la celda estudiada. El valor recomendable para el número de
Courant es de 0.25 según la documentación de programas computacionales de CFD como STAR-
CCM+ y ANSYS Fluent. La Ecuación (2) describe matemáticamente este número adimensional.
𝐶 = 𝑢𝐺Δ𝑡
Δ𝑥 (2)
Donde 𝐶 es el número de Courant, Δ𝑡 es el paso de tiempo y Δ𝑥 es la altura de una celda. La velocidad
𝑢𝐺 es calculada mediante el modelo Drift-Flux (Ujang, 2008) descrito en la Ecuación (3).
𝑢𝐺 = (1.2 +0.8
1+10𝑒−8𝑅𝑒𝑠2.55) 𝑢𝑀 + 0.35√𝑔𝐷 (3)
Para la ecuación anterior 𝑔 es la gravedad, 𝑅𝑒𝑠 es el número de Reynolds de la fase líquida y D es el
diámetro de la tubería.
8
2.4 Modelos matemáticos
Los modelos de flujo multifásico VOF y euleriano sólo considerarán la transferencia de momento y
masa. La transferencia de calor se omitirá asumiendo que la temperatura del sistema permanece
constante y uniforme a lo largo de la tubería. En esta sección se muestra la formulación matemática
de los modelos físicos estudiados en este proyecto (STAR-CCM+, 2014).
2.4.1 Modelo euleriano
Las Ecuaciones (4) y (5) representan la conservación de masa y la conservación de momento
respectivamente. En este modelo se presenta un set de estas ecuaciones por fase existente.
𝜕
𝜕𝑡(𝛼𝑖𝜌𝑖) + ∇(𝛼𝑖𝜌𝑖𝑢𝑖) = 0 (4)
𝜕
𝜕𝑡(𝛼𝑖𝜌𝑖𝑢𝑖) + ∇(𝛼𝑖𝜌𝑖𝑢𝑖𝑢𝑖) = −𝛼𝑖∇𝑃 + 𝛼𝑖𝜌𝑖𝑔 + ∇[𝛼𝑖(𝜏𝑖 + 𝜏𝑖
𝑡)] + 𝑀𝑖 (5)
Adicionalmente, debe cumplirse la relación descrita en la Ecuación (6).
∑ 𝑀𝑖 = 0 (6)
Para las ecuaciones anteriores 𝑖 representa la fase estudiada, 𝛼 es la fracción de volumen, 𝑢 es la
velocidad superficial, 𝑃 es la presión, 𝜏 es el esfuerzo molecular, 𝜏𝑡 es el esfuerzo turbulento, 𝜌 es la
densidad y 𝑀𝑖 representa la transferencia de momento en la interfase. Este último término puede
calcularse mediante la Ecuación (7).
𝑀𝑖 = ∑ (𝐹𝑖𝑗𝐷 + 𝐹𝑖𝑗
𝐿 + 𝐹𝑖𝑗𝑊𝐿 + 𝐹𝑖𝑗
𝑇𝐷)𝑗≠𝑖 (7)
Cumpliéndose que 𝐹𝑗𝑖 = −𝐹𝑖𝑗. En este proyecto sólo se consideraron cuatro fuerzas de interacción de
fase. En primer lugar, la fuerza de dispersión turbulenta 𝐹𝑖𝑗𝑇𝐷 hace referencia al efecto de la
turbulencia en la redistribución de la concentración de fases en zonas no uniformes. El segundo
término corresponde a la fuerza de lubricación en la pared 𝐹𝑖𝑗𝑊𝐿. Ésta representa la fuerza de repulsión
ejercida por la fase líquida hacia las burbujas evitando el contacto de estas con la pared. La tercera
fuerza a considerar es la de sustentación 𝐹𝑖𝑗𝐿 , la cual es la fuerza ejercida perpendicularmente a la
velocidad de una burbuja cuando el flujo de la fase continua no es uniforme. Finalmente, la fuerza de
arrastre 𝐹𝑖𝑗𝐷 corresponde a la fuerza de fricción ejercida en la dirección opuesta al movimiento de las
burbujas. Esta última se calcularon utilizando el modelo propuesto por Burns et al. (2004) descrito en
la Ecuación (8).
𝐹𝑖𝑗𝐷 =
3
4𝐶𝑖𝑗
𝐷 𝛼𝑗𝜌𝑖
𝑑𝑗(𝑢𝑗 − 𝑢𝑖)|𝑢𝑗 − 𝑢𝑖| (8)
El coeficiente de arrastre 𝐶𝑖𝑗𝐷 es calculado utilizando la correlación de Schiller-Naumann mostrada en
la Ecuación (9).
𝐶𝑖𝑗𝐷 = {
24
𝑅𝑒𝑑(1 + 0.15𝑅𝑒𝑑
0.687) 0 < 𝑅𝑒𝑑 ≤ 1000
0.44 𝑅𝑒𝑑 > 1000 (9)
9
Donde 𝑅𝑒𝑑 es el número de Reynolds de la fase dispersa. Sin embargo, otra forma de calcular este
coeficiente es mediante la Ecuación (10), expresión propuesta por Krishna et al. (1999). Esta ecuación
permite acoplar el modelo de Drift-Flux, el cual ha podido representar con gran exactitud la realidad.
𝐶𝑖𝑗𝐷 =
4
3
𝜌𝑖−𝜌𝑗
𝜌𝑖 𝑔𝑑𝑗
1
𝑢𝐷2 (10)
La velocidad de deriva 𝑢𝐷 se calcularon mediante las mejores correlaciones de Drift-Flux presentadas
en la base de datos de Pineda et al. (2014) y Coddington & Macian (2002). La selección de éste
dependió del error obtenido al compararse con datos experimentales.
2.4.2 Modelo VOF
Las Ecuaciones (11) y (12) representan la conservación de masa y la conservación de momento
respectivamente. En este modelo se presenta un único set de estas ecuaciones para todas las fases.
𝜕𝜌
𝜕𝑡+ ∇(𝜌𝑢) = 0 (11)
𝜕
𝜕𝑡(𝜌𝑢) + ∇(𝜌𝑢𝑢) = −∇𝑃 + 𝜌𝑔 + ∇(𝜏 + 𝜏𝑡) (12)
La densidad y la viscosidad se calculan en relación a la fracción de volumen como se evidencia en
las Ecuaciones (13) y (14) respectivamente.
𝜌 = ∑ 𝜌𝑖𝛼𝑖𝑖 (13)
𝜇 = ∑ 𝜇𝑖𝛼𝑖𝑖 (14)
Adicionalmente, las interfases presentes se modelan a partir de la ecuación de continuidad respecto a
la fracción de volumen como se expresa en la Ecuación (15).
𝜕𝛼𝑖
𝜕𝑡+ 𝑢∇(𝛼𝑖) = 0 (15)
3. RESULTADOS Y DISCUSIÓN
En este capítulo se reporta los resultados obtenidos para cumplir cada uno de los objetivos planteados.
En primer lugar, se muestra el análisis llevado a cabo para escoger el mejor modelo de Drift-Flux en
tuberías verticales. En segundo lugar, se realiza el test de independencia de mallado. El tercer análisis
consiste en reportar los resultados obtenidos para cada uno de los casos de estudio. Por último, se
determina el desempeño del acople entre el modelo de Drift-Flux y el modelo euleriano.
3.1 Selección de modelos Drift-Flux
Utilizando el programa computacional de Matlab v2014a llevado a cabo por el trabajo de Pineda et
al. (2014) se analizó la exactitud de 19 correlaciones empíricas de Drift-Flux. Este análisis fue
realizado junto a la base de datos que se resume en la Tabla 2. En la Tabla A1 se encuentra la
formulación matemática de cada uno de los modelos estudiados y su respectivo error porcentual al
intentar predecir los datos experimentales. El error fue calculado utilizando la Ecuación (16).
𝐸𝑟𝑟𝑜𝑟 = [1
𝑁∑ (
𝛼𝑐𝑎𝑙𝑐𝑢𝑙𝑎𝑑𝑎𝑖−𝛼𝑒𝑥𝑝𝑒𝑟𝑖𝑚𝑒𝑛𝑡𝑎𝑙𝑖
𝛼𝑒𝑥𝑝𝑒𝑟𝑖𝑚𝑒𝑛𝑡𝑎𝑙𝑖
)𝑁𝑖=1 ] ∗ 100% (16)
10
Tabla 2. Resumen de la base de datos utilizada.
Autor Puntos D (m) z/D ul (m/s) ug (m/s) Fluidos
Magrini
(2009) 20 0.0762 229.66
0.0034-
0.04
41.5-
79.3
Agua -
Aire
Felizola
(1992) 8 0.0510 294.11
0.05-
0.5
0.39-
1.88
Queroseno -
Aire
Schmidt
(2008) 87 0.0545 56.88
0.01-
3.42
0.05-
29.58
Providona y Agua -
Nitrógeno
Rosa
(2010) 73 0.0260 306.00
0.22-
3.08
0.12-
28.8
Agua-
Aire
Ghajar
(2012) 105 0.0127 173.00
0.08-
1.17
0.444-
20.059
Agua -
Aire
Majumder
(2010) 99 0.0191 178.48
0.585-
2.339
0.084-
0.699
Agua, Alcohol Amílico,
Glicerina A, Glicerina B -
Aire
La desviación de la realidad fue calculada mediante el error cuadrático medio entre los datos
predichos por la correlación y datos experimentales. En la Figura 4 se muestra la predicción de la
correlación de Maier (1997) que obtuvo el menor error porcentual.
Figura 4. Predicción de la correlación de Maier (1997) con respecto a datos experimentales.
Las cuatro correlaciones con mejor exactitud se mencionan en la Tabla 3. Adicionalmente, se
describen las ecuaciones (17) a (20) que representan sus respectivas expresiones matemáticas para
calcular la velocidad de deriva. Estas ecuaciones se acoplaron al modelo euleriano para llevar a cabo
la simulación. Por otro lado, las predicciones de las tres correlaciones restantes se muestran en las
Figura A1, A2 y A3.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
exp
corr
Maier
0
-30%
+30%
Error[%]27.07
11
Tabla 3. Correlaciones de Drift-Flux seleccionadas.
Autor Velocidad de deriva, 𝒖𝑫 (m/s) Ecuación Error (%)
Maier et al.
(1997) 𝑢𝐷 = (𝑣1𝑃2 + 𝑣2𝑃 + 𝑣3)𝐺 + (𝑣4𝑃2 + 𝑣5𝑃 + 𝑣6) (17) 27.07
Sonnenburg
(1989) 𝑢𝐷 =
𝐶𝑜(1 − 𝐶𝑜𝛼)
(𝐶𝑜𝛼 √𝑔∆𝜌 𝜌𝐺⁄⁄ ) + (1 − (𝐶𝑜𝛼 √𝑔∆𝜌 𝜌𝐿⁄⁄ )) (18) 29.25
Mattar et al.
(1974) 𝑢𝐷 = 0.7 (19) 33.37
Bestion
(1985) 𝑢𝐷 = 0.188√
𝑔𝐷(𝜌𝐿 − 𝜌𝑔)
𝜌𝑔
(20) 33.56
Los errores obtenidos con las correlaciones mostradas en la Tabla 3 se encuentran alrededor del 30%.
Este promedio resulta aceptable al estudiar la fracción de vacío, la cual es una variable que toma
valores pequeños (0 a 1). Sin embargo, existen herramientas computacionales que tienen modelos de
Drift-Flux más eficientes. Este análisis se realizó utilizando los programas OLGA v7.3 y LedaFlow
v1.2 de las compañías Schlumberger y Kongsberg respectivamente. Estas herramientas permiten
modelar el flujo bifásico para grandes trayectorias en una dimensión. En la Figura 5 se muestra el
resultado de estos modelos al predecir la base de datos resumida en laTabla 2. La predicción de estas
herramientas computacionales resulta más exacta que las correlaciones reportadas en la literatura al
obtener un error menor al 20%. No obstante, la privacidad de las compañías no permite obtener la
totalidad de sus modelos para el uso en este proyecto.
Figura 5. Predicción de los software OLGA y LedaFlow con respecto a datos experimentales.
12
3.2 Mallado de las geometrías
La modelación del experimento Caso D fue utilizado para realizar el test de independencia de
mallado. Krepper et al. (2005) utilizaron un sensor ubicado en z/D=60 para medir la fracción de vacío
con una alimentación de 𝑢𝑔 = 0.34 𝑚/𝑠 y 𝑢𝑙 = 1.00 𝑚/𝑠. La fracción de vacío promedio fue de
0.2618 con una desviación estándar del 10%. Los resultados obtenidos mediante el modelo VOF y el
modelo euleriano se muestran en la Figura 6.
Figura 6. Test independencia de mallado - Validación experimental
En la Figura 6 es evidente para el modelo VOF que entre más fina la malla se obtiene menor error
entre la respuesta de la simulación y lo datos experimentales. Sin embargo, siguiendo el criterio que
la desviación estándar es del 10%, sólo la malla de 415140 celdas permite modelar correctamente el
sistema. Por el contrario, los resultados del modelo euleriano muestran que la malla no tiene un
impacto significativo en la exactitud de los resultados. Además, se observa que este modelo obtiene
un error igual a la desviación estándar de los resultados experimentales.
El segundo criterio en la selección de malla es el tiempo de simulación. Éste fue analizado usando
únicamente un núcleo del procesador. En la Figura 7 se muestran los resultados obtenidos en este
estudio. Es evidente que para los dos modelos físicos se requiere de un mayor tiempo computacional
a medida que la malla es más fina. En este orden de ideas, la malla elegida para el modelo euleriano
es de 43400 celdas con el fin de reducir el tiempo computacional sin afectar la exactitud de los
resultados. Por otro lado, la malla seleccionada para el modelo VOF es de 415140 celdas con el fin
de tener una buena exactitud en los resultados a pesar de presentar un alto tiempo de simulación.
En la Figura 7 se compara el tiempo computacional gastado por el modelo euleriano frente al modelo
VOF. La simulación realizada requiere de 62000 iteraciones para completar el tiempo físico del
problema. En este estudio se comprobó que el modelo euleriano siempre presenta mayor tiempo de
simulación que el modelo VOF. Esto se debe a que el modelo euleriano presenta una mayor cantidad
de ecuaciones que el modelo VOF. No obstante, el modelo euleriano logró caracterizar las variables
de interés en 40000 iteraciones, se puede observar que el tiempo gastado sigue siendo ligeramente
mayor al modelo VOF.
13
Figura 7. Test independencia de mallado - Tiempo de simulación
3.3 Casos de estudio
Los experimentos de flujo bifásico descritos en la Tabla 1 se simularon utilizando el modelo euleriano
y el modelo VOF. En la Tabla 4 se muestra los resultados para los casos A, B, C y D donde la variable
analizada es la fracción de vacío. Por otro lado, los casos E y F analizan la velocidad total del gas,
sus resultados están descritos en la Tabla 5. Adicionalmente, estas tablas muestran los datos
experimentales obtenidos por los autores y la incertidumbre de su experimentación. Estos resultados
permitieron identificar que los modelos VOF y euleriano describen correctamente flujos bifásicos con
bajas fracciones de vacío. Esto se puede evidenciar en los resultados de CFD de los Casos A, B, C y
D, los cuales se encuentran dentro de la incertidumbre de la experimentación. Por el contrario, para
flujos con altas fracciones de vacío estos modelos presentaron errores mayores a la incertidumbre.
Tabla 4. Resultados de los casos de estudio A, B, C y D utilizando el modelo euleriano y el modelo VOF
Euleriano VOF
Caso αexperimental Incertidumbre (%) αCFD Error (%) αCFD Error (%)
A 0.058 24.31 0.063 9.24 0.065 12.31
B 0.115 49.77 0.092 20.13 0.101 11.92
C 0.196 46.42 0.165 15.61 0.118 39.86
D 0.262 42.78 0.234 10.47 0.230 12.05
Tabla 5. Resultados de los casos de estudio E y F utilizando el modelo euleriano y el modelo VOF
Euleriano VOF
Caso 𝒖𝑮,𝐞𝐱𝐩𝐞𝐫𝐢𝐦𝐞𝐧𝐭𝐚𝐥 Incertidumbre (%) 𝒖𝑮,𝐂𝐅𝐃 Error (%) 𝒖𝑮,𝐂𝐅𝐃 Error (%)
E 17.038 8.72 12.203 28.38 12.223 28.26
F 26.081 8.45 21.205 18.70 21.360 18.10
14
En la Figura 8 se muestra la predicción de los modelos euleriano y VOF para los casos de estudio A,
B, C y D. A partir de esta gráfica se puede observar que solamente en un caso el modelo VOF presenta
un error mayor al 30%, el cual representa un error significativo para la magnitud de la fracción de
vacío. Con el fin de escoger el modelo más exacto se utilizó la Ecuación (21) para compilar los
errores, obteniendo que el modelo euleriano presenta un error cuadrático medio (13.86%) inferior al
modelo VOF (19.04%) para flujos con baja fracción de vacío. Por el contrario, para flujos con alta
fracción de vacío los dos modelos presentan el mismo error relativo (≈ 23%).
𝐸𝑟𝑟𝑜𝑟 𝑇𝑜𝑡𝑎𝑙 = [1
𝑁∑ (
𝑋𝐶𝐹𝐷𝑖−𝑋𝑒𝑥𝑝𝑒𝑟𝑖𝑚𝑒𝑛𝑡𝑎𝑙𝑖
𝑋𝑒𝑥𝑝𝑒𝑟𝑖𝑚𝑒𝑛𝑡𝑎𝑙𝑖
)𝑁𝑖=1 ] ∗ 100% (21)
Figura 8. Predicción de los casos de estudio a partir de los modelos VOF y euleriano.
Los modelos físicos estudiados en este proyecto difieren en su formulación matemática como se
describió en el Capítulo 2. Esta diferencia causa que la apariencia de la solución sea distinta para
ambos modelos a pesar de obtener una caracterización similar del fluido bifásico. En la Figura 9 se
puede observar como el modelo VOF detalla con mayor resolución las burbujas en el fluido que el
modelo euleriano. Lo anterior es debido a que el modelo VOF resuelve la interfaz a partir de la
ecuación de continuidad respecto a la fracción de volumen (Ecuación 15) para poder discriminar las
variables de las fases ya que sus otras ecuaciones no lo hacen. Por otro lado, el modelo euleriano
especifica el tamaño de las burbujas y resuelve sus ecuaciones de tal modo que las variables describan
correctamente la cantidad de gas que debe haber en cada posición. Por esta razón, el color de las
imágenes del modelo euleriano es uniforme ya que corresponde a una aglomeración de burbujas. En
este orden de ideas, la Figura 9 evidencia que el modelo VOF puede ser utilizado para predecir los
patrones de flujo a diferencia del modelo euleriano.
Recordando que los casos de estudio se encuentran ordenados de manera ascendente según la fracción
de vacío es factible mencionar que los resultados obtenidos tienen sentido físico. Esto se corrobora
mediante la Figura 9. Adicionalmente, se puede observar que los casos A, B, C y D presentan un
patrón de burbuja dispersa y los casos E y F un patrón anular como se predijo en la Figura 2. No
0 0.05 0.1 0.15 0.2 0.25 0.3 0.350
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
exp
C
FD
VOF
Euleriano
0%
-30%
+30%
15
obstante, cabe resaltar que aunque el modelo VOF discrimina las interfaces presenta problemas para
describir la película líquida entre la pared y el flujo de gas. Lo anterior se debe a que el flujo anular
presenta burbujas arrastradas por la corriente gaseosa y una interface de líquido que tienen un
comportamiento muy dinámico. Por este motivo, se requiere de modelos físicos complementarios que
permitan resolver estas características en CFD (Tkaczyk, 2011; Van Der Meulen, 2012).
Figura 9. Fracción de vacío de los casos de estudio a partir de los modelos VOF y euleriano (1.74 m de tubería mostrada).
3.4 Divergencia del acople del modelo euleriano y el modelo de Drift-Flux
El desempeño del modelo euleriano podría aumentar en su exactitud si se calculara el coeficiente de
arrastre a partir de la Ecuación (10) usando un término que describa correctamente la velocidad de
deriva. En este proyecto se estudió los modelos de Drift-flux para obtener una expresión que describa
este parámetro debido a que la literatura muestra una gran precisión de estas correlaciones para
predecir fracciones de vacío. No obstante, las simulaciones realizadas presentaron problemas de
divergencia usando todos los modelos de Drift-flux reportados en este documento. El modelo
empírico reportado por Sonnenburg (1989) mostró el mejor comportamiento entre los modelos Drift-
Flux; sin embargo no alcanzó un número significativo de iteraciones. En la Figura 10 se muestra los
residuales obtenidos con esta correlación. A partir de esta gráfica se puede observar como la
simulación mantiene constante los residuales hasta alcanzar un comportamiento inestable en la
iteración 65.
16
Figura 10. Residuales de la simulación del modelo euleriano acoplada con el modelo de Sonnenburg (1989).
Los modelos de Drift-flux resultan ser expresiones matemáticas con coeficientes que ajustan una serie
de datos experimentales. Con el fin de lograr una correcta predicción de las fracciones de vacío estos
modelos cuentan con dos términos como se explicó anteriormente, donde el segundo término
corresponde únicamente a la velocidad de deriva. Por esta razón, estos modelos empíricos permiten
predecir de manera correcta la fracción de vacío al utilizar el modelo completo, sin embargo pueden
presentar inconvenientes al describir solamente un término. Lo descrito anteriormente probablemente
describe el problema presentado por las simulaciones.
4. CONCLUSIONES
En este proyecto, se estudiaron los modelos físicos euleriano y VOF usados en la herramienta
computacional CFD para caracterizar el flujo bifásico. Los resultados obtenidos permitieron validar
los modelos al comparar sus respuestas con datos experimentales. Adicionalmente, se logró observar
diferencias relevantes entre los modelos estudiados. Este trabajo evidenció que el modelo euleriano
presenta mayor exactitud que el modelo VOF por una diferencia de 5% para la fracción de vacío en
flujos en patrones de burbuja dispersa. Además, estos modelos presentaron resultados dentro del
rango de incertidumbre de la experimentación demostrando un sentido físico. Sin embargo, estos
modelos presentan inconvenientes modelando correctamente el flujo anular. Por otro lado, se
demostró que el modelo euleriano es independiente de la rigurosidad de la malla a diferencia del
modelo VOF. Por este motivo, la malla seleccionada para el modelo euleriano es más gruesa que la
utilizada por el modelo VOF, requiriendo menor tiempo computacional. Sin embargo, una simulación
con la misma cantidad de celdas necesita mayor tiempo computacional al utilizar el modelo euleriano
que el modelo VOF, debido a que este último presenta menor número de ecuaciones en su
formulación matemática. Otro aspecto importante a considerar es la apariencia de la solución. El
modelo VOF resuelve la interfaz permitiendo discriminar detalladamente las burbujas en el flujo. Por
el contrario, el modelo euleriano especifica el tamaño de las burbujas obligando al sistema a
agruparlas para describir correctamente la cantidad de gas existente en la tubería. Por todo lo anterior,
es más recomendable utilizar el modelo euleriano para la caracterización del flujo bifásico. No
obstante, si se requiere predecir el patrón del flujo del sistema se debe usar el modelo VOF.
El desarrollo del proyecto permitió seleccionar los mejores modelos de Drift-Flux en la literatura para
tuberías verticales. Las correlaciones seleccionadas fueron Maier et al. (1997), Sonnenburg (1989),
Mattar et al. (1974) y Bestion (1985) que presentaron errores alrededor del 30%. De esta manera, se
pudo estudiar el desempeño del modelo euleriano utilizando el modelo de Drift-Flux para calcular el
coeficiente de arrastre. Estos modelos ajustan dos términos para predecir correctamente las fracciones
17
de vacío. Sin embargo, acoplar los modelos de Drift-Flux junto al modelo euleriano para aumentar la
exactitud resulta deficiente. Lo anterior se debe a que las simulaciones presentan un comportamiento
inestable llevando el sistema a divergir. La posible razón se debe a que la expresión de velocidad de
deriva es solo un término de todo un modelo que ajusta una serie de datos experimentales. Por esta
razón, se debe encontrar una correlación que describa únicamente la velocidad de deriva.
5. REFERENCIAS
Abdulkadir, M. (2011). Experimental and computational fluid dynamics (CFD) studies of gas-liquid
flow in bends. Tesis doctoral. University of Nottingham, Nottingham, Inglaterra.
Bratland, O. (2010). Pipe flow 2. Multi-phase flow assurance. Chonburi, Tailandia: Dr. Ove Bratland
Flow Assurance Consulting.
Burns, A.D., Frank, T., Hamill, I. & Shi, J. (2004). The favre averaged drag model for turbulent
dispersion in eulerian multi-phase flows. 5th International Conference on Multiphase Flow, 392,
1-17
Choi, J., Pereyra, E., Sarica, C., Park, C. & Kang, J.M. (2012). An efficient Drift-Flux closure
relationship to estimate liquid holdups of gas-liquid two-phase flow in pipes. Energies, 5, 5294-
5306.
Coddington, P. & Macian, R. (2002). A study of the performance of void fraction correlations used
in the context of drift-flux two-phase flow models. Nuclear engineering and design, 215, 199-
216.
Hernandez, V., Abdulkadir, M. & Azzopardi, B.J. (2010). Grid generation issues in the CFD
modeling of two-phase flow in a pipe. Journal of Computational Multiphase Flows, 3(1), 13-26.
Hewitt, G.F., Delhaye, J.M. & Zuber, N. (1986). Multiphase science and technology (Volume 2).
Springer-Verlag, Alemania: Berlin.
Horgue, P., Augier, F., Quintard., M. & Prat., M. (2012). A suitable parametrization to simulate slug
flows with the Volume-Of-Fluid method. Comptes Rendus Mécanique, 340(6), 411-419.
Ishii, M. & Hibiki, T. (2011). Thermo-fluid dynamics of two-phase flow (2da Ed.). West Lafayette,
EE.UU.: Springer.
Jeyachandra, B.C., Gokcal, B., Al-Sarkhi, A., Sarica, C. & Sharma, A.K. (2012). Drift-velocity
closure relationships for slug two-phase high-viscosity oil flow in pipes. SPE Journal, June, 593-
601.
Krepper, E., Lucas, D. & Prasser, H. (2005). On the modeling of bubbly flow in vertical pipes.
Nuclear engineering and design, 235, 597-611.
Krishna, R., Urseanu, M.I., van Baten, J.M. & Ellenberger, J. (1999) Influence of scale on the
hydrodynamics of bubble columns operating in the churn-turbulent regime: experiments vs.
Eulerian simulations. Chemical Engineering Science, 54, 4903-4911.
Pineda, H.A., López, F.A., Oliveira, B.D., Ratkovich, N. & Carvalho, R.D.M. (2014). Comparison
between LedaFlow® liquid holdup results and experimental data for two-phase flows in
pipelines. Artículo no publicado. Universidad de los Andes, Colombia. Universidade Federal de
Itajubá, Brazil.
Shang, Z. (2014). A novel drag force coefficient model based on drift flux model for gas-water two-
phase flow in vertical pipes under different flow patterns. Artículo no publicado. Institute of
High Performance Computing, Singapur.
Sharaf, S. & Luna-Ortiz, E. (2014). Comparison between two-phase models and wire mesh sensor
measurements in medium and large diameter pipes. 14th AIChE Spring Meeting. New Orleans,
EE.UU.: Xodus group.
STAR-CCM+. (2014). Documentation. CD-adapco.
Sun, X., Paranjape, S., Kim, S., Ozar, B. & Ishii, M. (2004). Liquid velocity in upward and downward
air-water flows. Annals of Nuclear Energy, 31, 357-373.
Thome, J.R. (2004). Engineering Data Book III. Lausana, Suiza: Wolverine Tube, Inc.
18
Tkaczyk, P. (2011). CFD simulation of annular flows through bends. Tesis doctoral. University of
Nottingham, Nottingham, Inglaterra.
Ujang, P.M., Pan, L., Manfield, P.D., Lawrence, C.J. & Hewitt, G.F. (2008). Prediction of the
translational velocity of liquid slugs in gas-liquid slug flow using computational fluid dynamics.
Multiphase Science and Technology, 20(1), 25-79.
Van Der Meulen, G.P. (2012). Churn-annular gas-liquid flows in large diameter vertical pipes. Tesis
doctoral. University of Nottingham, Nottingham, Inglaterra.
Westende, J.M.C. (2008). Droplets in annular-dispersed gas-liquid pipe-flows. Tesis doctoral. Delft
University of technology, Netherlands.
Zhang, H.Q., Wang, Q., Sarica, C. & Brill, J. (2003). Unified model for gas-liquid pipe flow via slug
dynamics-Part 1: Model development. Journal of Energy Resources Technology, 125, 266-273.
6. ANEXOS
6.1. Correlaciones de Drift-Flux
Tabla A1. Correlaciones Drift-Flux y error con respecto a datos experimentales.
Autor Ecuación Error [%]
Bestion
(1985)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶𝑜 = 1; 𝑢𝐷 = 0.188√𝑔𝐷(𝜌𝐿 − 𝜌𝑔)
𝜌𝑔 33.56
Bonnecaze et al.
(1971) 𝛼 = 𝑢𝑔 (1.2𝑢𝑀 + 0.35√𝑔𝐷(1 − 𝜌𝐺 𝜌𝐿⁄ ))⁄ 49.95
Choi et al.
(2012)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶0 = 2 (1 + (𝑅𝑒 1000⁄ )2 + (1.2 − 0.2√𝜌𝐺 𝜌𝐿⁄ (1 − exp(−18𝛼))) (1 + (1000 𝑅𝑒⁄ )2)⁄ )⁄
𝑅𝑒 = 𝜌𝐿𝑢𝑀𝐷 𝜇𝐿⁄ ; 𝑈𝐷 = 0.0246 cos(𝜃) + 1.606(𝑔𝜎(𝜌𝐿 − 𝜌𝐺) 𝜌𝐿⁄ )1/4sin (𝜃) 45.49
Dix
(1971)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶0 =𝑢𝑔
𝑢𝑀(1 + (
𝑢𝑀
𝑢𝐺− 1)
(𝜌𝐺/𝜌𝐿)0.1
)
𝑢𝐷 = 2.9 (𝜌𝜎(𝜌𝐿 − 𝜌𝐺)
𝜌𝐿2 )
0.25
74.51
Filimonov et al.
(1957)
𝛼 = 𝑢𝑔 (𝑢𝑀 + 𝑢𝐷)⁄
𝑢𝐷 = (0.65 − 0.038𝑃)(𝐷 0.063⁄ )0.25 𝑝𝑎𝑟𝑎 𝑃 < 12.7
𝑢𝐷 = (0.33 − 0.00133𝑃)(𝐷 0.063⁄ )0.25 𝑝𝑎𝑟𝑎 𝑃 ≥ 12.7 350
Gomez et al.
(2000)
𝛼 = 𝑢𝑔 (1.5𝑢𝑀 + 𝑢𝐷)⁄
𝑢𝐷 = 1.53[(𝑔𝜎 (𝜌𝐿 − 𝜌𝐺) 𝜌𝐿2⁄ )0.25(sin (𝜃)√1 − 𝛼)] 50.32
Hibiki et al.
(2002)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶𝑜 = 1.2 − 0.2√𝜌𝐺 𝜌𝐿⁄ (1 − exp (−18𝛼));
𝑢𝐷 = [4𝑔𝜎 (𝜌𝐿 − 𝜌𝐺) 𝜌𝐿2⁄ ]0.25(1 − 𝛼)1.75
54.87
Ishii
(1977)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶0 = 1.2 − 0.2√𝜌𝐺 𝜌𝐿⁄ (1 − exp (−18𝛼))
𝑢𝐷 = (𝐶0 − 1)𝑢𝑀 + √2(𝑔𝜎(𝜌𝐿 − 𝜌𝐺)2)0.25 38.90
Kataoka et al.
(1987)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶𝑜 = 1.2 − 0.2√𝜌𝐺 𝜌𝐿⁄ ; 𝑁𝜇𝐿 = 𝜇𝐿 (𝜌𝐿𝜎√𝜎 (𝑔∆𝜌)⁄ )0.5
⁄ ; 𝐷𝐻 = 𝐷 √𝜎 (𝑔∆𝜌)⁄⁄
𝑆í 𝑁𝜇𝐿 > 2.25 × 10−3 → 𝑢𝐷 = 0.92 (𝜌𝐺 𝜌𝐿⁄ )−0.157𝑝𝑎𝑟𝑎 𝐷𝐻 ≥ 30
𝑠𝑖𝑛𝑜 𝑢𝐷 = 0.0019𝐷𝐻0.809(𝜌𝐺 𝜌𝐿⁄ )−0.157𝑁𝜇𝐿
−0.562 𝑝𝑎𝑟𝑎 𝐷𝐻 < 30
𝑢𝐷 = 0.030𝐷𝐻0.809(𝜌𝐺 𝜌𝐿⁄ )−0.157𝑁𝜇𝐿
−0.562 𝑝𝑎𝑟𝑎 𝐷𝐻 ≥ 30
94.16
Kokal et al.
(1989) 𝛼 = 𝑢𝑔 (1.2𝑢𝑀 + 0.345√𝑔𝐷(𝜌𝐿 − 𝜌𝐺) 𝜌𝐿⁄ )⁄ 50.03
Liao et al.
(1985)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶𝑜 = 1.2 − 0.2√𝜌𝐺 𝜌𝐿⁄ (1 − exp(−18𝛼))
𝑢𝐷 = 0.33(𝑔𝜎∆𝜌 𝜌𝐿2⁄ )0.25
54.81
Maier et al.
(1997)
𝛼 = 𝑢𝑔 [(𝐶𝑀𝐶1𝑃 + 𝐶𝑀𝐶2)⁄ 𝑢𝑀 + {(𝑣1𝑃2 + 𝑣2𝑃 + 𝑣3)𝐺 + (𝑣4𝑃2 + 𝑣5𝑃 + 𝑣6)}]
𝐶𝑀𝐶1 = 0.00257, 𝐶𝑀𝐶2 = 1.0062, 𝑣1 = 6.73 × 10−7, 𝑣2 = −8.81 × 10−5, 𝑣3 = 0.00105, 𝑣4 = 0.00563, 𝑣5 = 0.123, 𝑣6 = 0.8, 𝐺 = 𝑈𝑆𝐿𝜌𝐿
27.07
Mattar et al.
(1974) 𝛼 = 𝑢𝑔 (1.3𝑢𝑀 + 0.7)⁄ 33.37
19
Morooka et al.
(1989) 𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶𝑜 = 1.08; 𝑢𝐷 = 0.45 45.81
Nicklin et al.
(1962) 𝛼 = 𝑢𝑔 (1.2𝑢𝑀 + 0.35√𝑔𝐷)⁄ 49.94
Shi et al.
(2005)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶0 =𝐴
1 + (𝐴 − 1)𝛾2 ; 𝛾 =𝛽 − 𝐵
1 − 𝐵; 𝛽 = max (𝛼𝑔, 𝐹𝑣
𝛼𝑢𝑀
𝑢𝑔𝑓) ; 𝑢𝑐 = (
𝜎𝑔𝑙𝑔(𝜌𝑙 − 𝜌𝑔)
𝜌𝑙2 )
0.25
𝑢𝑔𝑓 = 𝐾𝑢 (𝜌𝑙
𝜌𝑔)
0.5
𝑢𝑐; 𝑢𝐷 =((1 − 𝛼𝐶0)𝐶0𝐾(𝛼)𝑢𝑐)
𝛼𝐶0 (𝜌𝑔
𝜌𝑙)
0.5
+ 1 − 𝛼𝐶0
𝑆í 𝛼 ≤ 𝑎1 𝐾(𝛼) = 1.53/𝐶0 𝑠𝑖𝑛𝑜 𝛼 ≥ 𝑎2 𝐾(𝛼) = 𝐾𝑢(�̂�) 𝑜 𝐼𝑛𝑡𝑒𝑟𝑝𝑜𝑙𝑎𝑐𝑖ó𝑛 𝑙𝑖𝑛𝑒𝑎𝑙
𝐴 = 1.2; 𝐵 = 0.3; 𝐹𝑣 = 1.0; 𝑎1 = 0.2; 𝑎2 = 0.4
53.24
Sonnenburg
(1989)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶𝑜 = 1 + 0.32 − 0.32√𝜌𝐺
𝜌𝐿;
𝑢𝐷 =𝐶𝑜(1 − 𝐶𝑜𝛼)
(𝐶𝑜𝛼 √𝑔∆𝜌 𝜌𝐺⁄⁄ ) + (1 − (𝐶𝑜𝛼 √𝑔∆𝜌 𝜌𝐿⁄⁄ ))
29.25
Woldesmayat
(2007)
𝛼 =𝑢𝑔
𝑢𝑔 (1 + (𝑢𝑙
𝑢𝑔)
(𝜌𝐺𝜌𝐿
)0.1
) + 2.9 [𝑔𝐷𝜎(1 + cos 𝜃)(𝜌𝐿 − 𝜌𝐺)
𝜌𝐿2 ]
0.25
(1.22 + 1.22 sin 𝜃)𝑃𝑎𝑡𝑚
𝑃𝑠𝑦𝑠𝑡𝑒𝑚
75.67
Zuber et al.
(1965)
𝛼 = 𝑢𝑔 (𝐶𝑜(𝑢𝑙 + 𝑢𝑔) + 𝑢𝐷)⁄
𝐶𝑜 = 1.2
𝑢𝐷 = 1.53(𝑔𝜎∆𝜌 𝜌𝐿2⁄ )0.25
45.59
6.2 Predicción de los modelos Drift-Flux
Figura A1. Predicción de la correlación de Sonnenburg (1989) con respecto a datos experimentales.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
exp
co
rr
Sonnenburg
0
-30%
+30%
Error[%]29.25
20
Figura A2. Predicción de la correlación de Mattar (1974) con respecto a datos experimentales.
Figura A3. Predicción de la correlación de Bestion (1985) con respecto a datos experimentales.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
exp
co
rr
Mattar
0
-30%
+30%
Error[%]33.37
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
exp
corr
Bestion
0
-30%
+30%
Error[%]33.56