Estudio numérico sobre sedimentación en el lecho de ríos

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería Civil Estudio numérico sobre sedimentación en el lecho de ríos Autor: Francisco Javier Moreno Hidalgo Tutor: Javier Dávila Martín Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 iii Trabajo Fin de Grado Grado en Ingeniería Civil Intensificación en Hidrología Estudio numérico sobre sedimentación en el lecho de ríos Autor: Francisco Javier Moreno Hidalgo Tutor: Javier Dávila Martín Profesor titular Dep. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 v Trabajo de Fin de Grado: Estudio numérico sobre sedimentación en el lecho de ríos Autor: Francisco Javier Moreno Hidalgo Tutor: Javier Dávila Martín El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2017 El Secretario del Tribunal vii A mi familia y Belén A mis profesores y compañeros de promoción ix Agradecimientos Agradecer, en primer lugar, a mi familia y a Belén, por todo el apoyo recibido durante estos cuatro intensos años. También a mis amigos y compañeros de promoción por haber hecho más ameno el trabajo diario. Y, por supuesto, a Javier Dávila por darme la oportunidad de introducirme en este trabajo y línea de investigación, y a Alonso Fernández, por su ayuda y aporte de conocimientos en Ansys Fluent y modelos multifásicos. Francisco Javier Moreno Hidalgo Sevilla, 2017 ÍNDICE DE TABLAS Tabla 2-1. Método isoyetas 5 Tabla 3-1. Resumen estudio de mercado de agitadores 23 Tabla 5-1. Características de la sección del canal 37 Tabla 5-2. Propiedades de la fase líquida 46 Tabla 5-3. Propiedades de la fase sólida 46 Tabla 5-4. Características fase sólida en modelo UDF 58 Tabla 5-5. Parámetros condición de contorno entrada 61 Tabla 5-6. Parámetros condición de contorno salida 62 Tabla 5-7. Parámetros condición de contorno lecho 62 Tabla 5-8. Agitadores empleados 66 Tabla 6-1. Valores de rugosidad 69 Tabla 6-2. Pérdida de carga 69 Tabla 6-3. Partículas consideradas en estudio previo 70 Tabla 6-4. Posiciones probadas de los agitadores 76 Tabla 6-5. Cálculo de la efectividad de los agitadores probados 81 Tabla 6-6. Rendimiento agitador con tobera 89 xvii ÍNDICE DE FIGURAS Ilustración 2-1. Ubicación de la zona de estudio en Andalucía 3 Ilustración 2-2. Cuenca del río Arenoso 4 Ilustración 2-3. Ubicación de la cuenca respecto al resto de ríos de Andalucía 4 Ilustración 2-4. Mapa isoyetas de la cuenca del río Arenoso 5 Ilustración 2-5. Presa y embalse del Arenoso 6 Ilustración 2-6. Evolución del volumen embalsado en los últimos años 7 Ilustración 2-7. Vista en planta de los elementos de la zona de estudio 8 Ilustración 2-8. Sección del canal desde p.k. 0+000 a p.k. 0+340 8 Ilustración 2-9. Sección del canal dede p.k. 0+340 hasta el final. 9 Ilustración 2-10. Cotas para un caudal de 10 m3/s 9 Ilustración 2-11. Materiales de la cuenca del río Guadalquivir 10 Ilustración 2-12. Nombre de los materiales de la Ilustración 2-11 11 Ilustración 2-13. Estaciones Automáticas de Alerta 12 Ilustración 2-14. Ficha de la estación 507 13 Ilustración 2-15. Ficha de la estación 512 14 Ilustración 2-16. Precipitación anual en estación meteorológica cercana a EAA 507 15 Ilustración 2-17. Precipitación anual en estación meteorológica cercana a EAA 512 15 Ilustración 2-18. Medidas EAA 507 meses octubre a diciembre 16 Ilustración 2-19. Precipitación anual en estación meteorológica cercana a EAA 518 17 Ilustración 2-20. Medidas EAA 518 meses julio a septiembre 17 Ilustración 2-21. Medidas EAA 518 meses octubre a diciembre 18 Ilustración 3-1. Aceleradores de corriente de distintos tamaños 20 Ilustración 3-2. Colocación de agitador 22 Ilustración 3-3. Colocación de agitadores 22 Ilustración 3-4. Alzado, planta y perfil de agitador de 900 mm de diámetro 24 Ilustración 3-5. Alzado, planta y perfil de agitador de 1200 mm de diámetro 24 Ilustración 3-6. Curva H-Q bomba de flujo axial 25 Ilustración 3-7. Curva H-Q bomba de flujo mixto 25 Ilustración 3-8. Curva H-Q bomba de flujo radial 26 Ilustración 3-9. Simplificaciones consideradas en curva H-Q 26 Ilustración 3-10. Tipo de curva H-Q empleada para caracterizar agitadores 28 Ilustración 4-1. Transporte de sedimento en un río 29 Ilustración 4-2. Ábaco de Shields 31 Ilustración 4-3. Movimiento turbulento liso 32 Ilustración 4-4. Movimiento turbulento rugoso 32 Ilustración 5-1. Volumen total discretizado en volúmenes de control 36 Ilustración 5-2. Geometría del canal 36 Ilustración 5-3. Elementos que forman la geometría en estudio 37 Ilustración 5-4. Sección transversal del canal 38 Ilustración 5-5. Nomenclatura de las caras del lecho 38 Ilustración 5-6. Nomenclatura de la entrada y la salida 39 Ilustración 5-7. Nomenclatura de la superficie libre 39 Ilustración 5-8. Tipos de celda en mallas 2D 40 Ilustración 5-9. Tipos de celda en mallas 3D 40 Ilustración 5-10. Componentes de una malla 2D y 3D 40 Ilustración 5-11. Malla estructurada cuadrilátera 41 Ilustración 5-12. Malla no estructurada cuadrilátera 41 Ilustración 5-13. Malla híbrida 42 Ilustración 5-14. Representación de los vectores que aparecen en ecuaciones 5-2 y 5-3 43 Ilustración 5-15. Forma de calcular la relación de aspect 43 Ilustración 5-16. Malla estructurada 44 Ilustración 5-17. Malla híbrida generada 44 Ilustración 5-18. Malla conforme 45 Ilustración 5-19. Malla en la zona de los agitadores 45 Ilustración 5-20. Malla en la zona central del canal 46 Ilustración 5-21. Definición de la gravedad y tipo de análisis 47 Ilustración 5-22. Subdivisiones zona cercana a la pared 49 Ilustración 5-23. Elección del modelo de turbulencia 50 Ilustración 5-24. Interpretación del código de la UDF en Ansys Fluent 54 Ilustración 5-25. Activación UDS 55 Ilustración 5-26. Activación MACRO para cálculo de tensión en el lecho 55 Ilustración 5-27. Modelos reológicos de viscosidad 56 Ilustración 5-28. Propiedades del agua 57 Ilustración 5-29. Términos fuente 57 Ilustración 5-30. Condiciones de contorno entrada y salida 58 Ilustración 5-31. Condición de contorno sobre el lecho 58 Ilustración 5-32. Caracterización de las inyecciones 59 Ilustración 5-33. Parámetros de seguimiento de la trayectoria de las partículas 60 Ilustración 5-34. Condición de contorno entrada y salida 62 Ilustración 5-35. Condición de contorno 63 Ilustración 5-36. Representación de los parámetros geométricos del agitador 64 Ilustración 5-37. Diseño de la geometría de los agitadores 64 Ilustración 5-38. Parámetros de los agitadores 65 xix Ilustración 5-39. Nomenclatura de las caras del agitador 66 Ilustración 5-40. Curva H-Q para agitador de 0.9 m de diámetro 66 Ilustración 6-1. Diagrama de Moody 69 Ilustración 6-2. Trayectoria de partículas tipo 1 70 Ilustración 6-3. Trayectoria partículas tipo 2 71 Ilustración 6-4. Trayectoria partículas tipo 3 71 Ilustración 6-5. Trayectoria partículas tipo 4 71 Ilustración 6-6. Trayectoria partículas tipo 1 72 Ilustración 6-7. Trayectoria partículas tipo 2 72 Ilustración 6-8. Residuos escalados 73 Ilustración 6-9. Presión en la entrada del canal 73 Ilustración 6-10. Valor de la presión en la entrada 73 Ilustración 6-11. Velocidad en la salida 74 Ilustración 6-12. Caudal en la salida 74 Ilustración 6-13. Balance de caudales 74 Ilustración 6-14. Valores de y+ 75 Ilustración 6-15. Concentración de sedimento en el canal sin agitadores 75 Ilustración 6-16. Trayectoria partículas 1 y 2, prueba 1, agitador ϕ=0.9 m 77 Ilustración 6-17. Trayectoria partículas 1 y 2, prueba 1, agitador ϕ=1.2 m 77 Ilustración 6-18. Trayectoria partículas 1 y 2, prueba 2, agitador ϕ=1.2 m 78 Ilustración 6-19. Trayectoria partículas 1 y 2, prueba 3, agitador ϕ=1.2 m 78 Ilustración 6-20. Trayectoria partículas 1 y 2, prueba 4, agitador ϕ=1.2 m 78 Ilustración 6-21. Trayectoria partículas 1 y 2, prueba 5, agitador ϕ=1.2 m 79 Ilustración 6-22. Trayectoria partículas 1 y 2, prueba 6, agitador ϕ=1.2 m 79 Ilustración 6-23. Trayectoria partículas 1 y 2, prueba 7, agitador ϕ=1.2 m 80 Ilustración 6-24. Trayectoria partículas 1 y 2, prueba 8, agitador ϕ=1.2 m 80 Ilustración 6-25. Rendimiento en función de la separación longitudinal 81 Ilustración 6-26. Rendimiento en función de la separación respecto al lecho 82 Ilustración 6-27. Rendimiento en función de la separación transversal 82 Ilustración 6-28. Contorno de velocidades en plano y = 0.1 metros 83 Ilustración 6-29. Contorno de velocidades en plano y = 0.5 metros 83 Ilustración 6-30. Contorno de velocidades en plano y = 1.2 metros 84 Ilustración 6-31. Contorno de velocidades en el plano y = 0.5 metros 84 Ilustración 6-32. Trayectoria partículas en canal completo 85 Ilustración 6-33. Velocidad generada por agitadores 85 Ilustración 6-34. Velocidad según eje y en plano z = 0 86 Ilustración 6-35. Concentración de sólidos en canal completo con agitadores 86 Ilustración 6-36. Concentración de sólidos en canal completo con agitadores (distinta escala) 87 Ilustración 6-37. Trayectoria de las partículas para Q=6 m3/s 87 Ilustración 6-38. Simulación mediante UDF para Q = 6 m3/s 88 Ilustración 6-39. Trayectoria de las partículas para Q=3 m3/s 88 Ilustración 6-40. Geometría de las toberas 89 Ilustración 6-41. Trayectoria partículas 1 y 2 con uso de toberas 89 1 1 INTRODUCCIÓN La presencia de material sólido en los ríos es un fenómeno importante, el cual debería suponer un gran interés en el estudio previo a cualquier actuación sobre éstos. Este material genera deterioro y supone un problema para las instalaciones que se encuentra a su paso, por ello un estudio detallado de este fenómeno en el lugar de actuación puede suponer ventajas en el mantenimiento de la instalación durante su vida útil. Son abundantes los casos en la actualidad que se ven afectados por este problema, sirviendo de ejemplo el caso de muchos embalses en el mundo, los cuales se están viendo colmatados al no haber considerado la posibilidad de presencia de material sólido o cómo disminuir su efecto. En el presente estudio, el sedimento desempeña un efecto similar a lo ocurrido en los embalses, perteneciendo éste a una línea de investigación iniciada en 2016 por Esperanza Navarro en su Trabajo de Fin de Máster para la Universidad de Sevilla [1], realizado en el mismo departamento y bajo la misma tutela. El objetivo de la línea de investigación iniciada por el trabajo mencionado es buscar una solución al problema que dificulta el correcto funcionamiento del canal de bombeo de la presa del río Arenoso, situado en el municipio de Montoro, Córdoba. El origen de estos estudios es resolver la adversidad que genera la decantación de lodos en las instalaciones de bombeo y turbinado del canal mencionado. El procedimiento empleado para la elaboración de este estudio ha sido, en primer lugar, comprobar las conclusiones establecidas en el trabajo que inicia la línea a seguir para, a continuación, incluir algunas mejoras y profundizar principalmente en una de las soluciones propuestas al problema inicial, la cual se trata del uso de agitadores sumergibles, también conocidos como aceleradores de corriente, para generar zonas controladas de depósito de sedimentos. En esta última idea se refleja el objetivo principal de este estudio concreto, optimizar el uso de agitadores para resolver el problema que se plantea. Mientras que el trabajo inicial tenía como objetivo la descripción de la problemática, de la cual se incluirán aquí las nociones básicas para establecer la situación a estudiar, y comprobar si el uso de agitadores como solución podía considerarse efectivo, en este caso se busca profundizar en este aspecto, ya que se concluyó que sí eran efectivos para disminuir el problema sin considerar cómo se podría optimizar la solución variando distintos parámetros que caracterizan a estos elementos. Para llevar a cabo el estudio, es necesario apoyarse de un software comercial de cálculo de dinámica computacional (CFD, en inglés Computational Fluid Dynamics), Ansys Fluent 16. El grueso del trabajo se encuentra en las simulaciones realizadas mediante este programa del comportamiento de las partículas en el canal de estudio. 1.1. Problemática a resolver La línea de investigación a la que pertenece el trabajo pretende resolver el problema de la decantación de lodos en el canal del Arenoso, el cual se emplea como elemento de trasvase de agua entre la presa del río Arenoso y el río Guadalquivir. Para describir la problemática es necesario incluir brevemente cuál es el funcionamiento del canal de estudio, ya que posteriormente se describirá con mayor detalle. El embalse del río Arenoso, además de ser un elemento de regulación, pretende paliar la necesidad de agua en época seca que requiere el río Guadalquivir, mientras que, en época húmeda recoge agua sobrante del Guadalquivir. Para realizar ambas funciones consta de equipos de bombeo (para trasvasar agua del Guadalquivir al embalse) y de equipos de turbinado (para aprovechar energéticamente el trasvase de agua desde el embalse al Guadalquivir). La funcionalidad de este canal se ve afectada por la presencia de sólidos, pudiendo destacar los siguientes problemas: Introducción 2 - Pérdida de sección del canal dificultando la llegada de agua a la zona de aspiración de las bombas. - Acumulación de sólidos en las cántaras de bombeo. - Obstrucción de los elementos (rejas y compuertas) de las instalaciones previas a las bombas. 1.2. Antecedentes Se considera necesario recoger las conclusiones establecidas en el trabajo anterior, ya que serán el punto de partida de este trabajo. A pesar de la necesidad de dragado del canal completo para eliminar los sólidos actualmente depositados es necesario llegar a alguna solución que garantice el correcto funcionamiento de la instalación durante su uso corriente sin llegar al extremo de la problemática anteriormente descrita. Entre las alternativas consideradas se establece la colocación de agitadores como aquella a profundizar en el estudio, por ser de bajo consumo energético y bajo impacto ambiental. En el estudio anterior se llegaron a las siguientes conclusiones en cuanto al empleo de agitadores como posible solución: - La colocación de agitadores es un método efectivo para reducir la concentración de partículas que llegan a las cántaras de aspiración. Aunque la entrada de sólidos al canal se sigue produciendo, se tiene una zona de sedimentación controlada. - La colocación de agitadores es un método efectivo para partículas de tamaño entre 0.01 mm y 0.05 mm y densidades en torno a 1500 – 1900 kg/m3. Éstas fueron las principales conclusiones establecidas, las cuales se comprobarán en el apartado correspondiente a la simulación en Ansys Fluent. Varias cuestiones quedaron abiertas en el estudio previo: - Valor de la fracción volumétrica de sólidos en el canal. Se tomó un valor altamente conservador, siendo éste el límite superior que presenta el modelo elegido para la simulación (10%). Teniendo en este caso como objetivo buscar un valor más real. - Estudiar el problema teniendo en cuenta que existen en el mercado diversos tipos de agitadores. - Estudio del funcionamiento de los agitadores con caudales de aguas bajas. 1.3. Objetivo del estudio Aunque ya se han recogido algunas ideas acerca del objetivo del presente trabajo, se recogen de nuevo con mayor detalle. Una vez conocido que la colocación de agitadores es efectiva en el caso a estudiar, se hace necesario optimizar el proceso. Para llevar a cabo esto, se ha procedido a un estudio del mercado de agitadores consultando catálogos de comerciantes que se dedican a su fabricación para caracterizar sus principales parámetros. La optimización se ha realizado desde el punto de vista de consumo energético, dimensiones y colocación. 3 3 Estudio numérico sobre sedimentación en el lecho de ríos 2 DESCRIPCIÓN DEL CASO Para la realización del estudio es necesario conocer el caso y los elementos que se involucran en él. 2.1. Localización Cerca de la localidad de Montoro (Córdoba) se encuentra el embalse del Arenoso (situado en el río Arenoso), el cual tiene como funciones: - Regulación de la cuenca propia. - Laminación de las avenidas. - Regulación de la cuenca del Guadalquivir. Asociado a este último objetivo del embalse surge el concreto problema a estudiar. La regulación de la cuenca del Guadalquivir se realiza en dos etapas, según la época del año: - En la época húmeda recibe agua mediante bombeo del río Guadalquivir a través del canal del Arenoso. - En la época seca aporta agua al Guadalquivir, incluyendo un aprovechamiento hidroeléctrico mediante el turbinado del caudal desembalsado. La función de esta regulación del Guadalquivir es almacenar agua sobrante durante el invierno para poder emplearla en uso agrícola durante la época de escasez de agua de lluvia. La conexión entre el embalse y el río Guadalquivir se hace a través del canal del Arenoso. Las bajas velocidades del agua en este canal provocan la deposición del material transportado en el lecho, lo cual produce una pérdida de sección del canal repercutiendo en la imposibilidad, en algunas ocasiones, de llevar a cabo el bombeo necesario desde el canal al embalse. Esta deposición de material ocurre principalmente en la época de lluvia, momento en el que el Guadalquivir va más cargado de partículas en suspensión. Ilustración 2-1. Ubicación de la zona de estudio en Andalucía Descripción del caso 4 2.2. Río Arenoso El río Arenoso nace en el municipio cordobés de Cardeña hasta su desembocadura en el río Guadalquivir, a cuya cuenca pertenece. Tiene una longitud de 40 kilómetros. Su cuenca es de 404 km2, sobre la cual cae una precipitación media anual de 740 mm (dato corroborado con método de isoyetas). Ilustración 2-2. Cuenca del río Arenoso Ilustración 2-3. Ubicación de la cuenca respecto al resto de ríos de Andalucía Aunque no es objetivo del presente estudio, se ha comprobado el dato de la precipitación media anual de la cuenca mediante el método de las isoyetas. 5 5 Estudio numérico sobre sedimentación en el lecho de ríos Tabla 2-1. Método isoyetas Precipitación (mm) Área (km2) 600 700 900 1000 1200 37.02 290.79 54.13 19.28 2.80 Este método calcula la precipitación media de una zona a partir de un plano de isoyetas (Ilustración 2-4) y usando la expresión 2-1. 𝑃= ∑𝑃𝑖𝐴𝑖 𝑛 𝑖=1 𝐴𝑇 (2-1) Siendo i el subíndice que representa cada porción de precipitación Pi y de superficie Ai. Obteniendo una precipitación media anual en la cuenca de aproximadamente 740 mm/año. Ilustración 2-4. Mapa isoyetas de la cuenca del río Arenoso 2.3. Embalse y presa del Arenoso La presa del Arenoso se encuentra sobre el río Arenoso, aproximadamente a 1.2 km de su desembocadura en el río Guadalquivir, siendo esta distancia la del canal de estudio, la cual se verá más adelante. Es una presa de reciente construcción, año 2007, con una altura máxima de 80 metros, longitud de 1481.28 metros y anchura de coronación de 11.30 metros. Se trata de una presa de tipo gravedad de materiales sueltos, cuyo núcleo central está formado de arcillas procedentes de propio vaso del embalse y protegido mediante escollera. Otras características de la presa son: - Nivel de coronación: 215 m.s.n.m. Descripción del caso 12 Ilustración 2-13. Estaciones Automáticas de Alerta Fuente: Confederación Hidrográfica del Guadalquivir La ubicación de la zona en estudio se encuenta en el curso medio – bajo del Guadalquivir, entre las estaciones 507 y 512 mostradas en la Ilustración 2-13, por lo que serán sobre las que se centre el estudio de sólidos. 13 13 Estudio numérico sobre sedimentación en el lecho de ríos Ilustración 2-14. Ficha de la estación 507 Fuente: Confederación Hidrográfica del Guadalquivir Descripción del caso 14 Ilustración 2-15. Ficha de la estación 512 Fuente: Confederación Hidrográfica del Guadalquivir De todas las medidas realizadas por las EAA interesan las relacionadas con los sólidos: turbidez y sólidos en suspensión. La medida de la tubidez se realiza en NTU (unidad nefelométrica de turbidez), la cual puede ser empleada para obtener una aproximación de los sólidos totales en suspensión. Además de a través de la turbidez, existe una medida directa de los sólidos en suspensión en algunas estaciones. La medida de sólidos en suspensión se realiza mediante un equipo con método de medición de luz difusa con ángulo dual dispersa a 90º. Este equipo permite medir turbidez, según el método EN ISO 7027 y sólidos en suspensión según DIN 37414 en un rango de 0.001 a 50 g/L. El problema que presenta la turbidez es que solamente permite medir valores inferiores a 500 NTU, por lo que su medida será representativa de situaciones que no se corresponden con avenidas, lo cual puede dar una idea de los sólidos en condiciones normales de funcionamiento. Por otro lado, la medida de sólidos no presenta este límite, sin embargo, solamente existen medidas del año 2011. Aún así, según un estudio de evolución de turbidez en el Guadalquivir realizado por la Confederación Hidrográfica del Guadalquivir [20], el periodo de 2009 a 2011, desde el punto de vista hidrológico, se 15 15 Estudio numérico sobre sedimentación en el lecho de ríos considera húmedo, siendo el periodo más lluvioso de los últimos treinta años, teniendo una precipitación total registrada de 630 mm frente a los 550 mm de la media. Esta última idea aumenta la valoración sobre el empleo del año 2011 como referencia en cuanto a sólidos, ya que la mayor parte de los sólidos en el canal son aportados durante las avenidas del río Guadalquivir. Establecidas las limitaciones que presenta este método para determinar los sólidos en suspension y la ventaja de tener el año 2011 como año húmedo, se pasa al análisis de las estaciones para obtener la fracción volumétrica de sólidos en la zona de estudio. Ilustración 2-16. Precipitación anual en estación meteorológica cercana a EAA 507 Fuente: Confederación Hidrográfica del Guadalquivir Ilustración 2-17. Precipitación anual en estación meteorológica cercana a EAA 512 Fuente: Confederación Hidrográfica del Guadalquivir Teniendo en cuenta la limitación que presentan los datos de turbidez se procede a analizar los datos de sólidos en suspensión del año 2011. Para la estación 507 los meses de octubre a diciembre, aunque no sean los más lluviosos, presentan datos significativos. Durante estos meses se recogen los datos mostrados en la Ilustración 2-18. Descripción del caso 16 Ilustración 2-18. Medidas EAA 507 meses octubre a diciembre Fuente: Confederación Hidrográfica del Guadalquivir Siendo el valor máximo próximo a 250 mg SS / L. El valor de sólidos en suspensión durante el funcionamiento normal se encuentra en torno a los 100 – 150 mg SS / L. Otros meses interesantes según se aprecia en la Ilustración 2-16 son: - Abril de 2008. - Diciembre de 2009. - Enero – Marzo 2010. - Diciembre 2010. En la mayoría de estos meses la EAA se encuentra inundada, impidiendo la obtención de datos. Por lo que puede establecerse que para la estación 507 un valor de 250 mg SS / L son representativos de episodios de lluvia, pero no extrema. En la estación 512, como se aprecia en la Ilustración 2-17, se tienen los siguientes meses interesantes: - Abril 2008. - Diciembre 2009. - Enero – Marzo 2010. - Diciembre 2010. La mayor parte de los días de los meses indicados la EAA ha permanecido inundada, lo cual provoca que los datos de esta estación sean poco representativos. Para no tener una única estación como referencia, se ha analizado también la EAA 518, que se encuentra aguas arriba de la 512, aunque esté alejada de la zona de estudio. Los datos de lluvia son similares a los representados anteriormente de las estaciones 507 y 512, por lo que las fechas representativas de avenidas se mantienen. Durante el mes de septiembre de 2011 se alcanza en esta estación un valor de 750 mg SS / L y durante los 17 17 Estudio numérico sobre sedimentación en el lecho de ríos meses de octubre a diciembre del mismo año se llega hasta los 500 mg SS / L. Ilustración 2-19. Precipitación anual en estación meteorológica cercana a EAA 518 Fuente: Confederación Hidrográfica del Guadalquivir Ilustración 2-20. Medidas EAA 518 meses julio a septiembre Fuente: Confederación Hidrográfica del Guadalquivir Descripción del caso 18 Ilustración 2-21. Medidas EAA 518 meses octubre a diciembre Fuente: Confederación Hidrográfica del Guadalquivir Teniendo en cuenta los datos analizados, se decide, para estar del lado de la seguridad, considerar en el canal una concentración de 1000 mg SS / L, ya que, como se ha podido comprobar, en algunos casos la concentración llega a valores cercanos a 750 mg SS / L. El siguiente paso será obtener la fracción volumétrica (V.F., del inglés volume fraction) de material que se tiene en el canal. 𝑉.𝐹.= 100𝑄𝑠 𝑄𝑇 (2-2) Siendo: - Qs: caudal sólido. - QT: caudal total. Teniendo en cuenta que el caudal de partida en el canal de estudio va a ser de 10 m3/s, se tiene una entrada de sedimentos de 10 kg/s. Como se verá más adelante, se van a considerar dos tipos de partículas (ρs = 1500 kg/m3 y ρs = 1900 kg/m3), asignándole a cada una de ellas 5 kg/s. De tal forma que se obtiene el siguiente caudal sólido: 𝑄𝑠=5(𝑘𝑔 𝑠)∗ 1 1500(𝑚3 𝑘𝑔)+5(𝑘𝑔 𝑠)∗ 1 1900(𝑚3 𝑘𝑔)≈6𝑒−3 (𝑚3 𝑠) (2-3) Con este valor de caudal sólido y un caudal total de 10 m3/s se obtiene una fracción volumétrica de 0.06%. 19 19 Estudio numérico sobre sedimentación en el lecho de ríos 3 DESCRIPCIÓN DE LA SOLUCIÓN A EMPLEAR La solución que se pretende implementar consiste en el uso de agitadores o aceleradores de corriente, cuya eficacia fue propuesta y verificada en el estudio anterior al presente. En el trabajo mencionado se nombran otras alternativas que se podrían haber considerado, sin embargo, el objetivo de esta línea de investigación tiende al empleo de agitadores, debido entre otras cosas a su efectividad y al bajo consumo energético. Aún así, a continuación, se incluyen otro tipo de medidas que se pueden emplear. 3.1. Métodos de extracción de sedimentos Dentro de los métodos de extracción de sedimentos podrían considerarse dos tipos: - Solución directa: métodos de extracción. - Solución indirecta: métodos para impedir o reducir la entrada de sólidos. 3.1.1 Soluciones directas Las soluciones directas consisten, como su propio nombre indica, en la extracción directa del material. Se tratan de medidas correctivas, puesto que corrigen el funcionamiento del sistema, pero no dan una solución definitiva. El aporte constante de sedimento por parte del río Guadalquivir volverá a aterrar el canal con el paso del tiempo. - Dragado. Se trata de un método de extracción de acción directa en el que se extrae el material depositado en el fondo del canal mediante medios mecánicos (dragas mecánicas: similar a excavación terrestre) o equipos hidráulicos (dragas de succión). - Remoción hidráulica de sedimentos. Esta solución consiste en emplear la fuerza erosiva del agua para limpiar el fondo del canal. Aprovechando la existencia de un embalse en uno de los extremos del canal, podría descargarse agua que permita el arrastre de partículas desde el canal al Guadalquivir. El problema que presenta es la pérdida de volumen almacenado en las descargas de limpieza. - Retiro de sedimentos con sifón. Gobernado por el efecto de vacío debido a la diferencia de niveles entre la toma de descarga y la extractora. Este sistema emplea la energía hidráulica disponible por diferencia de niveles, por lo que se necesita diferencia de cota considerable para que el agua recorra el sifón completo dejando los sólidos en la parte baja. 3.1.2 Soluciones indirectas La medida indirecta más eficaz sería ubicar las instalaciones en zonas donde el aporte de material sea muy bajo. Sin embargo, esto solo es posible en nuevos proyectos, y necesitan un estudio muy detallado de la cuenca en la que se va a situar. La presencia de la fase sólida en una corriente genera erosión (generada por las partículas), transporte y deposición. La solución de estos problemas podría estar en ir directamente a su origen, es decir, la erosión generada por otros factores sobre el material sólido para provocar el desprendimiento de partículas dando lugar a los sedimentos. Detener la erosión al completo es una tarea físicamente imposible, además de económicamente, pero sí que pueden llevarse a cabo medidas que lo frenen. En la actualidad, tanto en España como en otros países, la presencia de sedimentos está generando grandes problemas, por ejemplo, el aterramiento de algunos embalses o el daño de instalaciones de bombeo o turbinado, como es el caso del Descripción de la solución a emplear 20 presente estudio. Llevar a cabo medidas para intentar paliar este efecto en su origen puede suponer un gasto extra en la confección de los proyectos puesto que estos problemas varían en función de la situación (para una misma actuación, lo que en una zona genera una ventaja, en otra ocasiona un problema) y habría que hacer un estudio al máximo detalle, pero un alivio, incluso económico, desde el punto de vista de la conservación y mantenimiento de las instalaciones una vez se pone en marcha la actividad. Teóricamente existen actuaciones posibles a llevar a cabo para resolver el problema de la erosión. A continuación, se exponen brevemente, puesto que no es objetivo del estudio, algunas de las que se consideran que generan menor impacto ambiental: - Tratamiento vegetal. La presencia de vegetación sobre el terreno genera una protección del suelo que dificulta el arrastre y desprendimiento de material sólido. - Agricultura. Cultivar el campo siguiendo las líneas de nivel puede reducir la erosión generada por la escorrentía superficial. Estas medidas además suponen una protección para el terreno. Otra medida indirecta puede ser la coordinación con presas aguas arriba y aguas abajo del Guadalquivir en la zona de estudio. 3.2. Empleo de agitadores sumergibles Dentro de la clasificación anteriormente expuesta, este método se encontraría en un paso intermedio, puesto que el objetivo es crear una zona controlada de depósito de sedimentos que posteriormente hay que extraer de forma directa. Los agitadores sumergibles son elementos que aportan velocidad al flujo. Normalmente se emplean en plantas de tratamiento de aguas para remover el fluido y permitir la mezcla con reactivos añadidos. En este caso, el objetivo que persiguen es crear una recirculación del flujo que genere zonas de baja velocidad donde se favorezca la deposición de sedimentos. Para su modelado en Ansys Fluent ha sido necesario realizar un estudio de los catálogos de los principales fabricantes de estos dispositivos para caracterizarlos según sus parámetros de diseño. Ilustración 3-1. Aceleradores de corriente de distintos tamaños Fuente: fabricante Sulzer 21 21 Estudio numérico sobre sedimentación en el lecho de ríos Los agitadores pueden ser tanto verticales como horizontales, si nos referimos a su eje de giro, siendo interesantes para este caso los de eje horizontal. Dentro del grupo de agitadores de eje horizontal existen dos grandes grupos: - Agitadores de alta velocidad. Son de menor diámetro y generan en el flujo velocidades superiores a 3 m/s. Esta última característica conlleva mayor gasto energético. Estos agitadores no son objeto del presente estudio. - Agitadores de baja velocidad. Imprimen al fluido una velocidad máxima entre 1.5 y 3 m/s. Es el grupo en el que se centra el estudio. Sus diámetros los hacen competentes en la búsqueda de los objetivo de recirculación que se persigue. Un parámetro característico de un agitador es su eficiencia, la cual se define como: 𝐸𝑓𝑖𝑐𝑖𝑒𝑛𝑐𝑖𝑎= 𝐸𝑚𝑝𝑢𝑗𝑒 𝑃𝑜𝑡𝑒𝑛𝑐𝑖𝑎 𝑐𝑜𝑛𝑠𝑢𝑚𝑖𝑑𝑎 (3-1) Como se aprecia en la fórmula 3-1, un agitador será más eficiente mientras mayor cantidad de flujo mueva consumiendo la menor potencia. Esta fórmula se ha obtenido de los catálogos de fabricantes [11] y [12], los cuales se basan en la norma ISO 21630, específica para agitadores sumergibles. Sin embargo, para este trabajo no pasa de un concepto teórico, puesto que se buscará la eficiencia desde un punto de vista que se verá en el apartado correspondiente a los resultados. Realizado un estudio de varios fabricantes, se extraen los siguientes parámetros como característicos de su diseño: - Velocidad de giro. - Potencia máxima. - Diámetro. Estos tres parámetros, junto con su colocación en el canal serán los parámetros a considerar en el proceso. La colocación de los agitadores, para que se genere la recirculación deseada, ha de ser diagonal, habiendo que determinar: - Separación longitudinal entre ellos. - Altura de colocación. - Separación transversal. Algunas recomendaciones realizadas por los fabricantes en cuanto a su posición son [11] [12]: - Es recomendable que la altura de su eje sea igual a su diámetro (Ilustración 3-2). - La distancia máxima a la que llega el efecto de un agitador es de 2.5 veces la altura de lámina de agua. Estas recomendaciones serán objeto de revisión en las simulaciones a realizar. Descripción de la solución a emplear 28 Con las ecuaciones 3-7 y 3-8 se obtienen las constantes: 𝐵0=43𝐻𝑑 (3-9) 𝐵2=𝐻𝑑 3𝑄𝑑2 (3-10) Finalmente queda la siguiente curva altura – caudal: 𝐻 = 43𝐻𝑑−𝐻𝑑 3𝑄𝑑2𝑄2 (3-11) Donde los valores de Q oscilan entre 0.5Qd (o 0.25Qmax) y 1.5Qd (o 0.75Qmax). A continuación, se incluye una figura solamente para visualizar la forma de la curva, puesto que en el apartado correspondiente a la simulación en Ansys Fluent se incluirán las curvas de los agitadores empleados. Ilustración 3-10. Tipo de curva H-Q empleada para caracterizar agitadores 3.3. Solución esperada La respuesta al uso de agitadores esperada es la siguiente: - Recirculación del flujo en la zona entre agitadores. - Creación de zonas de remanso en las que se depositarán las partículas de foma que se controle la zona en la que se depositan. - Aumento de la concentración entre los agitadores, debido al hecho de la creación de zonas favorables para la deposición de partículas. 0.000 0.020 0.040 0.060 0.080 0.100 0.120 0.140 0.160 0.180 0.000 0.200 0.400 0.600 0.800 1.000 1.200 1.400 1.600 H (m) Q (m3/s) Curva H-Q 29 29 Estudio numérico sobre sedimentación en el lecho de ríos 4. TEORÍA DE SEDIMENTOS Antes de comenzar con la simulación del problema usando un software de cálculo es necesario plantear la teoría que hay detrás del problema. Con este apartado se pretende dar una visión teórica del posible origen de las partículas para que se genere la sedimentación y de la dinámica de pequeñas partículas en flujos no uniformes. 4.1. Origen de las partículas El origen del problema del presente estudio se encuentra en la presencia de partículas sólidas en el canal. Nos encontramos con 2 fases en el problema: - Fase líquida: agua. Caracterizada en el problema por su densidad, viscosidad y velocidad. Este último será uno de los responsables del movimiento de las partículas sólidas. A mayor velocidad, mayor posibilidad de movimiento del material sólido. - Fase sólida: material sólido. Sus características se dividen en las propias de cada partícula aislada y las del conjunto de partículas. Dentro de las primeras encontramos el peso específico, la forma, el tamaño y la velocidad de sedimentación entre otras. Dentro de las características de un grupo de partículas tenemos la distribución granulométrica y parámetros relacionados con esta como son la desviación típica granular y el coeficiente de graduación. En el problema a tratar se caracterizarán mediante su diámetro y densidad. La presencia de esta fase sólida junto con la fase líquida presenta consecuencias como pueden ser la erosión, transporte y deposición. La erosión es la responsable natural de la presencia de sedimentos en los ríos, ya sea ocasionada por la acción del viento o por la propia acción del agua. Además, existen fuentes artificiales como la explotación de minas y canteras. El transporte de sedimento en un río se clasifica según dos criterios: - Según su procedencia u origen: cuenca (también llamado material de lavado de la cuenca) o lecho. - Según el modo de transporte: en suspensión o de fondo. - Ilustración 4-1. Transporte de sedimento en un río Teoría de sedimentos 30 El material cuya procedencia es el lecho admite ambos tipos de transporte, mientras que el procedente de la cuenca solo admite el transporte en suspensión. Esto tiene que ver con el tamaño de las partículas. Un criterio para separar ambos orígenes es considerar que las partículas de tamaño inferior a 0.065 mm proceden de la cuenca y las mayores proceden del lecho [3]. La procedencia de las partículas será la que determine sus características, lo cual está relacionado con la zona por la que discurre el río. En cuanto al modo de transporte puede decirse que el material en suspensión es sostenido por la turbulencia del flujo, mientras que el material de fondo es arrastrado por la corriente rodando, deslizando o saltando. El corte de tamaño entre un modo de transporte y el otro está influenciado por la velocidad del flujo y la tensión que genere sobre el lecho. Por último, se presenta la deposición. Este fenómeno ocurre cuando la partícula no es capaz de seguir en movimiento junto al flujo. Suele ocurrir en zonas donde la velocidad de la corriente es baja provocando la acumulación de partículas en el lecho. Una consecuencia puede ser la pérdida de sección del cauce, aunque puede ser compensado por erosión en otra zona. En el caso concreto de este estudio esta deposición impide que las bombas tengan altura de lámina de agua suficiente para aspirar. 4.2. Teoría de umbral de movimiento Al tratar una fase granular que se mueve en una corriente de agua es necesario incluir el concepto de umbral de movimiento. El hecho de que un lecho granular soporte la circulación de una corriente de agua provoca que en algún momento este movimiento genere el desplazamiento de las partículas que se encuentran depositadas en el fondo. Esto quiere decir que el canal en estudio no presenta únicamente partículas que entran desde la confluencia con el Guadalquivir, sino que su lecho está formado por partículas que se han ido depositando y que, en función de las características de la corriente líquida puede entrar en movimiento. Conocer las condiciones en las que se produce el inicio del movimiento de una partícula que se encuentra sobre el lecho es el objetivo que persigue la teoría del umbral de movimiento. Añadiendo el concepto de umbral de movimiento se está incluyendo la posibilidad de resuspensión de partículas desde el lecho. Las variables que influyen en esta teoría son la tensión de corte o acción del agua sobre el lecho y la velocidad. Estas variables se conocen como tensión crítica y velocidad crítica en el punto exacto en el que se inicia el movimiento, punto entre el movimiento y el reposo. Analizando el problema teniendo en cuenta una única partícula, ésta ejerce una resistencia a las fuerzas que tratan de moverla, de tal forma que se puede establecer un equilibrio entre las fuerzas que tratan de desestabilizar y las que estabilizan. - Fuerza estabilizadora. Establece resistencia al movimiento. o Peso sumergido de la partícula. - Fuerzas desestabilizadoras. Tienden a mover la partícula. o Fuerzas de contacto. o Fuerzas de arrastre: turbulencia, viscosidad y flotación. Se define la tensión cortante adimensional (expresión 4-1) como un parámetro que relaciona la fuerza desestabilizadora con la estabilizadora. 𝜏= 𝜏0 (𝛾𝑝−𝛾𝑤)𝐷𝑝 (4-1) Donde los subíncides p y w representan a la partícula sólida y al agua respectivamente, y los demás términos son: - τ0: tensión cortante en el fondo. 31 31 Estudio numérico sobre sedimentación en el lecho de ríos - γ: peso específico. - D: diámetro. Además de la tensión cortante adimensional, esta teoría hace uso del número de Reynolds granular, Re*. 𝑅𝑒∗= 𝑣∗𝐷𝑝 𝜐 (4-2) Donde, además de lo ya definido: - v*: velocidad de corte. 𝑣∗= √𝜏0 𝜌𝑤 - υ: viscosidad cinemática. 𝜐= 𝜇𝑤 𝜌𝑤 En las expresiones anteriores, ρ y μ representan la densidad y viscosidad dinámica, respectivamente. La relación entre la expression 4-1 y la expression 4-2 se refleja en el ábaco de Shields (Ilustración 4-2). Ilustración 4-2. Ábaco de Shields Fuente: [3] La línea marcada en la anterior ilustración limita la zona de movimiento (zona superior) de la de reposo (zona inferior). El número de Reynolds granular refleja el grado de turbulencia mediante la relación entre el valor relativo de las fuerzas de inervia y las viscosas en el entorno de la partícula. Se distinguen tres zonas en función del número de Reynolds granular: - Re* < 5: movimiento turbulento liso. La subcapa límite granular cubre la altura del grano. No es habitual este caso. - 5 < Re* < 70: movimiento turbulento intermedio. Se conoce como zona de transición. - Re* > 70: movimiento turbulento rugoso. Cuando Re* alcanza el valor de 400, 𝜏 tiende a un valor constante de 0.056 sin depender de Re*. Teoría de sedimentos 32 Ilustración 4-3. Movimiento turbulento liso Fuente: [3] Ilustración 4-4. Movimiento turbulento rugoso Fuente: [3] En la Ilustración 4-3 e Ilustración 4-4, δ representa la subcapa límite. 4.3. Dinámica de pequeñas partículas en flujos no uniformes Para determinar la posición de una pequeña partícula que se mueve en un fluido se puede hacer uso de la segunda Ley de Newton (ecuación 4-3). ∑𝐹=𝑚𝑎 (4-3) Siendo las fuerzas que actúan sobre la partícula (representado en la ecuación 4-3 como Σ F): la gravedad, la resultante de las fuerzas de presión y la viscosidad que el fluido aplica sobre su superficie. Si se integra esta fuerza para el flujo alrededor de una esfera a bajos números de Reynolds se llega a (Maxey y Riley, Phys. Fluids 1983): (𝜌𝑝+𝐶𝑀𝜌𝑓)𝛺𝑑𝑣𝑝 𝑑𝑡=(𝜌𝑝−𝜌𝑓)𝛺𝑔−3𝜋𝜇𝑑𝑝(𝑣𝑝−𝑣𝑓)+𝜌𝑓(1+𝐶𝑀)𝛺𝐷𝑣𝑓 𝐷𝑡 (4-4) En la ecuación 4-4, los subíncides p y f se refieren a la partícula y al fluido respectivamente. Además: - ρ: densidad del fluido - μ: viscosidad dinámica del fluido - 𝛺= 𝜋𝑑𝑝 3 6 : volumen de la partícula - 𝑣: vector velocidad - 𝑑 𝑑𝑡: derivada temporal siguiendo a la partícula - 𝐷 𝐷𝑡: derivada temporal siguiendo al fluido La parte correspondiente a la izquierda de la ecuación 4-4 se refiere al producto de la masa por la aceleración de partícula al que se le añade un término debido a la inercia del fluido que se mueve alrededor de la partícula 33 33 Estudio numérico sobre sedimentación en el lecho de ríos y que se llama masa virtual (CM). El valor de CM es una constante, que para una esfera rígida tiene un valor de 0.5. A la derecha de la igualdad en la ecuación 4-4 se encuentran, por orden, los términos debidos a la fuerza de la gravedad y el empuje de Arquímedes, la resistencia de Stokes para una esfera (flujo uniforme a bajos números de Reynolds) y la aceleración del fluido por ser no uniforme o no estacionario. La validez de la expresión se reduce a ocasiones en las que el diámetro de la partícula es mucho menor que la escala característica más pequeña del flujo y la concentración de partículas es suficientemente pequeña como para poder despreciar su interacción con el fluido y entre las propias partículas. Como se ha visto anteriormente, se tiene una concentración sólida de 0.06%, inferior a 0.1%, por lo que dichas interacciones pueden ser despreciadas en los modelos. Conocido el campo fluido vf (x, t) puede obtenerse la trayectoria de una partícula integrando la ecuación 4-4 haciendo uso también de la relación 4-5. 𝑑𝑥𝑝 𝑑𝑡= 𝑣𝑝 (4-5) Donde xp representa la posición de la partícula. El problema habitual se encuentra en la obtención de una solución para el campo fluido, e incluso conociéndolo se llega a un sistema de ecuaciones diferenciales de gran complejidad. Esta complejidad está asociada a la no linealidad introducida por la función vf (x, t). Este es el motivo por el que se recurre a técnicas de simulación computacional para la resolución de este tipo de problemas. Adimensionalizando la expresión 4-4 se puede analizar algunos términos para comprobar su importancia en el sistema. La adimensionalización se va a hacer respecto de la longitud y velocidad características del flujo, L y U respectivamente. Se definen las siguientes variables adimensionales: 𝑉= 𝑣𝑝 𝑈 (4-6) 𝑢= 𝑣𝑓 𝑈 (4-7) 𝑋= 𝑥𝐿 (4-8) 𝜏= 𝑡 𝐿/𝑈 (4-9) Sustiyendo estas variables en la expresión 4-4 se obtiene: 𝑑𝑉 𝑑𝜏= 𝛽−1 𝛽+𝐶𝑀1 𝑆𝑡(𝑉𝑇+𝑢−𝑉)+1+𝐶𝑀 𝛽+𝐶𝑀𝐷𝑉 𝐷𝜏 (4-10) Donde se ha usado: 𝛽= 𝜌𝑝 𝜌𝑓 (4-11) 𝑆𝑡= 𝑡𝑝 𝑡𝑟 (4-12) Teoría de sedimentos 34 𝑉𝑇= 𝑡𝑝𝑔 𝑈 (4-13) 𝑡𝑝= 𝛽−1 𝑘𝑇𝑑𝑝2 18𝜐 (4-14) 𝑡𝑟= 𝐿𝑈 (4-15) El término representado en la ecuación 4-12 es el número de Stokes y mide la importancia relativa de la inercia frente a la resistencia aerodinámica a bajos números de Reynolds. El término representado en la ecuación 4-13 es la velocidad terminal de la partícula. La constante kT se emplea para corregir la desviación de la ley de Stokes para valores finitos del número de Reynolds relativo a la velocidad terminal de la partícula. De estas relaciones introducidas se puede deducir que, al disminuir el tamaño de la partícula, lo cual supone una reducción de los términos St y VT, los términos de inercia y gravitatorio ven reducida su importancia en relación con el término de fricción, por lo que la velocidad de la partícula se va aproximando a la del fluido. 35 35 Estudio numérico sobre sedimentación en el lecho de ríos 5. DESARROLLO NUMÉRICO USANDO TÉCNICAS CFD Este capítulo corresponde con la simulación del problema usando un software de dinámica de fluidos computacional como es Ansys Fluent. Puesto que la simulación no es una concreta, se ha realizado, en primer lugar, una calibración sin partículas para, posteriormente, pasar a realizar diversas pruebas cambiando la ubicación de los agitadores. 5.1. Introducción “La Dinámina de Fluidos Computacional o CFD (del inglés, Computational Fluids Dynamics) es el área de conocimiento que trata sobre la simulación numérica de flujos fluidos, transferencia de calor y fenómenos relacionados tales como reacciones químicas, combustión, aeroacústica, etc.” (KESSLER, M., 2016). La dinámina de fluidos computacional trata de resolver ecuaciones cuyo origen es la mecánica de fluidos haciendo uso de cálculo numérico. Mediante el uso de CFD se resuelven ecuaciones que de forma analítica no se podrían resolver, sino que es necesario hacer uso de la discretización del dominio. La principal ecuación que se resuelven mediante el uso de CFD es del tipo ecuación de transporte de una variable denominada ϕ. Esta ecuación tiene la siguiente forma: 𝜕 𝜕𝑡∫𝜌 𝜙 𝑑𝑉+ ∮ 𝜌 𝜙 𝑉 𝑑𝐴= ∮ 𝛤𝜙 ∇𝜙 𝑑𝐴 𝐴𝐴 𝑉+ ∫ 𝑆𝜙 𝑑𝑉 𝑉 (5-1) Donde los términos representan, por orden: - Término temporal. - Término convectivo. - Término difusivo. - Término fuente. Los métodos de discretización más conocidos son el método de las diferencias finitas, el método de los elementos finitos y el método de los volúmenes finitos, siendo el último el más empleado y el que usa Ansys Fluent. El método de volúmenes finitos permite la discretización del volumen de control para resolver numéricamente ecuaciones diferenciales. Divide el espacio fluido en una malla de discretización. En torno a cada punto de esta malla se construye un volumen de control de tal forma que el volumen total del fluido es igual a la suma de todos los volúmenes de control. Desarrollo numérico usando técnicas cfd 36 Ilustración 5-1. Volumen total discretizado en volúmenes de control Fuente: [22] Las ecuaciones que gobiernan la resolución del flujo media CDF se basan en leyes de conservación: - Conservación de la masa. - Conservación de la cantidad de movimiento. - Conservación de la energía. 5.2. Geometría La geometría del problema es la descrita en el apartado 2.7. Canal Arenoso, en el que se describe cómo es el canal de estudio. Se trata de un canal de 1200 metros de longitud con una sección transversal trapezoidal. Aunque la longitud sea de 1200 metros, el modelado de la geometría de Fluent permite un máximo de 1000 metros, por lo que las simulaciones se han tenido que efectuar para esta longitud. Ilustración 5-2. Geometría del canal El origen de coordenadas de la geometría se ha tomado en el centro del canal, lo cual significa que ambos extremos se encuentan en las coordenadas -500 metros (confluencia del canal con el Guadalquivir) y +500 metros (ubicación de las cántaras del bombeo). 37 37 Estudio numérico sobre sedimentación en el lecho de ríos La longitud del canal se divide en dos tramos, diferenciados por la rugosidad del material de la sección: - Tramo 1: desde la confluencia con el Guadalquivir hasta 660 metros. Recubrimiento de escollera en la margen derecha. - Tramo 2: resto del canal. Recubrimiento de escollera en ambas márgenes. A su vez, en la geometría, el tramo 1 se encuentra dividido en tres zonas, lo cual es debido al mallado con la presencia de agitadores. Ilustración 5-3. Elementos que forman la geometría en estudio Las características de la sección se muestran en la Tabla 5-1. Tabla 5-1. Características de la sección del canal Base menor (m) 8 Calado (m) 6 Talud 3H : 2V Base mayor (m) 26 Área (m2) 102 Perímetro mojado (m) 29.63 Desarrollo numérico usando técnicas cfd 44 Lo ideal para el tipo de problema que se presenta podría ser una malla formada por celdas tipo hexagonales, sin embargo, la presencia de los agitadores y el tipo de flujo que genera provoca que la calidad de una malla hexagonal sea baja. Por tanto, la posibilidad del uso de elementos hexagonales pasaría por generar una malla híbrida distinguiendo la zona de mayor agitación de la zona sin agitación. La malla ideal en caso de no tener agitadores sería una estructurada y formada por elementos hexagonales, como la que se presenta a continuación (Ilustración 5-16). En esta malla el valor de y+ (parámetro que se presenta en el apartado correspondiente al modelo de turbulencia) se controla manualmente mediante el uso de la herramienta sizing. Como adelanto, el valor de este parámetro debe mantenerse dentro de un rango y depende del espesor de los elementos en la zona cercana a las paredes del canal. Ilustración 5-16. Malla estructurada Sin embargo, los resultados que aporta al incluir los agitadores no son de calidad. El empleo de elementos tetraédricos a lo largo de todo el canal, además de generar un gasto computacional importante, no ofrece buenos resultados al tratarse de un flujo con una trayectoria regular en las zonas alejadas a los agitadores. Siguiendo estas ideas se genera una malla híbrida formada por elementos tetraédricos en la zona de los agitadores y elementos hexagonales en el resto del canal. Además, es una malla conforme, es decir, los nodos de las celdas se unen entre sí en los límites de las distintas zonas de la malla. Ilustración 5-17. Malla híbrida generada 45 45 Estudio numérico sobre sedimentación en el lecho de ríos Ilustración 5-18. Malla conforme La elección del tamaño de los elementos se ha hecho refinando la malla hasta que los resultados obtenidos sean similares usando una malla y otra. Ilustración 5-19. Malla en la zona de los agitadores Como se muestra en la Ilustración 5-18 y en la Ilustración 5-19, el refinamiento de la malla es mayor en la zona que rodea a los agitadores. Esto se ha conseguido generando una esfera de influencia alrededor de estos elementos en la que el tamaño de la malla se ha bajado respecto al resto del canal. El control del parámetro y+ se consigue haciendo uso de la herramienta inflation, la cual permite generar varias capas de celdas con el valor deseado. Desarrollo numérico usando técnicas cfd 46 Ilustración 5-20. Malla en la zona central del canal 5.4. Materiales Otro de los apartados a definir en la configuración del modelo de Ansys Fluent son los materiales. Las dos fases a considerar son el agua y las partículas sólidas. Las características de los materiales son: - Fase líquida: agua Tabla 5-2. Propiedades de la fase líquida ρw (kg/m3) 1000 μw (kg/(ms)) 0.001 υ (m2/s) 1.003e-6 - Fase sólida: partículas Tabla 5-3. Propiedades de la fase sólida Tipo de partícula ρs (kg/m3) d50 (mm) Tipo 1 Tipo 2 1500 1900 0.01 0.05 5.5. Modelos Se detallan a continuación los distintos modelos empleados para la resolución del problema. Previamente a la definición de los modelos, es necesario comentar que el problema a resolver se va a considerar estacionario, ya que se pretende obtener la posición final de las partículas, para así comprobar dónde se depositan según el efecto de los agitadores. El objetivo es obtener los puntos en los que se encuentran las partículas para impedir que avancen en el canal. 47 47 Estudio numérico sobre sedimentación en el lecho de ríos Ilustración 5-21. Definición de la gravedad y tipo de análisis También hay que indicar la gravedad, que en este caso se define según el eje y. 5.5.1. Turbulencia La turbulencia es un fenómeno que se da en la mayoría de flujos de la naturaleza e industriales. Se observa en flujos con un número de Reynolds alto y consiste en un movimiento del flujo irregular y aleatorio en tres dimensiones, del que no es posible obtener una solución analítica. La aparición de la Mecánica de Fluidos Computacional ha permitido calcular la turbulencia mediante las ecuaciones de Navier-Stokes. Algunos de los modelos existentes para resolver la turbulencia son los siguientes: - Direct Numerical Simulation (DNS). Tiene como fundamento la resolución de todas las escalas espaciales y temporales del flujo turbulento sin aproximaciones. El coste computacional es muy alto y su empleo se reduce a geometrías sencillas donde el número de Reynolds no es tan alto. - Reynolds-Averaged Navier-Stokes (RANS). Es el modelo más utilizado y el más económico desde el punto de vista de coste computacional. Es apropiado para la mayoría de los problemas de ingeniería ya que ofrece un nivel de precisión adecuado. Se resuelven las ecuaciones de Navier-Stokes usando aproximaciones. - Scale-Resolving Simulation (SRS). Es un modelo intermedio de los dos anteriores. En este caso se resuelve la turbulencia para al menos una parte del dominio del fluido. Además de estos modelos existen otros tantos que toman características de los presentados, conocidos como híbridos. Es importante la elección de un buen modelo de turbulencia, para lo que hay que considerar aspectos como la física del flujo, el nivel de precisión requerido, los recursos computacionales disponibles y el tiempo disponible para el cálculo. Reynolds-Averaged Navier-Stokes (RANS) Teniendo en cuenta lo anterior, es el modelo que se va a emplear en este caso, puesto que se adapta a las exigencias y recursos computacionales. En este modelo las variables de las ecuaciones de Navier-Stokes son descompuestas en una componente media Desarrollo numérico usando técnicas cfd 48 y otra de fluctuación debido a la turbulencia. ∅=∅+∅′ (5-4) Donde: - ∅: componente media - ∅′: componente de fluctuación Sustituyendo expresiones de esta forma (ecuación 5-4) en las ecuaciones de continuidad y de conservación de la cantidad de movimiento de Navier-Stokes, se obtienen las ecuaciones RANS: 𝜕𝜌 𝜕𝑡+𝜕 𝜕𝑥𝑖(𝜌𝑢𝑖)=0 (5-5) 𝜕 𝜕𝑡(𝜌𝑢𝑖)+𝜕 𝜕𝑥𝑗(𝜌𝑢𝑖𝑢𝑗)=−𝜕𝑝 𝜕𝑥𝑖+𝜕 𝜕𝑥𝑗[𝜇(𝜕𝑢𝑖 𝜕𝑥𝑗+𝜕𝑢𝑗 𝜕𝑥𝑖−23𝛿𝑖𝑗𝜕𝑢𝑙 𝜕𝑥𝑙)]+𝜕 𝜕𝑥𝑗(−𝜌𝑢′𝑖𝑢′𝑗        ) (5-6) Donde el término −𝜌𝑢′𝑖𝑢′𝑗        se conoce como tensor de esfuerzos turbulentos de Reynolds y se representa mediante 𝑅𝑖𝑗. Dentro del modelo de ecuaciones promediadas de Reynolds hay varios submodelos según la forma de modelar 𝑅𝑖𝑗. - Spalart-Allmaras - k-ε - k-ω - RSM - Transición laminar-turbulencia La lista anterior se presenta en orden ascendente de coste computacional. A medida que disminuye el coste computacional, disminuye también la precisión de la solución que ofrecen. Existen algunos submodelos más que no aparecen reflejados. Los más comunes de usar son el modelo k-ε y el modelo k-ω. Estos modelos resuelven el problema de la turbulencia usando dos ecuaciones adicionales e introducen un término llamado viscosidad turbulenta para resolver el tensor de esfuerzos turbulentos de Reynolds. La robustez, economía y razonable precisión del modelo k-ε lo hacen apropiado para resolver una gran cantidad de problemas. Dentro del submodelo kε se distingue entre Realizable, Standard y RNG. En Ansys Fluent es recomendable el empleo del modelo kε Realizable. Este modelo kε Realizable presenta las siguientes ecuaciones para k y ε: 𝜕 𝜕𝑡(𝜌𝑘)+𝜕 𝜕𝑥𝑗(𝜌𝑘𝑢𝑗)=𝜕 𝜕𝑥𝑗[(𝜇+𝜇𝑡 𝜎𝑘)𝜕𝑘 𝜕𝑥𝑗]+𝐺𝑘+𝐺𝑏−𝜌𝜀−𝑌𝑀+𝑆𝑘 (5-7) 𝜕 𝜕𝑡(𝜌𝜀)+𝜕 𝜕𝑥𝑗(𝜌𝜀𝑢𝑗)=𝜕 𝜕𝑥𝑗[(𝜇+𝜇𝑡 𝜎𝜀)𝜕𝜀 𝜕𝑥𝑗]+𝜌𝐶1𝑆𝜀−𝜌𝐶2𝜀2 𝑘+√𝜈𝜀+𝐶1𝜀𝜀𝑘𝐶3𝜀𝐺𝑏+𝑆𝜀 (5-8) Donde: - 𝐶1=𝑚𝑎𝑥[0.43,𝜂 𝜂+5] 49 49 Estudio numérico sobre sedimentación en el lecho de ríos - 𝜂=𝑆𝑘𝜀 - 𝑆=√2𝑆𝑖𝑗𝑆𝑖𝑗 - 𝐺𝑘 representa la generación de energía cinética turbulenta debido a los gradientes de la velocidad media 𝐺𝑘=−𝜌𝑢′𝑖𝑢′𝑗        𝜕𝑢𝑗 𝜕𝑥𝑖 - 𝐺𝑏 representa la generación de energía cinética turbulenta debido a la flotabilidad En este modelo la viscosidad turbulenta se modela de la siguiente forma: 𝜇𝑡=𝜌𝐶𝜇𝑘2 𝜀 (5-9) La diferencia de este modelo con otros modelos k-ε está en el valor de 𝐶𝜇, el cual, en este caso, no presenta un valor constante. Por ejemplo, para k-ε estándar tiene un valor de 0.09. El modelo emplea las siguientes constantes: - 𝐶1𝜀=1.44 - 𝐶2=1.9 - 𝜎𝑘=1 - 𝜎𝜀=1.2 5.5.1.1. Tratamiento cercano a las paredes (Near-Wall Treatment) Junto al modelo de turbulencia hay que elegir un modelo de resolución en las zonas cercanas a la pared. La presencia de una condición de contorno tipo pared afecta al comportamiento del flujo por ser una zona con cambio brusco de la velocidad. La zona próxima a las paredes puede dividirse en dos capas, una interior y otra exterior. A su vez, la interior se divide en tres subcapas: subcapa viscosa (viscous sublayer) en la que predominan los esfuerzos viscosos; subcapa de transición (buffer layer) donde la viscosidad y la turbulencia son importantes; subcapa turbulenta (turbulent region), donde la turbulencia tiene mayor importancia. Ilustración 5-22. Subdivisiones zona cercana a la pared Fuente: [14] Desarrollo numérico usando técnicas cfd 50 Existen dos enfoques para resolver la zona próxima a la pared. Por un lado, se encuentra la posibilidad de emplear un modelo de pared. Esta forma emplea relaciones semiempíricas para resolver la subcapa viscosa y de transición. Por otro lado, existe la posibilidad de emplear otro tipo de modelado el cual modifica el modelo de turbulencia empleado para resolver la zona entre la pared y la subcapa turbulenta. Esta forma de resolución se hace necesaria cuando la zona afectada por este fenómeno es la zona en la que se centra el estudio. Por este motivo, en este caso se va a emplear un modelo de pared, ya que la importancia de este estudio no reside esencialmente en el comportamiento del flujo en la zona próxima a las paredes, sino en todo el canal. Como se ve en la Ilustración 5-22, la región en la que nos encontramos depende del valor de y+, el cual se define como: 𝑦+=𝑢∞∆𝑦1 𝜈 (5-10) El valor de y+ se puede controlar mediante el parámetro Δy1, el cual se refiere a la distancia de la pared a la primera celda del mallado. El objetivo será mantener el valor de y+ en el rango de 30 a 300. Ilustración 5-23. Elección del modelo de turbulencia 5.5.2. Flujo multifásico Se entiende por flujo multifásico aquél que involucre a más de una fase. En este caso, se tienen dos fases: una fase líquida y una fase sólida. Existen distintos modelos para simular un flujo multifásico. Los avances en la mecánica de fluidos computacional han proporcionado la base para una mayor comprensión de la dinámica de flujos multifase. Existen actualmente dos enfoques para el cálculo numérico de flujos multifase: el enfoque Euler-Lagrange y el enfoque Euler-Euler. En el enfoque Euler-Lagrange, conocido como modelo discreto o lagrangiano (en inglés, Discrete Phase Model), la fase fluida (el flujo de agua en este caso) es tratada como una fase continua y se resuelve con las ecuaciones de Navier-Stokes, mientras que la fase dispersa (partículas en este caso) se calcula a partir del 51 51 Estudio numérico sobre sedimentación en el lecho de ríos campo fluido. La fase dispersa puede intercambiar masa, cantidad de movimiento y energía con la fase fluida. El empleo de este modelo va a permitir hacer un seguimiento de la trayectoria de las partículas en la zona de los agitadores, de tal forma que se controle si la recirculación generada por estos es la esperada o no lo es. En el enfoque Euler-Euler, se considera que el volumen de una fase no puede ser ocupado por otras fases, introduciendo el concepto de fracción volumétrica. Cada una de las fracciones se suponen funciones continuas en el tiempo y en el espacio y su suma es igual a la unidad. Las ecuaciones de conservación para cada fase se derivan para obtener un conjunto de ecuaciones de estructura similar para cada una de las fases. Estas ecuaciones se completan con relaciones empíricas o, en el caso de flujos granulares, aplicando teoría cinética. Además de estos dos enfoques, existe la posibilidad de implementar un modelo de transporte de sedimento haciendo uso de una UDS, el cual se va a definir a continuación. El empleo de este modelo va a permitir controlar las zonas en las que se aumenta la concentración debido al efecto de los agitadores. 5.5.2.1. Implementación mediante UDS El transporte de sedimento se puede modelar añadiendo a las ecuaciones del flujo principal (continuidad y cantidad de movimiento) una ecuación para un escalar. Si se define la concentración (kg sedimento / kg agua) como el escalar mencionado, la ecuación a implementar podría ser la ecuación de transporte genérica por convección-difusión para el escalar c, ecuación del tipo 5-1. 𝜕[𝜌𝑤𝑐] 𝜕𝑡 +𝜕[𝜌𝑤𝑢𝑖𝑐] 𝜕𝑥𝑖−𝜕[𝐽𝑐𝑖] 𝜕𝑥𝑖=𝑆𝑐 (5-11) Los términos representan, respectivamente: variación local de la concentración; flujo convectivo neto de la concentración; flujo difusivo de la concentración; generación volumétrica. Las unidades del escalar a definir son: 𝑐(𝑘𝑔 𝑠𝑒𝑑𝑖𝑚𝑒𝑛𝑡𝑜 𝑘𝑔 𝑎𝑔𝑢𝑎 ) El flujo difusivo se define como: 𝐽𝑐𝑖=Γ𝑐𝜕𝑐 𝜕𝑥𝑖 (5-12) Donde: - Γ𝑐: difusividad de la magnitud transportada Al tratar un flujo con sedimento hay que tener en cuenta los siguientes puntos para implementar el modelo: - Las partículas poseen movimiento propio, por lo que es necesario añadir la velocidad de sedimentación. Si se define la gravedad respecto al eje y, el término convectivo queda: 𝜕[𝜌𝑤𝑢𝑥𝑐] 𝜕𝑥 +𝜕[𝜌𝑤(𝑢𝑦−𝑤𝑦)𝑐] 𝜕𝑦 +𝜕[𝜌𝑤𝑢𝑧𝑐] 𝜕𝑧 (5-13) Donde wy es la velocidad de sedimentación, la cual depende del tamaño de grano adimensional 𝐷∗, la densidad y viscosidad cinemática del agua. 𝑤𝑗=10𝜇𝑤 𝑑50𝜌𝑤[(1+0.01𝐷∗3)1/2−1] (5-14) Desarrollo numérico usando técnicas cfd 52 𝐷∗=𝑑50[𝑔 𝜈2(𝜌𝑠 𝜌𝑤−1)]1/3 (5-15) Como la densidad del agua y la velocidad de sedimentación se consideran constantes, el término convectivo se puede expresar de la siguiente forma: 𝜕[𝜌𝑤𝑢𝑥𝑐] 𝜕𝑥 +𝜕[𝜌𝑤𝑢𝑦𝑐] 𝜕𝑦 +𝜕[𝜌𝑤𝑢𝑧𝑐] 𝜕𝑧 −𝜌𝑤𝑤𝑦𝜕𝑐 𝜕𝑦 (5-16) Así es más sencillo de implementar en el modelo UDF, puesto que queda la ecuación sin sedimento por un lado y el término fuente convectivo (𝑆𝑐) por otro. 𝑆𝑐=𝜌𝑤𝑤𝑦𝜕𝑐 𝜕𝑦 (5-17) - El movimiento es, generalmente, turbulento. Esto tiene efecto sobre la difusividad, que será la correspondiente al transporte difusivo turbulento en lugar del molecular. 𝐽𝑐𝑖=μt σt𝜕𝑐 𝜕𝑥𝑖 (5-18) Desarrollado queda: 𝜕𝐽𝑐𝑖 𝜕𝑥𝑖=μt σt𝜕2𝑐 𝜕𝑥2+ μt σt𝜕2𝑐 𝜕𝑦2+ μt σt𝜕2𝑐 𝜕𝑧2 (5-19) En esta expresión aparece el término de viscosidad turbulenta, del que se ha hablado en el apartado de modelado de la turbulencia. Este parámetro depende del modelo de turbulencia elegido. Para la UDF se va a emplear el modelo k-ε estándar, puesto que tiene un coste computacional menor. Para este modelo, la viscosidad turbulenta se define como: 𝜇𝑡 𝜌=𝐶𝜇𝑘2 𝜀 (5-20) Donde 𝐶𝜇 (constante de calibración del modelo de turbulencia empleado) vale 0.09 y el valor de σt (número de Prandtl) se tomará igual a 1. Una vez expuestas las particularidades para el problema de transporte de sedimento, la ecuación a resolver junto a las del flujo principal queda: 𝜕[𝜌𝑤𝑐] 𝜕𝑡 +𝜕[𝜌𝑤𝑢𝑥𝑐] 𝜕𝑥 + 𝜕[𝜌𝑤𝑢𝑦𝑐] 𝜕𝑦 +𝜕[𝜌𝑤𝑢𝑧𝑐] 𝜕𝑧 −μt σt[𝜕2𝑐 𝜕𝑥2+ 𝜕2𝑐 𝜕𝑦2+ 𝜕2𝑐 𝜕𝑧2]=𝜌𝑤𝑤𝑦𝜕𝑐 𝜕𝑦 (5-21) Para que el programa de cálculo tenga en cuenta la presencia de partículas, es necesario acoplar el modelo de transporte de sedimento con la resolución general del flujo. Para ello, hay que tener en cuenta el efecto del sedimento sobre las ecuaciones de continuidad, cantidad de movimiento y variables turbulentas. Este efecto se va a introducir haciendo uso de términos fuente (en inglés, source terms), los cuales se introducen en el apartado de condiciones de las celdas (en inglés, cell zone conditions). Los términos fuente a introducir son: 53 53 Estudio numérico sobre sedimentación en el lecho de ríos - Ecuación de cantidad de movimiento. La velocidad de sedimentación de las partículas genera una fuente de cantidad de movimiento en la ecuación que se corresponde con la componente de la velocidad alineada con la gravedad. En este caso la gravedad actúa según el eje OY, por lo que el término queda: 𝑆𝑢𝑦=(𝜌𝑠−𝜌𝑤)𝑔𝑐 (5-22) - Ecuación de energía cinética turbulenta. Al emplear el modelo de turbulencia k-ε, hay que considerar la generación de turbulencia del movimiento de sedimentación en la ecuación de la energía cinética. El término queda: 𝑆𝑘=−𝜌𝑠−𝜌𝑤 𝜌𝑤𝑔𝜕𝑐 𝜕𝑦 (5-23) - Ecuación para la tasa de disipación. Al igual que hay que introducir un término por el efecto de la sedimentación en la ecuación de energía cinética turbulenta, hay que considerar el mismo efecto sobre la tasa de disipación. 𝑆𝜀=−𝐶3𝜀𝜌𝑠−𝜌𝑤 𝜌𝑤𝑔𝜕𝑐 𝜕𝑦 (5-24) Donde: o 𝐶3𝜀=tanh{√𝑣𝑦2 (𝑣𝑥2+𝑣𝑧2)} La constante C3ε tiene en cuenta la verticalidad del vector velocidad en cada punto del dominio. En movimientos verticales (eje OY), el efecto de la sedimentación sobre la ecuación de la tasa de disipación es máxima, mientras que en movimientos horizontales no tiene ningún efecto. Sobre la ecuación de continuidad no es necesario actuar. Además del efecto del sedimento sobre la ecuación de transporte, es necesario definir unas condiciones de contorno. Se consideran dos tipos de condiciones de contorno sobre la concentración o escalar definido: - Valor impuesto. Se conoce el valor exacto o la distribución de valores de la concentración. 𝑐𝑥,𝑦,𝑧 ∈ 𝑠𝑢𝑝𝑒𝑟𝑓𝑖𝑐𝑖𝑒 =𝑐 (𝑡,𝑥,𝑦,𝑧) (5-25) Esta es la condición a imponer sobre la superficie de entrada y el lecho del canal. - Flujo impuesto. Se impone la derivada de la concentración sobre la superficie. 𝜕𝑐 𝜕𝑛 󰇍 |𝑥,𝑦,𝑧 ∈ 𝑠𝑢𝑝𝑒𝑟𝑓𝑖𝑐𝑖𝑒= 𝑞𝑐 (𝑡,𝑥,𝑦,𝑧) (5-26) Esta es la condición a imponer sobre las superficies no confinadas establecidas como condiciones de simetría y en la superficie de salida, de tal forma que se indique que la descarga de sedimento se Desarrollo numérico usando técnicas cfd 60 Siendo L la distancia que las partículas tienen que viajar, que en este caso se corresponde con la longitud del canal. Distintas pruebas han sido realizadas para ajustar el valor de estos parámetros, habiendo resultado como definitivas: - Número máximo de pasos: 50000 - Longitud de escala: 0.02 Ilustración 5-33. Parámetros de seguimiento de la trayectoria de las partículas Activando el modelo lagrangiano, hay que establecer unas condiciones de contorno para la trayectoria de las partículas sobre la superficie de entrada, salida y el lecho. - Condición de contorno tipo escape. Se establece cuando las partículas pueden atravesar el contorno debido al paso del flujo. Es la condición a definir en la entrada y en la salida. - Condición de contorno tipo reflect. Se establece cuando las partículas pueden rebotar contra el contorno al tratarse de una colisión sólido – sólido. Es el caso del lecho. 5.5.2.3. Modelo euleriano El enfoque Euler-Euler permite usar distintos modelos: - VOF. Cuando se tienen dos fluidos inmiscibles y la posición de la interfase entre ellos es de interés. Suele usarle para resolver flujos en canales donde los dos fluidos a definir son el agua y el aire. - Modelo mixto. Se emplea para resolver flujos con dos o más fases fluidas o de partículas. - Modelo de Euler. Es el modelo más complejo de los multifásicos. Resuelve un conjunto de n ecuaciones de continuidad y cantidad de movimiento para cada fase. Suele emplearse en fluidos con partículas en suspensión, en columnas de burbujas. De estos modelos, el VOF es apropiado para flujos en superficie libre, y el modelo mixto y euleriano para flujos en los que la fracción volumétrica de la fase de partículas exceda el 10%. Para flujos en los que la fracción volumétrica esté por debajo de este valor, es adecuado el uso del modelo lagrangiano o modelo disperso. Este motivo, añadido al gran coste computacional del modelo y la imposibilidad de tratar los agitadores usando la condición tipo fan cuando se activa el modelo euleriano, han hecho que se opte por no usar este modelo para la realización del estudio. 61 61 Estudio numérico sobre sedimentación en el lecho de ríos 5.6. Condiciones de contorno En todo estudio numérico es necesario definir las condiciones del problema a resolver en cada uno de los contornos existentes. Se establecen aquí las condiciones de contorno generales, habiendo especificado en cada modelo las particulares. En el modelo en estudio existen 4 contornos: - Entrada, llamada inlet. Se establece condición del tipo velocity inlet. Se emplea para definir la velocidad del flujo en la entrada, puesto que se conoce la sección y el caudal de estudio. Los parámetros necesarios para definir esta condición se especifican en la Tabla 5-5. Tabla 5-5. Parámetros condición de contorno entrada Parámetro Valor Velocidad (m/s) Re Intensidad turbulenta (%) Diámetro hidráulico (m) 0.1 1349832.44 2.74 13.77 El valor de estos parámetros se obtiene de las siguientes expresiones: o Intensidad turbulenta (I): 𝐼=𝑢′ 𝑢=0.16𝑅𝑒−1/8 (5-31) 𝑅𝑒=𝜌𝑤𝑣𝐷𝐻 𝜇𝑤 (5-32) o Diámetro hidráulico (DH): 𝐷𝐻=Á𝑟𝑒𝑎 𝑃𝑒𝑟í𝑚𝑒𝑡𝑟𝑜 𝑚𝑜𝑗𝑎𝑑𝑜 (5-33) Además de estos parámetros hay que definir la dirección del flujo: normal a la superficie. - Salida, llamada outlet. Se establece condición de contorno tipo pressure outlet. Ésta permite establecer la presión de salida en 0, con el objetivo de obtener la pérdida de carga en el canal mediante el conocimiento de la presión en la entrada. Los parámetros necesarios para definir esta condición de contorno se especifican en la Tabla 5-6. Desarrollo numérico usando técnicas cfd 62 Tabla 5-6. Parámetros condición de contorno salida Parámetro Valor Presión (Pa) Intensidad turbulenta (%) Diámetro hidráulico (m) 0 2.74 13.77 Además de estos parámetros hay que definir la dirección del flujo: normal a la superficie. - Lecho del canal, llamada wall. Se establece condición de contorno tipo wall. Se emplea en contornos que representan partes sólidas. Los parámetros necesarios para definir esta condición de contorno se especifican en la Tabla 5-7. Tabla 5-7. Parámetros condición de contorno lecho Parámetro Valor Altura de rugosidades del lecho (m) Altura de rugosidades de la escollera (m) 0.007 0.02 En este caso los parámetros a definir son los valores de la altura de rugosidades de los materiales que forman el lecho. El lecho se ha considerado formado por material arenoso pudiendo presentar alguna irregularidad debido a la presencia de vegetación. - Superficie libre, llamada symmetry. Para modelar el contorno de la superficie libre se ha empleado una condición de contorno de simetría, symmetry. Ilustración 5-34. Condición de contorno entrada y salida 63 63 Estudio numérico sobre sedimentación en el lecho de ríos Ilustración 5-35. Condición de contorno 5.7. Validación de los resultados El empleo de técnicas de dinámica de fluidos computacional proporciona una solución teórica de las simulaciones realizadas, teniendo que en cuenta las distintas simplificaciones que se han empleado. No todos los resultados que proporcionan las técnicas CFD son fiables, puesto que depende, entre otras cosas, de la precisión de la malla generada, de los modelos introducidos, etc. Para dar como válidos los resultados obtenidos en cada una de las simulaciones realizadas se ha tenido en cuenta lo siguiente: - Valor de los residuos escalados inferior a 10e-3. Los residuos escalados son la suma de los residuos de todas las ecuaciones, los cuales se van almacenando al final de cada iteración. - Control del valor de variables conocidas. Las variables empleadas han sido: presión en la entrada, velocidad del flujo en la salida, balance entre la cantidad de flujo a la entrada y a la salida. Hasta que dichas variables no se estabilicen, la solución no se considerará válida. - Control del valor y+. - Calibración de la malla sin uso de agitadores mediante el valor de la pérdida de carga en el canal. Se han realizado numerosas simulaciones, por lo que en el capítulo correspondiente a los resultados solo se mostrará la validación una única vez, entendiéndose que todos los resultados que se van a mostrar han seguido los mismos pasos de validación. 5.8. Agitadores Realizado el estudio de mercado de agitadores en el capítulo correspondiente a la solución a emplear, se pasa a describir el modelado de estos elementos en Ansys Fluent. Los agitadores pueden modelarse en Ansys Fluent de distintas formas encontrándose todas ellas en el apartado correspondiente a cell zone conditions que ofrece el software. - Fijar un valor (Fixed Value). Los agitadores se simulan fijando una velocidad. Es el caso que se empleó en el estudio anterior a este. - Término fuente (Source Term). Permite introducir un empuje para caracterizarlos. En este caso también se están tratando como un valor constante. - Condición tipo fan (3D Fan Zone). Es el método que permite introducir parámetros más reales del agitador. Desarrollo numérico usando técnicas cfd 64 Uno de los objetivos del estudio es introducir los agitadores con esta condición tipo fan. Ofrece como principal ventaja que el resultado se acerca a lo que se podría obtener mediante simulaciones en las que se incluyen los álabes del ventilador, sin ser necesario introducir geometrías complejas más allá de la del propio agitador que se describirá a continuación. Uno de los requisitos a cumplir es que la forma de los agitadores sea toroidal. El radio de los álabes (tip radius) y del cubo del que parten estos (hub radius) se reflejan mediante el radio exterior e interior, respectivamente. El espesor (thickness), las coordenadas del centro (fan origin), la cara de entrada del agitador (inlet fan zone) la cual debe ser de tipo interior y el punto de inflexión (inflection point) son parámetros geométricos más a introducir, representando el último de ellos la relación entre la longitud del álabe (b) y la distancia desde el radio interior al punto de mayor velocidad tangencial (a). De este último parámetro no se tiene dato, por lo que se dejará en valor por defecto de 0.83. Requisitos a cumplir para su uso: - Forma toroidal de los agitadores. - Todas las caras de la geometría del agitador han de definirse de tipo interior en las condiciones de contorno. Ilustración 5-36. Representación de los parámetros geométricos del agitador Fuente: [15] Ilustración 5-37. Diseño de la geometría de los agitadores 65 65 Estudio numérico sobre sedimentación en el lecho de ríos Además de la geometría, hay que incluir otros parámetros característicos entre los que se encuentran la velocidad de rotación (operating angular velocity) y su dirección (rotational direction), los límites de caudal (opcional) y la presión del agitador. Este último puede darse mediante un valor constante o introducir una curva característica (curva H-Q). También se piden unos valores correspondientes a valores de prueba (test angular velocity), considerando éstos del mismo valor que los de operación. La condición tipo fan tiene las siguientes limitaciones: - No puede emplearse con el modelo multifásico euleriano. - Es necesario activar la ecuación de energía al introducir una curva H-Q. Ilustración 5-38. Parámetros de los agitadores Para poder cumplir que todas las caras del agitador sean de tipo interior, es necesario nombrarlas, al igual que se hizo con las superficies de las condiciones de contorno. Desarrollo numérico usando técnicas cfd 66 Ilustración 5-39. Nomenclatura de las caras del agitador Debido al coste computacional que supone cada una de las simulaciones, la elección del agitador óptimo se va a realizar probando dos agitadores de los obtenidos en el estudio de fabricantes, siendo estos los de menor tamaño al corresponderse con los que menos potencia consumen. De este modo, se busca también la eficiencia económica. Realmente, durante el desarrollo del estudio se han realizado pruebas con más agitadores, pero al haber resultado efectivo uno de los mostrados en la Tabla 5-8, se ha decidido profundizar más en él. Las características de los agitadores empleados son las siguientes: Tabla 5-8. Agitadores empleados Agitador Diámetro (m) Potencia (kW) Velocidad de giro (rpm) 1 2 0.9 1.2 2.4 3 134 102 El tipo de curva H-Q empleada se muestra en el anexo 1. Ilustración 5-40. Curva H-Q para agitador de 0.9 m de diámetro 0.000 0.020 0.040 0.060 0.080 0.100 0.120 0.140 0.000 0.500 1.000 1.500 2.000 2.500 H (m) Q (m3/s) 67 67 Estudio numérico sobre sedimentación en el lecho de ríos 5.9. Opciones del solver Estas opciones también controlan la precisión y convergencia de los resultados. Hay dos tipos de solvers disponibles en Fluent: - Basado en presión. - Basado en densidad. En este caso se va a utilizar el basado en presión, aplicable a un amplio rango de regímenes de flujo, desde flujos incompresibles a baja velocidad hasta flujos compresibles de alta velocidad. Utiliza la cantidad de movimiento y las presiones como variables principales. Además, requiere menor memoria que el basado en densidad y ofrece mayor flexibilidad en el proceso de resolución. Cada una de las variables para resolver el modelo se calcula en el centro de las celdas, y deben ser interpoladas a las caras de los volúmenes de control. La precisión de esta interpolación viene determinada por los métodos de solución (en inglés, solution methods). El método empleado para cada una de las variables se presenta a continuación: - Acoplamiento presión-velocidad: SIMPLE. Es la opción por defecto. Se trata de una opción robusta y que se puede emplear en la mayoría de flujos. Al estar tratando un flujo estacionario no muy complejo se ha elegido esta opción. - Gradientes: Least-Squares Cell-Based. Es un método preciso. - Presión: Standard. Es el método por defecto. - Resto de variables: Second-Order Upwind. La convergencia es más lenta que si se usa el método de primer orden, sin embargo, genera un aumento de precisión. El proceso a seguir ha consistido en usar el primer orden para las primera iteraciones, para posteriormente usar el segundo orden. Resultados obtenidos 68 6. RESULTADOS OBTENIDOS Se muestra en el presente capítulo los resultados obtenidos. El procedimiento de resolución ha sido el siguiente: - Pruebas con el programa en un canal de sección cuadrada y 100 metros de longitud para probar el funcionamiento del software y los distintos modelos. - Pruebas en la geometría del estudio sin agitadores, con el objetivo de calibrar la malla. - Comprobación de las conclusiones alcanzadas en el anterior estudio. - Pruebas con agitadores en canal de 200 metros, debido al coste computacional. Estas pruebas han servido para obtener qué disposición de los agitadores puede considerarse efectiva. - Prueba en la geometría en estudio con las conclusiones obtenidas en las simulaciones de 200 metros. 6.1. Calibración del modelo sin agitadores Para comprobar si el modelo es adecuado, una de las validaciones a realizar es comparar la pérdida de carga calculada en la simulación con la pérdida de carga teórica. Otro de los aspectos a considerar en este apartado es la comprobación de los resultados obtenidos con los del estudio anterior. Esto viene motivado por el hecho de que la malla no es la misma, además de haber considerado otras diferencias en el modelo. 6.1.1. Pérdida de carga teórica del canal La pérdida de carga teórica se calcula mediante la siguiente expresión: ∆𝑃=𝑓𝐿 𝐷𝐻𝜌𝑤𝑣2 2 (6-1) Donde, además de los parámetros ya conocidos: - f: factor de fricción (depende del material) El valor del factor de fricción se obtiene a partir del Ábaco de Moody. Para su cálculo se necesita el valor de la rugosidad relativa del material y del número de Reynolds. 𝑅𝑢𝑔𝑜𝑠𝑖𝑑𝑎𝑑 𝑟𝑒𝑙𝑎𝑡𝑖𝑣𝑎= 𝑅𝑢𝑔𝑜𝑠𝑖𝑑𝑎𝑑 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑎 𝐷𝐻 (6-2) Donde, el valor de la rugosidad absoluta depende del material, y aparece definido en el apartado correspondiente a las condiciones de contorno (Tabla 5-7). 69 69 Estudio numérico sobre sedimentación en el lecho de ríos Tabla 6-1. Valores de rugosidad Material Rugosidad absoluta (m) Diámetro hidráulico (m) Rugosidad relativa Escollera Fondo 0.02 0.07 13.77 13.77 0.00145 0.00051 Teniendo en cuenta un valor del número de Reynolds de 1349832.44 y el diagrama de Moody (Ilustración 6-1) se obtiene los siguientes coeficientes de fricción y pérdidas de carga. Ilustración 6-1. Diagrama de Moody Tabla 6-2. Pérdida de carga Material f ΔP (Pa) Escollera Fondo 0.021 0.017 7.3 5.9 Obtenidos estos valores de pérdida de carga, se considerará válido el modelo si se encuentra dentro del rango establecido por ambos materiales, puesto que la sección del canal está compuesta por los dos. Un valor más concreto de la pérdida de carga podría obtenerse considerando el primer tramo de 660 metros con un factor de fricción correspondiente a 0.017, por ser predominante en la sección la zona sin escollera y, un segundo tramo de 340 metros con un factor de fricción de 0.021, por ser predominante la zona con escollera. Teniendo en cuenta esto, se obtiene un valor de pérdida de carga teórica de 6.67 Pa. Resultados obtenidos 76 6.2.1. Pruebas para encontrar ubicación óptima Debido al coste computacional de cada simulación en el canal completo, se ha buscado la posición más adecuada de los agitadores en pruebas con un canal de 200 metros. Estas pruebas se han realizado para dos tipos de agitadores. Teniendo en cuenta que el origen de coordenadas se ha establecido justamente en el centro del canal, como aparece reflejado en el apartado correspondiente a la geometría, las ubicaciones de los agitadores probadas han sido las siguientes: Tabla 6-4. Posiciones probadas de los agitadores Prueba Agitador X (m) Y (m) Z (m) 1 2 3 4 5 6 7 8 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 -10 10 -10 10 -10 10 -10 10 -10 10 -10 10 -15 15 -7.5 7.5 ϕ ϕ ϕ ϕ ϕ ϕ 2 2 2 2 2 2 ϕ ϕ ϕ ϕ -3 3 -4 4 -2.5 2.5 -4.75 4.75 -4 4 -3 3 -4 4 -4 4 Siendo ϕ el diámetro del agitador probado. Estas pruebas se han llevado a cabo para los agitadores mostrados en la Tabla 5-8, y, a continuación, se muestran los resultados de las trayectorias obtenidas. En la posición longitudinal se han probado a una distancia entre agitadores de 15, 20 y 30 metros (Valor X). En la separación transversal se han probado a distancias entre agitadores desde 5 a 9.5 metros (Valor Z). En la posición altitudinal se han probado a una distancia igual a su diámetro y a 2 metros. (Valor Y). - Para el agitador de 0.9 metros de diámetro. Los resultados obtenidos para este diámetro no han sido exitosos, por lo que la combinación de la potencia y diámetro probados no dan buen rendimiento. Se decide probar con diámetros superiores, puesto que en ellos se observa buen rendimiento, en lugar de estudiar un aumento de la potencia de éste, ya que este aumento la igualaría con la del siguiente agitador probado. Se incluye una imagen del resultado obtenido para la prueba 1 de la Tabla 6-4 para este diámetro de 77 77 Estudio numérico sobre sedimentación en el lecho de ríos agitador. Ilustración 6-16. Trayectoria partículas 1 y 2, prueba 1, agitador ϕ=0.9 m Solamente se produce la recirculación de algunas partículas, lo cual impide que se reduzca en gran medida la concentración aguas abajo de los agitadores. - Para el agitador de 1.2 metros de diámetro. Para este agitador se han obtenido mejores resultados. Si el rendimiento obtenido no hubiera sido bueno, se tendría que haber probado con otro tipo de agitador. Se incluyen a continuación la trayectoria de las partículas de las pruebas que mejores resultados han dado. Ilustración 6-17. Trayectoria partículas 1 y 2, prueba 1, agitador ϕ=1.2 m Resultados obtenidos 78 Ilustración 6-18. Trayectoria partículas 1 y 2, prueba 2, agitador ϕ=1.2 m Ilustración 6-19. Trayectoria partículas 1 y 2, prueba 3, agitador ϕ=1.2 m Ilustración 6-20. Trayectoria partículas 1 y 2, prueba 4, agitador ϕ=1.2 m 79 79 Estudio numérico sobre sedimentación en el lecho de ríos Ilustración 6-21. Trayectoria partículas 1 y 2, prueba 5, agitador ϕ=1.2 m Ilustración 6-22. Trayectoria partículas 1 y 2, prueba 6, agitador ϕ=1.2 m Resultados obtenidos 80 Ilustración 6-23. Trayectoria partículas 1 y 2, prueba 7, agitador ϕ=1.2 m Ilustración 6-24. Trayectoria partículas 1 y 2, prueba 8, agitador ϕ=1.2 m La efectividad de cada una de las pruebas realizadas se ha llevado a cabo mediante un balance entre las partículas que entran a la zona entre agitadores y las que salen. De esta forma, se considerará más efectiva aquella prueba en la que más partículas terminen su trayectoria entre los agitadores. 𝐸𝑓𝑒𝑐𝑡𝑖𝑣𝑖𝑑𝑎𝑑= 𝑃𝑎𝑟𝑡í𝑐𝑢𝑙𝑎𝑠 𝑞𝑢𝑒 𝑒𝑛𝑡𝑟𝑎𝑛−𝑃𝑎𝑟𝑡í𝑐𝑢𝑙𝑎𝑠 𝑞𝑢𝑒 𝑠𝑎𝑙𝑒𝑛 𝑃𝑎𝑟𝑡í𝑐𝑢𝑙𝑎𝑠 𝑞𝑢𝑒 𝑒𝑛𝑡𝑟𝑎𝑛 (6-3) Siendo la zona de agitación el volumen de control del que entran y salen las partículas. 81 81 Estudio numérico sobre sedimentación en el lecho de ríos Tabla 6-5. Cálculo de la efectividad de los agitadores probados Diámetro (m) Prueba Total partículas Entrada Salida Efectividad No agit 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 0.9 0 1 2 3 4 5 6 7 8 1 2270 2460 2460 2270 2270 2270 2270 2236 2270 2270 2270 2460 2460 2270 2270 2270 2270 2236 2270 2270 2270 1673 1082 1907 1725 1634 1544 1431 1634 1998 0 0.32 0.56 0.16 0.24 0.28 0.32 0.36 0.28 0.12 Puede observarse de la tabla anterior que la prueba más efectiva es la prueba 2 del agitador de 1.2 metros, lo cual es fácilmente apreciable en las ilustraciones mostradas. Para el agitador de diámetro 0.9 metros no se ha calculado ningún rendimiento además del de la prueba 1, ya que durante la realización de las pruebas se apreciaba su poca efectividad. Para comprobar la modificación del rendimiento del agitador de 1.2 metros en función de la variación de la posición se han realizado algunas pruebas más además de las mostradas. Un resumen de dichas pruebas y las anteriores se muestra en las tablas siguientes, donde X es la separación longitudinal, Y la separación del lecho y Z la separación transversal. Ilustración 6-25. Rendimiento en función de la separación longitudinal 0 0.1 0.2 0.3 0.4 0.5 0.6 0 5 10 15 20 25 30 35 Rendimiento Distancia x (m) Valor de X variable Para y = 1.2 m y z = 8 m Resultados obtenidos 82 Ilustración 6-26. Rendimiento en función de la separación respecto al lecho Ilustración 6-27. Rendimiento en función de la separación transversal Tras la realización de estas pruebas se concluye que una buena localización de los dispositivos aceleradores de corriente puede ser la de una separación longitudinal de unos 20 metros, una separación transversal de 8 metros y una localización vertical igual a su diámetro, teniendo en cuenta que se ha considerado como efectivo el agitador de 1.2 metros de diámetro. Para comprender lo que ocurre en la zona entre los agitadores, se va a hacer uso de la representación del contorno de la velocidad absoluta en distintos planos paralelos al lecho. Estas velocidades se corresponden con la simulación que se ha considerado como más efectiva. 0 0.1 0.2 0.3 0.4 0.5 0.6 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 Rendimiento Distancia y (m) Valor de Y variable Para x = 20 m y z = 8 m 0 0.1 0.2 0.3 0.4 0.5 0.6 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 Rendimiento Distancia z (m) Valor de Z variable Para x = 20 m e y = 1.2 m 83 83 Estudio numérico sobre sedimentación en el lecho de ríos Ilustración 6-28. Contorno de velocidades en plano y = 0.1 metros Ilustración 6-29. Contorno de velocidades en plano y = 0.5 metros Resultados obtenidos 84 Ilustración 6-30. Contorno de velocidades en plano y = 1.2 metros En las ilustraciones anteriores se observa la generación de zonas favorables a la decantación de partículas en las zonas internas entre los dos agitadores donde el valor de la velocidad es muy próximo a cero. Estas zonas favorables para la sedimentación son consecuencia de la recirculación generada por los agitadores, la cual produce zonas de baja velocidad en el centro de dicha recirculación. Otro aspecto que se visualiza en estas imágenes de contorno de velocidades es que tras la zona de agitación, buena parte del canal se mantiene con una velocidad menor a la velocidad de entrada de partículas, por lo que las partículas que logren superar la zona de agitación se mantendrán en el lado del canal con mayor velocidad. A modo de comparación, se incluye a continuación el contorno de velocidades para la prueba 3 del mismo agitador, la cual se diferencia de la elegida en que lo agitadores se separan transversalmente 5 metros en lugar de 8. Ilustración 6-31. Contorno de velocidades en el plano y = 0.5 metros Se puede observar en la ilustración anterior que no son tan claras las zonas favorables de sedimentación como 85 85 Estudio numérico sobre sedimentación en el lecho de ríos en el caso óptimo elegido, lo cual se traduce en una efectividad menor. No se alcanzan velocidades tan bajas, lo cual permite que un mayor número de partículas salgan de la zona de agitación. 6.2.2. Pruebas en canal completo De las simulaciones expuestas en el apartado anterior se puede concluir que la más efectiva en cuanto a recirculación de partículas está entre la prueba 1 y prueba 2 para el diámetro de 1.2 metros. Para la simulación del canal completo se ha considerado la ubicación especificada en la prueba 2. 6.2.2.1. Modelo lagrangiano Ilustración 6-32. Trayectoria partículas en canal completo Ilustración 6-33. Velocidad generada por agitadores Bibliografía 92 BIBLIOGRAFÍA [1] NAVARRO ORTIZ, E., Estudio de soluciones frente a la decantación de lodos en las instalaciones de bombeo y turbinado en el curso medio y bajo del Guadalquivir, Universidad de Sevilla, 2016. [2] GRAF, W. H., Fluvial hydraulics: flow and transport processes in channels of simple geometry, Chichester, 1998. [3] MARTÍN VIDE, J. P., Ingeniería de ríos, Barcelona, 2006. [4] MEYER, R. E., Introduction to mathematical fluid dynamics, Nueva York: Dover, 1982. [5] DÁVILA MARTÍN, J., Apuntes sobre dinámica de pequeñas partículas en flujos no uniformes, Universidad de Sevilla, 2017. [6] FERNÁNDEZ MORALES, A., Apuntes sobre modelo de transporte de sedimento, Universidad de Sevilla, 2017. [7] MORENO AREVALOS, V.M., Optimización de la succión de sólidos sumergidos utilizando bombas centrífugas con velocidad variable, Universidad Politécnica de Madrid, 2012. [8] CHAMORRO SOSA, G., Apuntes de curso de Obras Marítimas, Universidad de Sevilla, 2015. [9] Datos Espaciales de Referencia de Andalucía (DERA). Junta de Andalucía: Instituto de Estadística y Cartografía de Andalucía [consulta: 27 mayo 2017]. Disponible en: https://www.juntadeandalucia.es/institutodeestadisticaycartografia/DERA/ [10] Agitadores sumergibles Sulzer. Sulzer: agitadores sumergibles [consulta: 26 abril 2017]. Disponible en: http://www.sulzer.com/es/Products-and-Services/Agitators-Dynamic-and-Submersible- Mixers/Submersible-Mixers [11] Agitadores sumergibles Wilo. Wilo: agitadores sumergibles [consulta: 26 abril 2017]. Disponible en: http://www.sdm-sistemas.com/es/agitar/agitadores-sumergibles [12] Agitadores sumergibles Flygt. Flygt: agitadores sumergibles [consulta: 26 abril 2017]. Disponible en: http://www.xylemwatersolutions.com/scs/spain/eses/prensa/Eventos/Eficienciaenergetica/Documents/AGITADORES%20SUMERGIBLES.pdf [13] Red Saica Guadalquivir. Confederación Hidrográfica del Guadalquivir [consulta: 27 abril 2017]. Disponible en: http://www.chguadalquivir.es/calidad-de-las-aguas [14] ANSYS Inc, Fluent 15.0 Theory Guide, U.S.A., 2013. [15] ANSYS Inc, Fluent 15.0 Users Guide, U.S.A., 2013. [16] ANSYS Inc, Fluent 15.0 UDF Manual, U.S.A., 2013. [17] Roughness length. Wikipedia: roughness length [consulta: 6 mayo 2017]. Disponible en: https://en.wikipedia.org/wiki/Roughness_length [18] Embalse del Arenoso. Datos de embalses [consulta: 1 mayo 2017] [19] Modelos de viscosidad del agua en función de la concentración [consulta: 8 mayo 2017]. Disponible en: https://dspace.library.uu.nl/bitstream/handle/1874/553/c5.pdf;jsessionid=C85F8F3A84DFD267AB0FB4 75AC774065?sequence=4 [20] Confederación Hidrográfica del Guadalquivir, Estudio de la Evolución de la turbidez en el 93 93 Estudio numérico sobre sedimentación en el lecho de ríos Guadalquivir a través de las estaciones SAICA desde 2007 hasta 2012, 2012. [21] Bombas centrífugas. Aplicación, sistemas, principios fundamentales y selección. Monografías. [consulta: 4 mayo 2017] . Disponible en: http://www.monografias.com/trabajos36/bombascentrifugas/bombas-centrifugas2.shtml [22] KESSLER, M. Dinámica de Fluidos Computacional: ¿qué es?. [consulta: 1 agosto 2017]. Disponible en: http://www.esss.com.br/blog/es/2016/06/dinamica-de-fluidos-computacional-que-es/ Anexo I. Modelo de curva H-Q 94 ANEXO I. MODELO DE CURVA H-Q cell zone id 18 !CELL_ID Pressure(Pa) Flowrate(m3/s) 1236.299 0.758 1218.034 0.838 1197.943 0.918 1176.025 0.998 1152.281 1.078 1126.710 1.157 1099.313 1.237 1070.089 1.317 1039.039 1.397 1006.162 1.477 971.459 1.557 934.930 1.636 896.574 1.716 856.391 1.796 814.382 1.876 770.546 1.956 724.884 2.035 677.396 2.115 628.081 2.195 576.940 2.275 95 95 Estudio numérico sobre sedimentación en el lecho de ríos ANEXO II. MODELO UDF #include "udf.h" /*Hay que definir UDS para la concentración [kg sedimento/kg agua] */ /*Se van a guardar algunos términos creados en memoria mediante UDM*/ /* Se definen las constantes: G ( m/s2), wj (m/s), gs (t/m3), d50 (mm)*/ /* En función de la injection que estemos considerando, hay que cambiar wj y ps */ #define G 9.81 #define wj 0.000307 #define gs 1.70 #define d50 0.03 /*Fuente para la ecuación de sedimento*/ DEFINE_SOURCE(C_source,c,t,ds,eqn) { real x[ND_ND]; real source; C_CENTROID(x,c,t); source = C_R(c,t)*wj*C_UDSI_G(c,t,0)[1]; ds[eqn] = 0; C_UDMI(c,t,0) = C_R(c,t)*wj*C_UDSI_G(c,t,0)[1]; return source; } /*Fuente para la ecuación del momento en y*/ DEFINE_SOURCE(Ymon_source,c,t,dS,eqn) { real x[ND_ND]; real source; C_CENTROID(x,c,t); source = -(gs-1.00)*C_R(c,t)*G*C_UDSI(c,t,0); dS[eqn] = 0; C_UDMI(c,t,1) = -(gs-1.00)*C_R(c,t)*G*C_UDSI(c,t,0); return source; } Anexo II. Modelo UDF 96 /*Fuente de la energía cinética turbulenta (k)*/ DEFINE_SOURCE(k_source,c,t,dS,eqn) { real x[ND_ND]; real source; C_CENTROID(x,c,t); source = (gs-1.0)*(C_MU_T(c,t)/1.0)*G*C_UDSI_G(c,t,0)[1]; dS[eqn] = 0; C_UDMI(c,t,2) = (gs-1.0)*(C_MU_T(c,t)/1.0)*G*C_UDSI_G(c,t,0)[1]; return source; } /*Fuente para la ecuación de la tasa de disipación (epsilon)*/ DEFINE_SOURCE(epsilon_source,c,t,dS,eqn) { real x[ND_ND]; real source; real k1; real k2; real k3; real C3eps; C_CENTROID(x,c,t); k1 = (gs-1.0)*(C_MU_T(c,t)/1.0)*G*C_UDSI_G(c,t,0)[1]; k2 = (C_V(c,t)*C_V(c,t))/(C_U(c,t)*C_U(c,t)+C_W(c,t)*C_W(c,t)+0.00001); k3 = pow(k2,0.5); C3eps = tanh(k3); source = C3eps*k1; dS[eqn]=0; C_UDMI(c,t,3) = C3eps*k1; return source; } /*Difusividad*/ DEFINE_DIFFUSIVITY(Cdiff,c,t,i) { real mu_t = C_MU_T(c,t); 97 97 Estudio numérico sobre sedimentación en el lecho de ríos real sigma_t = 1; return mu_t / sigma_t; } /*Viscosidad*/ DEFINE_PROPERTY(viscosidad_concentracion,c,t) { real mu_r; real mu_1; real mu_2; real mu_3; real conc; conc = C_UDSI(c,t,0); mu_1 = 0.001; mu_2 = 0.00404527222+0.00304527222*(conc-0.70)/0.69; mu_3 = 0.00327*pow(10.0,0.132*conc); if (conc<=0.01) mu_r = mu_1; else if (conc<=0.70) mu_r = mu_2; else mu_r = mu_3; C_UDMI(c,t,4) = mu_r; return mu_r; } /*Cálculo del cortante*/ DEFINE_ADJUST(wall_shear_calc_fondo_inicio,dom) { Domain *domain = Get_Domain(1); Thread *ft, *ct; face_t f; cell_t c; real ts_force[ND_ND]; real area[ND_ND]; real tau; Anexo II. Modelo UDF 98 int wall_id=33; ft = Lookup_Thread(domain,wall_id); begin_f_loop(f,ft) { c=F_C0(f,ft); ct=THREAD_T0(ft); F_AREA(area,f,ft); N3V_V(ts_force,=,F_STORAGE_R_N3V(f,ft,SV_WALL_SHEAR)); tau=NV_MAG(ts_force)/NV_MAG(area); C_UDMI(c,ct,5)=tau; } end_f_loop(f,ft); } /*Perfil de concentracion de sedimento sobre el lecho (condición de contorno)*/ DEFINE_PROFILE(perfil_conc_lecho_fondo_inicio,thread,position) { Domain *domain = Get_Domain(1); Thread *ft, *ct; face_t f; cell_t c; real tension_shields; real tau_bed; real k4; real v_corte; real Re_granular; real mu; int wall_id = 33; ft = Lookup_Thread(domain,wall_id); begin_f_loop(f,ft) { c = F_C0(f,ft); ct = THREAD_T0(ft); tau_bed = C_UDMI(c,ct,5); tension_shields = tau_bed/(G*(gs-1.0)*d50); /*0.01 = 1000 (densidad agua)*0.00001m (d50)*/ C_UDMI(c,ct,6) = tension_shields; k4 = tau_bed/1000; mu = C_UDMI(c,ct,4); 99 99 Estudio numérico sobre sedimentación en el lecho de ríos v_corte = pow(k4,0.5); Re_granular = v_corte*d50/mu; C_UDMI(c,ct,7) = Re_granular; if (tension_shields < 0.25) F_PROFILE(f,thread,position) = 0.0; else if (tension_shields < 0.75) F_PROFILE(f,thread,position) = 0.60*(tension_shields-0.25); else F_PROFILE(f,thread,position) = 0.30; } end_f_loop(f,t) } /*Cálculo del cortante*/ DEFINE_ADJUST(wall_shear_calc_fondo_centro,dom) { Domain *domain = Get_Domain(1); Thread *ft, *ct; face_t f; cell_t c; real ts_force[ND_ND]; real area[ND_ND]; real tau; int wall_id=35; ft = Lookup_Thread(domain,wall_id); begin_f_loop(f,ft) { c=F_C0(f,ft); ct=THREAD_T0(ft); F_AREA(area,f,ft); N3V_V(ts_force,=,F_STORAGE_R_N3V(f,ft,SV_WALL_SHEAR)); tau=NV_MAG(ts_force)/NV_MAG(area); C_UDMI(c,ct,8)=tau; } end_f_loop(f,ft); } Anexo II. Modelo UDF 100 /*Perfil de concentracion de sedimento sobre el lecho (condición de contorno)*/ DEFINE_PROFILE(perfil_conc_lecho_fondo_centro,thread,position) { Domain *domain = Get_Domain(1); Thread *ft, *ct; face_t f; cell_t c; real tension_shields; real tau_bed; real k4; real v_corte; real Re_granular; real mu; int wall_id = 35; ft = Lookup_Thread(domain,wall_id); begin_f_loop(f,ft) { c = F_C0(f,ft); ct = THREAD_T0(ft); tau_bed = C_UDMI(c,ct,8); tension_shields = tau_bed/(G*(gs-1.0)*d50); /*0.01 = 1000 (densidad agua)*0.00001m (d50)*/ C_UDMI(c,ct,9) = tension_shields; k4 = tau_bed/1000; mu = C_UDMI(c,ct,4); v_corte = pow(k4,0.5); Re_granular = v_corte*d50/mu; C_UDMI(c,ct,10) = Re_granular; if (tension_shields < 0.25) F_PROFILE(f,thread,position) = 0.0; else if (tension_shields < 0.75) F_PROFILE(f,thread,position) = 0.60*(tension_shields-0.25); else F_PROFILE(f,thread,position) = 0.30; } end_f_loop(f,t) } /*Cálculo del cortante*/ 101 101 Estudio numérico sobre sedimentación en el lecho de ríos DEFINE_ADJUST(wall_shear_calc_final_1,dom) { Domain *domain = Get_Domain(1); Thread *ft, *ct; face_t f; cell_t c; real ts_force[ND_ND]; real area[ND_ND]; real tau; int wall_id=32; ft = Lookup_Thread(domain,wall_id); begin_f_loop(f,ft) { c=F_C0(f,ft); ct=THREAD_T0(ft); F_AREA(area,f,ft); N3V_V(ts_force,=,F_STORAGE_R_N3V(f,ft,SV_WALL_SHEAR)); tau=NV_MAG(ts_force)/NV_MAG(area); C_UDMI(c,ct,11)=tau; } end_f_loop(f,ft); } /*Perfil de concentracion de sedimento sobre el lecho (condición de contorno)*/ DEFINE_PROFILE(perfil_conc_lecho_final_1,thread,position) { Domain *domain = Get_Domain(1); Thread *ft, *ct; face_t f; cell_t c; real tension_shields; real tau_bed; real k4; real v_corte; real Re_granular; real mu; int wall_id = 32;