Numerical simulation of turbulent air flow with dispersed droplets in cooling tower eliminators

Resumen

Se estudia numéricamente el flujo turbulento de aire con gotas dispersas de agua a través de separadores de gotas en torres de refrigeración de tiro mecánico. Se resuelven numéricamente las ecuaciones de conservación promediadas de Navier-Stokes mediante un método de discretización de volúmenes finitos, utilizando los códigos Fluent y Phoenics, y empleando alternativamente los modelos de turbulencia k  − ϵ , k  − ω y de esfuerzos de Reynolds, con refinamientos «low-Re» y tratamiento mejorado del flujo junto a las paredes. Se comparan y validan los resultados obtenidos con otros numéricos y experimentales de la bibliografía. Se estudia el grado de aproximación obtenido con cada uno de los modelos de turbulencia considerados. Se investiga la influencia de considerar o no la simulación de la dispersión turbulenta de las gotas, así como los efectos de otros parámetros relevantes, sobre la eficiencia de captura y sobre el coeficiente de pérdida de carga. Centrando la atención sobre 4 tipos concretos de separador («Belgian wave», «H1-V», «L-shaped» y «Zig-zag»), se consideran rangos de velocidad entrante del flujo 1 ≤ U_e  ≤ 5 m/s, de diámetro de gota 2 ≤ D_p  ≤ 50 μm, de números de Reynolds 650 ≤ Re ≤ 8.500 y del parámetro inercial 0, 05 ≤ P_i  ≤ 5. Se comparan además los resultados numéricos obtenidos con Fluent y con Phoenics. Los mejores resultados obtenidos corresponden a las simulaciones llevadas a cabo con Fluent, empleando el modelo de turbulencia SST k  − ω , con valores de la distancia adimensionalizada a la pared y⁺ en el rango 0.2 − 0.5. Por último, se presentan correlaciones para predecir las condiciones de eficiencia de captura máxima (100 %), en función del parámetro geométrico de eficiencia de cada uno de los separadores, introducido en este trabajo.

Abstract

This work presents a numerical study on the turbulent flow of air with dispersed water droplets in separators of mechanical cooling towers. The averaged Navier-Stokes equations are discretised through a finite volume method, using the Fluent and Phoenics codes, and alternatively employing the turbulence models k  − ϵ , k  − ω and the Reynolds stress model, with low-Re version and wall enhanced treatment refinements. The results obtained are compared with numerical and experimental results taken from the literature. The degree of accuracy obtained with each of the considered models of turbulence is stated. The influence of considering whether or not the simulation of the turbulent dispersion of droplets is analyzed, as well as the effects of other relevant parameters on the collection efficiency and the coefficient of pressure drop. Focusing on four specific eliminators (‘Belgian wave’, ‘H1-V’, ‘L-shaped’ and ‘Zig-zag’), the following ranges of parameters are outlined: 1 ≤ U_e  ≤ 5 m/s for the entrance velocity, 2 ≤ D_p  ≤ 50 μm for the droplet diameter, 650 ≤ Re ≤ 8.500 for Reynolds number, and 0.05 ≤ P_i  ≤ 5 for the inertial parameter. Results reached alternately with Fluent and Phoenics codes are compared. The best results correspond to the simulations performed with Fluent, using the SST k  − ω turbulence model, with values of the dimensionless scaled distance to wall y⁺ in the range 0.2 to 0.5. Finally, correlations are presented to predict the conditions for maximum collection efficiency (100 %), depending on the geometric parameter of removal efficiency of each of the separators, which is introduced in this work.

Palabras clave

Eliminadores de gotas ; Eficiencia de captura ; Modelo de turbulencia ; CFD

Keywords

Drift eliminator ; Collection efficiency ; Turbulence model ; CFD

1. Introducción

En el interior de una torre de refrigeración de tiro forzado, se consigue el enfriamiento de un caudal de agua pulverizada a través del intercambio de masa y energía con una corriente de aire no saturado a contraflujo (Kloppers y Kröger [1] , Gan y Riffat [2] , y Gan et al. [3] ). Se evapora una cierta cantidad del agua pulverizada, mientras que el resto se enfría como consecuencia del calor sensible transferido por convección al aire y por el calor latente de cambio de fase. Una pequeña cantidad de gotas puede verse arrastrada hacia el exterior por la corriente de aire. Este arrastre puede resultar perjudicial por un lado porque puede generar corrosión y fallo de los equipos eléctricos, y por el otro porque da lugar a un aerosol que puede desarrollar la bacteria de la legionela. Existen otros efectos también importantes como la pérdida de eficiencia del sistema, o la degradación de los cultivos del entorno por la deposición de sales. Se han publicado distintos trabajos sobre la materia (véanse los trabajos numéricos de referencia de Meroney [4]  and [5] ).

Con objeto de reducir la cantidad de agua arrastrada se instalan en la sección de salida del aire, en el interior de la torre de refrigeración, unas estructuras estáticas que deflectan la corriente (separadores o eliminadores de gotas), y que actúan por impacto inercial. El comportamiento de los separadores puede cuantificarse mediante 2 factores. El primero de ellos es la eficiencia de captura de gotas, definida como la relación entre la masa de agua capturada y la masa total que entra en el separador en el seno de la corriente de aire. El segundo factor es la pérdida de presión producida. La existencia de una caída importante de carga en el interior del cuerpo principal de la torre puede bajar de forma drástica el gasto de aire y reducir por tanto la eficiencia térmica del sistema. En las torres de tiro mecánico puede producirse un incremento considerable en el consumo energético de los ventiladores para un mismo gasto de aire. Obviamente, el mejor diseño debería aunar por un lado una eficiencia de captura alta, y por el otro, una pérdida de carga inducida no demasiado importante.

Caben destacar varias contribuciones importantes al estudio del comportamiento de los separadores. Foster et al. [6] llevaron a cabo estudios analíticos y experimentales sobre 2 tipos de separadores de gotas: de doble fila de lamas de madera, y de perfiles senoidales de fibrocemento. Gardner y Lowe [7] incluyeron en sus estudios los 2 mismos tipos descritos. Chan y Golay [8] desarrollaron un modelo numérico para calcular la eficiencia de captura y la pérdida de presión de diferentes separadores. Becker y Burdik [9] han investigado experimentalmente la pérdida de presión para 2 tipos comerciales de separadores de gotas de geometría tridimensional. Más recientemente, la influencia de considerar canales de drenaje en eliminadores horizontales de humedad en forma de zig-zag ha sido estudiada por James et al. [10] , y Galletti et al. [11] . James et al. [12] han incluido en su trabajo estudios sobre la separación de la capa líquida pegada a las paredes para eliminadores también en forma de zig-zag.

En cuanto a la simulación numérica del flujo planteado, pueden encontrarse en la bibliografía en la mayor parte de los casos estudios con planteamiento euleriano-lagrangiano, que han proporcionado en general buenos resultados (Wang y James [13]  and [14] , Galletti et al. [11] ). Un enfoque alternativo consiste en el estudio del comportamiento detallado de una gota utilizando el concepto «Volume of Fluid, VOF», como han hecho Reis et al. [15] o Strotos et al. [16] , por ejemplo. En cualquier caso, la descripción lagrangiana, o siguiendo la trayectoria de la partícula fluida, conduce a predecir el comportamiento individual de cada gota, lo que resulta adecuado si lo que se persigue es hallar magnitudes tales como la eficiencia de captura. La interacción entre las partículas y las paredes sólidas es un fenómeno complejo (Mundo et al. [17] ). Esta interacción se puede caracterizar en función de la energía de impacto de las gotas (dependiente del número de Weber). No obstante, en este trabajo se supondrá que la energía de impacto no es significativa, de modo que las gotas se quedan pegadas a las paredes en el momento en que entran en contacto con ellas.

Con respecto a la simulación numérica de la turbulencia del flujo considerado, pueden encontrarse distintos estudios que no incluyen la dispersión turbulenta de gotas ejercida por el flujo medio. Wang y James [13]  and [14] , utilizando el código CFX con el modelo de turbulencia k  − ϵ , han encontrado que resulta fundamental, para obtener resultados ajustados a la experimentación, incluir la dispersión turbulenta de las partículas. Han utilizado una mejora del modelo de interacción «Eddy Interaction Model, EIM», llamado por ellos «varied EIM», basado en las sugerencias de Sommerfeld et al. [18] ) y Kallio y Reeks [19] . El mismo enfoque ha sido adoptado por Galletti et al. [11] , extendiendo sus resultados a simulaciones con el modelo SST k  − ω . Puede decirse que códigos CFD tales como Fluent y CFX son quizá los más utilizados en la simulación numérica planteada, siendo muy escasos los presentados por ejemplo con Phoenics (Verlaan [20] , aunque no incluyó la dispersión turbulenta de las gotas en sus simulaciones). Recientemente, Rafee et al. [21] ha presentando un estudio comparativo empleando el modelo de esfuerzos de Reynolds («Reynolds Stress Model, RSM»).

Sriveerakul et al. [22] han presentado un estudio comparativo de 3 tipos de separador, incidiendo en la conveniencia de expresar la eficiencia de captura en función de un cierto parámetro inercial, que será presentado después. En este trabajo se seguirá este enfoque, de modo que como se mostrará más adelante, se centrará la atención en 4 tipos básicos de separador, comparando sus comportamientos a través del cálculo de parámetros de funcionamiento relevantes.

Este trabajo trata del análisis de la influencia que tienen las distintas posibilidades de simulación numérica sobre los resultados obtenidos para la eficiencia de captura y para el coeficiente de pérdida de carga, en separadores de gotas en torres de refrigeración. Se comparán los resultados obtenidos con Fluent y con Phoenics. Se valorará la importancia de incluir o no la simulación de la dispersión turbulenta de las partículas, y se discutirá sobre la bondad de los resultados obtenidos con distintos modelos de turbulencia, comparando los resultados con otros experimentales y numéricos tomados de la bibliografía. Se realiza a continuación la descripción del problema físico a estudiar, para presentar después el modelo matemático y las condiciones de contornos empleadas; a continuación, se desarrollan algunos detalles relacionados con la simulación numérica de los casos. Después de la discusión de resultados (que incluye como aspecto novedoso la presentación de correlaciones para predecir la condición de eficiencia plena de los separadores), se concluye el trabajo relacionando las conclusiones más relevantes alcanzadas.

2. Definición del problema

2.1. Descripción geométrica

Con el fin de validar los resultados obtenidos, se han desarrollado simulaciones numéricas sobre la misma geometría considerada por Wang y James [13]  and [14] , y Phillips y Deakin [23] , consistente en un eliminador de humedad horizontal en forma de zig-zag, de longitud total 0,353 m y separación entre placas de 0,008 m, formado por 7 cambios de dirección principales de 66,2^o cada uno. También se considera un separador similar, con estructuras internas en forma de gancho, sobre el que Galletti et al. [11] ha obtenido resultados numéricos y Ghetti [24] datos experimentales (longitud total de 0,2925 m, anchura 0,0025 m, y 4 cambios de dirección principales de 45^o ).

Se consideran además 4 tipos de separador comercial, con morfología similar y siguiendo la idea expresada por Sriveerakul et al. [22] ,

• «Belgian wave» (que describe una forma senoidal).
• «H1-V» (también llamado de 3 segmentos).
• «L-shaped» (en forma de L inclinada).
• «Zig-zag».

Para comparar el comportamiento de cada uno de los separadores anteriores, se fijan las dimensiones principales de longitud L y de separación de placas b (véase la figura 1 ). Las dimensiones de referencia son b  = 0, 0125 m y L  = 0, 076 m, y por tanto la relación de aspecto de referencia es b /L  = 0, 1645 ≅ 1/6. Para valorar la influencia de la relación de aspecto, se consideran además casos con b /L  = 0, 00625/0, 076 ≅ 1/12 y b /L  = 0, 0250/0, 076 ≅ 1/3. A partir de ahora, los valores considerados de b /L se denominarán simplemente como 1/12, 1/6 y 1/3. Los rangos de parámetros considerados son para el diámetro de gota 2 × 10⁻⁶  ≤ D_p  ≤ 50 × 10⁻⁶ m y para la velocidad de entrada del flujo de aire 1 ≤ U_e  ≤ 5 m/s, lo que lleva a obtener un rango de números de Reynolds Re = ρU_eb /μ  = 650 − 8.500 aproximadamente, y un rango del parámetro de inercia P_i , que se define a continuación, P_i  = 0, 05 − 5.

Figura 1. Configuración geométrica de los 4 separadores estudiados: «Belgian wave», «H1-V», «L-shaped» y «Zig-zag» (fila superior, de izquierda a derecha). Detalle del mallado utilizado y dimensiones principales del canal de un eliminador (fila inferior).

2.2. Parámetros relevantes

El coeficiente de pérdida de presión a través del separador se define por

$K_{P}=-{\frac {\Delta P}{\rho U_{e}^{2}/2}},$

( 1)

siendo ΔP el salto total de presión. A partir de las ecuaciones de conservación y de las condiciones de contorno, la aplicación del análisis dimensional predice que

$K_{P}=f\left(Re,{\frac {b}{L}},forma\right),$

( 2)

donde la forma representa cualquier parámetro geométrico necesario para definir completamente la morfología considerada.

La eficiencia de captura se define por

$\eta =\left({\frac {flujo\quad m{\acute {a}}sico\quad de\quad agua\quad atrapada}{flujo\quad m{\acute {a}}sico\quad total\quad de\quad agua\quad entrante}}\right),$

( 3)

que puede calcularse a su vez como la relación entre el número de gotas que quedan atrapadas en las paredes y el número total de gotas inyectadas en la entrada, suponiendo que los efectos de evaporación son despreciables (más adelante se discutirá esto último). Como antes, analizando las ecuaciones de conservación y las condiciones de contorno, resulta apropiado poner por análisis dimensional que

$\eta =f\left(St,{\frac {L_{w}}{L}},{\frac {b}{L}},{\Theta }_{i},forma\right),$

( 4)

siendo St el número de Stokes; una forma típica de este número adimensional es ${\textstyle St={\rho }_{p}D_{p}^{2}U_{e}/(18\mu L)}$ . Sriveerakul et al. [22] ha denominado a la forma particular

$P_{i}={\frac {{\rho }_{p}D_{p}^{2}U_{e}}{\mu L}},$

( 5)

como el parámetro inercial   , que será el utilizado en el presente trabajo. Las relaciones ${\textstyle L_{w}/L}$ , b /L y los cambios de dirección de la línea media del separador Θ_i caracterizan formológicamente cada eliminador; ${\textstyle L_{w}}$ es la longitud total de las distintos tramos que forman una de las placas del canal (lógicamente, ${\textstyle L_{w}>L}$ ; ${\textstyle L_{w}}$ podría denominarse longitud mojada   ). Cuando los separadores considerados tienen una morfología parecida, que pueda caracterizarse fundamentalmente por ${\textstyle L_{w}/L}$ y Θ_i , es posible eliminar la «forma» de la dependencia funcional anterior, como se verá posteriormente.

3. Modelización matemática

3.1. Modelización del flujo turbulento

El flujo debe considerarse en general turbulento, de modo que las ecuaciones de conservación para el problema considerado son las promediadas de Navier-Stokes, particularizadas para un caso bidimensional y estacionario, reteniendo las fuerzas de flotación (que podrían llegar a ser relevantes para muy bajas velocidades del flujo forzado de aire), y despreciando la potencia de deformación de los esfuerzos viscosos y turbulentos en la ecuación de la energía,

${\frac {\partial (\rho U_{j})}{\partial x_{j}}}=0,$

( 6)

${\frac {\partial (\rho U_{i}U_{j})}{\partial x_{j}}}=-{\frac {\partial P}{\partial x_{i}}}+$$({\rho }_{\infty }-\rho )g_{i}$ $+{\frac {\partial }{\partial x_{j}}}\left(\mu \left[{\frac {\partial U_{i}}{\partial x_{j}}}+\right.\right.$$\left.\left.{\frac {\partial U_{j}}{\partial x_{i}}}-{\frac {2}{3}}{\delta }_{ij}{\frac {\partial U_{j}}{\partial x_{j}}}\right]-\right.$$\left.\rho {\overline {u_{i}u_{j}}}\right),$

( 7)

${\frac {\partial (\rho c_{p}TU_{j})}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\left(\kappa {\frac {\partial T}{\partial x_{j}}}-\right.$$\left.\rho c_{p}{\overline {T^{'}u_{j}}}\right),$

( 8)

donde U_j , T y P son respectivamente las componentes promediadas de velocidad, la temperatura promediada y la diferencia promediada entre la presión y la presión de un cierto estado de referencia (P_∞ ) en equilibrio hidrostático. Además, ${\textstyle -\rho {\overline {u_{i}u_{j}}}}$ es el tensor de esfuerzos turbulentos debido a las fluctuaciones de la velocidad u_j , y ${\textstyle -\rho c_{p}{\overline {T^{'}u_{j}}}}$ es el vector de transporte turbulento de calor debido a las fluctuaciones de T ′ y u_j . Estos 2 últimos términos vienen dados por el modelo de turbulencia empleado. Para los modelos de 2 ecuaciones de transporte,

$-{\overline {u_{i}u_{j}}}=2{\nu }_{t}S_{ij}-{\frac {2}{3}}k{\delta }_{ij},\quad \quad \quad \quad \quad \quad -$${\overline {T^{'}u_{j}}}={\frac {{\nu }_{t}}{{Pr}_{t}}}{\frac {\partial T}{\partial x_{j}}},$

( 9)

siendo ν_t y Pr_t la viscosidad cinemática y el número de Prandtl turbulentos, respectivamente; S_ij es el tensor de esfuerzos medios, S_ij  = [(∂U_i /∂x_j ) + (∂U_j /∂x_i )]/2, δ_ij la función delta de Kronecker, y k   la energía cinética turbulenta, ${\textstyle k={\sum }_{j=1}^{j=3}({\overline {u_{j}^{2}}})/2}$ .

3.1.1. Modelo k  − ϵ

El problema de cierre del sistema de ecuaciones planteado puede solventarse utilizando el modelo de 2 ecuaciones de transporte k  − ϵ (ϵ es el ritmo de disipación de energía cinética turbulenta, ϵ  ≈ dk /dt ), en su forma estándar en los casos en los que la turbulencia del flujo está desarrollada (altos números de Reynolds, con funciones o leyes de pared logarítmicas), o con una forma apropiada para flujo transicional (bajos números de Reynolds, «low-Re»). Las ecuaciones de la forma estándar son

${\frac {\partial k}{\partial t}}+U_{j}{\frac {\partial k}{\partial x_{j}}}=$${\frac {\partial }{\partial x_{j}}}\left[{\frac {{\Gamma }_{k}}{\rho }}{\frac {\partial k}{\partial x_{j}}}\right]+$${\frac {G_{k}}{\rho }}-\rho \epsilon ,$

( 10)

${\frac {\partial \epsilon }{\partial t}}+U_{j}{\frac {\partial \epsilon }{\partial x_{j}}}=$${\frac {\partial }{\partial x_{j}}}\left[{\frac {{\Gamma }_{\epsilon }}{\rho }}{\frac {\partial \epsilon }{\partial x_{j}}}\right]+$$C_{1\epsilon }{\frac {\epsilon }{\rho \quad k}}\quad G_{k}-$$C_{2\epsilon }{\frac {{\epsilon }^{2}}{k}},$

( 11)

donde G_k es la generación de energía cinética turbulenta debida a los gradientes de velocidad media, ${\textstyle G_{k}=-\rho {\overline {u_{i}u_{j}}}(\partial u_{j}/\partial x_{i})}$ , y Γ_k y Γ_ϵ vienen dadas por

${\frac {{\Gamma }_{k}}{\rho }}=\nu +{\frac {{\nu }_{t}}{{\sigma }_{k}}},\quad \quad \quad {\frac {{\Gamma }_{\epsilon }}{\rho }}=$$\nu +{\frac {{\nu }_{t}}{{\sigma }_{\epsilon }}},$

( 12)

siendo C_1ϵ , C_2ϵ , C_μ , σ_k y σ_ϵ las constantes experimentales del modelo, dadas por

$C_{1\epsilon }=1,44,C_{2\epsilon }=1,92,C_{\mu }=0,09,{\sigma }_{k}=$$1,{\sigma }_{\epsilon }=1,3,$

( 13)

y cumpliéndose que ν_t  = C_μ (k² /ϵ ).

En condiciones de turbulencia desarrollada y homogénea, se introducen leyes de pared logarítmicas para los perfiles de velocidad y de temperatura cerca de las paredes, de forma que el flujo turbulento no se resuelve en estas zonas. Cuando el número de Reynolds local no es lo suficientemente elevado, la forma estándar del modelo k  − ϵ se modifica para resolver con mayor precisión el flujo junto a las paredes (forma «low-Re» del modelo, véase por ejemplo Versteeg y Malalasekera [25] ). En concreto, para valores de la coordenada transversal adimensionalizada y⁺  = ρu_τy₁ /μ en el rango 30 < y⁺  < 300, se adopta el primer enfoque; y₁ es la distancia mínima desde el primer punto de cálculo a la pared y u_τ la velocidad característica turbulenta, u_τ  = (τ_p /ρ )^1/2 , siendo τ_p el esfuerzo cortante en la pared. En el caso de adoptar el segundo enfoque, es necesario un mallado muy denso junto a la pared para obtener valores de y⁺ en torno a la unidad, llegando incluso a escalas del orden de la subcapa límite laminar.

3.1.2. Modelo k  − ω

Puede obtenerse una descripción más detallada del comportamiento del flujo turbulento influenciado por paredes a través del modelo k  − ω . Este modelo fue el primero de 2 ecuaciones de transporte turbulentas, siendo propuesto por Kolmogorov [26] . En efecto, presenta sendas ecuaciones para la energía cinética turbulenta k y para su ritmo específico de disipación, ω  ≈ ϵ /k , y se formula a menudo incluyendo una mejora «low-Re», como puede verse por ejemplo en Wilcox [27] . A diferencia del modelo k  − ϵ , en el que las ecuaciones de k y de ϵ se derivan de las de Navier-Stokes, en el modelo k  − ω la ecuación para ω fue propuesta originalmente por Kolmogorov [26] a partir de razonamientos físicos y de análisis dimensional.

Las mismas conclusiones acerca de la influencia de los resultados obtenidos con respecto de y⁺ pueden indicarse para este modelo. Los detalles de la mejora SST («Shear-Stress Transport») del modelo k  − ω , propuesta por Menter [28] , pueden verse en [29] , por ejemplo. En este trabajo, se ha utilizado la forma SST del modelo k  − ω , que lógicamente incluye el tratamiento especial para bajos números de Reynolds, y además parece simular más correctamente los flujos de transición, como ha sido comprobado por los autores en trabajos previos ( [30] y [31] ).

3.1.3. Modelo de esfuerzos de Reynolds

En el modelo de esfuerzos de Reynolds («Reynolds Stress Model, RSM») se abandona la hipótesis de viscosidad turbulenta isotrópica en la que están basados los modelos de 2 ecuaciones de transporte tales como el k  − ϵ o el k  − ω ; el modelo RSM cierra las ecuaciones de Navier-Stokes promediadas resolviendo ecuaciones de transporte para los esfuerzos de Reynolds, junto con una ecuación para el ritmo de disipación turbulenta ϵ . Esto significa que es necesario resolver un número mayor de ecuaciones de transporte para simular el flujo turbulento. Lógicamente, los tiempos de computación son considerablemente mayores que los necesarios para los modelos k  − ϵ o k  − ω , y en ocasiones las dificultades de convergencia numérica son apreciables.

Por sus significativas diferencias con los modelos de 2 ecuaciones de transporte, se describe a continuación este modelo con algo más de detalle. Las ecuaciones exactas de transporte, para los esfuerzos de Reynolds, ${\textstyle R_{ij}={\overline {u_{i}^{'}u_{j}^{'}}}}$ son complejas. Una forma estándar compacta puede encontrarse en Versteeg y Malalasekera [25] ,

${\frac {{DR}_{ij}}{Dt}}=P_{ij}+D_{ij}-{\epsilon }_{ij}+$${\Pi }_{ij}+{\Omega }_{ij},$

( 14)

en la que el término de producción de esfuerzos turbulentos P_ij es

$P_{ij}=-\left(R_{im}\quad {\frac {\partial U_{j}}{\partial x_{m}}}+\right.$$\left.R_{jm}\quad {\frac {\partial U_{i}}{\partial x_{m}}}\right),$

( 15)

siendo necesario modelizar el resto de términos de forma conveniente. Existen distintas aproximaciones, propuestas por distintos autores. Por ejemplo, el término de transporte por difusión puede ponerse por

$D_{ij}={\frac {\partial }{\partial x_{m}}}\quad \left({\frac {{\nu }_{t}}{{\sigma }_{k}}}\quad {\frac {\partial R_{ij}}{\partial x_{m}}}\right),$

( 16)

el ritmo de disipación de R_ij por ϵ_ij  = (2/3)ϵδ_ij y el transporte de R_ij debido a la rotación,

${\Omega }_{ij}=2{\omega }_{k}(R_{jm}e_{ikm}+R_{im}e_{jkm}).$

( 17)

Ahora bien, el término de transporte de R_ij debido a la interacciones turbulentas de presión-tensión, Π_ij , es más complejo de modelizar. Una de las modelizaciones más sencillas es la descrita en [25] ,

${\Pi }_{ij}=-C_{1}\quad {\frac {\epsilon }{k}}\quad \left(R_{ij}-\right.$$\left.{\frac {2}{3}}\quad k{\delta }_{ij}\right)-C_{2}\quad \left(P_{ij}-\right.$$\left.{\frac {2}{3}}P{\delta }_{ij}\right).$

( 18)

Sin embargo, Fluent ofrece varias alternativas, más completas que esa última, habiéndose escogido una modelización lineal («linear pressure-strain model»), que puede verse con detalle por ejemplo en [29] . Se trata del mismo algoritmo utilizado por Rafee et al. [21] en sus simulaciones.

En las expresiones anteriores, ν_t  = C_μk² /ϵ , C_μ  = 0, 09, σ_k  = 1, 0, C₁  = 1, 8 y C₂  = 0, 6. En la ecuación (17) , se ha utilizado el operador e_ijk , que cumple lo siguiente: e_ijk  = 1 si i , j y k son diferentes y están en orden cíclico; e_ijk  = − 1 si i , j y k son diferentes, y están en orden no cíclico, y e_ijk  = 0 si 2 de los índices son iguales.

Con objeto de calcular de forma más precisa el flujo en el entorno de las paredes (enfoque «low-Re»), puede utilizarse el algoritmo «enhanced wall treatment», que introduce unas funciones de pared tipo «blended» o compuestas (laminares o lineales + turbulentas o logarítmicas). Se expondrá con algo más de detalle posteriormente.

3.2. Modelización euleriana-lagrangiana del flujo de aire con gotas dispersas

Se distinguirá entre el conjunto de ecuaciones que gobiernan la fase continua (aire seco mezclado con vapor de agua), de las que gobiernan la fase discreta (agua pulverizada). Se adopta pues un enfoque euleriano para la fase continua y otro lagrangiano para la discreta. El modelo que se describe a continuación es el que se utiliza en las simulaciones realizadas con Fluent [29] . Se ha comprobado que Phoenics utiliza básicamente el mismo modelo a través de un subprograma llamado «Gentra» [32] , aunque ciertas hipótesis parecen más simplificadas, según se puede comprobar haciendo una inspección de los manuales del código.

3.2.1. Modelización de la fase continua

Se incorporan al sistema de ecuaciones ((6)  and (7) ,8 ) distintos términos fuente,

${\frac {\partial (\rho U_{j})}{\partial x_{j}}}=S_{i^{'}},$

( 19)

${\frac {\partial (\rho U_{i}U_{j})}{\partial x_{j}}}=-{\frac {\partial P}{\partial x_{i}}}+$$({\rho }_{\infty }-\rho )g_{i}+{\frac {\partial }{\partial x_{j}}}\left(\mu \left[{\frac {\partial U_{i}}{\partial x_{j}}}+\right.\right.$$\left.\left.{\frac {\partial U_{j}}{\partial x_{i}}}-{\frac {2}{3}}{\delta }_{ij}{\frac {\partial U_{j}}{\partial x_{j}}}\right]-\right.$$\left.\rho {\overline {u_{i}u_{j}}}\right)+F_{i},$

( 20)

${\frac {\partial (\rho c_{p}TU_{j})}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\left(\kappa {\frac {\partial T}{\partial x_{j}}}-\right.$$\left.\rho c_{p}{\overline {T^{'}u_{j}}}\right)+{\frac {\partial }{\partial x_{i}}}\left(\sum _{i^{'}=1}^{n}h_{i^{'}}J_{i^{'}}\right)+$$S_{h},$

( 21)

siendo ${\textstyle S_{i^{'}}}$ las fuentes de masa existentes en el elemento diferencial considerado; F_i las fuentes de cantidad de movimiento; ${\textstyle \partial ({\sum }_{i^{'}=1}^{n}h_{i^{'}}J_{i^{'}})/\partial x_{i}}$ el transporte de entalpía debido a la difusión de las especies i ′, y S_h las fuentes de energía.

La ecuación de conservación de las especies (Ley de Fick) puede escribirse

$\rho u_{i}{\frac {\partial m_{i^{'}}}{\partial x_{i}}}=$$-{\frac {\partial J_{i^{'},i}}{\partial x_{i}}}+S_{i^{'}},$

( 22)

donde ${\textstyle J_{i^{'},i}}$ representa el flujo difusivo de la especie i ′ y viene dado por

$J_{i^{'},i}=-\rho D_{i^{'},m}{\frac {\partial m_{i^{'}}}{\partial x_{i}}},$

( 23)

siendo ${\textstyle m_{i^{'}}}$ la fracción másica de la especie i   ′, y ${\textstyle D_{i^{'},m}}$ el coeficiente de difusión de i ′.

3.2.2. Modelización de la fase discreta

Las ecuaciones para la fase discreta (gotas de agua líquida de diámetro D_p y velocidad U_p ) planteadas en un sistema de referencia lagrangiano, son las de balance de fuerzas sobre la partícula, de trayectoria r_p y de conservación de la energía en la partícula [29] ,

${\frac {dU_{p}}{dt}}={\frac {18\mu }{{\rho }_{p}D_{p}^{2}}}{\frac {C_{D}{Re}_{p}}{24}}\left(U-\right.$$\left.U_{p}\right)+{\frac {g({\rho }_{p}-\rho )}{{\rho }_{p}}}+$${\frac {\rho }{{\rho }_{p}}}U_{p}{\frac {\partial U}{\partial r_{p}}},$

( 24)

${\frac {dr_{p}}{dt}}=U_{p},$

( 25)

$m_{p}c_{p}{\frac {dT_{p}}{dt}}={hA}_{p}\left(T-T_{p}\right)+{\frac {dm_{p}}{dt}}h_{f},$

( 26)

siendo Re_p  = (ρD_p  ∣ U_p  − U  ∣)/μ el número de Reynolds de la partícula basado en su velocidad relativa, y

$C_{D}=a_{1}+{\frac {a_{2}}{{Re}_{p}}}+{\frac {a_{3}}{{Re}_{p}^{2}}},\quad F_{D}=$${\frac {18\mu }{{\rho }_{p}D_{p}^{2}}}{\frac {C_{D}{Re}_{p}}{24}},$

( 27)

donde los coeficientes a₁ , a₂ y a₃ son constantes utilizadas para el caso de partículas esféricas para diferentes rangos de valores de Re; ${\textstyle F_{D}\quad \left(U-U_{p}\right)}$ es la fuerza de resistencia por unidad de masa de la partícula, g (ρ_p  − ρ )/ρ_p la fuerza de gravedad por unidad de masa, (ρ /ρ_p )U_p (∂U /∂r_p ) la fuerza causada por el gradiente de presiones en el fluido, y dm_p /dt el promedio de masa evaporada en la partícula. Por último, h es el coeficiente de transmisión de calor por convección y h_f el calor latente de vaporización del agua.

3.2.3. Acoplamiento entre las fases continua y discreta

El acoplamiento entre las fases discreta y continua se efectúa a través de términos fuente, de una forma iterativa. Inicialmente se calcula la trayectoria de la partícula sin interaccionar con la fase continua, almacenando el calor, la masa y la cantidad de movimiento ganados o perdidos por la partícula y posteriormente estos valores se incorporan a través de los términos fuente al modelo de la fase continua, repitiéndose el proceso hasta converger. El término fuente de masa viene dado por:

$S_{i}^{'}={\frac {\Delta m_{p}{\dot {m}}_{p_{o}}}{m_{p_{o}}d\forall }},$

( 28)

donde Δm_p es el cambio de masa de la partícula en cada volumen diferencial d∀ en un dt   , ${\textstyle {\dot {m}}_{po}}$ el flujo másico inicial de partículas inyectadas en cada elemento diferencial y m_po la masa inicial de la partícula. El cambio de masa en cada elemento diferencial d∀ es

$\Delta m_{p}\left(d\forall \right)=m_{p}\left(t\right)-$$m_{p}\left(t-dt\right)=N_{v}M_{v}A_{p}dt,$

( 29)

donde dt  = ds /|U_p  + U |, y ds   es la fracción de trayectoria dentro de cada volumen diferencial d∀, ${\textstyle M_{v}}$ el peso molecular del vapor, A_p el área de la gota y ${\textstyle N_{v}}$ el flujo molar de vapor ${\textstyle N_{v}=K_{c}\left(C_{v,s}-C_{v}\right)}$ , con ${\textstyle C_{v,s}}$ la concentración de vapor en la superficie de la gota y ${\textstyle C_{v}}$ la concentración de vapor en el gas. El coeficiente de transferencia de masa K_c se obtiene de la siguiente correlación para el número de Nusselt

$K_{c}={\frac {(Nu)D_{v,m}}{D_{p}}}={\frac {\left[2+0,65({Re}_{p})^{1/2}(S_{c})^{1/3}\right]D_{v,m}}{D_{p}}},$

( 30)

siendo ${\textstyle D_{v,m}}$ el coeficiente de difusión del vapor en la mezcla.

Los términos fuente de cantidad de movimiento F_i y energía S_h vienen dados por:

$F_{i}=\left({\frac {18\mu }{{\rho }_{p}D_{p}^{2}}}{\frac {C_{D}Re}{24}}(U-\right.$$\left.U_{p})+{\frac {g({\rho }_{p}-\rho )}{{\rho }_{p}}}\right){\frac {{\dot {m}}_{p_{o}}dt}{d\forall }}$ $+\left({\frac {\rho }{{\rho }_{p}}}U_{p}{\frac {\partial U}{\partial r_{p}}}\right){\frac {{\dot {m}}_{p_{o}}dt}{d\forall }}$

( 31)

$S_{h}=\left[{\frac {\overline {m_{p}}}{m_{p_{o}}}}\quad c_{p}\quad \Delta T_{p}\right]{\frac {{\dot {m}}_{p_{o}}}{d\forall }}$ $+\left[{\frac {\Delta m_{p}}{m_{p_{o}}}}\left(-h_{f}+{\int }_{T_{ref}}^{T_{p}}c_{p,v}\quad dT\right)\right]{\frac {{\dot {m}}_{p_{o}}}{d\forall }},$

( 32)

donde ${\textstyle {\overline {m_{p}}}}$ es el promediado de masa de la partícula en el diferencial de volumen d∀ y ΔT_p el cambio de temperatura de la partícula en dicho diferencial de volumen.

Se ha preferido describir el modelo completo que tendría en cuenta los efectos de evaporación de las gotas, y por tanto la incorporación del vapor de agua a la corriente de aire. Se ha comprobado para todas las simulaciones que se presentan que los tiempos de residencia de las partículas líquidas en el seno de la corriente son tan bajos que los términos correspondientes a la evaporación del líquido en las ecuaciones de conservación resultan despreciables frente al resto de términos. De esta manera, los resultados numéricos que se presentan en este trabajo no se ven prácticamente afectados por los efectos de evaporación. De hecho, es posible renunciar a la solución de la ecuació de la energía sin menoscabo en la precisión de los resultados.

3.2.4. Dispersión turbulenta de las gotas

En Fluent, los efectos de dispersión turbulenta sobre las gotas se simulan utilizando un modelo de trayectorias aleatorias («Discrete Random Walk model, DRW»). La dispersión turbulenta puede predecirse a través de las ecuaciones de trayectoria para cada partícula individual, utilizando a lo largo de la integración sobre el camino seguido por cada partícula el concepto de velocidad instantánea del fluido, U_j  + u_j . La velocidad de fluctuación turbulenta o aleatoria u_j se supone constante sobre un cierto intervalo de tiempo t_e , que depende de la escala característica de tiempo o tiempo de existencia de los remolinos turbulentos; es decir, cada remolino turbulento se caracteriza por u_j y t_e . Calculando la trayectoria de las gotas a través de un número representativo de partículas («tries»), se puede llegar a simular los efectos de aleatoreidad que sobre la trayectoria de las partículas produce la turbulencia del flujo.

En primer lugar, se discutirá el tiempo de interacción t_i entre la partícula y el flujo medio. Por un lado, la escala de tiempo del flujo Lagrangiano puede aproximarse por t_L  = C_L (k /ϵ ). Es preciso tener en cuenta que C_L no está totalmente definido; puede tomarse igual a 0.15 para modelos tales como el k  − ϵ y similares; posteriormente se estudiará su influencia. El tiempo de existencia de los remolinos turbulentos puede tomarse como t_e  = − t_L  log(r ), siendo r un número uniformemente aleatorio comprendido entre 0 y 1. Por otro lado, el tiempo de cruce de remolino   por parte de la partícula is ${\textstyle t_{cross}=-\tau \quad ln[1-}$$(L_{e}/(\tau \left|U-U_{p}\right|)]$ , con τ un cierto tiempo de relajación para la partícula y L_e la escala de longitud de remolino   , que a su vez puede calcularse a partir del modelo de turbulencia por ${\textstyle L_{e}=C_{\mu }^{3/4}(k^{3/2}/\epsilon )}$ , con C_μ  = 0.09. Se decide entonces que el tiempo de interacción gota-remolino turbulento t_i sea el menor entre el tiempo de existencia de los remolinos t_e , y el tiempo de cruce partícula-remolino t_cross .

En segundo lugar, se discutirá la velocidad de fluctuación turbulenta u_j . Las componentes promediadas de velocidad turbulenta calculadas en modelos de turbulencia tales como k  − ϵ o k  − ω   son iguales en cualquier dirección, ${\textstyle {\overline {u_{j}^{2}}}=(2/3)k}$ , suponiendo lógicamente isotropía del flujo. Así, las componentes turbulentas de velocidad que actúan durante el período t_i pueden ponerse por ${\textstyle u_{j}=\zeta ({\overline {u_{j}^{2}}})^{1/2}}$ , siendo ζ un número distribuido normalmente, que conduce a una distribución Gaussiana de velocidades turbulentas fluctuantes. En el caso del modelo de turbulencia RSM, u_j se calcula directamente a partir de los esfuerzos de Reynolds, con un procedimento similar al que se acaba de describir.

El mismo procedimiento que se ha explicado ha sido empleado por Wang y James [13]  and [14] y Galletti et al. [11] para simular la dispersión turbulenta de las gotas (el «Eddy Interaction Model, EIM»). No obstante, alguna diferencia ha quedado de manifiesto, como puede ser la determinación de t_e (los autores mencionados suponen t_e  = 0.201 [k /ϵ ]). Además, ellos usan una mejora del EIM, el llamado «varied EIM»; el principal refinamiento consiste en que tanto U_j como u_j se recalculan cuando la partícula atraviesa cada volumen de control de la malla. Wang y James [13]  and [14] han utilizado incluso otra corrección del EIM, denominada «GQ EIM», en la que se establece además una cierta mejora del tratamiento de la región próxima a las paredes.

En Phoenics, el modelo de dispersión turbulenta de las partículas de Gentra [32] parte de la misma base expuesta más arriba, pero no parece ser tan completo. Por ejemplo, el tiempo de interacción t_i se elige también como el menor entre el tiempo de existencia medio de los remolinos turbulentos t_e y el tiempo de cruce t_cross , haciéndose directamente ${\textstyle t_{e}=L_{e}/\vert u_{j}^{'}\vert }$ y t_cross  = L_e /|U  − U_p |, de forma que aunque L_e se calcula del mismo modo (${\textstyle L_{e}=C_{\mu }^{3/4}(k^{3/2}/\epsilon )}$ ), no se introduce por ejemplo de forma explícita un cierto coeficiente del tipo de C_L para actuar sobre la escala de tiempo lagrangiano, ni por tanto sobre el tiempo de interacción t_i , como puede efectuarse en Fluent. Posteriormente se verá cómo el valor escogido de C_L es determinante en el cálculo de la eficiencia de captura.

3.3. Condiciones de contorno

– Sección de entrada . En el dominio computacional de la figura 1 , el aire entra por la sección inferior con una cierta velocidad uniforme U_e . Los valores iniciales de k y ϵ (o ω ) en la entrada se imponen a través del concepto de intensidad de turbulencia , definida por I  = [(2/3)k ]^1/2 /U_e . De este modo,

$k_{ent}={\frac {3}{2}}I^{2}U_{e}^{2},\quad \quad \quad {\epsilon }_{ent}=$${\frac {k_{ent}^{2}}{{\nu }_{t}}},\quad \quad \quad {\omega }_{ent}=$${\frac {k_{ent}}{{\nu }_{t}}},$

( 33)

con un valor inicial ν_t  ≈ 40 ν . Con objeto de obtener resultados sistemáticos, todos las simulaciones se han obtenido con una intensidad de turbulencia del 5 %. Se han realizado no obstante algunas simulaciones con 2 ≤ I  ≤ 20 %, llegando a la conclusión de que I no afecta a η ni a K_P de forma significativa. La entrada de gotas dispersas se efectúa mediante un número de inyecciones suficientemente elevado en la sección de entrada del dominio. La mayoría de casos se realizaron con 100 inyecciones en el caso de Phoenics y con 500 inyecciones en el caso de Fluent.

• Sección de salida . Se fija la presión igual a la ambiente y se suponen despreciables las variaciones de las propiedades fluidas en la dirección de la corriente.
• Paredes . Todas las paredes se suponen adiabáticas. La condición de adherencia se impone haciendo nulos los valores de U_j y u_j . Cuando se utiliza la versión «low-Re» del modelo k  − ϵ ,

$k_{walls}=0;\quad \quad \quad {\epsilon }_{walls}=$$\nu \quad {\left({\frac {{\partial }^{2}k}{\partial y^{2}}}\right)}_{y=0},$

( 34)

mientras que para el k  − ω ,

$k_{walls}=0,\quad \quad \quad {\omega }_{walls}={\frac {K\nu }{{\beta }_{\infty }^{\ast }y_{1}^{2}}},$

( 35)

con K  = 2 y ${\textstyle {\beta }_{\infty }^{\ast }=0.09}$ , una constante del modelo de turbulencia.

Con respecto al modelo RSM, se ha impuesto la condición «enhanced wall treatment» en las paredes, en función del valor calculado de y⁺ ,

$u^{+}={\frac {U}{u_{\tau }}}=e^{\Gamma }u_{lam}^{+}+e^{1/\Gamma }u_{turb}^{+},$

( 36)

siendo ${\textstyle u_{lam}^{+}}$ y ${\textstyle u_{tur}^{+}}$ las funciones propuestas para u⁺ en la subcapa límite laminar y en la capa turbulenta logarítmica, respectivamente; el factor de composición es Γ = − [a (y⁺ )⁴ ]/[1 + by⁺ ], con a  = 0, 01 y b  = 5. Aunque los mejores resultados se alcanzan para valores de y⁺  < 1, este tratamiento especial de los contornos sólidos debería proporcionar buenos resultados incluso con valores de y⁺  ≈ 5. Una descripción detallada se puede encontrar en [29] .

Con respecto a la determinación numérica de η , siendo la energía de impacto de las gotas sobre las paredes relativamente pequeña, con ausencia de fenómenos de ebullición, se considera como apropiada para los objetivos del presente trabajo una condición de contorno consistente en la deposición tras el primer contacto de la gota con la pared. Lógicamente, con una condición de contorno de este tipo no es posible simular los fenómenos de formación de película y de separación de la misma sobre las paredes, que merecerían una atención más detallada.

4. Detalles sobre la modelización numérica

Como se ha indicado anteriormente, se emplea un procedimiento de discretización de volúmenes finitos para resolver los problemas planteados, alternativamente mediante los códigos Fluent y Phoenics. La simulación ha sido bidimensional y estacionaria. En la mayoría de los casos, se ha utilizado un procedimiento «upwind» de segundo orden para los términos convectivos de las ecuaciones. Se han realizado simulaciones adicionales empleando de forma alternativa esquemas tales como el «quick»y el «muscl», para asegurar que los efectos de difusión numérica no afectaban a los resultados obtenidos. No se han registrado diferencias significativas al emplear estos esquemas numéricos alternativos en ninguno de los 2 códigos.

Para resolver el acoplamiento entre las ecuaciones de cantidad de movimiento y de continuidad a través de la presión, se utiliza un algoritmo del tipo «Simple», tanto en Fluent como en Phoenics. El criterio de convergencia se ha establecido como ∣(ϕ^{i +1}  − ϕⁱ )/ϕⁱ  ∣ ≤ 10⁻⁴ , donde ϕ denota cualquier variable e i el número de iteración; además, los residuos normalizados de masa, cantidad de movimiento, energía y variables turbulentas han caído en todos los casos a niveles del orden de 10⁻³  − 10⁻⁵ .

Se han utilizado mallados estructurados, afinando la malla junto a las paredes, y en los entornos de las secciones de entrada y de salida en la mayoría de casos. Algunas de las simulaciones planteadas pueden considerarse como transicionales, por lo que se ha hecho especial hincapié en el tratamiento del flujo junto a las paredes. Para los 4 tipos de separadores estudiados, las mallas fueron de 50 − 80 × 400 − 700 elementos en la mayoría de los casos. Con la condición de obtener valores máximos de y⁺ en torno a la unidad junto a las paredes, se ha comprobado tanto en Phoenics como en Fluent que la independencia de los resultados obtenidos con respecto a la densidad de la malla se alcanza para valores de y⁺ en el rango 0,2–0,5. Para valores inferiores de y⁺ , las dificultades de convergencia numérica aumentaban de forma considerable, no obteniéndose una ganancia clara en la precisión de los resultados. Los tiempos de computación han variado en función de la complejidad del separador. En el caso de las morfologías más complejas, incluso fueron necesarias unas 10.000 iteraciones para alcanzar resultados satisfactorios. En el caso de Phoenics, para conseguir la convergencia de los resultados al emplear perentoriamente el sistema de mallado ajustado al cuerpo BFC («Body Fitted Coordinates»), ha sido necesario activar el algoritmo llamado GCV («General Colocated Velocity») [33] para las componentes de velocidad, lo que ha producido en ocasiones tiempos de cálculo algo más elevados.

5. Discusión de resultados

5.1. Validación con resultados numéricos y experimentales

En las figuras Figura 2 , Figura 3 , Figura 4 , Figura 5 , Figura 6  and Figura 7 se muestran las distintas comparaciones que se han realizado entre los resultados numéricos obtenidos con Fluent y los resultados numéricos y experimentales presentados por otros autores.

En las figuras 2 y 3 pueden observarse las validaciones llevadas a cabo con los resultados numéricos de Wang y James [13]  and [14] , que utilizaron un modelo «low-Re» k  − ϵ para un eliminador horizontal de tipo zig-zag (el mismo que el descrito y estudiado experimentalmente por Phillips y Deakin [23] ). La figura 2 , para U_e  = 2 m/s, muestra que existen diferencias significativas entre los resultados obtenidos incluyendo o no la modelización de la dispersión turbulenta de las partículas. La eficiencia de captura resulta mayor, para un mismo valor de D_p , cuando se incluye el modelo EIM de dispersión turbulenta, siendo los resultados más ajustados a los experimentales con la mejora «varied EIM». En cuanto a los resultados numéricos presentados en este trabajo, puede verse que los que más se aproximan a los experimentales de Phillips y Deakin [23] son los correspondientes al modelo SST k  − ω , con dispersión turbulenta.

Figura 2. Comparación con los resultados numéricos para η de Wang y James [13]  and [14] y experimentales de Phillips y Deakin [23] , para un separador horizontal en zig-zag. Ue  = 2 m/s.

Figura 3. Comparación con los resultados numéricos para η de Wang y James [13]  and [14] y experimentales de Phillips y Deakin [23] . Influencia del factor CL . Ue  = 1 m/s.

Como se ha indicado anteriormente, Rafee et al. [21] han encontrado que el modelo de turbulencia RSM, con tratamiento mejorado en las paredes, conducía a obtener resultados muy próximos a los de Wang y James [13]  and [14] y Phillips y Deakin [23] . Sin embargo, la misma tendencia no se ha podido reproducir en este trabajo. En la figura 3 puede comprobarse lo que se acaba de exponer. Se ha estudiado además la influencia del factor C_L en los resultados obtenidos. Aunque los mejores resultados deberían obtenerse con C_L ≈ 0,3 (según la recomendación leída en [29] ), probablemente la utilización de la opción «enhanced wall treatment» lleva a obtener valores excesivamente altos de η , de modo que se precisa de un valor bajo de C_L para conseguir resultados de η cercanos a los experimentales. Para bajos valores de C_L , la tendencia es a obtener valores de η bajos para valores de D_p altos, de modo que el valor predicho para el D_p que consigue la eficiencia plena resulta excesivamente alto. Para valores de C_L altos, la tendencia es inversa. Los mejores resultados se obtuvieron con C_L  = 0,075, la mitad del empleado por Rafee et al.[21] . Imponiendo C_L  = 0,15, como Rafee et al. [21] , los resultados obtenidos resultaron considerablemente más altos que los de Rafee et al. [21] , aunque la condición de eficiencia plena no varió demasiado. El mismo comportamiento se observa en la figura 4 , en esta ocasión para U_e  = 2 m/s. En esta misma figura se ve que la influencia de C_L sobre η parece ser mayor cuando se emplea el modelo k  − ω en lugar del RSM. Los mejores resultados se siguen produciendo con el modelo k  − ω y C_L  = 0, 15.

Figura 4. Comparación con los resultados numéricos para η de Wang y James [13]  and [14] y experimentales de Phillips y Deakin [23] . Influencia del factor CL . Ue  = 2 m/s.

Las discrepancias registradas con los resultados de Rafee et al. [21] utilizando el modelo RSM de turbulencia merecen una atención algo más detallada. En la figura 5 se han representado los resultados de Rafee et al. [21] para la pérdida de carga sobre el separador estudiado por Wang y James [13]  and [14] , y se han comparado con los obtenidos en este trabajo. Pueden extraerse las 2 siguientes conclusiones: por un lado, la pérdida de carga obtenida resulta mayor que la de Rafee et al. [21] , y por el otro, la tendencia de ΔP con la velocidad resultó parabólica (turbulenta), en contraposición a la de Rafee et al. [21] , fundamentalmente lineal (laminar). Teniendo en cuenta que también se han comparado los resultados de la pérdida de carga con los numéricos de Galletti et al. [11] (a su vez validados con los experimentales de Ghetti [24] ), obteniendo una variación media solo del 7 % y con una misma tendencia cuadrática, podría conjeturarse con la existencia de una discrepancia en las condiciones de contorno o en la geometría reproducida por Rafee et al. [21] . También es preciso apuntar que estos autores no especifican el código con el que han obtenido sus resultados.

Figura 5. Comparación con los resultados numéricos para la caída de presión ΔP de Rafee et al. [21] .

En las figuras 6 y 7 se muestran las comparaciones realizadas con los resultados numéricos de Galletti et al. [11] y experimentales de Ghetti [24] , respectivamente para U_e  = 3 y 5 m/s. Tal como se ha indicado anteriormente, la utilización del modelo SST k  − ω de turbulencia con valores de y⁺ suficientemente bajos en las paredes e incluyendo la dispersión turbulenta de las gotas, deviene en obtener resultados muy ajustados a los experimentales, probablemente por la descripción detallada del flujo que intenta conseguir el modelo en la proximidad de las paredes (lógicamente, la eficiencia de captura debe de depender de forma muy importante de la veracidad del flujo simulado numéricamente junto a las paredes). Las discrepancias más grandes registradas en η , en función de utilizar o no la dispersión turbulenta de las gotas tienen lugar para valores pequeños de D_p , probablemente por el hecho de que para gotas con poca masa, el efecto aleatorio de la turbulencia predomina sobre la inercia de las partículas a seguir las líneas de corriente del flujo medio de aire. Finalmente, puede concluirse que la modelización de la dispersión turbulenta de las gotas conseguida a través del código Fluent es similar a la obtenida con el «varied EIM» por parte de autores como Wang y James [13]  and [14] y Galletti et al. [11] .

Figura 6. Validación con los resultados numéricos para η de Galletti et al. [11] y experimentales de Ghetti [24] , para un separador horizontal en forma de zig-zag dotado de ganchos. Ue  = 3 m/s.

Figura 7. Validación con los resultados numéricos para η de Galletti et al. [11] y experimentales de Ghetti [24] , para un separador horizontal en forma de zig-zag dotado de ganchos. Ue  = 5 m/s.

5.2. Influencia de considerar la dispersión turbulenta sobre los separadores considerados

En el apartado anterior se ha expuesto que los efectos de la dispersión turbulenta de las gotas son más significativos cuanto más pequeño es el tamaño de la gota. Por otro lado, cabe esperar que la morfología del separador sea un factor igualmente influyente. Con tal fin, en las figuras 8 a, b, c y d se han representado los resultados obtenidos para η en función de D_p , respectivamente para 2 tipos considerados de separador, y para las 3 relaciones de aspecto en estudio. Puede decirse que en general la diferencia en los resultados obtenidos crece al disminuir D_p (siguiendo el comportamiento expuesto antes) y al aumentar la relación de aspecto b /L . Cuando la geometría del separador es sencilla, como en el caso del «Belgian wave» (fig. 8 a), la diferencia entre los resultados es relativamente pequeña, llegando a hacerse despreciable para tamaños de gota suficientemente grandes. Si la morfología del separador es más compleja, como en el caso del «L-shaped» (fig. 8 b,c y d), las diferencias son más significativas (esto puede confirmarse de forma cualitativa observando que el patrón de flujo es más complejo, y que la aleatoreidad del camino seguido por las gotas es mayor), no registrándose una dependencia importante con respecto de la velocidad del flujo en el rango 1.5 ≤ U_e  ≤ 5 m/s.

Figura 8. Evaluación de la influencia de considerar o no la dispersión turbulenta de las gotas sobre la eficiencia de captura η , para diferentes relaciones de aspecto b /L . a) «Belgian wave», Ue  = 2, 5 m/s. b) «L-shaped», Ue  = 1, 5 m/s. c) «L-shaped», Ue  = 2, 5 m/s. d) «L-shaped», Ue  = 5 m/s.

5.3. Influencia del parámetro inercial y de la relación de aspecto

El análisis dimensional muestra que la eficiencia de captura puede ponerse en función del parámetro inercial P_i , introducido por autores como Sriveerakul et al. [22] . Manteniendo constantes el resto de variables, P_i recoge el efecto combinado de U_e y de D_p . De alguna manera, distintos valores de la velocidad y del tamaño de las partículas pueden producir el mismo efecto inercial, de forma que la deposición sea equivalente, y por tanto también lo sea la eficiencia de captura. En las figuras 9 y 10 , para los separadores «Belgian wave» y «L-shaped» respectivamente, puede observarse cómo las distintas series de puntos (cada uno de los 3 tipos corresponde a velocidades U_e  = 1, 5, 2,5 y 5 m/s)), convergen básicamente a una misma tendencia salvo pequeñas diferecias, para cada una de las relaciones de aspecto consideradas. Es importante reseñar que aunque el separador «L-shaped» es geométricamente más complicado, tiene una dependencia similar al «Belgian wave» con respecto a P_i .

Figura 9. Eficiencia de captura en función del parámetro inercial Pi . Para cada relación de aspecto b /L , los diferentes símbolos corresponden a velocidades Ue  = 1, 5, 2,5 y 5 m/s. Separador «Belgian wave».

Figura 10. Eficiencia de captura en función del parámetro inercial Pi . Para cada relación de aspecto b /L , los diferentes símbolos corresponden a velocidades Ue  = 1, 5, 2,5 y 5 m/s. Separador «L-shaped».

La eficiencia de captura es tanto mayor cuanto menor es la relación de aspecto, para un mismo valor de P_i , lo que puede explicarse simplemente por el hecho de que el fluido tiene un espacio más reducido para atravesar el separador (aunque como inconveniente, la pérdida de carga puede aumentar de forma significativa). La tendencia asintótica que aparece para b /L  = 1/3 y el separador «Belgian wave» se debe a que la relativamente alta separación de las placas permite la formación de una corriente de aire con gotas que apenas se ven perturbadas por la presencia de las paredes, lo que claramente invalida el objetivo de obtener una eficiencia de captura plena, y llevaría a descartar este tipo de separador. Lógicamente, para geometrías más complicadas, este efecto llegaría a producirse para mayores valores de b /L .

5.4. Comparación de los resultados numéricos Fluent-Phoenics

Como se ha indicado anteriormente, uno de los objetivos de este trabajo consiste en comparar la bondad de la simulación numérica obtenida a través del código Phoenics. Hasta el momento, se ha comprobado que la simulación numérica conseguida con Fluent es satisfactoria y próxima a resultados tomados de la bibliografía. Se procederá a continuación a comparar los resultados obtenidos alternativamente con Fluent y Phoenics. En la figura 11 se muestran ejemplos de la simulación de las trayectorias lagrangianas obtenidas con Phoenics, para los 4 separadores en estudio. Como cabía esperar, las mayores eficiencias de captura se consiguen en general en el orden «Belgian wave» (menor η )–«H1-V»–«L-shaped»–Zig-zag (mayor η ).

Figura 11. Trayectorias de las gotas obtenidas para los 4 separadores en estudio, utilizando Phoenics (modelo k  − ω , incluyendo dispersión turbulenta de las partículas). Dp  = 15 × 10−6 m, Ue  = 2, 5 m/s y b /L  = 1/6.

En las figuras 12 a, b, c y d se muestran los resultados obtenidos para el coeficiente de pérdida de presión K_P en función del número de Reynolds, para los 4 tipos de separador en estudio, respectivamente. Las tendencias obtenidas parecen ser en cualquier caso lógicas, siendo fundamentalmente descendentes con Re para el caso del «Belgian wave« y del «H1-V», y más constantes en el caso de los tipos de forma más compleja, como el «L-shaped» y el «Zig-zag». En estos 2 últimos perfiles, llegan a invertirse las tendencias seguidas para las diferentes relaciones de aspecto, lo que puede vincularse con los cambios drásticos que se producen en los patrones de flujo, sobre todo por el desarrollo de zonas de recirculación para valores de Re suficientemente altos, y morfologías más difíciles. Es de reseñar que tanto Fluent como Phoenics reproducen el comportamiento descrito, siendo los patrones de flujo simulado muy parecidos en cualquier caso. Sin embargo, existe una diferencia mantenida en las 4 figuras, siendo en general los valores de K_P siempre más altos con Fluent que con Phoenics, a pesar de que los valores medios de y⁺ se mantuvieron en un mismo rango en todas las simulaciones, y que los mallados fueron similares. Esto puede llevar a pensar que las pequeñas diferencias que puedan existir en la implementación de los modelos en cada código produzca ciertas discrepancias en flujos transicionales, que esta ocasión no pueden considerarse despreciables.

Figura 12. Coeficiente de pérdida de presión KP en función de Re. Comparación con los resultados obtenidos con Fluent y Phoenics. a) «Belgian wave». b) «H1-V». c) «L-shaped». d) «Zig-zag».

En cuanto a la eficiencia de captura, de la figura 13 pueden extraerse varias consideraciones de interés. En primer lugar, para el tipo «Belgian wave», los valores obtenidos para η son, salvo alguna diferencia que se comenta despúes, básicamente los mismos en ambos códigos. Además, el comportamiento registrado para los distintos valores de b /L es también coincidente, incluyendo la tendencia asintótica para b /L  = 1/3. Puede afirmarse entonces que la modelización lagrangiana del comportamiento de las partículas de Phoenics, incluyendo las dispersión turbulenta, consigue una aproximación correcta al problema, pero con algunas restricciones. Por ejemplo, no ha sido posible conseguir una buena aproximación para valores muy pequeños de P_i , para los que la dispersión turbulenta parecía ser preponderante (sobre todo para b /L pequeños), resultado valores de η muy por encima de los conseguidos con Fluent. El efecto comentado no se produce si no se incluye la dispersión turbulenta de partículas, como puede verse en las figuras 14 a y b, para los separadores «H1-V» y «L-shaped», respectivamente. En general, las diferencias entre los resultados obtenidos crecen con la complejidad de la morfología del separador, agravándose la falta de precisión en el cálculo de la deposición de las gotas para bajos valores de D_p (y de P_i . En cualquier caso, resulta curioso destacar que la eficiencia de captura plena se ha simulado en ambos códigos prácticamente para el mismo valor de P_i , siempre que se incluya en la simulación la dispersión turbulenta (incluso para el tipo «Zig-zag», más complejo).

Figura 13. Eficiencia de captura η en función del parámetro inercial Pi . Comparación de los resultados obtenidos con Fluent y con Phoenics. Separador «Belgian wave».

Figura 14. Eficiencia de captura η en función del parámetro inercial Pi . Comparación de los resultados obtenidos con Fluent y con Phoenics. a) Separador «H1-V». b) Separador «L-shaped».

Por todo lo explicado anteriormente, se conviene en que la modelización numérica obtenida con Fluent para el problema propuesto parece ser más aproximada y completa que la de Phoenics, aunque naturalmente siempre es posible incorporar más variables y chequear más posibilidades de los códigos, de forma que no puede interpretarse como una verdad absoluta. La necesidad de utilizar la herramienta BFC en Phoenics para reproducir la geometría de forma fideligna en casos como los estudiados provoca la presencia de más dificultades de convergencia numérica que en Fluent (en este trabajo, las mallas fueron construidas con el mallador Gambit), de forma que los tiempos de computación han sido más largos en general en Phoenics. Es posible que las dificultades comentadas hayan tenido alguna influencia en la bondad de los resultados finales obtenidos en el caso de Phoenics para η y para K_P .

5.5. Predicción del punto de eficiencia plena

Se ha visto hasta el momento que la evolución de la eficiencia de captura η en función de P_i sigue una tendencia creciente, normalmente con forma de S, que lleva a pensar en una curva sigmoide (es decir, como el inverso de una exponencial negativa). Aunque correlacionar este comportamiento parece complicado para todas las situaciones estudiadas (se deja planteado para investigaciones futuras), sí ha sido posible encontrar funciones exponenciales simples para predecir el valor de la eficiencia de captura máxima (η_max  = 100 %) en función del parámetro inercial P_i , para los 4 tipos de separador estudiados. Para ello, los separadores se caracterizan mediante el parámetro geométrico de eficiencia Λ , que se introduce a continuación como

$\Lambda ={\left({\frac {L_{w}}{L}}\right)}^{2}\quad \sum _{i}^{n}{\Theta }_{i}.$

( 37)

Este parámetro cuantifica de alguna forma la influencia de la relación entre la longitud de las paredes mojadas del separador ${\textstyle L_{w}}$ y la longitud neta del mismo L , así como del número n y de la magnitud Θ_i de los cambios de dirección de la línea media del separador. El valor de Λ se incrementa con la complejidad de la morfología del separador, de modo que cabe esperar que valores altos de Λ se correspondan con valores altos de η . Los valores de Λ calculados fueron de 2,651, 4,806, 6,945 y 7,857, respectivamente para los separadores «Belgian wave», «H1-V», «L-shaped» y «Zig-zag».

De esta manera, la expresión

$P_{i,max}=1,62\quad e^{-0,191\Lambda },$

( 38)

se ajusta con un 5 % de error medio a los resultados numéricos obtenidos para b /L  = 1/12, mientras que la expresión

$P_{i,max}=4,09\quad e^{-0,259\Lambda },$

( 39)

lo hace con un 8 %, para b /L  = 1/6.

En cuanto a la relación de aspecto b /L  = 1/3, efectos tales como la tendencia asintótica que se ha reproducido para valores altos de P_i imposibitan la obtención de una correlación similar. De forma alternativa, si consideramos solo los valores de P_{i ,max} correspondientes a los separadores «H1-V» y «Zig-zag», la correlación obtenida sería la siguiente:

$P_{i,max}=5,13\quad e^{-0.263\Lambda },$

( 40)

y que lógicamente no se ajusta a la simulación numérica obtenida para los tipos «Belgian wave» y «L-shaped».

Los valores prospectivados por las expresiones (38) , (39) y (40) se han representado mediante los círculos de mayor tamaño en las figura 15 a, b y c. Puede afirmarse que el ajuste realizado para las relaciones de aspecto b /L  = 1/12 y 1/6 es satisfactorio, mientras que para b /L  = 1/3 (figura 15 c), se revelan efectos adicionales significativos; por un lado, un cambio de tendencia en el comportamiento del separador «L-shaped», y por el otro, el comportamiento asintótico explicado antes para el «Belgian wave». En la figura 15 c, se han escrito los valores de Λ característicos de cada separador que corresponderían a cada condición predicha por las correlaciones propuestas.

Figura 15. Predicción del punto de eficiencia de captura plena (η  = 100%), a través de las correlaciones propuestas en este trabajo. Los resultados presentados se han obtenido con Fluent. a) b /L  = 1/12. b) b /L  = 1/6. c) b /L  = 1/3.

6. Conclusiones

• Se han obtenido resultados numéricos satisfactorios del flujo de aire con gotas dispersas que tiene lugar en los separadores de gotas de las torres de refrigeración de tiro mecánico. Se han obtenido resultados para parámetros relevantes tales como la eficiencia de captura y el coeficiente de pérdida de carga para 4 tipos distintos de separadores comerciales.
• Se han comparado los resultados con otros de la bibliografía, tanto numéricos como experimentales. Se han obtenido altos grados de acuerdo con Wang y James [13]  and [14] , Phillips y Deakin [23] , Galletti et al. [11] y Ghetti [24] , y alguna discrepancia significativa con Rafee et al. [21] .
• Se ha probado la bondad de los resultados obtenidos con 3 modelos de turbulencia: k  − ϵ , k  − ω y RSM, empleando opciones para bajos números de Reynolds, y tratamiento mejorado cerca de las paredes. Los resultados más fiables se obtuvieron con el modelo SST k  − ω , incluyendo la dispersión turbulenta de las gotas. La restricción de valores de la distancia adimensionalizada a la pared y⁺ a un rango 0, 2 − 0, 5 ha sido determinante para la consecución de resultados aproximados a los experimentales, e independientes de la densidad de malla.
• Se ha determinado la influencia del coeficiente C_L para el cálculo del tiempo característico lagrangiano en la aproximación de los resultados numéricos cuando se simula la dispersión turbulenta de las gotas. En cualquier caso, es necesario la consideración de la dispersión turbulenta, y en mayor medida, cuando la complejidad de la morfología del separador aumenta.
• Aunque se han obtenido resultados aceptables con el código Phoenics, los resultados obtenidos con Fluent parecen ser más consistentes y ajustados a los experimentales disponibles en la bibliografía.
• Se han presentado correlaciones para calcular el valor del parámetro inercial P_i para el que se produce la eficiencia máxima de captura (η_max  = 100 %), en función de cierto parámetro geométrico Λ que caracteriza la morfología de cada separador, y que ha sido introducido en este trabajo.
• En futuras investigaciones, se espera poder correlacionar la eficiencia de captura η en función de P_i , e incluso plantear un análisis técnico-económico para poder escoger el mejor separador posible de entre los 4 estudiados, para una situación o requerimiento dado. Desde el punto de vista numérico, puede ser interesante evaluar los efectos de la tridimensionalidad del flujo.

Agradecimientos

Esta investigación ha sido subvencionada a través de los Proyectos de Investigación ENE2007-68106-C02-02 y ENE2010-21679-C02-01 del Plan Nacional de I+D+I (Ministerio de Ciencia e Innovación).

References

1. [1] J.C. Kloppers, D.G. Kröger; Cooling tower performance evaluation: Merkel, Poppe, and e-NTU Methods of Analysis; ASME J. Eng. Gas Turbines Power, 127 (2005), pp. 1–7
2. [2] G. Gan, S.B. Riffat; Numerical simulation of closed wet cooling towers for chilled ceiling systems; Appl. Thermal Eng., 19 (1999), pp. 1279–1296
3. [3] G. Gan, S.B. Riffat, L. Shao, P. Doherty; Application of CFD to closed-wet cooling towers; Appl. Thermal Eng., 21 (2001), pp. 79–92
4. [4] R.N. Meroney; CFD prediction of cooling tower drift; J. Wind End. Ind. Aerodyn., 94 (2006), pp. 463–490
5. [5] R.N. Meroney; Protocol por CFD prediction of cooling-tower drift in an urban environment; J. Wind End. Ind. Aerodyn., 96 (2008), pp. 1789–1804
6. [6] P.M. Foster, M.I. Williams, R.J. Winter; Droplet behaviour and collection by counterflow cooling tower eliminators; Atmos. Environ., 8 (1974), pp. 349–360
7. [7] B.R. Gardner, H.J. Lowe; The research and development background to the environmental problems of natural draught cooling towers; Atmos. Environ., 8 (1974), pp. 313–320
8. [8] J. Chan, M.W. Golay; Comparative performace evaluation of current design evaporative cooling tower drift eliminators; Atmos. Environ., 11 (1977), pp. 775–781
9. [9] B.R. Becker, L.F. Burdick; Effect of drift eliminator design on cooling tower performance; ASME J. Eng. Gas Turbines Power, 114 (1992), pp. 632–642
10. [10] P.W. James, Y. Wang, B.J. Azzopardi, J.P. Hughes; The role of drainage channels in the performance of wave-plate mist eliminators; Trans. IChemE, 81 (A) (2003), pp. 639–648
11. [11] C. Galletti, E. Brunazzi, L. Tognotti; A numerical model for gas flow and droplet motion in wave-plate mist eliminators with drainage channels; Chem. Eng. Sci., 63 (2008), pp. 5639–5652
12. [12] P.W. James, Y. Wang, B.J. Azzopardi, J.P. Hughes; A model for liquid film flow and separation in a wave-plate mist eliminator; Chem. Eng. Res. Des., 83 (A5) (2005), pp. 469–477
13. [13] Y. Wang, P.W. James; The calculation of wave-plate demister efficiencies using numerical simulation of the flow field and droplet motion; Trans. IChemE, 76 (A) (1998), pp. 980–985
14. [14] Y. Wang, P.W. James; Assessment of an eddy-interaction model and its refinements usig predictions of droplet deposition in a wave-plate demister; Trans. IChemE, 77 (A) (1999), pp. 692–698
15. [15] N.C. Reis, R.F. Griffiths, J.M. Santos; Numerical simulation of the impact of liquid droplets on porous surfaces; J. Comput. Phys., 198 (2004), pp. 747–770
16. [16] G. Strotos, M. Gavaises, A. Theodorakakos, G. Bergeles; Numerical investigation of the cooling effectiveness of a droplet impinging on a heated surface; Int. J. Heat Mass Transfer, 51 (2008), pp. 4728–4742
17. [17] C. Mundo, M. Sommerfeld, C. Tropea; Droplet-wall collitions: experimental studies of the deformation and breakup process; Int. J. Multiphase Flow, 21 (1995), pp. 151–173
18. [18] M. Sommerfeld, G. Kohnen, M. Ruger; Some open questions and inconsistencies of Lagrangian particle dispersion models; En Ninth Symposium Turbulent Shear Flows, 2005, Kyoto, Japan (1993), pp. 16–18
19. [19] G.A. Kallio, M.W. Reeks; A numerical simulation of particle deposition in turbulent boundary layers; Int. J. Multiphase Flow, 15 (3) (1989), pp. 433–446
20. [20] C. Verlann, Performance of novel mist eliminators. Ph.D. Thesis, Delf University of Technology, Delf, Holanda (1999).
21. [21] R. Rafee, H. Rahimzadeh, G. Ahmadi, Numerical simulations of airflow and droplet transport in a wave-plate mist eliminator. Chem. Eng. Res. Des., in press (2010).
22. [22] T. Sriveerakul, I.L. McLaine-cross, M. Behnia; Performance prediction for cooling tower eliminators; En The 17th Annual Conference of Mechanical Engineering Network, Prajinburi, Tailandia (2003)
23. [23] H. Phillips, A.W. Deakin; Measurements of the collection efficiency of various demister devices; En Proceedings of the 4th Annual Meeting of The Aerosol Society, 1990, Loughborough, Reino Unido (1990)
24. [24] S. Ghetti, Investigation of entrainment phenomena in inertial separators. Ph.D. Thesis, University of Pisa, Pisa, Italia (2003).
25. [25] H.K. Versteeg, W. Malalasekera; An introduction to Computational Fluid Dynamics; Prentice-Hall (1995)
26. [26] A.N. Kolmogorov; Equations of turbulent motion of an incompressible fluid; Izv. Acad. Sci., USSR; Physics, 6 (1,2) (1942), pp. 56–58
27. [27] D.C. Wilcox; Turbulence modeling for CFD; (2nd ed.)DCW Industries, USA (2003)
28. [28] F.R. Menter; Two-equation eddy-viscosity turbulent models for engineering applications; AIAA J., 32 (8) (1994), pp. 1598–1605
29. [29] Fluent INC, Fluent 6.3 Users Guide (2006). http://www.fluent.com [consultado 3 Feb 2011].
30. [30] B. Zamora, A.S. Kaiser, A. Viedma; On the effects of Rayleigh number and inlet turbulence intensity upon the buoyancy-induced mass flow rate in sloping and convergent channels; Int. J. Heat Mass Transfer, 51 (2008), pp. 4985–5000
31. [31] A.S. Kaiser, B. Zamora, A. Viedma; Numerical correlation for natural convective flows in isothermal heated, inclined and convergent channels, for high Rayleigh numbers; Comp. Fluids, 38 (2009), pp. 1–15
32. [32] Gentra User Guide (Phoenics), TR211 (2006). http://www.cham.co.uk [consultado 3 Feb 2011].
33. [33] Phoenics support: Encyclopaedia. Disponible en: http://www.cham.co.uk [consultado 3 Feb 2011].