Estudio preliminar de un sistema de tratamiento de salmuera para vertido líquido cero (Zero Liquid Discharge) alimentado por gases de escape de micro turbina de gas

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales El proyecto SOLMIDEFF está financiado por el Plan Estatal Retos de Investigación del Ministerio de Ciencia, Innovación y Universidades (RTI2018-102196-B-I00) Estudio preliminar de un sistema de tratamiento de salmuera para vertido líquido cero (Zero Liquid Discharge) alimentado por gases de escape de micro turbina de gas Autor: Francisco Béjar Hernández Tutores: David Tomás Sánchez Martínez Lourdes García Rodríguez Departamento de Ingeniería Energética Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 iii Trabajo Fin de Grado Grado en Ingeniería de tecnologías industriales Estudio preliminar de un sistema de tratamiento de salmuera para vertido líquido cero (Zero Liquid Discharge) alimentado por gases de escape de micro turbina de gas Autor: Francisco Béjar Hernández Tutores: David Tomás Sánchez Martínez Lourdes García Rodríguez Profesores titulares Dpto. de Ingeniería Energética Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 v Trabajo de Fin de Grado: Estudio preliminar de un sistema de tratamiento de salmuera para vertido líquido cero (Zero Liquid Discharge) alimentado por gases de escape de micro turbina de gas Autor: Francisco Béjar Hernández Tutores: David Tomás Sánchez Martínez, Lourdes García Rodríguez El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2020 El Secretario del Tribunal vii A mis maestros A mi familia que siempre estuvo ahí A mi novia que me apoyó en los peores momentos ix Resumen El proyecto SOLMIDEFF (Solar micro gas turbine-driven desalination for environmental off-grid applications), financiado por el plan Nacional de Investigación del Gobierno de España tiene como objetivo el desarrollo de un sistema combinado de desalación de agua alimentado por energía solar, apto para instalaciones en en lugares aislados de la red eléctrica. En este trabajo se realiza un estudio teórico preliminar de un sistema de vertido cero o ZLD (Zero Liquid Discharge), instalado aguas debajo del escape de una microturbina de gas solar, el cual consiste en un proceso térmico para recuperar agua de la corriente de salmuera que, de otro modo, sería vertida al ambiente. Esta característica de SOLMIDEFF incorpora ventajas muy importantes respecto de los sistemas convencionales basados en tecnología fotovoltaica. El sistema ZLD consiste en un regenerador al cual entra un gas caliente en flujo horizontal y, tras pasar por un codo, asciende y se encuentra con la salmuera procedente del proceso de ósmosis, la cual se inyecta a contracorriente y por gravedad a cierta altura dentro del canal, evaporando el agua de la salmuera para separarla de la sal. En este estudio se incluye el análisis de la geometría del regenerador, el mallado empleado y unas condiciones de contorno suficientes para que el sistema alcance el objetivo de evaporar el agua de la salmuera introducida (caso base). Tras un estudio completo sobre el efecto de cada variable en el sistema se demuestra que el del tamaño de las gotas introducidas y su caudal toman un papel fundamental en la evaporación de las gotas. Finalmente, empleando la información recabada a lo largo del trabajo, se desarrollan mejoras para el caso base, consiguiendo duplicar la cantidad de salmuera evaporada en relación a un cierto caudal de gases de escape de la microturbina de gas. ÍNDICE DE FIGURAS Figura 1. Objetivos principales de SOLMIDEFF [2] 1 Figura 2. Sistema de disco solar de OMSoP [4] 3 Figura 3. Diagrama de un proceso ZLD [6] 4 Figura 4. Diagrama de flujo de un proceso ZLD con sistemas de membrana y de tratamiento térmico [7] 5 Figura 5. Tecnologías de pre-tratamiento de salmuera: (a) electrodiálisis, (b) osmosis inversa, (c) destilación de membrana [5] 5 Figura 6. Flujo a través de un módulo de membrana [13] (citado en [11, p.5]) 6 Figura 7. Vasija a presión [16] (citado en [11, p.5]) 6 Figura 8. Fase de evaporación/cristalización en un dispositivo convencional [5] 8 Figura 9. Mapa de los regímenes de interacción entre partículas y turbulencia [20] 11 Figura 10. Entalpía de cambio de fase de disoluciones de sales marinas, en kJ/kg 14 Figura 11. Elevación del punto de ebullición de disoluciones de sales marinas, en K 14 Figura 12. Geometría del caso base 17 Figura 13. Dimensiones de la geometría del caso base 17 Figura 14. Geometría alternativa con más anchura en el canal vertical 18 Figura 15. Geometría alternativa con entrada brusca en la zona de intercambio 18 Figura 16. Función 𝑦+ [37] 20 Figura 17: Valor de 𝑦+ a lo largo de las paredes del dominio 20 Figura 18. Mallado completo del sistema 22 Figura 19. Malla aumentada 22 Figura 20. Codo aumentado con la función Inflation cerca de la pared 23 Figura 21. Conexión del codo con el canal vertical 23 Figura 22. Canal horizontal 24 Figura 23. Canal vertical 24 Figura 24. Inyecciones de salmuera en la fase continua. 26 Figura 25. Contorno de la velocidad a lo largo del dominio 27 Figura 26. Diámetro de las partículas 28 Figura 27. Humedad relativa 29 Figura 28. Temperatura 29 Figura 29. Gráfica de energía cinética turbulenta 30 Figura 30. Fracción másica de H2O 30 Figura 31. Presión 31 Figura 32. Evaporación de DPM 31 xvii Figura 33. Gráfica de concentración de DPM 32 Figura 34. Gráfica de evaporación de DPM 32 Figura 35. Gráfica de fracción másica de O2 33 Figura 36. Gráfica de fracción másica de H2O 33 Figura 37. Gráfica de temperatura 34 Figura 38. Gráfica de velocidad 34 Figura 39. Gráfica de velocidad radial 35 Figura 40. Gráfica de velocidad tangencial 35 Figura 41. Concentración de DPM a distintas alturas del eje y. 38 Figura 42. Energía cinética turbulenta a distintas alturas del eje y 38 Figura 43. Fracción másica de H2O a distintas alturas del eje y 39 Figura 44. Humedad relativa a distintas alturas del eje y 39 Figura 45. Temperatura del flujo a distintas alturas del eje y 40 Figura 46. Evaporación de DPM a distintas alturas del eje y 40 Figura 47. Velocidad a distintas alturas del eje y 41 Figura 48. Representación de la temperatura y la humedad relativa para varios caudales de salmuera. Las barras muestran el parámetro en cuestión (izquierda) y la línea el CV (derecha) 43 Figura 49. Diámetro de las partículas para distintos caudales de salmuera 44 Figura 50. Diámetro de las partículas para distintos tamaños de partícula 46 Figura 51. Representación de la temperatura y la humedad relativa para varios tamaños de partícula. Las barras muestran el parámetro en cuestión (izquierda) y la línea el CV (derecha) 46 Figura 52. Representación de la temperatura y la humedad relativa para varias temperaturas de entrada de salmuera. Las barras muestran el parámetro en cuestión (izquierda) y la línea el CV (derecha) 49 Figura 53. Representación de la temperatura y la humedad relativa para varias velocidades del aire a la entrada. Las barras muestran el parámetro en cuestión (izquierda) y la línea el CV (derecha) 50 Figura 54. Gráfica de perfiles de la fracción másica de H2O para v=7 m/s a distintas alturas dentro del canal vertical 50 Figura 55. Representación de la temperatura y la humedad relativa para varias temperaturas del aire a la entrada. Las barras muestran el parámetro en cuestión (izquierda) y la línea el CV (derecha) 51 Figura 56. Diámetro de partículas en el caso 1 55 Figura 57. Gráfica de la fracción másica de H2O a distintas alturas del eje de ordenadas para el caso 1 55 Figura 58. Gráfica de la temperatura a distintas alturas del eje de ordenadas para el caso 1 56 Figura 59. Posición de las líneas de inyección de salmuera en el caso 2. El valor de las cotas se encuentra en la Figura 13 57 Figura 60. Gráfica de la fracción másica de H2O a distintas alturas del eje de ordenadas para el caso 2 57 Figura 61. Diámetro de partículas en el caso 2 58 Figura 62. Temperatura en el caso 2 58 Figura 63. Diámetro de partículas en el caso 3 59 Figura 64. Gráfica del número de elementos de la malla en los distintos casos 62 Figura 65. Gráficas de variación de parámetros índice según el caso 64 Figura 66. Gráficas de variación relativa al caso anterior de los parámetros índice según el caso 65 xix Notación ZLD Zero Liquid Discharge SOLMIDEFF Solar Micro gas turbine-driven Fesalination for Environmental Off-grid applications SMGT Solar Micro Gas Turbine OMSoP Optimised Microturbine Solar Power system FO Ósmosis convencional RO Ósmosis inversa ED Electrodialisis EDR Electrodialisis reversible MD Destilación de membrana PV Tubo a presión CFD Computational Fluid Dynamics DNS Simulaciones numéricas directasB CV Coeficiente de variación DPM Lagrangian discrete phase model EES Engineering Equation Solver MIT Instituto tecnológico de Massachusetts HR Humedad relativa UDF Función definida por el usuario INTRODUCCIÓN 1.1. Motivación del proyecto l objetivo de este trabajo es el diseño, modelado y análisis de un sistema ZLD (Zero Liquid Discharge), uno de los procesos que forman parte de SOLMIDEFF (Solar Micro gas turbine-driven Desalination for Environmental Off-grid applications), un proyecto multidisciplinar cuyo objetivo principal es desarrollar y demostrar una tecnología innovadora, a pequeña escala que posibilite la producción de agua en ubicaciones remotas y sin acceso a electricidad [1]. En este punto se hace una introducción a dicho proyecto, de qué partes se compone, qué áreas de la ingeniería implica y qué aportación ofrece este trabajo. SOLMIDEFF es un proyecto financiado por el Ministerio de ciencia e innovación y el Fondo Europeo de Desarrollo Reginal que consiste, esencialmente, en el desarrollo conceptual de un sistema de producción de energía eléctrica y agua potable a partir de energía solar, empleando para ello una microturbina de gas solar (SMGT). Este motor, con un diseño optimizado, produce energía eléctrica y térmica que es empleada por una unidad de ósmosis inversa en combinación con un sistema de vertido líquido nulo para producir agua potable con un impacto ambiental prácticamente nulo [2]. La integración de estos sistemas puede verse en la Figura 1 de manera conceptual. Figura 1. Objetivos principales de SOLMIDEFF [2] Como se detalla en la propia página web de SOLMIDEFF [1], estos sistemas podrían constituir una solución ideal y rentable para el suministro eléctrico en comunidades aisladas, ya que permite un único sistema generador de energía con mucho menos mantenimiento y costes de operación que otras tecnologías combinadas como fotovoltaica + baterías o fotovoltaica + diésel. Las microturbinas solares en el rango de 5 a 30 kWe han sido probadas recientemente en Europa en proyectos como OMSoP [3], demostrando el potencial de producir electricidad (con un coste de 10-15 c€/kWh) y energía térmica a alta temperatura a 250 ºC a partir de la energía solar, y permitiendo una fácil hibridización con combustible líquido o gas para regímenes de operación prolongados (o incluso continuos). E Introducción 2 El proyecto SOLMIDEFF se compone de activididades teóricas y experimentales, y este trabajo pretende contribuir a este desarrollo mediante la simulación y diseño del proceso óptimo de desalación del sistema ZLD que encaje con las condiciones y características que requiere el proyecto. La motivación principal de desarrollar un proceso de desalación óptimo y a pequeña escala es contribuir al progreso de comunidades en ubicaciones remotas produciendo agua que pueda servir tanto para el consumo humano como para servicios agrícolas. Ya existen tecnologías disponibles para solventar este problema, pero o requieren una potencia eléctrica no disponible (ósmosis inversa dependiente de la energía) o se apoyan en sistemas que no permiten una producción continua de electricidad. Además, el coste de esta tecnología podría ser inadmisible en ciertas ubicaciones [1]. 1.2. Organización del proyecto Este trabajo comienza con una descripción del proyecto SOLMIDEFF en el capítulo 2, incluyendo de qué partes se compone y cuál es su fin último, haciendo especial hincapié en el sistema ZLD, marco en el que se desarrolla este trabajo. A lo largo del la subsección 2.1 se desarrolla cómo funciona este dispositivo y los usos que se le puede dar en distintas industrias. El capítulo 3 trata los modelos numéricos empleados según la literatura para proyectos de la misma tipología, describiendo los fenómenos físicos que ocurren en sistemas similares en la subsección 3.1. En las subsecciones 3.2 y 3.3 se centra la atención en la simulación en CFD mediante ANSYS Fluent, detallando cómo se aplican los modelos teóricos a las simulaciones. El capítulo 4 constituye el eje central del proyecto, tratando el caso base propuesto, así como analizando variaciones del mismo. En primer lugar se hace una descripción del sistema y sus características, comentándose las características de la geometría en la subsección 4.1, de la malla diseñada para la simulación en la subsección 4.2 y de las condiciones fluidodinámicas de contorno en ella subsección 4.3. Tras esto, se exponen los resultados obtenidos y se comenta su sentido y fiabilidad en la subsección 4.4. En este último se analizan varias aspectos de los resultados obtenidos: se discute sobre el caso base (término que se empleará a lo largo de este trabajo) en la subsección 4.4.1; y se analiza la influencia que cada una de las variables de contorno tiene sobre el sistema, dejando en evidencia los parámetros que mejor definen los fenómenos que se observan (subsecciones 4.3.1 y 4.3.2). Una vez conocidos todos estos datos, se proponen nuevas condiciones de contorno con el objetivo de mejorar la actuación del caso base en la subsección 4.5. En el capítulo 5 se ejecuta un análisis de sensibilidad de la malla empleada en el caso base, validandola y analizando la influencia que una mayor o menor bondad del mallado puede tener sobre el sistema. El proyecto finaliza con un apartado de conclusiones en el que se hace un resumen de todas las lecciones aprendidas durante la realización del mismo. PROYECTO SOLMIDEFF l funcionamiento del sistema, según la propia documentación de SOLMIDEFF [2] es el siguiente. Un disco parabólico recoge y concentra la energía solar en un plano focal donde se instala la unidad de conversión de potencia (Power Conversion Unit o PCU), la cual se compone de un receptor solar y una microturbina de gas. El compresor de la microturbina aspira aire atmosférico, el cual se comprime y, tras su calentamiento en el intercambiador regenerativo, se introduce en el receptor solar. El receptor también es un intercambiador de calor que transfiere la energía solar concentrada recibida del colector a la corriente entrante de gas presurizado, devolviéndola a alta temperatura y presión. Este aire, típicamente alrededor de los 800ºC es expandido en la turbina, donde se produce la potencia para accionar tanto el compresor como el generador eléctrico. Finalmente, antes de ser descargados a la atmósfera, los gases de salida de la turbina son conducidos a través del regenerador para precalentar el aire de entrada en el compresor. Dentro del regenerador, las corrientes fría y caliente suelen disponerse a contracorriente, disposición empleada también en el proyecto OMSoP [3]. Figura 2. Sistema de disco solar de OMSoP [4] El segundo subsistema del concepto SOLMIDEFF se compone de dos elementos [2]: un sistema avanzado de desalación basado en tecnología de osmosis inversa y una unidad de ZLD (Zero Liquid Discharge). El sistema de osmosis inversa produce agua dulce y salmuera con una alta concentración de sales, la cual es tratada en un dispositivo de recuperación de agua de la corriente de salmuera por medios térmicos. Este sistema, conocido generalmente como Zero Liquid Discharge (ZLD), es el objeto de este trabajo. El ZLD es el elemento en el que se centra este trabajo y su funcionamiento, prestaciones y características serán descritas a lo largo del mismo. Ambos sistemas que componen el sistema inferior incorporan profundas innovaciones con respecto a la tecnología actual. 2.1. Sistema Zero Liquid Discharge (ZLD) 2.1.1 Objetivos, aplicaciones y funcionamiento En este subapartado se explica brevemente en qué procesos consiste y de qué partes se compone el sistema ZLD convencional, dejando de lado consideraciones económicas Como ya se adelantó anteriormente, el Zero Liquid Discharge (ZLD) es realmente proceso de tratamiento de efluentes que se puede emplear en una gran variedad de industrias y cuyo objetivo es eliminar, en la medida de lo posible, el residuo líquido de un vertido. Estos sistemas se pueden emplear para [5]: E Proyecto SOLMIDEFF 4  Reducir el gasto de agua residual de una industria, disminuyendo los gastos de gestión de la misma.  Reciclar agua y, por tanto, disminuir la cantidad de agua demandada por una industria  Disminuir la cantidad de residuos que vierte al ambiente  Producir agua limpia apta para su reutilización. Estos sistemas se emplean en una gran cantidad de campos: industrias químicas, textiles, alimentarias, en plantas de potencia, para potabilización de aguas municipales, reutilización de agua para usos agrícolas, etc. Además, más allá de los beneficios que ofrece, el incremento de las regulaciones medioambientales por parte de los gobiernos alrededor del mundo hace necesaria la implantación de ZLD en un futuro cercano. En la Figura 3 se representa un diagrama que destaca el lugar que generalmente ocupa un ZLD en un proceso industrial: Figura 3. Diagrama de un proceso ZLD [6] En el caso que nos ocupa, el motivo principal del empleo de un ZLD, como ya se adelantó previamente, es emplear la energía del gas residual procedente de la microturbina de gas para aumentar la salinidad de una salmuera, extrayendo de ella agua que se pueda emplear en otros procesos. Esto tiene el beneficio múltiple de reducir el impacto ambiental, mejorar la manejabilidad del residuo y aumentar el aprovechamiento energético de la energía solar capturada. El diseño de un ZLD depende de tres factores principalmente: los contaminantes específicos que se puedan encontrar en el flujo que se vaya a tratar, la cantidad de flujo y el volumen del material disuelto. En sistemas convencionales se suelen considerar dos etapas [5]: una primera etapa de pre-concentración, la cual se consigue generalmente con concentradores de salmuera de membrana o electrodiálisis con el objetivo de incrementar la salinidad y recuperar hasta el 60-80% del agua; y una segunda etapa de evaporación/cristalización, la cual se lleva a cabo mediante procesos térmicos, evaporando el agua y recogiéndola para su reutilización, además de cristalizando las impurezas para poder tratarlas como sólidos. 2.1.2 Tratamientos de membrana En realidad, un ZLD puede ser diseñado únicamente mediante procesos térmicos con concentradores y cristalizadores, tecnologías ampliamente desarrolladas en grandes industrias tras años de desarrollo tecnológico. Sin embargo, los procesos térmicos requieren una gran inversión de capital y consumen enormes cantidades de energía. Por tanto, el coste de los ZLD puede ser reducido significativamente empleando membranas previamente al tratamiento térmico, como se puede apreciar en la Figura 4. Figura 4. Diagrama de flujo de un proceso ZLD con sistemas de membrana y de tratamiento térmico [7] En la literatura se han ido proponiendo distintas tecnologías de membrana [7] como la osmnosis convencional (FO), la osmosis inversa (RO), la electrodiálisis (ED), la electrodiálisis reversible (EDR) y la destilación de membrana (MD). En la Figura 5 se presenta un esquema del funcionamiento de estas tecnologías. Figura 5. Tecnologías de pre-tratamiento de salmuera: (a) electrodiálisis, (b) osmosis inversa, (c) destilación de membrana [5] a) La ED es un proceso de membrana que emplea electrodos para crear un campo eléctrico que empuja los iones positivos y negativos a través de membranas semipermeables que atrapan los iones positivos y negativos respectivamente b) La RO coinsiste en una membrana osmótica semipermeable a través de la cual fluye el agua cuando es sometida a una presión superior a la presión osmótica de la disolución en sentido contrario al flujo natural de agua desde la disolución de menor concentración a la más concentrada. De esta forma, se logra separar la sustancia concentrada de la permeada, esta última baja en contenidos disueltos [8]. Ya que este es el sistema empleado en SOLMIDEFF, se profundiza en el funcionamiento de la osmosis inversa en el siguiente subapartado. c) La MD es un proceso de transporte accionado térmicamente en el que la diferencia de presión de vapor entre los dos lados de la membrana hidrófoba permite un intercambio de calor y masa a través de la misma. 2.1.3 Ósmosis inversa En el proyecto SOLMIDEFF el proceso empleado es el de ósmosis inversa (RO), el cual está bien establecido Modelo de cálculo 12 12 Para tener una idea superficial sobre el comportamiento del sistema, se observará en qué categoría de las tres descritas se podría clasificar a nuestro sistema y se mencionará junto con los resultados en el subapartado correspondiente 3.2. Métodos de resolución en CFD de flujos turbulentos conductores de partículas Una forma completa de determinar el comportamiento de un fluido turbulento es el método de simulación numérica directa (DNS). Este consiste en la resolución completa de las ecuaciones de Navier-Stokes para obtener el campo de velocidades completo del flujo turbulento en todas las escalas (sin emplear modelos de aproximación para las escalas de menor longitud características). Aunque este podría ser el medio más deseable para tener la mayor precisión del comportamiento de los fluidos, en la práctica no es posible su aplicación debido a la gran complejidad y coste computacional que supondría, pues este método conlleva la resolución exacta de las ecuaciones de Navier-Stokes de cada partícula del fluido, computando todas las escalas de torbellinos que puedan existir en el problema. Este método, debido a las limitaciones anteriormente descritas, está restringido a dominios pequeños. Otra opción mucho menos compleja que la anterior es emplear aproximaciones adecuadas de las ecuaciones de Navier-Stokes incorporando modelos de turbulencia para caracterizar el comportamiento de los fenómenos turbulentos de menor longitud característica. Esta técnica de simulación se conoce como ‘Reynolds-Averaged Navier-Stokes (RANS) Equations’y, aunque no proporciona la solución exacta del dominio fluido, sí permite abordar problemas complejos y de gran carga computacional. Existen diversas clases de modelos de turbulencias para la simulaciones mediante RANS [30]: (i) Zero-equation models. En esta categoría solo se resuelve un sistema diferencial de ecuaciones parciales (PDEs) para el flujo principal. (ii) One-equation models. Con respecto al (i), esta categoría implica una ecuación de transporte adicional para el cálculo de la velocidad en la escala turbulenta, normalmente expresada en términos de la energía cinética turbulenta media (k). (iii) Two-equation models. Con respecto al (ii), esta categoría implica una ecuación de transporte adicional para el cálculo de la longitud característica de la escala turbulenta, normalmente en términos de la disipación de energía cinética turbulenta (ε). (iv) Stress-equation models. Con respecto a (i), esta categoría implica un número de ecuaciones adicionales de transporte para los componentes del tensor de tensiones de Reynolds (𝜏𝑖𝑗) y uno por el índice escalar de la disipación de energía cinética turbulenta (ε). Por este motivo, los modelos de la categoría (iv) son también llamados los modelos (𝜏𝑖𝑗−𝜀) Un método intermedio entre los anteriores es el conocido como Large Eddy Simulations (LES) o simulación de grandes vórtices [20] en el cual se resuelven las ecuanciones del flujo para los fenómenos turbulentos de mayor longitud características mientras los de menor escala son resueltos mediante modelos simplificados. Esta solución constituye un punto intermedio de complejidad entre la resolución completa de las ecuaciones de Navier-Stokes (DNS) y la simplificación que implica emplear modelos para conocer el comportamiento de los fluidos (RANS). En el apartado 164 se especifica qué método se emplea en este problema y por qué motivos. 3.3. Propiedades termodinámicas y estimaciones preliminares Una estimación de la cantidad de agua que puede evaporar una cierta cantidad de gas caliente (proveniente en la aplicación real del escape de una turbina de gas) podría ser realmente útil para afrontar la simulación del sistema ZLD mediante CFD, teniendo unas nociones básicas de cómo se comporta el sistema. En primer lugar, se presenta un cálculo aproximado de la relación de masas entre fases, de manera simplificada, suponiendo una entalpía de cambio de fase del agua salada y una temperatura de vaporización que no se ven afectados por la salinidad de la misma. A fin de realizar este cálculo, se considera el sistema ZLD como si fuera un intercambiador de calor en el cual se aplica el principio de conservación de energía. Para ello se considera igualmente que el intercambio de calor entre fases (fluidos) sólo se produce en el proceso de evaporación y que la diferencia mínima de temperaturas entre fluidos es nula (𝛿𝑇=0): 𝑚𝑔·𝐶𝑝𝑔·(𝑇𝐸−𝑇𝑠𝑎𝑡)≈𝑚𝑠·(∆ℎ𝑣𝑎𝑝) (1) 𝑚𝑠 𝑚𝑔=0,1969≈0,2 (2) Siendo:  𝐶𝑝𝑔=1,05 𝑘𝐽/(𝑘𝑔· 𝐾) la capacidad térmica isobárica másica del gas  ∆ℎ𝑣𝑎𝑝=2257 𝑘𝑔/𝑘𝐽 la entalpía de cambio de fase del agua  𝑇𝐸=250 °𝐶 la temperatura de entrada del gas  𝑇𝑠𝑎𝑡=100 °𝐶 la temperatura de saturación del agua  𝑚𝑔 el caudal de gas  𝑚𝑠 el caudal de salmuera No obstante lo anterior, la existencia de salmuera (en lugar de agua dulce) podría cambiar estos resultados, pues la temperatura de vaporización y la entalpía de cambio de fase cambiarían con la concentración de la misma. Estas variaciones se pueden observar en las figurasFigura 10. y Figura 11. y en las tablasTabla 1 yTabla 2, documentación procedente de la librería de propiedades termofísicas del agua salada del Instituto Tecnológico de Massachusetts MIT [31] referidas a disoluciones de sales marinas. Estos datos están recogidos a una presión de 1 atm para temperaturas menores de 100°C y a la temperatura de la presión de vapor de agua de mar para temperaturas mayores de 100°C. Tabla 1. Entalpía de cambio de fase de disoluciones de sales marinas, en kJ/kg Modelo de cálculo 14 14 Figura 10. Entalpía de cambio de fase de disoluciones de sales marinas, en kJ/kg Tabla 2. Elevación del punto de ebullición de disoluciones de sales marinas, en K Figura 11. Elevación del punto de ebullición de disoluciones de sales marinas, en K Este efecto se estudiará en la siguiente sección, sabiendo que el orden de magnitud de la salinidad de la corriente de salmuera que sale del proceso de osmosis inversa se estudiará en el rango entre 40000 y 50000 ppm que, en mg/kg, implica un rango entre 40 y 50 mg/kg. Aunque no se ha tratado aún la temperatura de la salmuera entrante en el caso que nos ocupa, en estos datos aportados vemos una elevación máxima de la temperatura de saturación que ronda el 2% y una disminución máxima del calor latente que ronda el 12%. Estos datos se deben tener en cuenta en las simulaciones en CFD pues los cambios en estos parámetros y su influencia sobre el sistema pueden llegar a ser muy significativos según la concentración y temperatura de la corriente de salmuera entrante. Simulación con CFD 16 16 SIMULACIÓN CON CFD l flujo bifásico en sistemas de enfriamiento evaporativo, como el existente en nuestro sistema, se ve influenciado por gran cantidad de parámetros: velocidad, temperatura y humedad relativa de la fase continua; distribución de tamaños de gota, velocidad y temperatura de las gotas; y las interacciones que puedan existir entre gotas o entre gotas y fase continua. Para conocer el funcionamiento de nuestro sistema de manera fiable y lo suficientemente precisa, es importante saber el impacto de cada parámetro. Los programas CFD (Computational Fluid Dynamics) son herramientas valiosas para analizar las variaciones y comportamiento de estos parámetros. Como se menciona en [32], existen tres métodos principales para evaluar la actuación del enfriamiento evaporativo: medidas a gran escala, medidas en túnel de viento y simulaciones numéricas con CFD. Las medidas a gran escala ofrecen la ventaja de que se estudia la situación real, teniendo en cuenta toda la complejidad del problema. Sin embargo, no se puede tener control total sobre las condiciones de contorno por lo que la experimentación queda limitada. Las medidas a pequeña escala en túnel de viento permiten un gran control sobre las condiciones de contorno, aunque se requiere garantizar un buen número de condiciones de semejanza entre el experimento y el sistema real, algo que frecuentemente no es posible. Además, las medidas en el túnel de viento son habitualmente tomadas en un número limitado de puntos en el espacio de estudio (sección de ensayo) mientras que el CFD nos ofrece datos de los parámetros relevantes sobre el dominio completo. A diferencia del test en el túnel de viento, las herramientas CFD permiten recrear con exactitud todos los detalles y parámetros del sistema que queremos medir. Las simulaciones en CFD permiten evaluar distintas configuraciones de diseño, especialmente cuando las distintas configuraciones son, a priori, en base al mismo dominio computacional y mallado. Sin embargo, es necesario asegurar la precisión y fiabilidad del CFD por lo que es imperativa la verificación y validación de los estudios realizados. Con estas consideraciones, el uso del CFD se está incrementando en el estudio de un enorme rango de procesos atmosféricos y medioambientales como, por ejemplo, la dispersión de la contaminación en el aire, el efecto del viento sobre la lluvia o la transferencia de calor convectiva. Dadas las ventajas que esta herramienta proporciona, podría proporcionar una aproximación muy buena al funcionamiento fluidodinámico del ZLD. El software empleado en este trabajo ha sido ANSYS Fluent. Según la propia documentación de ANSYS [33], este programa es capaz de analizar los siguientes datos del modelo: cálculo de la trayectoria de la fase discreta usando formulación que incluya la inercia de la misma, arrastre hidrodinámico y fuerza de la gravedad; predicción de los efectos de la turbulencia en la dispersión de las partículas debido a torbellinos turbulentos presentes en la fase continua; calentamiento o enfriamiento de la fase discreta; vaporización y ebullición de gotas líquidas; combustión de partículas; y ruptura y coalescencia de gotas. Todas estas opciones de modelado de las que dispone el programa permiten hacer una simulación precisa, fiel y completa del proceso de introducción de las gotas y evaporación en la zona de intercambio, incluyendo todos los aspectos físicos y fluidodinámicos presentes en el problema. Los apartados sucesivos se centran en el sistema que nos ocupa, haciendo un análisis profundo de las peculiaridades del sistema, los valores de los parámetros que lo definen y de los resultados obtenidos. En primer lugar se describe el caso base para, en el apartado de resultados, tratar casos modificados con el objetivo de compararlos. 4.1. Geometría La geometría del sistema será bidimensional, lo cual representa suficientemente el problema para poder ser estudiado con fiabilidad. La geometría del evaporador consiste en un canal vertical que se conecta mediante un codo con una entrada horizontal en el que el intercambio de energía se produce a contracorriente: el gas caliente E entra horizontalmente y, tras pasar por un codo, asciende y se encuentra con la salmuera, la cual se inyecta a contracorriente y por gravedad a cierta altura dentro del canal. Al final del canal, en la parte superior, se encuentra la salida (outlet) por donde sale el gas humedecido por el agua de la salmuera que se ha evaporado. El resto de superficies de la geometría son paredes adiabáticas (las propiedades de las paredes se comentarán más adelante). Se ha adoptado un modelo bidimensional con la función Design Modeler dentro de ANSYS, una herramienta lo suficientemente versátil como para obtener una geometría sencilla como es el caso que nos ocupa. Concretamente, la geometría base a partir de la cual se ha simulado el problema queda representada en la Figura 12, con sus dimensiones dadas en la Figura 13. El canal tiene una anchura de 1 metro en toda su longitud y está dividido en tres partes para afinar el posterior mallado: canal horizontal, codo y canal vertical. Esta división se puede ver en la figura con las tres partes divididas por línas azules. La geometría tiene una longitud suficiente como para garantizar la no influencia de las condiciones de contorno. Además, como se verá en siguientes apartados, con esta geometría se satisfacen los objetivos que se persiguen en este proyecto. Figura 12. Geometría del caso base Figura 13. Dimensiones de la geometría del caso base Además de esta geometría, que resultará ser definitiva, también se ha evaluado el comportamiento del sistema con otras geometrías. Las dos más significativas han sido la geometría de la Figura 14. Geometría alternativa con más anchura en el canal vertical, con un canal con el doble de anchura en la zona vertical; y la geometría de la Figura 15. Geometría alternativa con , con una conexión brusca entre la entrada horizontal y el canal vertical. Estas geometrías se descartaron durante el estudio debido a su peor comportamiento. Simulación con CFD 18 18 Figura 14. Geometría alternativa con más anchura en el canal vertical Figura 15. Geometría alternativa con entrada brusca en la zona de intercambio 4.2. Mallado Como se comenta en [34], generar el mallado óptimo en el dominio consume un tiempo importante. Esto hace que haya que establecer un equilibrio entre la precisión buscada y el coste computacional (aspectos más profundamente detallados en el apartado 5). En este trabajo, la bondad del refinamiento de la malla se valora de dos maneras: mediante un análisis de sensibilidad, el cual se lleva a cabo detenidamente en el apartado 5 y medante el parámetro 𝑦+, descrito a continuación. Para ilustrar y comentar la importancia del parámetro 𝑦+ nos apoyaremos en la Figura 16, en la que se presenta la velocidad tangencial adimensional del flujo 𝑈+ frente a la distancia normal a la pared 𝑦+, ambos parámetros definidos en las ecuaciones (3) y (4), siendo y la distancia medida desde la pared, U la velocidad del flujo en función de dicha distancia y 𝑢𝜏 una velocidad de referencia basada en el esfuerzo cortante en zona próxima a la pared. 𝑢𝜏 queda definido en la ecuación (5), donde 𝜏𝜔 es dicho esfuerzo y 𝜌 la densidad del flujo. Tal y como se comenta en [35], numerosos experimentos han demostrado que la forma de la función resultante corresponde a la línea roja. La región cercana a la pared puede ser dividida, por sus diferentes características en tres capas. En la capa más cercana a la pared, llamada subcapa viscosa (viscous sublayer), el flujo es casi laminar, y la viscosidad toma un papel fundamental en la cantidad de movimiento y la transferencia de masa o calor. En la capa más alejada de la pared, llamada capa completamente turbulenta (fully-turbulent layer), la energía turbulenta es el parámetro que predomina. Finalmente, entre ambas capas existe la capa amortiguada (buffer layer), donde los efectos de la viscosidad y la energía turbulenta son igualmente importantes. Estas tres zonas se muestran también en la Figura 16. 𝑦+=𝑦𝑢𝜏 𝜈 (3) 𝑈+=𝑈 𝑢𝜏 (4) 𝑢𝜏=√𝜏𝜔 𝜌 (5) Para entender el interés de conocer el valor de 𝑦+ y su importancia, es necesario conocer cómo modela ANSYS Fluent el comportamiento del flujo cerca de la pared. ANSYS es capaz de predecir el comportamiento del flujo aproximando la curva roja de dos formas: haciendo un ajuste lineal si 𝑦+ es menor o igual a 11,067 o un ajuste logarítmico si 𝑦+ es mayor que 11,067. Cabe aclarar que el ajuste lineal aparece representado con una línea curva y viceversa porque la escala del eje horizontal es logarítmica, característica que puede llevar a confusión. Si el centroide de la celda adyacente a la pared tiene un 𝑦+<10, ANSYS resolverá el flujo de manera precisa en la subcapa viscosa mediante el ajuste lineal mencionado. Por otro lado, si el centroide de la célda adyacente a la pared se encuentra en el intervalo 30<𝑦+<200, ANSYS también podría predecir el comportamiento del flujo con cierta precisión empleando funciones de pared. Así pues, es recomendable evitar que el 𝑦+ sea menor que 10 o menor que 30, zona en la que los resultados serían calculados con menos precisión, ya que ANSYS no puede ajustar correctamente ninguno de los dos modelos [36]. En la Figura 16 puede comprobarse cómo la capa amortiguada es donde los ajustes propuestos por ANSYS (líneas azules) están más lejos del comportamiento real del sistema (línea roja). La mayor precisión, como se comenta en la literatura, se consigue para bajos 𝑦+, e idealmente se busca que 𝑦+<1, un valor que asegura la precisión total del modelo. Las regiones cercanas a la pared tienen mayores gradientes de los parámetros característicos del flujo y para ello se requiere medir con precisión estos bruscos cambios cercanos a la pared, para así poder obtener un resultado fiable en estas zonas más comprometidas. Por último, cabe destacar que este parámetro viene afectado por el refinamiento de la malla en las zonas a la pared y el número de Reynolds en dichas zonas, por lo que conseguir un 𝑦+ aceptable también viene determinado por las condiciones del fluido y no únicamente por las características del mallado. Sin embargo, dadas las condiciones del problema y la bondad de la malla creada, limitada por el poder de computación disponible durante la realización de este proyecto, no ha sido posible alcanzar el objetivo de 𝑦+< 1. En la Figura 17 se representa el valor de 𝑦+obtenido en el caso base a lo largo del dominio en función de la posición en el eje x. En dicha figura no se quiere hacer hincapié en la forma de la función representada, sino más bien en el rango de valores en el que se encuentra y la precisión de la malla en distintos puntos de esta. Simulación con CFD 20 20 Figura 16. Función 𝑦+ [37] Figura 17: Valor de 𝑦+ a lo largo de las paredes del dominio En cuanto al tipo de malla empleada en el sistema, esta puede ser estructurada o no estructurada. La empleada en este sistema es una malla estructurada que, como se explica en [38], se caracteriza por tener los vértices interiores topológicamente similares, pues forman una red de cuadriláteros (en nuestro caso bidimensional; en el caso tridimensional formarían hexaedros) a lo largo del dominio. En general, a diferencia de las mallas no estructuradas, las mallas estructuradas ofrecen una mayor simplicidad y mejor convergencia. Una malla no estructurada se adaptaría mejor a un dominio de geometría compleja por lo que no es necesaria en nuestro caso, donde la geometría se podría dividir únicamente se compone de dos rectángulos y un codo, formas sencillas de mallar con cierto grado de precisión. Una opción más precisa es hacer un mallado híbrido, un enfoque que disfruta de las ventajas de ambos métodos. El mallado se ha hecho a partir de algunas funciones sencillas: (i) Face Meshing. Con esta función, ANSYS hace un mapeado del dominio y crea una malla estructurada con cuadriláteros. Es útil en dominios rectangulares o aproximadamente rectangulares, por lo que es especialmente útil en los canales vertical y horizontal, aunque también es capaz de crear una malla ordenada y fiable en el codo. (ii) MultiZoneQuad/Tri Method. El objetivo principal del uso de esta función es poder emplear la función Inflation, además de descomponer el mallado en las tres partes mencionadas anteriormente (iii) Inflation. Esta función permite hacer un mallado en el que las celdas van disminuyendo de tamaño conforme nos acercamos a la pared. De esta forma, podemos conseguir una mayor precisión en estas zonas conflictivas. Se emplea una función inflation en cada una de las tres partes del dominio: En el codo el número de capas es 40 (se requiere una gran precisión en esta zona, donde el fluido cambia de dirección y hay una fuerte influencia de la pared); en el canal horizontal el número de capas es 6; y en el canal vertical el número de capas es 25. En las tres zonas el índice de crecimiento entre capa y capa es de 1,2. Configurado para tener un tamaño de elemento máximo de 0,008 metros, ejecutar el mallado con todas las restricciones impuestas da como resultado una malla formada por 355760 elementos y 357900 nodos. A continuación se presentan algunas imágenes de la malla resultante. La Figura 18 es el mallado completo del sistema, pero debido a la densidad de la malla no se puede apreciar cómo es esta en detalle, por lo que se representan las siguientes imágenes enfocando únicamenta a partes concretas de la misma. Cabe destacar un par de peculiaridades:  Se puede observar en la Figura 19 una mayor densidad del mallado entre el codo y los canales horizontal y vertical. Para aplicar el Face Meshing en el dominio (y que, por tanto, la malla sea estructurada) ANSYS debe detectar que la cara donde se aplica dicha función es un cuadrilátero, o aproximadamente un cuadrilátero. De esta forma, para conseguir la malla estructurada fue necesario dividir el dominio en tres partes, canal vertical, horizontal y codo. Además, luego se aplicó la función Inflation, la cual necesita ser aplicada, al menos, en tres bordes de cada cuadrilátero. Todo esto da como resultado que cada una de las tres partes tenga un mallado más fino en, al menos, tres de sus bordes, apareciendo esa mayor densidad en zonas donde no hay pared.  En la Figura 23 hay una mayor densidad de malla a la salida para obtener un resultado más fiable en esa zona, pues es de gran interés tener una buena precisión sobre las condiciones de temperatura y saturación del gas a la salida. Simulación con CFD 28 28 concreta (la salmuera) se reduce a cero a partir de cierto punto, por lo que podemos confirmar que se evapora completamente. De hecho, se puede afirmar que la evaporación se produce de forma bastante eficaz, pues ocurre en apenas los primeros centímetros del canal tras el punto de inyección. Nuevamente está presente el efecto de la pared en el tiempo de residencia de las partículas, pues las partículas más cercanas a esta son transportadas por el fluido a menor velocidad. Como ya se adelantaba en apartados anteriores, se produce el enfriamiento evaporativo de la fase continua. Esto se pone de manifiesto en el aumento de la humedad relativa en la corriente y el descenso de la temperatura. Estos aspectos se pueden observar, cualitativamente, en las figuras Figura 27 y Figura 28. También se puede observar el efecto del codo en las turbulencias producidas a lo largo de todo el canal vertical. Esto se puede ver gráficamente y de manera cuantitativa en la Figura 29. Esta figura muestra varias líneas que representan la variación de energía turbulenta a lo largo del eje horizontal a varias alturas del eje vertical, incluyendo la salida. La mayor perturbación del flujo se produce en la cara interna del canal y esta va disminuyendo a medida que aumenta la distancia desde el codo. Para tener un conocimiento completo y preciso sobre el comportamiento del sistema es necesario el análisis de los datos desde varios puntos de vista. Las figuras Figura 30,Figura 31 yFigura 32 se presentan algunos diagramas de contorno adicionales a los ya expuestos. Por otro lado, las figuras comprendidas entre la Figura 33 y la Figura 40 representan algunos gráficos que muestran información adicional sobre la variación de algunos de los parámetros característicos que definen el sistema a lo largo del eje horizontal para varias alturas del eje vertical. Figura 26. Diámetro de las partículas Figura 27. Humedad relativa Figura 28. Temperatura Simulación con CFD 30 30 Figura 29. Gráfica de energía cinética turbulenta Figura 30. Fracción másica de H2O Figura 31. Presión Figura 32. Evaporación de DPM Simulación con CFD 32 32 Figura 33. Gráfica de concentración de DPM Figura 34. Gráfica de evaporación de DPM Figura 35. Gráfica de fracción másica de O2 Figura 36. Gráfica de fracción másica de H2O Simulación con CFD 34 34 Figura 37. Gráfica de temperatura Figura 38. Gráfica de velocidad Figura 39. Gráfica de velocidad radial Figura 40. Gráfica de velocidad tangencial Siguiendo el mismo criterio de medición que con la energía turbulenta se pueden obtener los resultados de varios Simulación con CFD 36 36 parámetros significativos que definen la variación de las condiciones del flujo a lo largo del canal vertical. Estos resultados se pueden observar de forma numérica en la Tabla 3. Parámetros significativos en el canal vertical., donde se presenta los valores medios de ciertos parámetros para varios valores del eje y en el canal vertical. Algunos de los parámetros presentan variaciones apreciables en un corte transversal al eje del conducto vertical por lo que el valor medio puede no representar adecuadamente los valores máximo y mínimo del parámetro en cuestión en la localización correspondiente, pero la tabla es una aproximación que puede ser intuitiva en la mayoría de casos. Adicionalmente, para poner de manifiesto esta variabilidad, estos valores vienen representados junto a su coeficiente de variación (CV), el cual representa el porcentaje de dispersión de los datos medios presentados. Los datos con un CV menor del 25% pueden considerarse homogéneos, por lo que la magnitud media representada tiene cierto grado de fiabiliad; mientras que si el CV es mayor que 25% los datos comienzan a ser más heterogéneos y estos datyos medios van perdiendo fiabilidad conforme se acercan a 100%. Los datos de la Tabla 3 vienen representados gráficamente en las figuras comprendidas entre entre la Figura 41 y la Figura 47, en las que se representa en azul el parámetro en cuestión (izquierda) y en rojo el CV (derecha). Existen datos medios cuyo CV superan el 100%, lo cual significa que el dato no es fiable. Todas estas medidas han sido tomadas en ANSYS fluent, al igual que las anteriores. Las desviaciones típicas empleadas para el cálculo del CV se pueden consultar en la Tabla 5. Tabla 3. Parámetros significativos en el canal vertical. Posición Parámetros y=4,5m y=4,9m y=5m y=5,5m y=6m y=6,5m y=7m Outlet Concentración DPM (kg/m³) 0 0,0115 0,0077 0 0 0 0 0 Energía cinética de las turbulencias (m²/s²) 0,5815 0,5762 0,5913 0,5774 0,5595 0,5467 0,5356 0,5264 Fracción másica de H2O 0 0,00056 0,00736 0,01829 0,18270 0,01823 0,01821 0,01819 Humedad relativa (%) 0 0,00206 0,02231 0,08069 0,07837 0,07737 0,07680 0,07642 Temperatura (K) 573 570,2 550,7 521,7 521,7 521,8 521,9 521,9 Velocidad (m/s) 16 15,93 15,57 15,02 15,02 15,02 15,02 15,02 Evaporación de DPM (kg/s) 0 3,6∙10⁻⁵ 6,9∙10⁻⁵ 7,1∙10⁻⁶ 0 0 0 0 Tabla 4. Desviación típica de parámetros significativos en el canal vertical Posición Parámetros y=4,5m y=4,9m y=5m y=5,5m y=6m y=6,5m y=7m Outlet Concentración DPM (kg/m³) 0 0,0658 0,0225 0 0 0 0 0 Energía cinética de las turbulencias (m²/s²) 0,5318 0,5124 0,5194 0,4884 0,4530 0,4253 0,4067 0,3081 Fracción másica de H2O 0 0,00361 0,00544 0,00508 0,00395 0,00340 0,00300 0,00201 Humedad relativa (%) 0 0,02240 0,04123 0,06846 0,04300 0,03415 0,02805 0,01807 Temperatura (K) 0 13,32 16,44 14,42 11,35 9,84 8,56 5,77 Velocidad (m/s) 2,36 2,39 2,36 2,32 2,28 2,35 2,19 1,49 Evaporación de DPM (kg/s) 0 1,2∙10⁻⁴ 8,9∙10⁻⁵ 1,0∙10⁻⁵ 0 0 0 0 A continuación, las representaciones gráficas junto con algunas observaciones y conclusiones que se pueden extraer de la evolución de los distintos parámetros.  Concentración de DPM: Solo está presente en las zonas cercanas a las inyecciones, pues el agua se evapora mayoritariamente entre las cotas y=4,9m e y=5,5m. Además, existe un gran efecto de pared que hace que la concentración sea mucho mayor en las paredes, haciendo que el CV se dispare, como Simulación con CFD 44 44 Figura 49. Diámetro de las partículas para distintos caudales de salmuera En segundo lugar, se evalúa el impacto de variar el tamaño de las gotas. Como ya se ha comentado en apartados anteriores, aunque en la mayoría de sistemas de inyección se obtiene un cierto rango de gotas de distinto tamaño, en este subapartado también se suponen todas las gotas del mismo tamaño. Estas comparaciones se abordan de la misma forma que la comparación anterior de caudales de entrada de la fase discreta, y se pueden consultar en la Tabla 6 y la Figura 51. Nuevamente esta comparación se presenta de forma gráfica sobre diagramas de contorno del diámetro de las partículas, pues nos da una visión muy intuitiva de lo que le ocurre a las partículas. Tabla 6. Evaporación de las partículas aumentando el diámetro de las gotas introducidas Diámetro de las gotas introducidas (μm) Temperatura media de salida (K) Humedad relativa media de salida (%) Evaporación de las partículas 20 521,89 0,07646 30/30 40 (base) 521,89 0,07643 30/30 60 521,90 0,07643 30/30 80 521,88 0,07650 30/30 100 526,45 0,06484 29/30 El único cambio observado está en el caso de 100 μm, escenario en el cuál no se evaporan todas las partículas introducidas, por lo que deja de ser interesante seguir aumentando el tamaño de las mismas. En este caso, debido a que el enfriamiento evaporativo es menor, la temperatura de salida media es mayor y la humedad relativa m=0,018 kg/s m=0,024 kg/s media es menor. Como se puede observar, el sistema es bastante sensible a ligeros cambios del tamaño de las gotas para la configuración planteada en el caso base, por lo que convendría tener un gran control sobre esta variable. Nuevamente, gracias a este estudio se puede comprobar que el sistema admite variaciones y mejoras, pues es posible introducir en el sistema un tamaño mayor de partículas, que siempre es la posibilidad más barata y sencilla de llevar a cabo en la mayoría de procesos de creación de gotas. D = 20 μm D = 40 μm D = 60 μm D = 80 μm Simulación con CFD 46 46 Figura 50. Diámetro de las partículas para distintos tamaños de partícula Figura 51. Representación de la temperatura y la humedad relativa para varios tamaños de partícula. Las barras muestran el parámetro en cuestión (izquierda) y la línea el CV (derecha) Como ya se ha comentado, en la práctica, se ha empleado agua en lugar de salmuera en el caso base por simplicidad. Cabe esperar que, si las propiedades del fluido inyectado cambian, también lo haga el comportamiento del sistema y, por tanto, las características del flujo de salida. Para el cálculo de las nuevas características de la fase discreta, considerando la salmuera, se ha empleado EES (Engineering Equation Solver). En dicho software se ha empleado una librería de funciones desarrollada por MIT [42] que permite calcular distintas propiedades del agua marina en función de la temperatura, la presión y la salinidad. De esta forma, a 0 10 20 30 40 50 60 70 80 90 100 519 520 521 522 523 524 525 526 527 20 μmBase 60 μm80 μm100 μm Temperatura media outlet (K) 0 10 20 30 40 50 60 70 80 90 100 0,058 0,06 0,062 0,064 0,066 0,068 0,07 0,072 0,074 0,076 0,078 20 μmBase 60 μm80 μm100 μm Humedad relativa media outlet (%) D = 100 μm continuación, se observa el impacto de cambiar los parámetros que definen el agua por los que definen la salmuera, cambiando su densidad, su calor específico y su capacidad calorífica para hacer que el fluido inyectado se parezca a agua salada lo máximo posible La salinidad de la corriente de salmuera que entra en el ZLD es, concretamente, de 46183 ppm, por lo que es razonable un estudio del efecto de la salinidad en un rango lo suficientemente amplio como para poder apreciar el efecto de la salinidad (entre 40000 ppm y 80000 ppm). En la Tabla 7 se observan las características del flujo inyectado según el grado de salinidad; parámetros calculados a partir de la función de EES mencionada. Tomando el calor latente constante con la temperatura, como se ha venido haciendo hasta ahora en todo este proyecto, el cambio de comportamiento del flujo a lo largo del canal y los cambios los parámetros del flujo a la salida del mismo apenas son apreciables. Esto se puede comprobar en la Simulación con CFD 48 48 Tabla 8. En ambas tablas se incluye una fila en referencia al caso base (salinidad 0) para facilitar la comparación. Tabla 7. Características del flujo inyectado según la concentración de sales Salinidad (g/kg) Densidad (kg/m3) Capacidad térmica isobárica (J/(kg·K)) Entalpía de cambio de fase (kJ/kg) Punto de ebullición (K) 0 998,2 4182 2263 373 40 1029 3973 2335 373,3 46,18 1034 3942 2340 373,4 50 1036 3924 2331 373,4 60 1044 3875 2306 373,5 70 1052 3829 2282 373,7 80 1060 3784 2257 373,8 Tabla 8. Parámetros en el outlet según la concentración de sales Salinidad (g/kg) Temperatura media a la salida (K) Humedad relativa a la salida (%) CV temperatura media a la salida (%) CV humedad relativa a la salida (%) 0 521,9 0,07642 1,106 23,65 40 520,6 0,07819 1,136 24,031 46,18 520,9 0,07779 1,129 23,94 50 521,1 0,07754 1,113 23,98 60 521,6 0,07687 1,131 23,60 70 522,0 0,07624 1,102 24,51 80 522,5 0,07559 1,091 23,44 Simular un fluido inyectado con estas características conlleva una mayor aproximación a la realidad que la inyección de una fase discreta compuesta únicamente por agua pura. Como se ha comprobado, esta nueva aproximación no supone grandes cambios en el comportamiento del sistema con respecto del caso base por lo que se podría tomar el fluido de salinidad 40,18 g/kg como referencia de cara a futuros trabajos que desarrollen este proyecto. Finalmente, cabe mencionar que el cambio de velocidad o temperatura en la corriente inyectada no introduce apenas cambios significativos en el comportamiento del sistema (al menos dentro de un margen lógico y acepable para la aplicación que se le quiere dar al recuperador). No obstante, se expone en la Tabla 9 el efecto que estas variables producen en la temperatura y la humedad relativa medias en el outlet, sin entrar en más detalle. Tabla 9. Variaciones de la temperatura y la humedad relativa medias en el outlet relativas al caso base si se modifican la velocidad o la temperatura de entrada de la fase discreta. Modificación Variación de la temperatura media de salida Variación de la humedad relativa media de salida Aumento de la velocidad de entrada de las gotas a 13 m/s a contracorriente + 0,0019 % + 0,2 % Aumento de la temperatura de entrada de las gotas a 353 K + 0,8 % - 7,2 % Las variaciones de la temperatura y la humedad cuando varía la velocidad de entrada de la salmuera son apenas apreciables, pudiendo ser causadas por la variación introducida en las variaciones del propio software. Sí se aprecia cierto cambio cuando se aumenta la temperatura de entrada de las gotas, produciéndose un aumento de la temperatura y una disminución de la humedad relativa. También se puede comparar dicho efecto en la Figura 52. Al subir la temperatura de la salmuera se consume menor calor sensible del gas y por lo tanto este se enfría Simulación con CFD 50 50 menos, dando como resultado un flujo más caliente a la salida. El aire húmedo, para la misma humedad absoluta y mayor temperatura, presenta una humedad relativa inferior. Además, la energía térmica disponible del gas de entrada independientemente de la temperatura de la salmuera, escenario que podría no darse si la energía térmica disponible del gas fuera menor. Figura 52. Representación de la temperatura y la humedad relativa para varias temperaturas de entrada de salmuera. Las barras muestran el parámetro en cuestión (izquierda) y la línea el CV (derecha) 4.3.2 Variaciones de los parámetros de la fase continua El análisis del sistema no solo se limita a observar el comportamiento del mismo ante cambios en la fase discreta. También son necesarias experimentaciones y mediciones con respecto a variaciones de los parámetros de la fase continua. Los parámetros de la fase continua cuya influencia se va a analizar son los siguientes:  Velocidad del aire a la entrada  Temperatura del aire a la entrada  Humedad del aire a la entrada El primer parámetro a analizar de la fase continua es la velocidad del aire a la entrada. Esta velocidad debe ser la suficientemente elevada como para evitar que la fase discreta se deposite en la base del canal por acción de la gravedad, pero tampoco demasiado alta para que arrastre las gotas sin que estas se evaporen. En un escenario ideal, las gotas se quedarían en suspensión con tiempos de residencia altos hasta que se evaporaran. Se estudian los cambios producidos en el sistema para casos en los que la velocidad disminuye, que es el terreno donde se pueden observar mejoras, pues empleando más velocidad que en el caso base (16 m/s) solo se consigue el arrastre de las gotas hacia el exterior, consiguiendo un enfriamiento evaporativo menos eficiente. En los casos aquí estudiados se reduce la velocidad del aire a 10 m/s y posteriormente a 7 m/s. La variación de la temperatura media y la humedad relativa media en el outlet se representa en la Figura 53, nuevamente representando también el coeficiente de variación para denotar la fiabilidad de estos datos medios. Como se puede observar en estas representaciones, el fenómeno del enfriamiento evaporativo es mucho más notorio conforme menor es la velocidad, haciendo que el sistema sea más eficiente. Esto conlleva, puede intuirse, un aumento en el tiempo de residencia: el tiempo de residencia máximo en el caso base es de 6,03e-02 segundos, mientras que en el caso en el que v=7 m/s, el tiempo de residencia máximo es 2,52e-01 segundos; más de un órden de magnitud mayor. También se puede comprobar el aumento del enfriamiento evaporativo con la disminución de la velocidad en la Figura 54 que, si se compara con la Figura 36, se puede observar que los valores de fracción másica de H2O en la corriente resultante son mayores para cualquier valor del eje vertical. 0 10 20 30 40 50 60 70 80 90 100 518 519 520 521 522 523 524 525 526 278 K Base 308 K 323 K 338 K Temperatura outlet media 0 10 20 30 40 50 60 70 80 90 100 0,069 0,07 0,071 0,072 0,073 0,074 0,075 0,076 0,077 0,078 0,079 278 K Base 308 K 323 K 338 K Humedad relativa Figura 53. Representación de la temperatura y la humedad relativa para varias velocidades del aire a la entrada. Las barras muestran el parámetro en cuestión (izquierda) y la línea el CV (derecha) Figura 54. Gráfica de perfiles de la fracción másica de H2O para v=7 m/s a distintas alturas dentro del canal vertical Pese a que hubiera sido interesante continuar bajando la velocidad para determinar el comportamiento del sistema hasta llegar al límite en el que las gotas caen por la fuerza de la gravedad, la simulación se vuelve inestable y los resultados obtenidos no son fiables. Por ello la mínima velocidad estudiada en este proyecto es 7 m/s. Por lo tanto, este trabajo puede ser mejorado y ampliado en el futuro con un estudio más exhaustivo y abarcando un mayor rango de velocidades. Otro parámetro clave en este estudio es la temperatura del aire a la entrada del canal. Si se observa con 0 10 20 30 40 50 60 70 80 90 100 420 430 440 450 460 470 480 490 500 510 520 530 Base 10 m/s 7 m/s Temperatura outlet media 0 10 20 30 40 50 60 70 80 90 100 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 Base 10 m/s 7 m/s Humedad relativa Simulación con CFD 52 52 detenimiento el caso base se puede observar que la cantidad de energía aportada a través de la fase continua es capaz de evaporar rápidamente las gotas, pues estas pasan a vapor en centésimas de segundo y recorriendo apenas unos centímetros. A causa de los procesos previos al ZLD, es razonable estimar que la corriente de aire estará a una temperatura perteneciente al intervalo 250-300 ºC. Este es el rango que se tiene en cuenta en este estudio. Disminuyendo la temperatura de la corriente de entrada, la Figura 55 muestra las variaciones de la corriente de salida si esta entra a 473 K ó a 574 K. Como se puede observar, conforme se disminuye la temperatura de entrada del flujo, mayor es el efecto del enfriamiento evaporativo. Considerando estos datos y los obtenidos en el análisis de la velocidad de entrada del flujo, se puede concluir que el flujo que entra en el canal en el caso base se caracteriza por tener demasiada temperatura y velocidad, por lo que para conseguir un enfriamiento evaporativo más eficaz sería más eficiente bajar estos valores disminuyendo el valor de ambos parámetros antes de entrar en el ZLD. Figura 55. Representación de la temperatura y la humedad relativa para varias temperaturas del aire a la entrada. Las barras muestran el parámetro en cuestión (izquierda) y la línea el CV (derecha) En el caso base se tomó la humedad del aire de entrada como nula, cuando dicho fluido no tiene por qué ser seco. De hecho, el aire empleado tendrá cierta humedad casi con toda seguridad. Como valores índice para evaluar el comportamiento del sistema ante distintos grados de humedad, se toman valores que se pueden encontrar en Sevilla, por ejemplo, a lo largo del año. Para ello, los datos empleados de la humedad relativa de la zona a lo largo del año se han tomado de [43]. En ANSYS la humedad del aire se introduce en términos de fracción másica, y para calcular este término se emplea una calculadora psicrométrica [44] en la que se introducen dos valores extremos de la humedad relativa en esta región, calculando ambas situaciones a 25ºC para estar del lado de la seguridad. De esta forma, las fracciones másicas empleadas en ambos casos son los siguientes:  HR = 25% → Fracción másica H2O: 0,0049  HR = 75% → Fracción másica H2O: 0,015 Como se puede intuir, para valores de humedad relativa mayores, mayor es la humedad relativa del flujo a su paso por el ZLD. De hecho, el aumento de este parámetro apenas afecta a otras variables del sistema, por lo que la aproximación del caso base (HR=0%) es suficiente como para obtener unos resultados fiables. Los cambios introducidos por la variación de este parámetro se recogen en la Tabla 10. 0 10 20 30 40 50 60 70 80 90 100 0 100 200 300 400 500 600 Base 473 K 423 K Temperatura outlet media 0 10 20 30 40 50 60 70 80 90 100 0 0,2 0,4 0,6 0,8 1 1,2 1,4 Base 473 K 423 K Humedad relativa Tabla 10. Parámetros en el outlet según la humedad relativa en el inlet Humedad relativa (%) Temperatura media outlet (K) Humedad relativa outlet (%) CV temperatura media outlet (%) CV humedad relativa outlet (%) 0 521,9 0,08 1,106 23,65 25 522,0 0,10 1,105 20,91 75 522,1 0,14 1,103 17,81 Simulación con CFD 60 60 priorizado el tratamiento de un mayor caudal de agua, manteniendo el mismo tamaño de gota que en el caso base (4·10-5 m) Se han introducido dos líneas de inyección con 30 entradas de gotas cada una, una inyección en y=2,5 m y otra en y=3,15 m, dejando una distancia prudente desde las inyecciones hasta la salida para aumentar el tiempo de residencia y que se produzca una evaporación eficaz del agua. Las distribuciones de temperatura y otros parámetros son similares al caso 2. En la Figura 63 se puede visualizar el comportamiento de la salmuera en este caso. Este caso constituye, para un diámetro de partículas de 40 micras, la variación más eficaz del caso base en este proyecto, consiguiendo saturar la corriente de agua hasta casi un 39% en la sección de salida (en el caso base la humedad relativa en el outlet era del 0,076%) y tratando un 50% más de salmuera, todo ello estableciendo un set-up menos simple, teniendo en cuenta unas condiciones de entrada de la salmuera y del aire de entrada más cercanas a las que podrían darse en la aplicación práctica final de este sistema. Figura 63. Diámetro de partículas en el caso 3 Finalmente, la Tabla 12 muestra de manera resumida las condiciones de salida en el outlet del flujo en cada uno de los casos, mostrando nuevamente la temperatura media y la humedad relativa media en esa zona, con sus respectivos coeficientes de variación. Tabla 12. Parámetros en el outlet según el caso estudiado Caso Temperatura media outlet (K) Humedad relativa outlet (%) CV temperatura media outlet (%) CV humedad relativa outlet (%) 1 461,6 0,75 2,831 51,49 2 385,3 8,16 1,140 15,54 3 353,6 38,85 5,860 49,63 Análisis de sensibilidad 62 62 ANÁLISIS DE SENSIBILIDAD na malla eficaz es vital para cualquier simulación en CFD. Ya que ningún software permite, a día de hoy, un mallado completamente automático, el desarrollo de una malla óptima suele considerarse el cuello de botella de cualquier proyecto, pues si esta no es capaz de recoger todas las características y singularidades del comportamiento de un sistema, la simulación se aleja de la realidad y el propio uso del CFD pierde el sentido. Aunque en el caso base o de referencia se ha propuesto una malla lo suficientemente eficaz como para que cumpla el propósito impuesto, una forma de evaluar la veracidad de los resultados y la fiabilidad de la solución es haciendo un análisis de sensibilidad sobre la malla. Aunque existen métodos para generar una malla correctamente e índices que nos proporcionan indicaciones sobre la fiabilidad de la misma, en este proyecto se hace un análisis de sensibilidad. Este consiste en simular el caso base sobre varias mallas de distinta complejidad, analizando algunos de los parámetros que se obtienen como resultado (a los cuales denominaremos parámetros índice) y observando si cambian de una malla a otra. Lo que se persigue es diseñar una malla lo suficientemente precisa como para que los resultados de la simulación no varíen si se aumenta la densidad de esta, lo que significaría que se habría conseguido la densidad mínima de celdas como para obtener unos resultados fiables. Los parámetros índice sobre los que se evalúa la eficacia de la malla son los siguientes:  Temperatura media a la salida del canal vertical  Humedad relativa media a la salida del canal vertical  Energía cinética turbulenta en y=5m Para afinar la malla no se han hecho grandes cambios con respecto a la estructura del mallado ni a los valores de los parámetros establecidos en el subapartado 4.2, sino únicamente ligeras variaciones de las mismas técnicas empleadas en el caso base, dando lugar a celdas más pequeñas y un mallado más fino para conseguir mayor precisión en las medidas. Los cambios efectuados para ello han sido: disminución del tamaño de los elementos para aumentar el número de celdas en todo el dominio; y mayor densidad en las zonas cercanas a la pared con la función Inflation, pues estas zonas son las que incluyen grandes gradientes en los distintos parámetros intervinientes y las que exigen una mayor precisión. En este proyecto el análisis se centra en cuatro casos numerados, correspondiendo el caso 1 al caso base y aumentando la densidad de la malla hasta el caso 4. Cabe destacar que se ha elegido como caso más complejo (caso 4) aquel que se encuentra en el límite del poder computacional del que se dispone en la realización de este proyecto para un tiempo de cálculo asumible, por lo que este apartado es susceptible de ser ampliado y mejorado con nuevos casos de mayor complejidad en trabajos posteriores. Los datos de los distintos casos se muestran en la Tabla 13, a partir de la cual se presenta una gráfica del número de elementos de los distintos casos para mostrar de forma más visual la diferencia entre ellos (Figura 64). Por otro lado, en la U Figura 64. Gráfica del número de elementos de la malla en los distintos casos 300000 350000 400000 450000 500000 550000 600000 1 2 3 4 Número de elementos en los distintos casos Análisis de sensibilidad 64 64 Tabla 14 se presentan los resultados de los parámetros índice anteriormente mencionados tras los cálculos en los diferentes casos. Aunque ésta refleja de manera precisa los cambios que se han producido en los parámetros, ya que la finalidad de este apartado es comparar la diferencia de los distintos casos con el de referencia, quizá sea más interesante aún conocer la diferencia relativa de estos resultados con respecto al caso base, lo cual puede comprobarse en la Tabla 15, donde se refleja la variación relativa de estos nuevos resultados en relación al caso base, con su posterior representación gráfica en la Figura 65. Además, para apreciar estas variaciones desde otro punto de vista, la Tabla 16Figura 66 expresa la variación relativa al caso anterior de los parámetros índice, mostrando cuánto cambian estos parámetros de un caso a otro, también con su posterior representación gráfica en la Figura 66. Tabla 13. Casos de afino del mallado Caso Opciones inflation codo Opciones inflation canal vertical Tamaño de elementos (m) Número de elementos 1 Number of layers: 40 Growth Rate: 1,2 Maximum Thickness: 0,007m Number of layers: 25 Growth Rate: 1 Maximum Thickness: 0,1m 0,08 355760 2 Number of layers: 50 Growth Rate: 1,05 Maximum Thickness: 0,1m Number of layers: 25 Growth Rate: 1 Maximum Thickness: 0,1m 0,007 471603 3 Number of layers: 65 Growth Rate: 1,04 Maximum Thickness: 0,1m Number of layers: 25 Growth Rate: 1 Maximum Thickness: 0,1m 0,007 515658 4 Number of layers: 65 Growth Rate: 1,04 Maximum Thickness: 0,1m Number of layers: 55 Growth Rate: 1,04 Maximum Thickness: 0,1m 0,065 586528 5 Number o layers: 65 Growth Rate: 1,04 Maximum Thickness: 0,1m Number of layers: 55 Growth Rate: 1,04 Maximum Thickness: 0,1m 0,005 950220 Figura 64. Gráfica del número de elementos de la malla en los distintos casos 300000 350000 400000 450000 500000 550000 600000 1 2 3 4 Número de elementos en los distintos casos Análisis de sensibilidad 66 66 Tabla 14. Tabla de variación de parámetros índice según el caso Caso Energía cinética turbulenta (m²/s²) Temperatura media de salida (K) Humedad relativa media de salida (%) Rango de 𝑦+ 1 0,6066 521,9 0,08 12-68 2 0,6066 521,9 0,08 2,5-37 3 0,6303 521,9 0,08 2-37 4 0,7507 521,8 0,08 2-37 Tabla 15. Tabla de variaciones relativas al caso inicial de parámetros índice en cada caso expresados en % Caso Energía cinética turbulenta Temperatura media de salida Humedad relativa media de salida 1 0 0 0 2 0 0,001916 0,0001308 3 3,76 0,003832 0,0007847 4 22,86 0,01341 0,003008 0 5 10 15 20 25 1 2 3 4 Variación de la energía cinética turbulenta (%) 0 0,002 0,004 0,006 0,008 0,01 0,012 0,014 0,016 1 2 3 4 Variación de la temperatura media en outlet (%) Figura 65. Gráficas de variación de parámetros índice según el caso Tabla 16. Tabla de variaciones relativas al caso anterior de parámetros índice en cada caso Caso Energía cinética turbulenta Temperatura media de salida Humedad relativa media de salida 1 0 0 0 2 0 0,001916 0,01308 3 3,097 0,001916 0,06539 4 22,86 0,009581 0,2222 0 0,0005 0,001 0,0015 0,002 0,0025 0,003 0,0035 1234 Variación de la humedad relativa media en outlet (%) 0 5 10 15 20 25 2 3 4 Variación de la energía cinética turbulenta respecto al caso anterior (%) 0 0,002 0,004 0,006 0,008 0,01 0,012 234 Variación de la temperatura media en outlet (%) respecto al caso anterior Análisis de sensibilidad 68 68 Figura 66. Gráficas de variación relativa al caso anterior de los parámetros índice según el caso Observando todos estos datos podemos obtener varias conclusiones: (i) La malla inicial es susceptible de mejora, aunque se puede considerar de una calidad aceptable, pues las variaciones relativas de los distintos parámetros no son desorbitados, siendo <0,2% en el caso de la temperatura y la humedad (ii) Pese a que el aumento en el número de elementos se aumenta de forma casi lineal a lo largo de los casos (Figura 64), las pendientes de las variaciones de los distintos parámetros cada vez es mayor, por lo que los aumentos en la complejidad de la malla implican un aumento de la precisión cada vez mayor (iii) El rango del parámetro 𝑦+ disminuye con respecto al caso base y se acerca más al rango de valores óptimo (más detalles sobre el parámetro 𝑦+ la subsección 4.2), por lo que los resultados cercanos a las paredes son más precisos, justo lo que se perseguía con la función Inflation. Como ya se ha mencionado anteriormente, este es el mayor nivel de precisión que se puede alcanzar con los medios disponibles para este proyecto, y durante la realización de este proyecto no se ha alcanzado el caso ideal que se buscaba, pues los resultados de las simulaciones sucesivas continúan alejándose de los resultados del caso base. Sin embargo, este trabajo (y más concretamente este subapartado) puede ser de gran importancia de cara a una ampliación de las ideas aquí expuestas y la continuación de su desarrollo. 0 0,05 0,1 0,15 0,2 0,25 2 3 4 Variación de la humedad relativa media en outlet (%) respecto al caso anterior CONCLUSIONES e ha realizado un estudio teórico preliminar sobre un sistema que logra la nula descarga líquida (ZLD) para el tratamiento de concentrados de procesos de ósmosis inversa. Este estudio es un primer paso de las investigaciones teóricas y experimentales del proyecto SOLMIDEF en lo que respecta al desarrollo de este tipo de sistemas. Con este trabajo se ha establecido una metodología válida para el diseño preliminar de los sistemas ZLD, que podrá aplicarse a diversas configuraciones y geometrías, así como a nuevos parámetros de diseño y operación con la finalidad de adaptarlo a la aplicación concreta con microturbinas de gas solares. Las principales aplicaciones contempledas en el projecto SOLMIDEF son la desalación de agua de mar o salobre y el tratamiento de aguas industriales. Se ha demostrado que las variables que más influyen sobre el comportamiento del sistema son el tamaño de las gotas y el caudal de salmuera introducido, variables que se deben tener especialmente en cuenta en futuros desarrollos de este proyecto. En esta primera aproximación se ha conseguido evaporar un máximo de 0,36 kg/s de agua (lo que corresponde a 0,053 kg de salmuer por kg de aire) de una salmuera de concentración de 46183 ppm (condiciones detalladas en el apartado 4.5), lo que demuestra la validez de esta tecnología innovadora de ZLD, la cual es susceptible de ser mejorada. Existe un amplio abanico de mejoras aplicables a este proyecto. Los nuevos estudios que se realicen sobre esta tecnología podrían enfocarse en la introducción de gotas desde distintos puntos, empleando distintas geometrías del regenerador o introduciendo álabes en el codo que favorezcan el cambio de dirección que sufre el gas al pasar a través de él. Estos cambios de configuración podrían permitir el tratamiento de un mayor caudal de salmuera o la introducción de gotas con un mayor tamaño, permitiendo hacer estudios en los que se concrete y desarrolle la tecnología empleada en la introducción de la fase discreta. En cualquier caso, este estudio contribuye al desarrollo del proyecto SOLMIDEFF empleando una tecnología prometedora. S