Socavación en lechos de ríos y fondos marinos

Socavación en lechos de ríos y fondos marinos i i Trabajo Fin de Máster Máster en Ingeniería Aeronáutica Socavación en lechos de ríos y fondos marinos Autor: Javier Rodríguez Suárez Tutor: Javier Dávila Martín Departamento de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 ii Resumen ii Socavación en lechos de ríos y fondos marinos iii iii Trabajo Fin de Máster Máster en Ingeniería Aeronáutica Socavación en lechos de ríos y fondos marinos Autor: Javier Rodríguez Suárez 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, 2019 iv Resumen iv Socavación en lechos de ríos y fondos marinos v v Trabajo de Fin de Máster: Socavación en lechos de ríos y fondos marinos Autor: Javier Rodríguez Suárez 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: El Secretario del Tribunal: Fecha: vi Resumen vi Socavación en lechos de ríos y fondos marinos vii vii RESUMEN Este documento versa sobre la creación de un modelo eficiente y robusto para el análisis de los procesos erosivos y de la formación de la socavación en torno a obstáculos que se encuentran en el fondo de los ríos. El interés de estos análisis reside en que la aparición de la socavación puede provocar problemas estructurales que pueden tener consecuencias catastróficas. Para la realización de estos modelos se ha empleado el Software ANSYS. Se emplean distintos tipos de modelos para estudiar la socavación, modelos monofásicos y bifásicos, analizando tanto casos bidimensionales como tridimensionales. Los resultados obtenidos son comparados entre ellos y con datos de la literatura para así poder determinar que herramienta obtiene mejores resultados tanto en términos cualitativos como en términos de tiempo necesario para su ejecución. viii Resumen viii Socavación en lechos de ríos y fondos marinos ix ix ABSTRACT This document deals with the creation of an efficient and robust model to analyse the erosive processes as well as the development of the scour around the obstacles which can be found on the bottom of the rivers. The benefit of these analysis is that the appearance of the scour may cause problems with catastrophic consequences. To develop these models ANSYS Software has been employed. Several types of models have been developed to analyse the scour processes such us single-phase models as well as multiphase models in order to analyse two-dimensional and three-dimensional cases. The obtained results have been compared among them and with data from several articles. Consequently, it could be determined which is the most powerful tool in terms of quality of the results and solving time. xvi Índice de tablas xvi Socavación en lechos de ríos y fondos marinos xvii xvii ÍNDICE DE ECUACIONES Ecuación 1: Parámetro de Shields .......................................................................................... 8 Ecuación 2: Fórmula de van Rijn ............................................................................................. 8 Ecuación 3: Ecuación de van Rijn. Aclaración ......................................................................... 8 Ecuación 4: Fórmula empleada para la velocidad de entrada................................................ 9 Ecuación 5: Ecuaciones de transporte para el método kε estándar .................................. 11 Ecuación 6: Modelado de la viscosidad turbulenta .............................................................. 11 Ecuación 7: Número de Courant ........................................................................................... 12 Ecuación 8: Condición de Courant ........................................................................................ 12 Ecuación 9: Aplicación de la condición de Courant al caso de malla dinámica .................... 13 Ecuación 10: Tensor de tensiones de corte sólidas .............................................................. 20 Ecuación 11: Definición de la fuerza de masa virtual ........................................................... 21 Ecuación 12: Definición del coeficiente de resistencia ......................................................... 21 Ecuación 13: Ecuaciones de conservación de la masa.......................................................... 21 Ecuación 14: Ecuaciones de conservación del momento ..................................................... 22 Ecuación 15: Definición de tensores de tensiones ............................................................... 22 Ecuación 16: Definición de las componentes de la viscosidad del sedimento ..................... 22 Ecuación 17: Función de distribución radial ......................................................................... 22 Ecuación 18: Aplicación de la condición de Courant al caso multifásico ............................. 23 Ecuación 19: Cálculo del número de Stokes ......................................................................... 23 Ecuación 20: Función de predicción de la socavación .......................................................... 24 Ecuación 21: Ecuaciones de transporte ................................................................................ 30 Ecuación 22: Definición de la velocidad de sedimentación .................................................. 30 Ecuación 23: Definición de la tensión de cizalladura turbulenta de Reynols ....................... 30 xviii Índice de ecuaciones xviii Ecuación 24: Definición del flujo escalar turbulento de Reynols ......................................... 30 Ecuación 25: Ecuaciones del método k-ε estándar .............................................................. 30 Ecuación 26: Definición del término del transporte de sedimentos .................................... 31 Ecuación 27: Parámetro de Shields crítico ........................................................................... 31 Ecuación 28: Definición de la concentración de sedimentos cerca del fondo ..................... 31 Ecuación 29: Definición del arrastre de sedimentos ............................................................ 31 Ecuación 30: Velocidad de desarrollo de la socavación ....................................................... 32 Ecuación 31: Definición del perfil de velocidades ................................................................ 32 Ecuación 32: Definición del perfil de la energía cinética ...................................................... 32 Ecuación 33: Definición del perfil de la escala de longitud .................................................. 32 Ecuación 34: Definición de la tasa de disipación de la turbulencia ...................................... 33 Ecuación 35: Cálculo del tiempo de paso a partir del número de Courant para el caso tridimensional ............................................................................................................................. 39 Ecuación 36: Definición del nuevo parámetro para definir el tiempo de paso máximo ...... 42 Socavación en lechos de ríos y fondos marinos xix xix NOMENCLATURA ACRÓNIMOS ❖ Computational Fluid Dynamics CFD ❖ Large Eddy Simulation LES ❖ Standard 𝑘−𝜀 SKE ❖ Realizable 𝑘−𝜀 RKE ❖ User Defined Function UDF ❖ User Defined Memory UDM ❖ User Defined Scalar UDS ❖ Volume of Fluid VOF xx Nomenclatura xx NOTACIÓN ❖ Altura de la zona rugosa 𝑘𝑠 ❖ Amplitud de la perturbación inicial 𝛿 ❖ Ángulo de fricción interna 𝜑 ❖ Ángulo de reposo del sedimento 𝜙 ❖ Área mínima de las caras del mallado 𝐴𝑚í𝑛 ❖ Coeficiente de resistencia CD ❖ Coeficiente de restitución 𝑒𝑠𝑠 ❖ Concentración de sedimentos 𝑐 ❖ Concentración de sedimentos cerca del fondo 𝑐0 ❖ Constante para el cálculo de la velocidad 𝜅 ❖ Constante para el cálculo del número de Stokes 𝑘𝑇 ❖ Contribución de la dilatación fluctuante a la tasa de disipación YM ❖ Densidad del fluido 𝜌𝑓 ❖ Densidad del sedimento 𝜌𝑠 ❖ Diámetro medio de la arena 𝐷50 ❖ Diámetro medio de la arena adimensional 𝐷∗ ❖ Energía cinética de la turbulencia k ❖ Esfuerzo de cortadura crítico 𝜏𝑐𝑟 ❖ Esfuerzo de cortadura 𝜏 ❖ Espesor aproximado de la capa límite 𝛿 ❖ Flujo escalar turbulento de Reynols 𝐽𝑐,𝑗 ❖ Fracción volumétrica del fluido 𝛼𝑓 ❖ Fracción volumétrica del sedimento 𝛼𝑠 ❖ Función de distribución radial 𝑔0,𝑠𝑠 ❖ Generación de energía cinética turbulenta debida a la flotabilidad Gb ❖ Generación de energía cinética turbulenta debida a los gradientes de velocidad Gk ❖ Gravedad g ❖ Longitud característica 𝐿 ❖ Nivel de turbulencia del flujo libre 𝑇 Socavación en lechos de ríos y fondos marinos xxi xxi ❖ Número de Reynols Re ❖ Parámetro de Shields 𝜃 ❖ Parámetro de Shields crítico 𝜃𝑐𝑟 ❖ Parámetro de Shields crítico en el fondo 𝜃𝑐𝑟,0 ❖ Pendiente del fondo 𝛼 ❖ Porosidad de la arena 𝛾 ❖ Presión de fricción 𝑝𝑓𝑟𝑖𝑐𝑡𝑖𝑜𝑛 ❖ Presión estática 𝑝+ ❖ Ratio de densidades 𝑠 ❖ Rugosidad del fondo 𝑧0 ❖ Segundo invariante del tensor de tensiones desviador 𝐼2𝐷 ❖ Tasa de disipación de la turbulencia ε ❖ Temperatura granular Θ𝑠 ❖ Tensión de cizalladura turbulenta de Reynols 𝜏𝑖𝑗 ❖ Tensor de tensiones de cizalladura del fluido 𝜏𝑓 ❖ Tensor de tensiones de cizalladura del sedimento 𝜏𝑠 ❖ Tensor identidad 𝐼 ❖ Tiempo de relajación 𝑡𝑝 ❖ Tiempo de residencia 𝑡𝑟 ❖ Velocidad 𝑢 ❖ Velocidad característica U ❖ Velocidad de desarrollo de la socavación 𝑉 ❖ Velocidad de sedimentación 𝑤𝑠 ❖ Velocidad del fluido 𝑣𝑓 ❖ Velocidad del flujo libre 𝑢0 ❖ Velocidad del sedimento 𝑣𝑠 ❖ Velocidad terminal para la fase sólida 𝑣𝑟,𝑠 ❖ Velocidad terminal (número de Stokes) 𝑉𝑇 ❖ Viscosidad cinemática 𝜈 ❖ Viscosidad turbulenta µt ❖ Viscosidad absoluta o dinámica del fluido µf ❖ Viscosidad absoluta o dinámica del sedimento µs ❖ Viscosidad cinética 𝜇𝑠,𝑘𝑖𝑛 ❖ Viscosidad colisional 𝜇𝑠,𝑐𝑜𝑙 xxii Nomenclatura xxii ❖ Viscosidad friccional 𝜇𝑠,𝑓𝑟 Socavación en lechos de ríos y fondos marinos xxiii xxiii xxiv Nomenclatura xxiv Socavación en lechos de ríos y fondos marinos 1 - 1 - 1. INTRODUCCIÓN 1.1 OBJETIVOS Y MOTIVACIÓN En la actualidad, los fondos de los ríos están repletos de tuberías para transporte de gas u otros elementos, así como de pilares empleados para cimentar grandes estructuras. La presencia de estos elementos provoca la aparición de unos fenómenos de erosión que pueden provocar, con el paso del tiempo, tanto problemas derivados del transporte de sedimentos en la zona cercana a las tuberías (Sumer & Fredsøe, 2002), como provocar la debilitación de los cimientos de grandes estructuras, especialmente en puentes (Whitehouse, 1998). Esta erosión puede tener graves consecuencias, ya no sólo en cuanto a la pérdida de tiempo y dinero requeridos para reparar los desperfectos provocados por la erosión, sino que también puede acarrear importantes problemas para la seguridad de las personas, especialmente en el caso de los puentes. Este fenómeno de erosivo se conoce como erosión fluvial (debido a que es provocada por la corriente de un río) que se produce cuando la fuerza de la corriente del agua consigue vencer la resistencia de los diferentes materiales presentes en el fondo del río, iniciándose así el proceso de socavación. La socavación es esa remoción de material producida por la corriente que, como se ha comentado antes, puede llegar a tener efectos devastadores si no se analiza correctamente su comportamiento y se tienen en cuenta sus consecuencias. Existen dos tipos principales de socavación: la socavación general que es la bajada generalizada del fondo del río, provocando así un aumento de su profundidad; y la socavación local que se produce por la presencia de una estructura o un estrechamiento en el cauce del río, y por lo tanto sólo afecta Figura 1: Esquema de la socavación en torno a pilares https://upload.wikimedia.org/wikipedia/commons/7/77/Local_scour.gif 8 Validación del modelo - 8 - 2.1 MODELO MONOFÁSICO 2.1.1 Definición del modelo El primer método para predecir la socavación es el empleo de una malla dinámica. Para validarlo, se van a intentar replicar los resultados obtenidos por Kazaminezhad y Yeganeh- Bakhtiary en sus dos artículos (Kazaminezhad & Yeganeh-Bakhtiary, 2011) (Yeganeh-Bakhtiary, Kazaminezhad, Etemad-Shahidi, Baas, & Chieng, 2011), aunque en estos dos artículos, el método empleado para su resolución es distinto al que se va a emplear en el caso bajo estudio. El método que se va a emplear es similar al utilizado en el artículo (Ganesh Elapolu, Majumdar, Lottes, & Kostic, 2012). El método consiste en el cálculo de la cantidad de sedimento que debe ser eliminado entre dos actualizaciones del mallado, y como eso afectaría a la forma de la interfase entre el agua y la arena (hay que recordar que, en este primer modelo, sólo se tendrá fase de agua). Esto se va a realizar a partir de la definición de la velocidad vertical a la que se desplaza el fondo del dominio fluido. Para ello, es utiliza el criterio de Shields que se basa en el empleo de un número adimensional que indica el inicio del movimiento de las partículas del sedimento que se encuentran en un campo fluido (Shields, 1936). 𝜃= 𝜏 𝑔(𝜌𝑠−𝜌𝑓)𝐷50 Ecuación 1: Parámetro de Shields Por otro lado, el cálculo de la velocidad vertical de descenso con la que se desplaza el fondo del dominio para simular la socavación se obtiene a partir de la ecuación de van Rijn (C. van Rijn, 1984). 𝑉𝑦=0.00033(Δ𝑔𝐷50)12𝐷∗0.3 1−𝛾 (𝜃 𝜃𝑐𝑟−1) Ecuación 2: Fórmula de van Rijn Δ=𝜌𝑠−𝜌𝑓 𝜌𝑓 𝐷∗=𝐷50[Δ𝑔 𝜈2] 𝜃𝑐𝑟=0.23 𝐷50+0.054[1−𝑒−𝐷∗0.85 23 ] Ecuación 3: Ecuación de van Rijn. Aclaración Los valores necesarios para resolver las ecuaciones anteriores son: • 𝜌𝑠=2650𝑘𝑔 𝑚3 • 𝜌𝑓=998.2𝑘𝑔 𝑚3 • 𝐷50=0.00036 𝑚→ 𝑑𝑖á𝑚𝑒𝑡𝑟𝑜 𝑚𝑒𝑑𝑖𝑜 𝑑𝑒𝑙 𝑠𝑒𝑑𝑖𝑚𝑒𝑛𝑡𝑜 • 𝑔=9.8𝑚 𝑠2 • 𝜈=0.001003𝑘𝑔 𝑚𝑠 • 𝛾=1−0.63=0.37→𝑝𝑜𝑟𝑜𝑠𝑖𝑑𝑎𝑑 𝑑𝑒𝑙 𝑠𝑒𝑑𝑖𝑚𝑒𝑛𝑡𝑜 Socavación en lechos de ríos y fondos marinos 9 - 9 - Introduciendo estos valores en las Ecuaciones 2 y 3 se obtiene el siguiente resultado 𝑉= 0.0000195(𝜏 𝜏𝑐𝑟−1). Una vez descrita la idea principal de este primer método, se pasa a definir la geometría y las características del mismo. En la Figura 5 se muestran las dimensiones del caso analizado. El diámetro que se ha empleado para la circunferencia es de D=0.1m. Por lo tanto, las dimensiones de la zona fluida son de 3x0.4 m, estando la circunferencia situada a un metro de la entrada de la zona fluida y por lo tanto a dos metros de la salida. Como se puede observar, la zona fluida que hay detrás de la tubería es mayor que la que hay delante. Esto se debe a que así se puede observar mejor el comportamiento de los torbellinos que se deben formar aguas abajo de la circunferencia debidos al flujo incidente. Como se ha comentado previamente, la zona fluida está compuesta únicamente por agua en este primer método. Las propiedades del agua de densidad y viscosidad dinámica 𝜌= 998.2 𝑘𝑔 𝑚3 𝑦 𝜈=0.001003𝑘𝑔 𝑚𝑠. Figura 5: Esquema modelo malla dinámica Las condiciones de contorno empleadas en el análisis son: • En la entrada: velocidad de entrada del agua. Para ello, se ha empleado una función definida por el usuario (UDF) a partir de la cual, se define el perfil de velocidades a la entrada del dominio fluido. La fórmula empleada para ello es la mostrada en la Ecuación 4, donde 𝑢𝑚á𝑥=0.5 𝑚/𝑠 e 𝑦𝑚á𝑥=0.4 𝑚. La UDF empleada para imponer la condición de velocidad de agua a la entrada se muestra en el Anexo A. • En la salida: se ha empleado la condición de presión de salida. • En la parte superior: se ha empleado la condición de simetría. • En la parte inferior: se ha empleado la condición de pared. Esta parte es la que simula la interfase entre el agua y la arena, y será la parte en la que se aplique la malla dinámica a partir de otra UDF como se explicará más adelante. • En el cilindro: se ha empleado la condición de pared. 𝑢(𝑦)=𝑢𝑚á𝑥(𝑦 𝑦𝑚á𝑥)17=0.5(𝑦 0.4)17 Ecuación 4: Fórmula empleada para la velocidad de entrada x y 10 Validación del modelo - 10 - Para poder obtener unos resultados correctos, se ha añadido una pequeña perturbación sinusoidal inicial bajo la circunferencia, facilitando así, la formación de la socavación. Esta perturbación inicial será de pequeñas dimensiones, en torno al 2-3% del diámetro de la tubería para intentar que esa perturbación afecte lo menos posible a los resultados que se obtengan. Figura 6: Detalle perturbación inicial caso malla dinámica El mallado empleado se muestra en las Figuras 7 y 8. Se ha realizado un mallado estructurado, disminuyendo el tamaño de las celdas en la zona cercana al cilindro, así como en la zona cercana a la interfase entre el agua y la arena. Esto se debe a que en esas zonas es donde es más importante captar con mayor precisión los valores de las diferentes magnitudes para poder obtener unos mejores resultados. Se ha empleado un mallado no uniforme, con 40 divisiones en la dirección vertical y 190 en la dirección horizontal. Se tiene con todo esto, un total de 7743 celdas. Figura 7: Mallado caso malla dinámica Figura 8: Detalle mallado caso malla dinámica 𝛿 Socavación en lechos de ríos y fondos marinos 11 - 11 - 2.1.2 Resolución del modelo Una vez descritas la geometría y las características del modelo empleadas para este primer método se van a mostrar en esta sección las principales ecuaciones empleadas por Fluent en la resolución del modelo, así como las principales variables y parámetros empleados por el mismo. En primer lugar, el modelo viscoso empleado para la resolución del caso bajo estudio es el k-ε estándar. Este mtodo fue propuesto por Launder y Spalding, y permite determinar tanto el tamaño como la escala temporal de la turbulencia (Launder & Spalding, 1972). Las ecuaciones empleadas por este modelo para la determinación de la energía cinética de la turbulencia, k, y la tasa de disipación de esta, ε, se muestran en la Ecuación 5. 𝜕 𝜕𝑡(𝜌𝑘)+𝜕 𝜕𝑥𝑖(𝜌𝑘𝑢𝑖)=𝜕 𝜕𝑥𝑖[(𝜇+𝜇𝑡 𝜎𝑘)𝜕𝑘 𝜕𝑥𝑗]+𝐺𝑘+𝐺𝑏−𝜌ε−𝑌𝑀+𝑆𝑘 𝜕 𝜕𝑡(𝜌ε)+𝜕 𝜕𝑥𝑖(𝜌ε𝑢𝑖)=𝜕 𝜕𝑥𝑖[(𝜇+𝜇𝑡 𝜎ε)𝜕ε 𝜕𝑥𝑗]+𝐶1εε𝑘(𝐺𝑘)−𝑐2ε𝜌ε2 𝑘+𝑆ε Ecuación 5: Ecuaciones de transporte para el método kε estándar En la Ecuación 5, Sk y Sε en el caso bajo estudio son nulas ya que no se ha empleado ninguna función definida por el usuario (UDF) para determinar la viscosidad turbulenta. En cuanto a las constantes que aparecen en las ecuaciones, se han empleado los siguientes valores: 𝐶𝜇=0.09,𝐶1ε=1.44,𝐶2ε=1.92,𝜎𝑘=1 𝑦 𝜎ε=1.3, que son los valores que proporciona Fluent por defecto. Por otro lado, la viscosidad turbulenta se ha modelado como se muestra en la Ecuación 6. 𝜇𝑡=𝜌𝐶𝜇𝑘2 ε Ecuación 6: Modelado de la viscosidad turbulenta En cuanto a la malla dinámica empleada en la resolución del modelo, se ha creado una función definida por el usuario (UDF) para modelar la zona de la interfase entre el agua y la arena. El propósito de esta función es determinar cuándo y dónde es necesario actualizar la malla para captar correcta y progresivamente la formación de la socavación bajo la circunferencia, la cual se está analizando. Esta UDF se muestra en el Anexo B. Unos de los principales problemas que se identificó durante la ejecución del modelo, fue la deformación de la parte delantera que se muestra en la Figura 9. Este error se debe a la condición de contorno impuesta en la entrada del dominio. Para resolver este problema, se ha impuesto que los cinco primeros milímetros del fondo del dominio no puedan deformarse. 12 Validación del modelo - 12 - Figura 9: Error en la deformación de la interfase Para llevar a cabo el mallado dinámico en Fluent, se han empleado dos de los esquemas que este software presenta para la resolución de la malla dinámica: layering y smoothering. El primero de ellos implica la creación y destrucción de las celdas cuando la zona crece o se encoge; mientras que el segundo método hace que las celdas se muevan como si estuvieran conectadas con resortes. Por otro lado, para que no aparezcan problemas en el resto del dominio, se han añadido los siguientes tipos de movimiento a las diferentes partes de este: • Parte inferior: función definida por el usuario (UDF), para simular el desplazamiento vertical del fondo del dominio a partir de la velocidad obtenida en la Ecuación 2. • Entrada y salida: deforming, para que puedan cambiar su tamaño en caso de que sea necesario debido a la variación del mallado. • Parte superior y circunferencia: stationary, ya que en esas zonas no está permitido el movimiento. Un parámetro importante para la resolución del modelo es el paso de tiempo que se utiliza para resolverlo. Este parámetro depende principalmente de la velocidad del fluido y del tamaño del mallado. Para determinar este parámetro se ha empleado en número de Courant (Courant, Friedrichs, & Lewy, 1967 [1928]). En las Ecuaciones 7 y 8, se muestra la definición del número de Courant y la condición de Courant, que es una condición de convergencia para ecuaciones diferenciales en derivadas parciales. 𝐶=𝑢Δ𝑡 Δ𝑥 Ecuación 7: Número de Courant 𝑢Δ𝑡 Δ𝑥<𝐶 𝐷𝑜𝑛𝑑𝑒 0<𝐶≤1 Ecuación 8: Condición de Courant En el caso bajo estudio, el máximo valor para el intervalo de tiempo (suponiendo C=1) es Δ𝑡𝑚á𝑥=0.0318 𝑠. Este valor se ha obtenido aproximando el tamaño de las celdas del mallado en la dirección del flujo. Para ello, se ha dividido el tamaño total del dominio fluido en la dirección del flujo, que son tres metros, entre el número de celdas impuesto en dicha dirección, que es 190. Este cálculo se muestra en la Ecuación 9. Una vez obtenido ese valor, se ha decido Socavación en lechos de ríos y fondos marinos 13 - 13 - emplear como intervalo de tiempo Δ𝑡=0.03 𝑠 que, como se verá más adelante, es un valor correcto ya que los resultados no divergen. Este valor es más pequeño que el valor límite, por lo que según el número de Courant los resultados deberían ser correctos. Pero para lograr unos resultados más precisos, se ha decido elegir ese valor. Los parámetros necesarios para el cálculo del número de Courant son: • El tamaño medio de cada celda en la dirección x (∆𝑥): este valor se ha calculado a partir del tamaño del dominio en la dirección x, que son tres metros, y el número de divisiones en esa misma dirección que, en este caso, son 190. • La velocidad de entrada al dominio fluido: en este caso 0.5 𝑚𝑠 ⁄. Δ𝑡𝑚á𝑥=Δx 𝑢≅3 190 0.5=0.0318 𝑠 → 𝐶= Δ𝑡 Δ𝑡𝑚á𝑥=0.03 0.0318=0.943<1 Ecuación 9: Aplicación de la condición de Courant al caso de malla dinámica 2.1.3 Resultados Una vez presentadas las principales características del problema, así como las ecuaciones y parámetros que definen el modelo, se está en disposición de mostrar los resultados que se obtienen con su resolución. Se van a mostrar en esta sección los resultados obtenidos con el método monofásico presentado previamente. La mayoría de los datos necesarios para su ejecución se han mostrado en las secciones previas, aunque es importante comentar que la amplitud que se le ha dado a la perturbación inicial es de 𝛿=2 𝑚𝑚. Es interesante comentar la importancia de esta pequeña perturbación inicial ya que, si no se impusiera, se producirían discontinuidades en los resultados, especialmente en la parte inicial del desarrollo de la socavación como se pueden observar en la Figura 10. El principal inconveniente que presenta esta situación es la aparición de ese ‘pico de socavación’ que provoca que lo resultados posteriores no sea correctos ya que el proceso de socavación debe ser más uniforme. Figura 10: Problemas por ausencia de perturbación inicial 14 Validación del modelo - 14 - Debido a que el objetivo es comprobar cuál es la profundidad máxima que alcanza la socavación, hay que analizar el proceso de socavación hasta que comienza a estabilizarse el valor de dicha profundidad. Se muestran en la Figura 11 los resultados obtenidos para la profundidad máxima durante los primeros 90 segundos de ejecución del modelo. El máximo que se alcanza a los 90 segundos es de 2.6 mm aproximadamente, es decir, un 26% del diámetro de la circunferencia. Figura 11: Profundidad máxima de la socavación Otro resultado interesante es el campo de velocidades, y su evolución en el tiempo. En la Figura 12 se muestra este resultado para el caso anteriormente descrito. El principal problema es que el modelo utilizado predice la socavación, pero no predice el transporte de sedimentos que hacen que se produzca una acumulación de arena detrás de la circunferencia. Esta acumulación de sedimentos se muestra en la Figura 13 para los análisis realizados por Kazaminezhad y Yeganeg-Bakhtiary. Figura 12: Evolución del campo de velocidades con el tiempo (2 mm) t=2 seg t=10 seg t=90 seg t=25 seg Socavación en lechos de ríos y fondos marinos 15 - 15 - En la ejecución del caso bajo estudio se pudo ver que los resultados, en términos de profundidad máxima alcanzada por la socavación, dependían claramente del tamaño de la perturbación inicial impuesta en el dominio. Por lo tanto, se resolvió el modelo con distintos valores de la perturbación inicial para comprobar como varía el desarrollo de la socavación bajo el cilindro. Los resultados obtenidos con este método comparan con los resultados obtenidos por en el estudio realizado por Kazaminezhad y Yeganeg-Bakhtiary (Kazaminezhad & Yeganeh- Bakhtiary, 2011). Este estudio consistió en el análisis del proceso de socavación análogo al aquí estudiado, aunque el método empleado fue distinto, basado en el modelado de las fases de agua y arena como fluidos independientes. Figura 13: Formación de la socavación según el artículo de Kazaminezhad y Yeganeg-Bakhtiary Para terminar con el análisis de los resultados obtenidos con este primer modelo, se van a mostrar algunos datos del tiempo de ejecución del mismo. Como ya se comentó previamente, se ha elegido Δ𝑡=0.01 𝑠 como intervalo de tiempo. Teniendo en cuenta esto, y que se han simulado los primeros 90 segundos para cada caso, se tiene un total de 9000 iteraciones. Estas 9000 iteraciones han necesitado de en torno a siete horas de tiempo real para ser ejecutado. Figura 14: Comparación de la profundidad máxima de la socavación para varios valores de la perturbación inicial y con los resultados del artículo (Kazaminezhad & Yeganeh-Bakhtiary, 2011) 16 Validación del modelo - 16 - Se muestra en la Figura 14 el desarrollo de la socavación para varios valores de la perturbación inicial (3, 2.5 y 2 mm) junto con los resultados presentados en el artículo (Kazaminezhad & Yeganeh-Bakhtiary, 2011). Se puede observar claramente la dependencia de los resultados con el valor de la amplitud de la perturbación inicial, cabe destacar lo siguiente: • Cuanto mayor es la amplitud de la perturbación inicial 𝛿, mayor es el tamaño final de la socavación. • Cuanto menor es 𝛿, menor es el tiempo que tarda la solución en estabilizarse, es decir, la socavación tarda menos en alcanzar un valor de la profundidad máxima estable. • De los distintos casos analizados, el que mejor aproxima sus resultados a los de la literatura es el caso con 𝛿=2.5 𝑚𝑚. La profundidad máxima esta sobreestimada en la parte inicial del experimento en el caso bajo estudio (los 25 primeros segundos aproximadamente), pero tras esto los resultados son muy similares a los obtenidos por Kazaminezhad y Yeganeg-Bakhtiary en su experimento. Se muestra en la Figura 15, al igual que se hizo para el caso en el que la perturbación inicial era de dos milímetros, la evolución del perfil de velocidades en función del tiempo para el caso en el que la perturbación inicial es de tres milímetros. Figura 15: Evolución del campo de velocidades con el tiempo (3 mm) Una vez analizados los resultados obtenidos con este modelo de análisis de socavación, se puede decir que los principales problemas que tiene son los siguientes: • Aunque el modelo predice la socavación e incluso llega a estabilizarse llegando a un máximo en términos de profundidad, no predice la formación de una acumulación de arena justo detrás del cilindro como, según la literatura, debería producirse. • El desarrollo de la socavación depende claramente de la perturbación inicial. Esto es un problema debido a que habría que decidir qué tipo de perturbación inicial debería ser impuesta para obtener buenos resultados. Por otro lado, es importante recordar que esta perturbación inicial es imprescindible cuando se emplea este modelo. t=45 seg t=90 seg t=10 seg t=2 seg Socavación en lechos de ríos y fondos marinos 17 - 17 - 2.2 MODELO MULTIFÁSICO 2.2.1 Definición del modelo Para intentar solucionar los problemas que presenta el modelo con malla dinámica presentado anteriormente, se presenta a continuación un segundo método para resolver el problema que se está estudiando. Este segundo método se basa en el uso de la aproximación Euler-Euler de dos fases, concretamente se ha utilizado el modelo euleriano. Además, se usarán las ecuaciones de Navier-Stokes para la definición de ambas fases, tanto la principal (agua) como la secundaria (arena). Por último, se empleará un modelo de cierre de turbulencia k-ε modificado. En cuanto a las fuerzas de interacción, se va a emplear un modelo de resistencia; mientras que la sustentación no se tendrá en cuenta. Existen artículos en la literatura que si emplean la sustentación en la resolución del modelo (Kazaminezhad & Yeganeh-Bakhtiary, 2011) (Yeganeh-Bakhtiary, Kazaminezhad, Etemad-Shahidi, Baas, & Chieng, 2011), habitualmente se emplea 𝐶𝐿=4/3. Sin embargo, el artículo en el que se basa esta sección para la validación del modelo no introduce la sustentación a las interacciones entre el agua y la arena (Zhao & Fernando, 2007), es por eso por lo que no se va a tener en cuenta en el presente documento. En primer lugar, se va a definir la geometría del modelo, así como las condiciones de contorno impuestas para su resolución. Se muestra en la Figura 16 un esquema de la geometría que se va a emplear para la resolución de este segundo modelo. Se tiene un dominio fluido de 2 x 0.4 m; mientras que, la zona la en la que se encuentra en el sedimento, tiene unas dimensiones de 2 x 0.1 m. En este segundo caso, es necesario definir tanto el dominio fluido (agua) como el sedimento (arena). Las condiciones de contorno que se han empleado son las siguientes: • En la entrada: se ha impuesto la condición de velocidad de entrada, pero teniendo en cuenta que el flujo de entrada debe ser exclusivamente de agua, por lo tanto, hay que imponer que la fracción volumétrica del sedimento sea nula, es decir, 𝛼𝑠=0. Para ello, se ha empleado de nuevo la fórmula mostrada en la Ecuación 4. La diferencia con el caso anterior es que ahora 𝑢𝑚á𝑥=0.31 𝑚/𝑠, por lo tanto 𝑢(𝑦)=0.31(𝑦 0.4)17. Este cambio se ha hecho para poder comparar los resultados del modelo con los obtenidos por Zhao y Fernando (Zhao & Fernando, 2007). Esta condición se ha impuesto, al igual que en el primer método, mediante una función definida por el usuario. Esta función es la presentada en el Anexo B. • En la salida: se ha empleado la condición de presión de salida. • En la zona de entrada de arena: se ha impuesto la condición de pared para evitar que exista un flujo tanto de agua como de arena en esa zona, el cual perturbaría los resultados. • En la zona de salida arena: se ha impuesto igualmente la condición de pared. • En la parte superior: se ha empleado la condición de simetría. • En la parte inferior: que es la interfase entre el fluido (agua) y la arena, se ha empleado la condición de pared. 24 Validación del modelo - 24 - Con todo lo anterior, se está en disposición de realizar el proceso iterativo para calcular el número de Stokes. Este proceso iterativo se ha realizado mediante su implementación en el software matemático MATLAB. Este proceso se muestra en el Anexo C. Una vez realizado el proceso iterativo, se obtiene un número de Stokes de 𝑆𝑡=0.0175. Por otro lado, el número de Reynols que se obtiene es de 𝑅𝑒𝑇=19.849, que está cercano a la unidad, por lo que el cálculo de 𝑘𝑇 es correcto. Teniendo en cuenta el valor obtenido, es decir, un valor bajo del número de Stokes (𝑅𝑒≪ 1), indica que los términos de inercia y gravitatorios son menos importantes en relación al término de fricción, de esta manera la velocidad de la partícula se va aproximando a la del fluido. El número de Stokes mide la relación de parada de una partícula frente a la dimensión característica del obstáculo, en este caso, el diámetro del cilindro. Por lo tanto, este bajo número de Stokes indica que las partículas siguen las líneas de corriente del fluido, es decir, no continúan en una línea recta impactando con el obstáculo mientras el fluido evita lo evita, que es lo que ocurriría si el valor del número de Stokes fuera bastante mayor que la unidad (𝑅𝑒≫1). La consecuencia de este bajo valor del número de Stokes en el caso bajo estudio es que seguramente es necesario emplear un valor más bajo para el paso de tiempo que el obtenido a partir del número de Courant. Se ha elegido teniendo en cuenta todo lo anterior un valor del intervalo de tiempo de Δ𝑡=0.001 𝑠. Este valor es cincuenta veces menor que el obtenido a partir del número de Courant. 2.2.3 Resultados Una vez presentadas la geometría del modelo y las ecuaciones que rigen su comportamiento, se está en disposición de obtener los resultados necesarios para la validación de este segundo modelo realizado en ANSYS. El artículo de Zhao y Fernando en el que se basa este método realiza una simulación de trescientos minutos. Para evitar la gran cantidad de tiempo que sería necesario para llevar a cabo esta simulación, se va ha emplear la siguiente fórmula, Ecuación 20, para la predicción de la socavación, y para así evitar tener que ejecutar los trescientos minutos de la simulación. 𝐹=𝐹𝑒𝑞[1−𝑒−𝑡/𝑇]=0.06[1−𝑒−𝑡/15] Ecuación 20: Función de predicción de la socavación Los parámetros empleados para el cálculo de esta función se han obtenido a partir de los resultados presentados en artículo de Zhao y Fernando (Zhao & Fernando, 2007). Los parámetros necesarios para su cálculo son los siguientes: • 𝐹𝑒𝑞: que es el valor de equilibrio máximo de la profundidad de la socavación, que en el artículo es aproximadamente 𝐹𝑒𝑞=0.06 𝑚. • 𝑇: que es el 63% del valor del tiempo para el cual se obtiene el equilibrio comentado en el punto anterior. El valor del tiempo de equilibrio es de aproximadamente 25 minutos, por lo que se va a emplear 𝑇=15 𝑚𝑖𝑛. Socavación en lechos de ríos y fondos marinos 25 - 25 - Se muestra en la Figura 19 una gráfica en la que se puede observar el proceso de determinación de cuál debería ser el tiempo total de ejecución: • Si se quiere obtener un 63% del valor de equilibrio, el tiempo debería ser aproximadamente 17 minutos. • Si se quiere obtener un 90% del valor de equilibrio, el tiempo debería ser aproximadamente 37 minutos. • Si se quiere obtener un 95% del valor de equilibrio, el tiempo debería ser aproximadamente 48 minutos. • Si se quiere obtener un 99% del valor de equilibrio, el tiempo debería ser aproximadamente 67 minutos. Figura 19: Predicción de la socavación Teniendo en cuenta lo anterior, la simulación durará 50 minutos para así tener una buena aproximación del valor de equilibrio de la socavación en el artículo de Zhao y Fernando. Por lo tanto, teniendo en cuenta que la simulación durará 50 minutos y que el paso de tiempo empleado es de Δ𝑡=0.001 𝑠, son necesarias tres millones de iteraciones para poder analizar el problema. Antes de ejecutar el problema que se ha descrito en la sección anterior, se ejecutó el problema suponiendo que sólo existe agua en el dominio y que la interfase entre la fase del agua y la fase en la que debería haber arena es una pared. El propósito de este proceso es tener un campo fluido totalmente desarrollado en términos de velocidad y turbulencia antes de añadir la fase de arena. Para ello se ejecutan veinte segundos de tiempo real con un paso de tiempo de Δ𝑡=0.1 𝑠. Estos veinte segundos son aproximadamente tres veces el tiempo que tarda el fluido en atravesar el dominio 𝑡=3𝐿𝑥 𝑢𝑥=3 2 𝑚 0.31 𝑚𝑠 ⁄≅20 𝑠𝑒𝑔. Se muestra en la Figura 20 el campo de velocidades tras estos veinte segundos. 0 0,01 0,02 0,03 0,04 0,05 0,06 0,07 020 40 60 80 100 120 140 160 180 200 220 240 260 280 300 Profundidad de la socavación [m] Tiempo [min] Resultados aproximados del artículo Error del 63% Error del 90% Error del 95% Error del 99% 26 Validación del modelo - 26 - Figura 20: Campo de velocidades inicial Una vez realizada esta primera parte de la resolución, se está en disposición de realizar el principal análisis en el cual se centra este segundo método. Como se muestra a continuación, los resultados obtenidos no son nada buenos. Se muestra en la Figura 21 la evolución de la fracción volumétrica de arena en el dominio fluido. El color rojo indica una concentración 0.63 es decir, que es la fase de arena; mientras que el color azul indica una concentración de arena nula y por lo tanto que la fase es agua. Como se puede observar, los resultados no tienen sentido físico ya que, aunque en un principio parece que la socavación comienza a simularse de una forma aparentemente correcta, llega un Figura 21: Evolución de la fracción volumétrica de arena en función del tiempo t=50 min t=5 min t=10 min Socavación en lechos de ríos y fondos marinos 27 - 27 - momento en el que la fase de arena comienza a desaparecer lo cual, teniendo en cuenta las condiciones impuestas en el dominio, no debería suceder. Para intentar solucionar este problema, se han intentado emplear otros modelos para simular tanto las interacciones entre ambas fases, como para definir la fase de sedimento. Se ha intentado, por ejemplo, emplear las funciones de Gidaspow (Gidaspow, Bezburuah, & Ding, 1992) para definir tanto las tensiones de corte sólidas como el modelo de resistencia entre las fases de agua y arena en vez de las funciones de Syamlal-Obrien empleadas en la obtención de los resultados presentados anteriormente. Los resultados obtenidos tras este cambio son muy similares a los obtenidos con Syamlal-Obrien y por lo tanto siguen siendo erróneos. Se ha decidido dejar este método y la corrección de sus errores para estudios posteriores. 2.3 COMPARACIÓN DE LOS MODELOS En esta sección se van a comparar los resultados obtenidos con ambos modelos, y se explicará que modelo se ha elegido para continuar realizando análisis de socavación. Debido a que el modelo multifásico no ha obtenido unos resultados con sentido físico, es evidente que se va a optar por continuar los análisis con un modelo monofásico. Pero es interesante comentar los aspectos que han llevado a esta decisión, así como las mejoras que se pueden realizar en el modelo monofásico para obtener unos mejores resultados. En primer lugar, se ha decidido optar por el modelo monofásico no sólo porque el modelo multifásico obtenga unos resultados sin sentido físico (ya que este problema seguramente se pudiera solucionar) sino por el tiempo de ejecución de ambos métodos: • El modelo monofásico emplea un paso de tiempo de Δ𝑡=0.03 𝑠𝑒𝑔 y necesita alrededor de dos minutos para comenzar a llegar a un valor de la socavación prácticamente constante, es decir, unas 4000 iteraciones. En tiempo real supondría alrededor de 3 o 4 horas de ejecución en un equipo con Intel Core i7-7700HQ con un procesador de 16.0 GB de RAM. • El modelo multifásico emplea un paso de tiempo mucho más pequeño (Δ𝑡= 0.001 𝑠𝑒𝑔) y necesita simular al menos 50 minutos. Como se comentó anteriormente esto supone 3 millones de iteraciones que, además de ser un número tres órdenes de magnitud superior al del primer modelo, es más lento ya que Fluent tarda más en analizar un caso multifásico que un monofásico. Con todo lo anterior, se tiene que le tiempo de ejecución de este modelo en el mismo equipo que el primero, necesitaría de varios meses para terminar su cálculo. Como es de esperar, los resultados mostrados en la sección anterior se obtuvieron realizando el análisis en un equipo con unas prestaciones muy superiores, aún así el modelo tardó más de una semana en simularse. Esta evidente diferencia en el tiempo de simulación de ambos modelos fue el motivo principal que llevó a elegir el modelo monofásico y no intentar resolver los problemas que el modelo multifásico tiene. 28 Validación del modelo - 28 - Por otro lado, aunque el modelo monofásico haya sido el elegido para continuar con los análisis tiene un gran problema, que no tiene en cuenta el transporte y la deposición de los sedimentos y, por lo tanto, no simula el montículo que se forma aguas abajo de la tubería. Por lo tanto, aunque se haya elegido el modelo monofásico se van a añadir las ecuaciones necesarias para que simule correctamente tanto el transporte como la deposición de los sedimentos. Socavación en lechos de ríos y fondos marinos 29 - 29 - 3. MODELO DEFINITIVO En esta sección se va a desarrollar el modelo definitivo que como se comentó anteriormente se basa al igual que el primer modelo presentado, en el uso de una malla dinámica. Además, se va a emplear un modelo monofásico en el que el único elemento es agua. Para mejorar los resultados obtenidos con el primer modelo, se van a añadir ecuaciones para simular el transporte, la deposición y el arrastre de los sedimentos para simular tanto la aparición de la socavación como la formación del montículo aguas abajo del obstáculo. Este último modelo, a diferencia de lo que se hizo en los dos casos anteriores, se va a realizar un modelo tridimensional para poder simular mejor el arrastre, la deposición y el transporte de los sedimentos. 3.1 ECUACIONES DEL MODELO Este modelo se basa en el uso de las ecuaciones de Navier-Stokes promediadas por el número de Reynols (RANS) similares a las empleadas en el modelo multifásico. Antes de presentar las ecuaciones es importante indicar el sentido de los ejes en referencia al dominio: • El eje x: va en la dirección del fluido. • El eje z: va en la perpendicular al plano que incluye el fondo del dominio. • El eje y: completando un sistema dextrógiro. 𝜕𝑢𝑗 𝜕𝑥𝑗=0 𝜕 𝜕𝑡(𝜌𝑢𝑖)+𝜕 𝜕𝑥𝑗(𝜌𝑢𝑖𝑢𝑗)=−𝜕𝑝+ 𝜕𝑥𝑖+𝜕𝜏𝑖𝑗 𝜕𝑥𝑗+𝜌(𝑠−1)𝑐𝑔𝑗 30 Modelo definitivo - 30 - 𝜕 𝜕𝑡(𝜌𝑐)+𝜕 𝜕𝑥𝑗(𝜌𝑐(𝑢𝑗−𝑤𝑗))=𝜕𝐽𝑐,𝑗 𝜕𝑥𝑗 Ecuación 21: Ecuaciones de transporte En la Ecuación 21 aparecen varios parámetros que necesitan ser definidos antes de continuar con el desarrollo del modelo: • La velocidad de sedimentación: 𝑤𝑗 cuyas componentes son nulas excepto aquella componente en la dirección perpendicular al plano que contiene el fondo del dominio que en este caso será, como se muestra más adelante, el eje z. Por lo tanto, la velocidad de sedimentación en la dirección del eje z se define como se muestra en la Ecuación 22. En dicha ecuación el parámetro 𝐷∗ es el diámetro medio del sedimento adimensional y se definió en la Ecuación 3. 𝑤𝑠=10𝜈 𝐷50[(1+0.01𝐷∗3)12−1] Ecuación 22: Definición de la velocidad de sedimentación • Tensión de cizalladura turbulenta de Reynols: 𝜏𝑖𝑗 cuya definición se muestra en la Ecuación 23, en esta ecuación aparece el térmico de viscosidad turbulenta que se define del mismo modo que en la Ecuación 6. 𝜏𝑖𝑗=−𝜌𝑢𝑖′𝑢𝑗′       =𝜇𝑇(𝜕𝑢𝑖 𝜕𝑥𝑗+𝜕𝑢𝑗 𝜕𝑥𝑖)−23𝛿𝑖𝑗𝜌𝑘 Ecuación 23: Definición de la tensión de cizalladura turbulenta de Reynols • Flujo escalar turbulento de Reynols: 𝐽𝑐,𝑗 cuya definición de muestra en la Ecuación 24. 𝐽𝑐,𝑗=−𝜌𝑐′𝑢𝑗′      =𝜇𝑇 𝜎𝑇 𝜕𝑐 𝜕𝑥𝑗 𝐷𝑜𝑛𝑑𝑒 𝜎𝑇=1 Ecuación 24: Definición del flujo escalar turbulento de Reynols El método de cierre turbulento que se va a emplear en la resolución del método es el k-ε estándar (Launder & Spalding, 1972). 𝜕 𝜕𝑡(𝜌𝑘)+𝜕 𝜕𝑥𝑗(𝜌𝑘𝑢𝑗)=𝜕𝐽𝑘𝑗 𝜕𝑥𝑗+𝑃+𝐺−𝜌ε 𝜕 𝜕𝑡(𝜌ε)+𝜕 𝜕𝑥𝑗(𝜌ε𝑢𝑗)=𝜕𝐽ε𝑗 𝜕𝑥𝑗+ε𝑘(𝐶ε,1𝑃+𝐶ε,3𝐺−𝐶ε,1𝜌ε) Ecuación 25: Ecuaciones del método k-ε estándar Los términos P y G de la Ecuación 25 son las dos partes de la energía cinética turbulenta. El primer término es la producción media debida a la tensión 𝑃=2𝜇𝑇(𝜕𝑢𝑖 𝜕𝑥𝑗+𝜕𝑢𝑗 𝜕𝑥𝑖)2 y el segundo es la debida a la flotabilidad 𝐺=−(𝜌𝑠𝑎𝑛𝑑 𝜌𝑤𝑎𝑡𝑒𝑟−1)𝐽𝑐,𝑗𝑔𝑗=−(𝑠−1)𝐽𝑐,𝑗𝑔𝑗. Por otro lado, la Socavación en lechos de ríos y fondos marinos 31 - 31 - constante de producción de energía cinética turbulenta debida a la flotabilidad es 𝐶ε,3= 𝑡𝑎𝑛ℎ(√𝑤2 𝑢2+𝑣2). Una vez presentadas las principales ecuaciones que rigen el comportamiento del modelo, se va a presentar a continuación la ecuación que define la formación de la socavación, la cual sustituye a la Ecuación 2 en el primer modelo monofásico, en el que simplemente se definía una velocidad vertical descendente del fondo del dominio a partir de una tensión crítica. En este caso, los movimientos de las partículas se basan en tres procesos (Reid, Bathurst, Carling, Walling, & Webb, 1997): • Transporte de sedimentos (𝑞𝑗): este primer término es el debido al movimiento de los sedimentos. Se define, al igual que el primer modelo, a partir de una tensión crítica basada en el parámetro de Shields. Se define en la Ecuación 26 el término que mide el transporte de sedimentos, en esta ecuación C es una constante del modelo cuyo valor es C=1.5. En dicha ecuación el parámetro 𝜃 es el parámetro de Shields que fue definido en la Ecuación 1 (Shields, 1936); mientras que 𝜃𝑐𝑟 es el parámetro de Shields crítico y se define en la Ecuación 27. 𝑞𝑗=𝑞0𝜏𝑗 |𝜏|−𝐶𝑞0𝜕ℎ 𝜕𝑥𝑗 𝑞0=12√𝑔(𝑠−1)𝐷50 3𝜃(𝜃−𝜃𝑐𝑟) 𝑠𝑖 𝜃>𝜃𝑐𝑟 Ecuación 26: Definición del término del transporte de sedimentos 𝜃𝑐𝑟=𝜃𝑐𝑟,0sin (𝛼+ϕ) sin (ϕ) Ecuación 27: Parámetro de Shields crítico • Deposición (D): este término simula el proceso de deposición de las partículas de sedimento transportadas que se produce en el proceso de socavación. Este proceso se modela como 𝐷=𝑤𝑠𝑐0, donde 𝑐0 es la concentración de sedimentos cerca del fondo y se define como se muestra en la Ecuación 28 (Whitehouse, 1998). 𝑐0={ 0.3 𝑠𝑖 𝜃>0.75 0.3 𝜃−0.25 0.75−0.25 𝑠𝑖 0.25< 𝜃<0.75 0 𝑠𝑖 𝜃<0.25 Ecuación 28: Definición de la concentración de sedimentos cerca del fondo • Arrastre (E): el último parámetro se muestra en la Ecuación 29 y es el que hace referencia a aquellas partículas del sedimento que se mueven debido al arrastre. El término 𝐽𝑐,3 hace referencia a la tercera componente del vector del flujo escalar turbulento de Reynols. 𝐸=𝐽𝑐,3 𝜌 Ecuación 29: Definición del arrastre de sedimentos 32 Modelo definitivo - 32 - Teniendo en cuenta esto tres términos, la ecuación que modela la formación de la socavación es la que se muestra en la Ecuación 30. En esta ecuación aparecen los tres términos presentados anteriormente. 𝑉=𝜕ℎ 𝜕𝑡=1 1−𝛾[𝜕𝑞𝑗 𝜕𝑥𝑗+𝐷−𝐸] Ecuación 30: Velocidad de desarrollo de la socavación Por último, se van a definir las condiciones de contorno que son necesarias imponer a la entrada del dominio fluido. A diferencia de los modelos anteriores no solo es necesario definir el perfil de velocidades a la entrada del dominio, también es necesario definir los perfiles de la energía cinética turbulenta y la escala de longitud de la turbulencia que ayudan a definir el modelo turbulento que en este caso será el modelo k-ε turbulento. En primer lugar, se va a mostrar el perfil de velocidades a la entrada del dominio. En la Ecuación 31 de muestra la definición de dicho perfil de velocidad donde aparece un parámetro el cual es necesario explicar antes de continuar con la explicación del modelo, la rugosidad del fondo del dominio fluido (𝑧0), este término se define como 𝑧0=𝑘𝑠 30, donde 𝑘𝑠 es la altura de la rugosidad, la cual se define como 𝑘𝑠=2.5𝐷50. Por otro lado, 𝑢0 es la velocidad del fluido sin perturbar. 𝑢(𝑧)=𝑚𝑖𝑛[𝑢∗ 𝜅ln(𝑧 𝑧0),𝑢0] 𝑢∗=𝜅𝑢0 𝑙𝑛(𝛿 𝑧0) Ecuación 31: Definición del perfil de velocidades En cuanto al perfil de la energía cinética a la entrada del dominio fluido, se muestra en la Ecuación 32 la definición de dicho perfil. En dicha ecuación 𝑇 es el nivel de turbulencia del flujo libre. 𝑘(𝑧)=𝑚𝑎𝑥[𝐶𝜇−1/2(1−𝑧𝛿)2𝑢∗2,1.5(𝑇𝑢0)2] Ecuación 32: Definición del perfil de la energía cinética Por último, se muestra en la Ecuación 33 la definición del perfil de la escala de longitud de la turbulencia. En esta ecuación la constante 𝐶𝜇 tiene el siguiente valor 𝐶𝜇=0.09. 𝑙(𝑧)=𝑚𝑎𝑥[𝜅𝑧(1+1.5𝑧𝛿)−1,𝐶𝜇𝛿] Ecuación 33: Definición del perfil de la escala de longitud Para definir la tasa de disipación de la turbulencia (ε), se emplea la definición que muestra en la Ecuación 34. ε(z)=Cμ3/4𝑘(𝑧) 𝑙(𝑧) Socavación en lechos de ríos y fondos marinos 33 - 33 - Ecuación 34: Definición de la tasa de disipación de la turbulencia Se tienen, por lo tanto, todas las ecuaciones necesarias para la resolución de este tercer modelo. Una vez mostradas dichas ecuaciones el siguiente paso es aplicarlas a distintos casos con distintas geometrías tridimensionales para analizar el desarrollo de la socavación en dichas situaciones. 3.2 APLICACIÓN DEL MODELO El primer caso que se va a analizar con este modelo tridimensional basado en el empleo de una malla dinámica es un campo fluido tridimensional en el que se encuentra un cilindro dispuesto verticalmente. Esta geometría puede simular un pilar (cilíndrico) de un puente que se encuentra situado en el cauce de un río, por lo que los resultados obtenidos serán muy interesantes desde el punto de vista del análisis de la estabilidad de los puentes y estructuras similares dispuestas en el cauce de un río, así como el estudio de la generación de la socavación en este tipo de estructuras que puede generar como se ha comentado previamente consecuencias catastróficas. En primer lugar, se va a definir la geometría del problema, así como las características de físicas del caso. Como se ha comentado previamente, el primer problema tridimensional que se va a analizar con este tercer modelo consiste en un prisma rectangular en cuyo lecho se encuentra un cilindro dispuesto verticalmente cuyo diámetro es 𝐷=0.2 𝑚. El dominio fluido, que se puede ver en la Figura 23, está compuesto por agua y sus dimensiones son de 4x2 metros en el fondo y de 0.5 metros de altura. Por otro lado, el centro del cilindro está situado a dos metros de la entrada del dominio, estando la salida del mismo, a cuatro metros del centro del cilindro; en cuanto a la dirección transversal al movimiento del fluido, el cilindro está situado justo en la mitad del mismo, es decir, a un metro tanto del lado derecho como del lado izquierdo. De forma análoga al primer modelo analizado, la zona del fondo del dominio simulará la interfase entre las fases de agua y arena. En cuanto a la altura del cilindro, Figura 22: Ejemplo puente https://pixnio.com/architecture/bridges/pillar-bridge- river-arch-construction-architecture-city# 40 Modelo definitivo - 40 - A la vista de estos resultados, se pensó en primer lugar que el error podría estar en el mallado empleado ya que podría ser demasiado simple, o no estar lo suficientemente refinado en la zona cercana al obstáculo. Para intentar solucionar dicho problema se definieron varios mallados con los que se pudiera evitar la aparición de dichas discontinuidades en los resultados del mallado dinámico. En primer lugar, se intentó obtener unos mejores resultados empleando un mallado mucho más estructurado que el empleado en la primera prueba. El mallado estructurado se muestra en las Figuras 30 y 31 y como se puede ver, es bastante más complejo y más refinado que el empleado previamente. Para realizar este mallado se ha dividido el dominio en tres partes para poder realizar un mallado estructurado. Este mallado está formado por 54000 elementos, cinco veces más que el mallado empleado en la primera ejecución, por lo que la ejecución del modelo será considerablemente más lenta. Por ejemplo, realizar 1000 iteraciones con el primer modelo requería unos 45 minutos; mientras que ejecutar 1000 iteraciones con este segundo mallado necesita aproximadamente 70 minutos, por lo que el tiempo total de ejecución aumenta notablemente. Además del problema evidente del tiempo, los resultados con este segundo mallado fueron muy similares a los obtenidos con el mallado anterior, es decir, llegaba un momento en el que la divergencia en los resultados provocaba un fallo en el mallado dinámico. Por lo tanto, se llegó a la conclusión de que el empleo de un mallado estructurado no era la solución del problema. Figura 29: Aparición de la divergencia en la malla dinámica Socavación en lechos de ríos y fondos marinos 41 - 41 - Figura 30: Mallado estructurado 3D Figura 31: Detalle mallado estructurado 3D 42 Modelo definitivo - 42 - A la vista de estos resultados, se optó por el empleo de un mallado formado por tetraedros. La intención del empleo de este mallado es poder utilizar el método remeshing para resolver el modelo. La ventaja de este método es que cada vez que se actualiza el mallado dinámico, la malla se adapta. Se muestra en la De nuevo, los resultados obtenidos empleando el mallado formado por tetraedros no fue el esperado ya que el mallado volvía a dar errores al ejecutarse. Teniendo en cuenta que el cambio de mallado no es la solución para el problema, se ha de afrontar la resolución de los problemas de resolución del mismo desde otro punto de vista. Se decidió que la solución al problema podría ser una reducción en el paso de tiempo. Para ello se emplea un criterio en el que se calcula un ‘tiempo difusivo’, este criterio es el que se conoce como ‘diffusive time stepping criteria’. Este criterio aparece recomendado en un artículo en el que se calcula un modelo similar (Solberg, Hjertager, & Bove, 2006). ∆𝑡𝑚á𝑥=(1−𝑛)𝑔(𝑠−1)𝐴𝑚í𝑛 24𝐶𝑢03=(1−0.4)∗9.81∗(2650 998.2−1)∗0.0001395 24∗1.5∗0.53=0.000302 𝑠 Ecuación 36: Definición del nuevo parámetro para definir el tiempo de paso máximo El valor obtenido para el paso de tiempo con esta segunda ecuación es ∆𝑡𝑚á𝑥=0.000302 𝑠. Teniendo en cuenta este valor, se ha decido utilizar ∆𝑡=0.00029 𝑠=0.29 𝑚𝑠 como paso de tiempo. Este paso de tiempo es de dos órdenes de magnitud menor que el obtenido a partir del número de Courant. Con el paso de tiempo obtenido y teniendo en cuenta que los procesos de socavación tardan horas, se puede prever que la ejecución de un modelo que necesite de varias horas de simulación y con un paso de tiempo tan pequeño, va a requerir mucho tiempo para obtener los resultados. 3.3 RESULTADOS En esta sección se muestran los resultados obtenidos con este último modelo. Figura 32: Mallado tetraédrico Socavación en lechos de ríos y fondos marinos 43 - 43 - En primer lugar, se muestran en la Figura 33 los resultados obtenidos para la coordenada z a los 0.5 segundos de simulación. En esta imagen se puede ver como la socavación se va formando de la forma esperada, es decir, en la zona del cilindro con la que choca la corriente incidente se va creando la socavación; mientras que en la zona lateral y aguas abajo del cilindro es donde comienzan a acumularse los sedimentos que han sido removidos de la zona delantera. Figura 33: Resultados para la coordenada z (vertical) a los 0.5 seg Figura 34: Resultados para los esfuerzos a los 0.5 seg 44 Modelo definitivo - 44 - En la Figura 34 se muestran, para el mismo instante de tiempo, los esfuerzos en el fondo del dominio. Los esfuerzos mínimos aparecen en las zonas delantera y trasera del cilindro, como era de esperar; mientras que los esfuerzos máximos aparecen en los laterales del mismo. Otro parámetro interesante para reforzar las conclusiones obtenidas de las figuras anteriores es la velocidad de desarrollo de la socavación definida en la Ecuación 30. En la Figura 35 se puede observar que en la zona delantera (zona azul) la velocidad es negativa, lo cual indica que se está produciendo una socavación en esa zona. Por otro lado, la parte del dominio en la que aparecen colores rojos y amarillos se está produciendo un aumento en la cota del fondo del dominio, es decir, en esas zonas se está produciendo una acumulación de los sedimentos que se van retirando de la zona delantera. Por último, y para poder ver de una forma más clara la formación de la socavación y de la acumulación de sedimentos en torno al cilindro, se muestran en la Figura 36 las cotas de la coordenado z del fondo del dominio. Al igual que se mostró previamente, es en la zona delantera del cilindro donde se produce la socavación; mientas que aguas abajo y en los laterales del obstáculo comienzan a acumularse los sedimentos. Debido a que el tiempo de simulación es muy pequeño, el proceso de socavación no es demasiado grande en términos numéricos, aunque si se puede ver que la simulación es correcta. El máximo de la coordenada z que se obtiene a los 0.5 segundos de socavación es de 𝑧𝑚á𝑥= −1.2 𝑚𝑚; mientras que el mínimo es 𝑧𝑚í𝑛=−2.38 𝑚𝑚. Figura 35: Resultados para la velocidad de desarrollo de la socavación a los 0.5 seg Socavación en lechos de ríos y fondos marinos 45 - 45 - Como conclusión a los resultados obtenidos con este último modelo se puede decir que, a pesar de que la simulación está todavía en una etapa primaria, los resultados que se obtienen tienen sentido físico y son lo que se podía esperar antes de realizar la simulación. Por lo tanto, se tiene un modelo robusto para la simulación de la socavación que era el principal objetivo del presente documento. El principal inconveniente que presenta este modelo es claramente el gran coste de tiempo que requiere su simulación, aunque este tiempo depende en gran medida del mallado elegido, de las condiciones del problema y sobre todo del equipo informático empleado para su resolución. Por lo que la gran cantidad de tiempo necesaria para su ejecución en este primer cálculo podría reducirse. Una vez se refine el funcionamiento de este modelo, se podría aplicar a distintos casos y con distintas geometrías realizando pequeños cambios en el mismo. El poder aplicar este mismo método para distintos casos realizando simplemente pequeños cambios, es decir, su flexibilidad, una de las principales ventajas de este modelo. Otra de las principales ventajas de este modelo es que, a diferencia del primer modelo desarrollado, simula la acumulación de sedimentos en torno al obstáculo debido a que tiene en cuenta el transporte, el arrastre y la deposición de los sedimentos. Figura 36: Cotas para la coordenada z en el fondo del dominio a los 0.5 seg 46 Conclusiones - 46 - 4. CONCLUSIONES Una vez realizados análisis de los procesos de socavación para distintos casos empleando para ellos distintos modelos, se ha llegado a una serie de conclusiones sobre el tema estudiado las cuales van a ser expuestas en esta sección. Tras todos los análisis realizados queda claro que la creación de un modelo robusto y eficiente para el estudio de los procesos erosivos y el análisis de la formación de la socavación es una tarea compleja debido a la dificultad a la hora de elegir un modelo adecuado para su resolución. Tanto si se decide resolver el problema con un modelo monofásico como con un modelo bifásico, es importante decidir que modelo de cierre se va a emplear o que paso de tiempo es el adecuado para obtener unos buenos resultados sin que el tiempo total de la simulación sea demasiado grande. En primer lugar, se van a comentar las principales conclusiones alcanzadas tras el desarrollo del modelo bifásico para la resolución de problemas bidimensionales. La principal ventaja de este método frente al resto de métodos estudiados es su velocidad, ya que la cantidad de tiempo requerida para la simulación es mucho menor que la requerida con los otros métodos. Por otro lado, los resultados obtenidos son bastante buenos en términos de predicción de la socavación que se forma bajo el obstáculo, ya que son bastante similares a los obtenidos en la literatura. La principal desventaja del método empleado es que no predice el transporte y la deposición de sedimentos y, por lo tanto, no predice la acumulación de arena que se produce aguas abajo des obstáculo. Esto se debe a que el modelo desarrollado no tiene en cuenta este tipo de efectos y por lo tanto no puede simularlos. El otro inconveniente de este modelo es la dependencia del resultado con el tamaño de la perturbación inicial que se impone en el problema. Esto provoca que cuando se analiza un caso hay que tener mucho cuidado en la elección de esa perturbación inicial ya que el resultado dependerá de ese valor. En cuanto al modelo bifásico desarrollado para el análisis de casos bidimensionales, no se ha conseguido llegar a un modelo que prediga correctamente la formación de la socavación. Aunque no se haya llegado a una solución correcta, es evidente que un problema importante de este tipo de modelos es el tiempo de ejecución debido principalmente a tres motivos: • El bajo paso de tiempo (Δ𝑡=0.001 𝑠) que hay que emplear para la resolución del modelo debido al bajo número de Stokes obtenido. Socavación en lechos de ríos y fondos marinos 47 - 47 - • El proceso de socavación es lento por lo que hay que analizar al menos 50 minutos para obtener una buena aproximación del máximo valor que alcanza el proceso. Si se quisiera analizar el proceso completo, hasta alcanzar el valor de equilibrio, habría que analizar varias horas del proceso erosivo. • Al emplear un modelo bifásico con un modelo viscoso por fase la resolución del problema se ralentiza. Por otro lado, los principales problemas del modelo están en la definición de las propiedades de la fase secundaria, la fase de arena, así como la elección las interacciones entre ambas fases, especialmente que métodos se emplean para definirlas. Por último, se van a comentar las principales conclusiones alcanzadas sobre el último modelo empleado. Este último modelo es similar al primero ya que basa su funcionamiento en el empleo de una malla dinámica para simular el proceso erosivo. A diferencia del primer modelo, este último se emplea para simular casos tridimensionales y tiene en cuenta los procesos de arrastre y deposición de sedimentos, por lo que simula la acumulación de sedimentos que se produce aguas abajo del obstáculo, en este caso, un cilindro. Aunque este modelo también necesita una gran cantidad de tiempo para su ejecución, los resultados son bastante buenos y simulan tanto la socavación como la acumulación de la arena. Para terminar con esta sección se puede decir que se han obtenido dos herramientas para analizar los procesos erosivos, ambas basadas en el empleo de una malla dinámica: • La primera de ellas sirve para el análisis de casos bidimensionales, su principal ventaja es el tiempo de resolución; mientras que el principal inconveniente es que no predice la acumulación de sedimentos aguas abajo del cilindro. • La segunda herramienta se emplea para el análisis de casos tridimensionales, su principal ventaja es la calidad de los resultados; mientras que el principal inconveniente es que necesita mucho más tiempo para ejecutarse. De forma general, podría ser interesante emplear el primer modelo para realizar una primera aproximación para un caso bidimensional ya que, aunque no se tenga en cuenta el arrastre y la deposición de sedimentos, es una herramienta bastante rápida. Una vez que se tengan estos resultados y se vea que es interesante seguir con los análisis, se podría realizar un estudio de un problema tridimensional con el último método desarrollado en el que ya si se prediga tanto la socavación en torno al obstáculo como la acumulación de sedimentos. 48 Líneas de estudio futuras - 48 - 5. LÍNEAS DE ESTUDIO FUTURAS Debido a la gran amplitud de este cambio de estudio, quedan muchas líneas de estudio abiertas para ser estudiadas posteriormente las cuales se van a tratar en esta sección. En cuanto al primer método analizado queda abierta la posibilidad de analizar casos con distintas geometrías, aunque quizás no se demasiado interesante ya que no tiene en cuenta el arrastre y la deposición de sedimentos. El segundo método empleado es la línea de estudio que evidentemente más queda por desarrollar en un futuro ya que no se ha conseguido en el actual estudio que este modelo obtenga unos resultados correctos. El error de este modelo podría estar en algunos de los métodos empleados para definir la fase de arena; o en alguno de los métodos empleados para definir las interacciones entre ambas fases. El error podría estar también en el paso de tiempo empleado, aunque esto parece menos probable porque el paso de tiempo que se ha utilizado ya es bastante pequeño y se han tenido en cuenta tanto el número de Courant como el número de Stokes para su cálculo. Además, es un tiempo análogo al empleado en estudios similares. En cuanto al tercer método sería interesante, ya que se ha obtenido una herramienta que obtiene buenos resultados, analizar casos para distintas geometrías del obstáculo (más complejas que la analizada en este documento) y para distintas condiciones de entrada al dominio. En definitiva, aunque se han conseguido dos buenos métodos para predecir la socavación, todavía queda bastante por analizar en este campo. Además de todo lo comentado previamente en esta sección, queda abierto para un análisis futuro, el estudio de nuevos métodos para predecir este tipo de procesos. Socavación en lechos de ríos y fondos marinos 49 - 49 - 6. REFERENCIAS Brennen, C. (2005). Fundamentals of Multiphase Flow. Chapter 13.6. Pasadena, CA, Estados Unidos: Cambridge University Press. doi:10.1017/CBO9780511807169 Brørs, B. (1999). Numerical modeling of flow and scour at. Journal of Hydraulic Engineering, 125(5), 511–523. C. van Rijn, L. (1984). Sediment Pick-Up Functions. Journal of Hydraulic Engineering. Holanda, 1494-1502. doi:10.1061/(ASCE)0733-9429(1984)110:10(1494) Cheng, L., Yeow, K., Zhang, Z., & Teng, B. (2009). Three-dimensional scour below pipelines in steady currents. Coastal Engineering, 56(5-6), 577-590. doi:10.1016/j.coastaleng.2008.12.004 Courant, R., Friedrichs, K., & Lewy, H. (1967 [1928]). On the Partial Difference Equations of Mathematical Physics. Berlin: IBM Journal of Research and Development. doi:10.1147/rd.112.0215 Dalla Valle, J. (1948). Micromeritics: the technology of fine particles. Londres, Inglaterra: Pitman Publishing Corporation. Ding, J., & Gidaspow, D. (1990). A Bubbling Fluidization Model Using Kinetic Theory of Granular Flow. AIChE Journal, 36 (4), 523-538. doi:10.1002/aic.690360404 Drew, D., & Lahey, R. (1993). In Particulate Two-Phase Flow. Boston, MA, Estados Unidos: Butterworth-Heinemann. Ganesh Elapolu, P., Majumdar, P., Lottes, S., & Kostic, M. (2012). Development of a Three- Dimensional Iterative Methodology Using a Commercial CFD Code for Flow Scouring Around Bridge Piers. ASME 2012 Heat Transfer Summer Conference (págs. 1073-1085). Rio Grande, Puerto Rico, Estados Undios: ASME. doi:10.1115/HT2012-58491 Gidaspow, D., Bezburuah, R., & Ding, J. (1992). Hydrodynamics of Circulating Fluidized Beds, Kinetic Theory Approach. In Fluidization VII, Proceedings of the 7th Engineering Foundation Conference on Fluidization, (págs. 75-82). Birsbane, Australia. doi:5896246 Guo, J. (2002). Advances in Hydraulics and Water Engineering (Vol. 2). Singapur: Hunter Rouse and Shields Diagram. doi:doi.org/10.1142/5041 56 Anexo C: Cálculo del número de Stokes - 56 - ANEXO C: CÁLCULO DEL NÚMERO DE STOKES En este anexo se muestra el proceso iterativo realizado en MATLAB para el cálculo del número de Stokes. En primer lugar, se definen las variables necesarias para realizar el cálculo. Tras esto se calcula el tiempo de residencia como la longitud característica entre la velocidad característica del problema, y ese valor es el que se toma como condición inicial para el proceso iterativo. Una vez hecho esto, se realiza el proceso iterativo en el que, en primer lugar, se calcula el número de Reynols a partir del tiempo de relajación, tras esto se calcula la constante 𝑘𝑇 y por último se calcula el nuevo valor del tiempo de relajación a partir de ese número de Reynols. El proceso iterativo sigue hasta que la tolerancia (definida como la diferencia entre el tiempo de relajación de una iteración y el de la anterior) sea menor que un cierto valor, concretamente, 𝑇𝑂𝐿=1𝑒−15. Una vez se consigue un valor que cumpla con la tolerancia especificada, se calcula el número de Stokes a partir de la definición mostrada en el documento. %Definición de las variables rho_water=998.2; rho_sand=2650; beta=rho_sand/rho_water; mu=0.001003; nu=mu/rho_water; g=9.81; dp=0.00036; U=0.31; L=0.1; %Cálculo del tiempo de residencia tr=L/U; tp_aux=tr; TOL=1; %Proceso iterativo while TOL>1e-15 VT=tp_aux*g; Socavación en lechos de ríos y fondos marinos 57 - 57 - ReT=VT*dp/nu; kT=1+0.15*ReT^(2/3); tp=(beta-1)*dp^2/kT/18/nu; TOL=abs(tp_aux-tp); tp_aux=tp; end tp; VT=tp*g; ReT=VT*dp/nu; %Cálculo del número de Stokes St=tp/tr 58 Anexo D: UDF empleada para el modelo 3D - 58 - ANEXO D: UDF EMPLEADA PARA EL MODELO 3D En este anexo se muestra la UDF empleada para la resolución del modelo tridimensional. Esta UDF se ha dividido en cuatro partes para simplificar la explicación de la misma. A continuación, se van a explicar cada una de las partes: • En la primera parte (la enmarcada en color azul) es en la que se calculan los perfiles necesarios para definir las condiciones de entrada al dominio a partir de las Ecuaciones 31, 32 y 34. En primer lugar, se calcula el perfil de velocidades a la entrada del dominio, tras lo cual se calculan los perfiles de energía cinética turbulenta y la tasa de disipación turbulenta necesarios para definir el modelo viscoso que se va a emplear. • En la segunda parte (la enmarcada en color rojo) se definen los cuatro términos fuentes necesarios para la definición de las condiciones del fluido que son: el término fuente para la UDS-0, el momento en torno al eje z, la energía cinética turbulenta y la tasa de disipación turbulenta. En esta parte también se define la difusividad de la fase de agua. • En la tercera parte (la enmarcada en color verde) se definen los dos DEFINE_ADJUST necesarios para la definición de los distintos parámetros necesarios para la correcta ejecución de Fluent. En esta parte, también se define (mediante un DEFINE_PROFILE) la concentración de sedimentos cerca del fondo del dominio. • Por último, en la cuarta parte (la enmarcada en color amarillo) se calcula la velocidad vertical del fondo del dominio (dhdt) la cual se utiliza posteriormente para la definición del movimiento de la malla dinámica (utilizando el DEFINE_GRID_MOTION). /* Perfil de velocidades a la entrada */ DEFINE_PROFILE(Ux_inlet,thread,position) { real x[ND_ND]; real z; real aux1v; Socavación en lechos de ríos y fondos marinos 59 - 59 - face_t f; begin_f_loop(f, thread) { F_CENTROID(x,f,thread); z = x[2]; aux1v = 0.0488*log(48008.0*z); if (aux1v <= 0.50) F_PROFILE(f,thread,position) = aux1v; else F_PROFILE(f,thread,position) = 0.50; } end_f_loop(f,thread) } /* Perfil de energía cinética turbulenta a la entrada */ DEFINE_PROFILE(K_inlet,thread,position) { real x[ND_ND]; real z; real aux1k; real aux2k; face_t f; begin_f_loop(f, thread) { F_CENTROID(x,f,thread); z = x[2]; aux1k = 0.00133*(1.0 - 2.0*z)*(1.0 - 2.0*z); aux2k = 0.00015; if (aux1k >= aux2k) F_PROFILE(f,thread,position) = aux1k; else F_PROFILE(f,thread,position) = aux2k; } end_f_loop(f,thread) } /* Perfil de tasa de disipación turbulenta a la entrada */ DEFINE_PROFILE(EPS_inlet,thread,position) { real x[ND_ND]; real z; real aux1l; real aux2l; real aux1k; 60 Anexo D: UDF empleada para el modelo 3D - 60 - real aux2k; real l; real k; face_t f; begin_f_loop(f, thread) { F_CENTROID(x,f,thread); z = x[2]; aux1l = 0.41*z/(1.0+3.0*z); aux2l = 0.045; aux1k = 0.00133*(1.0 - 2.0*z)*(1.0 - 2.0*z); aux2k = 0.00015; if (aux1k >= aux2k) k = aux1k; else k = aux2k; if (aux1l <= aux2l) l = aux1l; else l = aux2l; F_PROFILE(f,thread,position) = 0.164*pow(k,1.5)/l; } end_f_loop(f,thread) } /* Definición del término fuente para la UDS-0 */ DEFINE_SOURCE(Ceq_src,c,t,dS,eqn) { real source; real x[ND_ND]; C_CENTROID(x,c,t); source = 0.035*C_R(c,t)*C_UDSI_G(c,t,0)[2]; dS[eqn] = 0; /* UDM-0: fuente de la UDS-0 */ C_UDMI(c,t,0) = 0.035*C_R(c,t)*C_UDSI_G(c,t,0)[2]; return source; } /* Definición de la difusividad */ Socavación en lechos de ríos y fondos marinos 61 - 61 - DEFINE_DIFFUSIVITY(Ceq_diff,c,t,i) { return C_MU_T(c,t) / 1.0; } /* Definición del término fuente para el momento en torno al eje z */ DEFINE_SOURCE(Zmom_src,c,t,dS,eqn) { real source; real x[ND_ND]; C_CENTROID(x,c,t); source = -C_R(c,t)*C_UDSI(c,t,0)*9.81*(2.66-1.00); dS[eqn] = 0; /* UDM-1: fuente de para el momento en torno al eje z */ C_UDMI(c,t,1) = -C_R(c,t)*C_UDSI(c,t,0)*9.81*(2.66-1.00); return source; } /* Definición del término fuente para la energía cinética turbulenta */ DEFINE_SOURCE(Keq_src,c,t,dS,eqn) { real k_sc; real x[ND_ND]; C_CENTROID(x,c,t); k_sc = 9.81*(2.66 - 1.00)*(C_MU_T(c,t)/1.0)*C_UDSI_G(c,t,0)[2]; dS[eqn] = 0; /* UDM-2: fuente de para la energía cinética turbulenta */ C_UDMI(c,t,2) = 9.81*(2.66 - 1.00)*(C_MU_T(c,t)/1.0)*C_UDSI_G(c,t,0)[2]; return k_sc; } /* Definición del término fuente para la tasa de disipación turbulenta */ DEFINE_SOURCE(EPSeq_src,c,t,dS,eqn) { real G; real aux1; real aux2; real c_3eps; real eps_sc; real x[ND_ND]; 62 Anexo D: UDF empleada para el modelo 3D - 62 - C_CENTROID(x,c,t); G = 9.81*(2.66 - 1.00)*(C_MU_T(c,t)/1.0)*C_UDSI_G(c,t,0)[2]; aux1 = (C_W(c,t)*C_W(c,t))/(C_U(c,t)*C_U(c,t) + C_V(c,t)*C_V(c,t) + 0.00001); aux2 = pow(aux1,0.5); c_3eps = tanh(aux2); eps_sc = G*c_3eps; dS[eqn] = 0; /* UDM-3: fuente de para la tasa de disipación turbulenta */ C_UDMI(c,t,3) = G*c_3eps; /* UDM-4: 𝑪𝟑𝜺*/ C_UDMI(c,t,4) = c_3eps; return eps_sc; } DEFINE_ADJUST(parametros_adjust,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 PI=3.14159265; real tau, shields_sce, shields_be; real alpha, beta, gamma; real omega_x, omega_y, omega_z; real aux1, aux2, aux3, qb; real qb_x, qb_y, qb_z; int wall_id=9; ft = Lookup_Thread(domain,wall_id); begin_f_loop(f,ft) { c=F_C0(f,ft); ct=THREAD_T0(ft); Socavación en lechos de ríos y fondos marinos 63 - 63 - 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); alpha=acos(-area[0]/NV_MAG(area))*180.0/PI; beta=acos(-area[1]/NV_MAG(area))*180.0/PI; gamma=acos(-area[2]/NV_MAG(area))*180.0/PI; omega_x=acos(-ts_force[0]/NV_MAG(ts_force))*180.0/PI; omega_y=acos(-ts_force[1]/NV_MAG(ts_force))*180.0/PI; omega_z=acos(-ts_force[2]/NV_MAG(ts_force))*180.0/PI; shields_sce=0.2456*0.1890*cos(gamma*PI/180.0); /*En las aux_i aparece el ángulo de fricción interna*/ aux1=pow((0.1890/tau)*sin(gamma*PI/180.0)/tan(25.0*PI/180.0),2); aux2=2.0*0.1890*cos(omega_z*PI/180.0)/tau; aux3 = (1 + aux1 - aux2)/tan(25.0*PI/180.0); shields_be = 0.2456*tau*pow(aux3,0.50); if (shields_be < shields_sce) qb = 0; else qb=0.0001914*(shields_be-shields_sce)*pow(shields_be,0.5); /*qb_x=qb*cos(omega_x*PI/180.0) - qb_x*1.5*tan((90-alpha)*PI/180.0);*/ /*qb_y=qb*cos(omega_y*PI/180.0) - qb_y*1.5*tan((90-beta)*PI/180.0);*/ /*qb_z=qb*cos(omega_z*PI/180.0) - qb_z*1.5*tan(gamma*PI/180.0);*/ qb_x=qb*cos(omega_x*PI/180.0); qb_y=qb*cos(omega_y*PI/180.0); qb_z=qb*cos(omega_z*PI/180.0); /* UDM-5: esfuerzo cortante*/ C_UDMI(c,ct,5)=tau; /* UDM-6: ángulo entre la normal al suelo y el eje X*/ C_UDMI(c,ct,6)=alpha; /* UDM-7: ángulo entre la normal al suelo y el eje Y*/ C_UDMI(c,ct,7)=beta; /* UDM-8: ángulo entre la normal al suelo y el eje Z*/ C_UDMI(c,ct,8)=gamma; /* UDM-9: ángulo entre el cortante y el eje X*/ C_UDMI(c,ct,9)=omega_x; /* UDM-10: ángulo entre el cortante y el eje Y*/ C_UDMI(c,ct,10)=omega_y; /* UDM-11: ángulo entre el cortante y el eje Z*/ C_UDMI(c,ct,11)=omega_z; /* UDM-12: cortante crítico*/ C_UDMI(c,ct,12)=shields_sce; 64 Anexo D: UDF empleada para el modelo 3D - 64 - /* UDM-13: exceso de cortante sobre el crítico*/ C_UDMI(c,ct,13)=shields_be; /* UDM-14: flujo de sedimentos*/ C_UDMI(c,ct,14)=qb; /* UDM-15: flujo de sedimentos en la dirección del eje X*/ C_UDMI(c,ct,15)=qb_x; /* UDM-16: flujo de sedimentos en la dirección del eje Y*/ C_UDMI(c,ct,16)=qb_y; /* UDM-17: flujo de sedimentos en la dirección del eje Z*/ C_UDMI(c,ct,17)=qb_z; } end_f_loop(f,ft); } /* Concentración de sedimentos cerca del fondo*/ DEFINE_PROFILE(C0_bed,thread,position) { Domain *domain = Get_Domain(1); Thread *ft, *ct; face_t f; cell_t c; real tau_bed; real shields; int wall_id=9; 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); shields = tau_bed/(1000.0*9.81*1.66*0.00025); if (shields < 0.25) F_PROFILE(f,thread,position) = 0.0; else if (shields < 0.75) F_PROFILE(f,thread,position) = 0.60*(shields - 0.25); else F_PROFILE(f,thread,position) = 0.30; } end_f_loop(f,ft); } DEFINE_ADJUST(dqi_adjust,dom) { Domain *domain = Get_Domain(1); Thread *t; Socavación en lechos de ríos y fondos marinos 65 - 65 - cell_t c; face_t f; domain = Get_Domain(1); thread_loop_c (t,domain) { begin_c_loop (c,t) { /* UDS-1: flujo de sedimentos en la dirección del eje X*/ C_UDSI(c,t,1) = C_UDMI(c,t,15); /* UDS-2: flujo de sedimentos en la dirección del eje Y*/ C_UDSI(c,t,2) = C_UDMI(c,t,16); /* UDS-3: flujo de sedimentos en la dirección del eje Z*/ C_UDSI(c,t,3) = C_UDMI(c,t,17); } end_c_loop (c,t) } } /*Cálculo de la velocidad vertical del fondo*/ DEFINE_EXECUTE_AT_END(dhdt_calc) { Domain *domain = Get_Domain(1); Thread *ft, *ct; face_t f; cell_t c; real D; real E; real DmE; real dqxdx; real dqydy; real dqzdz; real dhdt; int wall_id=9; ft = Lookup_Thread(domain,wall_id); begin_f_loop(f,ft) { c=F_C0(f,ft); ct=THREAD_T0(ft); D = 0.035*C_UDSI(c,ct,0); E = (C_MU_T(c,ct)/1.0)*C_UDSI_G(c,ct,0)[2]/C_R(c,ct); DmE = D - E; /* UDM-18: deposición menos arrastre de sedimentos*/ C_UDMI(c,ct,18)=DmE; dqxdx = C_UDSI_G(c,ct,1)[0]; dqydy = C_UDSI_G(c,ct,2)[1];