XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - XLIV Jornadas de AUTOMÁTICA Z A R A G O Z A2 0 2 3 - Fuente: ZT Gatha Font (gratis para uso comercial) Logo vertical Logo horizontal Ejemplo en camisetas Propuesta de logo para las XLIV Jornadas de automática Autor: Ignacio Cuiral-Zueco https://jautomatica.es Métodos de Simulación de Fluidos aplicados a Lagunas y Embalses Ferrero-Losada, Samuel.a,∗, López-Orozco, José Antonio.a, Besada-Portas, Eva.a, Carazo-Barbero, Gonzalo.b, Risco-Martı́n, José Luis.b aDepartamento de Arquitectura de Computadores y Automática, Universidad Complutense de Madrid, Plaza Ciencias, nº1, 28040, Madrid España. bDepartamento de Arquitectura de Computadores y Automática, Universidad Complutense de Madrid, C/ Profesor José Garcı́a Santesmases, nº9, 1, 28040, Madrid España. To cite this article: Ferrero-Losada, S., López-Orozco, J.A., Besada-Portas, E., Carazo-Barbero, G., Risco-Martı́n, J. L. 2023. Fluid simulation methods applied to lakes and reservoirs. XLIV Jornadas de Automática, 393-398 https://doi.org/10.17979/spudc.9788497498609.393 Resumen Existe una gran variedad de métodos para simular un fluido incompresible como el agua. Estos se pueden clasificar en diferentes tipos según utilicen partı́culas, una malla o una combinación de ambos como soporte de la simulación; o según el enfoque en que se aborde la incompresibilidad del fluido (débilmente compresible o realmente incompresible). Con esta variedad de opciones, es necesario comparar los métodos para asegurar la utilización del más adecuado, considerando las ventajas y desventajas de cada uno. Por ello, se analizan de manera general algunos de estos métodos con un escenario especı́fico en mente: un lago con una entrada y una salida de agua en régimen estacionario. El objetivo es discernir cuál de estos métodos es mejor para llevar a cabo dicha simulación, o tiene el menor número de problemas en cuanto a las condiciones de contorno, aplicación de fuerzas externas, o inestabilidades numéricas. Finalmente, se presenta un caso de prueba sencillo empleando la opción considerada más adecuada. Los resultados de este artı́culo se emplearán en futuros trabajos en el estudio de blooms de cianobacterias en dichos cuerpos acuáticos. Palabras clave: Hidroinformática, Modelado e identificación de sistemas ambientales, Administración de Recursos Naturales. Fluid simulation methods applied to lakes and reservoirs. Abstract There are various methods to simulate an incompressible fluid such as water. These can be grouped into different types according to their simulation framework (grid-based, particle-based, or mixed) or how they approach the fluid’s incompressibility (weakly compressible or truly incompressible). Facing such variety to choose from, a comparison between them becomes a must to ensure the method used in the simulation scenario is adequate, considering their advantages and disadvantages in the different aspects of the simulation. A general analysis of some of these methods will be done with a specific scenario in mind: a stationary lake with one water entrance and one exit. The aim will be to discern which of the studied methods is better suited to carry out this simulation with the minor problems regarding boundary conditions, external forces treatment, or numerical instabilities. Finally, a simple test case for the more adequate method is presented. This work will be ultimately used in the study of lake cyanobacteria blooms. Keywords: Hydroinformatics, Modeling and identification of environmental systems, Natural Resources management. 1. Introducción Cuando abordamos un problema de dinámica de fluidos me- diante simulación, una de las decisiones más relevantes es la elección del método. Esto se debe a que definirá no solo la for- ma en que el fluido se representará en el estudio, sino también las condiciones de contorno a implementar, las magnitudes in- volucradas en la simulación y el tipo de análisis que se podrán llevar a cabo con los resultados. Debemos tener en cuenta que existen varios tipos de fluidos, con distintas propiedades fı́sicas, que deben ser precisamente plasmadas para emular adecuada- ∗Autor para correspondencia: saferr03@ucm.es Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0) 393 mente su comportamiento. Estas propiedades suelen depender de factores internos y ambientales. Por ello, se debe compren- der detalladamente dichas propiedades y comportamientos in- dispensables en la simulación para lograr el nivel de realismo deseado. Además, es importante considerar las posibles situa- ciones que puedan surgir durante la computación y las limita- ciones de ciertos métodos para simular adecuadamente. Actualmente existen varios programas comerciales, como COMSOL™ (Multiphysics 1998), y de código abierto, como OpenFOAM (Weller et al. 1998), con la capacidad para simu- lar flujos turbulentos. Por lo tanto, uno podrı́a preguntarse si es realmente necesario implementar un simulador propio. La ne- cesidad surge porque la mayorı́a de estos programas requieren una gran cantidad de parámetros de entrada para funcionar, lo que los hace difı́ciles de sintonizar. Además, en muchos casos proporcionan resultados más precisos de lo necesario, lo que aumenta el coste computacional y el tiempo de ejecución in- necesariamente. En el ámbito de los simuladores numéricos hi- drodinámicos de cuerpos acuáticos (p.e. Delft3D, Roelvink and Van Banning 1995), el problema es aún mayor. Estos tienden a incorporar efectos termodinámicos y atmosféricos, lo que hace que un método inicialmente sencillo sea demasiado complicado y dificulte mucho la obtención de información para otros aspec- tos del problema, como la optimización de rutas y la predicción de comportamientos. Debido a esta situación, se inicia este pro- ceso de análisis con el objetivo de desarrollar un método más fácil de parametrizar y con un menor tiempo de ejecución que los mencionados anteriormente. Nuestro objetivo es estudiar el comportamiento dinámico de cuerpos de agua semejantes a lagos o embalses, donde el régi- men turbulento es predominante, limitándonos a las corrientes estacionarias. Este tipo de movimiento se caracteriza por una variación caótica en la presión y velocidad de flujo, lo que com- plica significativamente la simulación y restringe los métodos que se pueden utilizar para llevarla a cabo. Cabe destacar que el agua es un fluido incompresible en la mayorı́a de los casos, lo que debe ser considerado al momento de seleccionar el méto- do apropiado para este estudio. Aunque la elección del método depende del régimen de movimiento, la incompresibilidad sólo establece una condición a mantener durante toda la simulación que modifica principalmente el modelo de ecuaciones a utilizar, introduciendo una relación aún más restrictiva entre presión y velocidad. El objetivo de obtener la velocidad del flujo estacionario del agua es para utilizarlo como base para determinar el mo- vimiento de las cianobacterias en el cuerpo de agua en traba- jos futuros. Las cianobacterias, también conocidas como algas verde-azuladas, son organismos unicelulares capaces de reali- zar la fotosı́ntesis, fijar CO2 y generar O2. Algunas especies producen cianotoxinas que pueden ser peligrosas para la salud cuando se reproducen en exceso creando afloraciones (llamadas Blooms en inglés). Estas cianobaterias se mueven mediante fac- tores biológicos (con unas vesı́culas que se hinchan y les pro- porcionan flotabilidad) y externos (debido a corrientes y vien- to). Como no siempre son visibles en la superficie, es impor- tante prever la aparición de un bloom en un punto determinado para anticipar posibles riesgos. Por lo tanto, este trabajo impli- ca la elección de un método adecuado para analizar, comparar y simular el flujo de un cuerpo de agua, con el objetivo de obtener resultados con el menor tiempo y parámetros necesarios aplica- bles en la predicción del desplazamiento de las cianobacterias. El esquema que se seguirá es el siguiente: en la sección 2 se presentan las principales opciones disponibles y los candida- tos seleccionados en cuanto al tipo de estructura que soporta la simulación, analizando de forma crı́tica las caracterı́sticas po- sitivas y negativas de los métodos mencionados; la sección 3 presenta un ejemplo con el método seleccionado. Por último, se comentan las conclusiones extraı́das de este estudio. 2. Modelos de simulación de fluidos A continuación se describen tres opciones de simulación: el método de Navier-Stokes promediado por Reynolds (RANS, Reynolds 1895) discretizado a través de volúmenes finitos (FVM), el método Lattice-Boltzmann (LBM, Chen and Doo- len 1998) y los métodos de Smoothed-Particle Hydrodinamics (SPH, Gingold and Monaghan 1977) incompresibles. 2.1. Navier-Stokes promediado por Reynolds y FMV El modelo de Navier-Stokes constituye el enfoque más an- tiguo y común para abordar la simulación de flujos turbulentos. El método RANS somete este modelo a una descomposición de Reynolds: en una componente media (promediada en el tiempo si el flujo de interés es estacionario) y en otra fluctuante. De esta forma, para calcular cantidades de interés tales como el campo de velocidades o presiones (con baja resolución), coeficientes de fuerzas sobre superficies, etc., serı́a suficiente con estudiar la primera de ellas. Sin embargo, esta descomposición añade un término a las ecuaciones relativo a los estreses que sufre el flujo medio debido a la componente fluctuante, conocido como estrés de Reynolds. Este término no es lineal y requiere un modela- do adicional para poder resolver el sistema de manera rápida y computacionalmente eficiente. A raı́z de ello, numerosos mo- delos de turbulencia han sido desarrollados con el objetivo de cubrir distintas casuı́sticas, clasificándose de acuerdo al número de ecuaciones adicionales requeridas para modelar el efecto de la turbulencia. Además, cada uno añade diferentes limitaciones a las soluciones debido a las aproximaciones y suposiciones en las que se basa para simplificar el término, de manera que para ciertas aplicaciones puede ser necesario emplear modelos más complejos que en otras. Si bien RANS solo hace referencia al tratamiento que se le da a las ecuaciones de Navier-Stokes para construir un modelo dinámico del fluido, este método siempre se apoya en mayas de puntos, elementos o volúmenes para llevar a cabo la simu- lación en el espacio discreto. De todas estas opciones la más frecuentemente utilizada en dinámica de fluidos es la discreti- zación en Volúmenes Finitos (FVM) (Kolditz, 2002), que con- siste en transformar las ecuaciones diferenciales (en términos de volúmenes infinitesimales) en relaciones algebraicas discre- tas que actúan sobre volúmenes finitos (o celdas). Para esto, pri- mero, es necesario discretizar el espacio en elementos no super- puestos, sobre los que se integrarán las ecuaciones diferenciales para realizar la transformación. De esta forma, algunos térmi- nos de las ecuaciones, como, por ejemplo, las velocidades, se convierten en flujos a través de las paredes de los elementos de volumen (de modo que la velocidad saliente de un elemento es 394 XLIV Jornadas de Automática 2023. Modelado, Simulación y Optimización. Ferrero-Losada, S. et al., pp. 393-398 la que entra en el contiguo). Mientras que otros, como la densi- dad, se evalúan en el centroide de cada volumen. Este esquema facilita en gran medida la implementación de condiciones de contorno sencillas de una manera no invasiva, ya que se puede alterar directamente el comportamiento de las fronteras entre elementos. Todas estas cualidades hacen que la discretización por volúmenes finitos funcione muy bien a la hora de simular movimientos de fluidos en una malla, ya que la conservación de transferencia de masa, calor o flujo está asegurada por su propia naturaleza. Por otro lado, el uso de este tipo de mallas cuenta con desventajas a la hora de implementar condiciones de contorno más avanzadas. En casos con geometrı́as complicadas, contornos o fronteras móviles, y/o superficies libres, la creación de una malla adecuada se convierte en un proceso difı́cil y cos- toso en cuanto a tiempo y recursos se refiere. Dado que este método es el más comúnmente utilizado exis- te una gran variedad de ejemplos de uso. Además, debido al gran número de modelos de turbulencia desarrollados, sobre los que se sigue trabajando y mejorando continuamente, las publi- caciones tienden a hacer referencia a los modelos empleados. Entre los principales se encuentran el modelo K-ϵ (Jones and Launder, 1972), de dos ecuaciones y el más común; el K-ω (Wilcox, 2008), similar al anterior pero empleando un paráme- tro distinto para la disipación de las turbulencias; y el Menter’s Shear Stress Transport o SST (Menter, 1994), que emplea K-ϵ en las zonas alejadas de los contornos y K-ω cerca de ellos. 2.2. Método de Lattice-Boltzmann El método de Lattice-Boltzmann trata el fluido como una malla de densidades sobre la que se simulan partı́culas ficticias sometidas a procesos de difusión y colisión alternativamente. Estos pasos ejecutados de forma recurrente hacen evolucionar las densidades de cada nodo de la malla. Estas densidades tie- nen distintas componentes dependiendo del tipo de malla que se emplee, pues, a diferencia de lo que sucede en el método anterior, son las conexiones entre nodos las que propagan la evolución, en lugar de las caras de cada celda. De manera que el número de conexiones entre una partı́cula y sus vecinas in- fluye en gran manera en el desarrollo de la simulación, ya que regula directamente la cantidad y el tipo de colisiones posibles en cada nodo, a la vez que puede favorecer un tipo de compor- tamiento difusivo frente a otros. Dada su importancia, existen diversos estudios acerca de los efectos de distintas mallas sobre la simulación (Espinoza-Andaluz et al., 2019). El proceso de difusión en este método consiste en el movi- miento instantáneo de las partı́culas a través de las conexiones de la malla, de forma que a cada paso de tiempo el desplaza- miento de cada partı́cula es fijo, determinado por la longitud de la conexión recorrida. Por otro lado, las colisiones de las partı́culas en los nodos se rigen por el modelo Bhatnagar-Gross- Krook (BGK, Bhatnagar et al. 1954), que define la relajación vı́a colisiones a un estado de equilibrio en densidad para cada una de las conexiones de la malla, teniendo en cuenta la densi- dad en los nodos conectados. Este proceso está regulado por un parámetro temporal, que directamente determina (o está deter- minado por) la viscosidad cinemática del fluido simulado. Con este esquema, se evita tener que lidiar con ecuaciones complicadas como Navier-Stokes, optando por emplear mode- los microscópicos y ecuaciones cinemáticas más sencillas que aglutinen la fı́sica de los procesos que se dan en la micro y me- soescalas, de manera que los comportamientos promedio ma- croscópicos obedezcan las ecuaciones deseadas. Además, este método evita también el tener que seguir a cada partı́cula por separado, como ocurre en las simulaciones de dinámica mole- cular. Por lo tanto, con una mezcla de mallas y partı́culas (ficti- cias) se consigue eliminar desventajas propias de métodos que emplean solo uno de estos elementos. Entre las principales ventajas de este método se encuentra una gran facilidad para implementar geometrı́as complejas (ta- les como medios porosos), debido a la estructura de nodos y conexiones sobre las cuales se pueden imponer condiciones di- rectamente. En cuanto a las fronteras multifase, este método permite realizar simulaciones con gran resolución, en escena- rios donde intervienen burbujas (gas en medio lı́quido) o gotı́cu- las (lı́quido en medio gaseoso), por ejemplo. Por otro lado, en lo que se refiere al rendimiento, desde sus inicios su diseño ha estado enfocado para ser ejecutado de manera lo más rápida posible en arquitecturas paralelas. El método de Lattice-Boltzmann está en constante desarro- llo desde su creación. A lo largo del tiempo se ha ido puliendo y mejorando para afianzarse como una técnica prometedora en el ámbito de la simulación de dinámica de fluidos. En el pasado se han presentado gran variedad de modelos con LBM para di- ferentes tipos de flujos y escenarios (p. ej. Amati et al. 1997 y Flekkøy and Herrmann 1993). Sin embargo, aún no están com- pletamente testados y verificados. En el caso de los modelos multifase, este método destaca a la hora de tratar problemas isotérmicos. Aún con ello, el desarrollo de modelos térmicos más fiables es una de las principales vı́as de estudio, ya que abrirı́a el camino para simular problemas con transferencias de calor además de los fenómenos de superficie.(Chen and Doo- len, 1998) 2.3. Smoothed Particle Hydrodinamcis (SPH) SPH es un método inicialmente desarrollado para proble- mas de astrofı́sica (Gingold and Monaghan, 1977), aunque ac- tualmente su uso abarca muchos más campos de investigación. En el ámbito de la simulación de fluidos ha ido ganando impor- tancia en los últimos años, presentándose como una alternativa útil a los métodos convencionales. A diferencia de estos, SPH no necesita discretizar el espacio con una malla para calcular derivadas espaciales. En su lugar, estas se obtienen por diferen- ciación analı́tica de fórmulas de interpolación, de manera que cualquier función puede ser expresada en términos de su va- lor en puntos desordenados del espacio, representados por las partı́culas. En este caso, éstas no son elementos puntuales, sino que interaccionan entre sı́ a través de lo que se denomina fun- ción núcleo (kernel function), cuyo radio de efecto viene dado por un parámetro tı́picamente representado por h. Estas funcio- nes deben estar normalizadas a la unidad y tender a la delta de Dirac cuando h tiende a 0. Con ellas se puede obtener el valor de cualquier magnitud fı́sica para una partı́cula dada. Basta con su- mar los valores de dicha magnitud en las partı́culas, A j, pesadas con la función núcleo, W(r⃗i j), y su masa, m j (1). Generalmente, suelen usarse como funciones núcleo la función gaussiana o un spline cúbico (o incluso de quinto orden). A(r⃗i) = ∑ j m jA j(r⃗ j)W(r⃗i j) (1) 395 XLIV Jornadas de Automática 2023. Modelado, Simulación y Optimización. Ferrero-Losada, S. et al., pp. 393-398 El hecho de que este método no emplee mallas le confie- re una serie de ventajas respecto al resto, comenzando por la simplificación en la implementación, al no tener que generar la malla. Además, el coste computacional por número de partı́cu- las es significativamente menor que el coste para un sistema con celdas; siempre y cuando el parámetro de interés esté relaciona- do directamente con la densidad del fluido. Esto se debe a que al emplear partı́culas, la resolución espacial está directamente relacionada con la densidad de fluido. Sin embargo, en los ca- sos en los que basta con una malla uniforme, los métodos de celdas tienen una eficiencia mayor, ya que harán falta muchas partı́culas para llenar volúmenes en los que no es necesaria tanta resolución. Por otro lado, en lo relativo a los escenarios preferi- dos para este método, SPH se presenta como una solución ideal para tratar problemas dominados por dinámicas complejas en las fronteras, tales como flujos con superficies libres o fronteras móviles, situaciones en las que la estaticidad de la malla juega en contra del interés del usuario. Sin embargo, cabe destacar que se ha llegado a afirmar que “la implementación de condi- ciones de contorno es ciertamente uno de los puntos técnicos más difı́ciles del método SPH” (Shadloo et al., 2016), debido a que las partı́culas en la vecindad de las fronteras cambian con el tiempo. Este método es relativamente nuevo en el campo de la dinámica de fluidos y aún existen ciertas cuestiones que es ne- cesario pulir (p. ej.: paso de tiempo muy pequeño, inestabili- dades numéricas,...). Sin embargo existen numerosas lı́neas de trabajo y publicaciones recientes que tratan de incrementar la eficiencia, precisión y áreas de aplicación de esta metodologı́a, p. ej. Mintu et al. 2021, Shi et al. 2019. 2.4. Análisis Comparativo En la Tabla 1 se sintetizan las caracterı́sticas señaladas en los métodos anteriores, indicando positivas con un +, y nega- tivas con un -. A continuación se discutirá cada uno de estos métodos teniendo en cuenta estas ventajas y desventajas, desde el punto de vista de la simulación de la dinámica de un cuerpo de agua natural. Como se puede ver en la tabla 1, las tres opciones presen- tadas tienen buenas cualidades para llevar a cabo la simulación deseada. Sin embargo, hay que recordar que ésta se sitúa en el contexto del análisis dinámico de embalses y lagunas, y que es un paso inicial en el desarrollo de un modelo de cuerpo de agua lenticular que permita al resto de modelos involucrados en el proyecto realizar predicciones relacionadas con la aparición de blooms de cianobacterias. Por tanto, de cara a la elección de una metodologı́a de simulación, se dará más importancia a unos as- pectos que a otros de acuerdo a los objetivos futuros para dicho modelo dentro del proyecto. En lo que respecta a la simulación de un lago, una de las problemáticas que pueden surgir es la geometrı́a de la cuenca y la implementación de entradas y salidas de caudal. Dado que en nuestro caso se trabaja con distintos embalses y lagunas, con perfiles batimétricos dispares, será de gran importancia garan- tizar que las condiciones de contorno no resultan un proble- ma a la hora de implementarse, y, en la medida de lo posible, que sean fáciles de implementar y cambiar entre configuracio- nes. Por ello, LBM y SPH se presentan, a priori, como mejores candidatas para la tarea, ya que tienden a ser más versátiles en cuanto a la geometrı́a del escenario que los métodos más clási- cos. Por otro lado, de cara a trabajos futuros, es de interés la capacidad del candidato de realizar simulaciones multifase, ya que la estratificación de temperaturas en este tipo de cuerpos de agua es de gran importancia a la hora de estudiar la posición de las cianobacterias. Analizando las capacidades de LBM y SPH en este aspecto, se llega a la conclusión de que, si bien LBM es un sistema válido para este tipo de escenarios, cuenta con cier- tas carencias en el apartado del modelado termodinámico; en el caso de SPH, el hecho de que la interfase quede determinada de manera orgánica por la naturaleza de la simulación, y que la estratificación esté estrechamente relacionada con la densidad hacen que este sea la opción elegida para ser desarrollada en la siguiente sección. 3. Caso de uso Este ejemplo sencillo muestra, como prueba de concepto, lo que se pretende conseguir finalmente con esta lı́nea de trabajo. Se emulará el efecto del viento sobre las corrientes de un cuerpo de agua con geometrı́a simple utilizando un método SPH. 3.1. Modelo matemático El conjunto de ecuaciones que se emplean a la hora de reali- zar el ejemplo en cuestión está inspirado por el trabajo de Rouz- bahani and Hejranfar 2017, donde se emplea el método desa- rrollado por Chorin (Chorin, 1967) para imponer la condición de incompresibilidad en combinación con la versión lagrangia- na adimensional de la ecuación de Navier-Stokes (ecuación 2). Las distintas magnitudes se han dividido por unas cantidades de referencia para hacerlas adimensionales (3). DV⃗ Dt = −∇⃗p + 1 Re ∇2V⃗ + 1 Fr2 f⃗ (2) t = t∗ Lre f /Vre f , xi = x∗i Lre f , Vi = V∗i Vre f (3) p = p∗ ρV2 re f , f⃗ = f⃗ ∗ gre f El método de Chorin introduce una compresibilidad artifi- cial a través de la ecuación de continuidad, a la que se le añade un término de variación de densidad para luego sustituirla por una presión a través de la relación (4), en la que β usualmente tiene un valor entre 1 y 10. Dando la expresión (5). Nótese que esta expresión usa un tiempo artificial, τ, en lugar del real, t. ∂ρ̃ ∂τ + ∇⃗ · V⃗ = 0, ρ̃ = β2 · p (4) 1 β2 ∂p ∂τ + ∇⃗ · V⃗ = 0 (5) Pasando esta expresión a forma lagrangiana se obtiene (6). Dp Dτ = (V⃗ · ∇)p − β2∇ · V⃗ (6) En el trabajo de Rouzbahani and Hejranfar se propone un esquema de paso de tiempo dual implı́cito, con el cual es posi- ble realizar simulaciones temporalmente precisas. Este esque- ma añade una derivada parcial de la velocidad respecto al tiem- po artificial en el primer miembro de la ecuación (2). De manera 396 XLIV Jornadas de Automática 2023. Modelado, Simulación y Optimización. Ferrero-Losada, S. et al., pp. 393-398 Tabla 1: Comparación de las ventajas (+) y desventajas (-) de RANS-FVM, LBM y SPH. RANS-FVM Lattice-Boltzmann SPH + El más común. + Algoritmo sencillo (2 pasos). + Sin malla (simple y paralelizable). +Versátil (modelos de turbulencias). + Variedad de mallas. + Conservación de m. + Fácil cond. contorno sencillas. + Geometrı́as complejas (estáticas). + Fronteras complejas, dinámicas. + Conservación de masa, v y calor. + Diseño eficiente para arq. paralelas. + Coste menor (métrica relacionada con ρ). + Simulaciones estables. +Multi-fase fácil (isotérmicos) + Interfases por construcción. +Meso- y micro-escalas. + Resolución espacial ∝ dens. de partı́culas. - Coste mayor (métrica ∝ ρ). - Elección de modelo de adecuado. - Mayor influencia de la malla. - Relativamente nuevo en el campo. - Estrés de Reynolds. - Modelos requieren más tests. - Peor en casos con geometrı́a simple. - Geometrı́as complejas. - Sin modelo térmico consistente. - Implementar cond. de contorno. que la expresión de las velocidades queda como: DV⃗ Dτ + DV⃗ Dt = (V⃗ · ∇⃗)V⃗ − ∇⃗p + 1 Re ∇2V⃗ + 1 Fr2 f⃗ (7) Sin embargo, dado que el estudio dinámico de los afloramientos se realiza a lo largo de varios dı́as, nosotros nos centraremos en la parte estacionaria, que permite obtener las corrientes princi- pales que afectan a dicha dinámica. Por ello, basta escoger un paso de tiempo real suficientemente grande para despreciar la derivada respecto a t de la expresión anterior. Una vez establecido el modelo dinámico, vemos cómo im- plementar la metodologı́a SPH. Partiendo de que las magnitu- des que intervienen en las expresiones anteriores pueden expre- sarse por (1), y que a la notación actual se traduce como (8), donde ∆V j es el volumen de cada partı́cula dentro del radio de influencia de la función W, se obtienen los gradientes por (9), y las divergencias de manera similar, cambiando la última multi- plicación por un producto escalar. Ai ≈ ∑ j ∆V j(Ai − A j)W(ri j) (8) ∇Ai ≈ ∑ j ∆V j(Ai ∓ A j)∇W(ri j) (9) En cuanto a la función núcleo, se ha empleado el siguiente spline cúbico para dos dimensiones: W(ri j, h) = 15 7πh2  2 3 − ( ri j h )2 + 1 2 ( ri j h )3 0 ≤ ri j < h 1 6 (2 − ri j h )3 h ≤ ri j < 2h 0 ri j ≤ 2h (10) Finalmente, con todos estos elementos, y dando el mismo valor de ∇V⃗ para todas las partı́culas, se obtienen las siguientes ecuaciones de evolución para la presión y la velocidad: Dpi Dτ = V⃗i · ∇V ∑ j (pi + p j)∇⃗W (⃗ri j) − β2(V⃗ j − V⃗i)∇⃗W (⃗ri j) (11) DV⃗i Dτ = ∇VV⃗i ∑ j (V⃗i + V⃗ j)T ∇⃗W (⃗ri j)− ∇V ∑ j ( pi + p j + R(|pi| + |p j|) W(ri j) W(ri j) ) ∇⃗W (⃗ri j)+ 8∇V Re ∑ j ∇⃗W(ri j) · r⃗i j ∥ri j∥ 2 (V⃗i − V⃗ j) + 1 Fr2 f⃗ (12) Estas ecuaciones tienen ciertas diferencias con el trabajo de referencia. Esto se debe a que se han detectado errores en cier- tos términos y que han sido subsanados a partir de la consulta en (Lee et al. 2008 y Shao and Lo 2003). Además, también se ha añadido una fuerza de repulsión por presiones en la ecuación 12 para evitar algunos problemas debidos a la función núcleo que hemos escogido (ver Monaghan 2000). La intensidad de esta fuerza se regula a través del parámetro R. Al ejecutar la simulación se espera que, a medida que los parámetros evolucionan, ciertos términos de variación de pre- sión y velocidad se hagan nulos, o casi nulos, alcanzado ası́ el estado pseudo-estacionario, al quedar solo los términos con el operador (∇ · V⃗) para regir el movimiento incompresible. En lo que se refiere a las condiciones de contorno, las fronteras sólidas del escenario estarán definidas por capas de partı́culas fijas (para imponer V⃗ = 0) y su presión vendrá de- terminada por las de la capa interior, que está en contacto con el fluido. Las presiones de esta capa se replicarán en el resto de la pared en dirección perpendicular a la superficie interior, de manera que se garantice la condición de presión constante en esa dirección (ver Lee et al. 2008, Figura 1). 3.2. Simulación y Resultados Sin pérdida de generalidad se puede realizar la simulación de fluidos en 2D, lo que permite observar más fácilmente el comportamiento del método escogido y reducir su tiempo de cómputo. El escenario a simular consta de una caja cuadrada de lado unidad que contiene 40x40 partı́culas equiespaciadas. La pared superior emulará un flujo de viento constante (Vpared) sobre la superficie del agua en el sentido positivo del eje X. Inicialmente el campo de presiones será constante en todo el espacio (incluidas las paredes) con un valor de p0. En la Tabla 2 se muestran todos los parámetros empleados. Tabla 2: Parámetros de la simulación. Lre f Vre f gre f δ dτ β 2 1 9.8 1/40 0.001 3 Fr Vpared p0 h Re R 0.226 0.5 0.8 2,5δ 400 3 · 10−5 En la Figura 1 se puede observar el campo de velocidades resultante de correr la simulación durante más de una hora. Con este simple ejemplo se puede observar que se obtiene un flujo 397 XLIV Jornadas de Automática 2023. Modelado, Simulación y Optimización. Ferrero-Losada, S. et al., pp. 393-398 rotatorio que puede hacer sumergirse y emerger a las posibles cianobacterias (sin tener en cuenta otros factores). Además, la aparición del vórtice central indica claramente el correcto fun- cionamiento de la dinámica del fluido. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 x adimensional 0.2 0.4 0.6 0.8 1 y a d im e n s io n a l Figura 1: Campo de velocidades. Azul: vel. adimensionales. Naranja: variación. Sin embargo, los resultados muestran también algunos fa- llos que presenta el modelo. Las flechas naranjas (DV⃗/Dτ) de la figura anterior indican que el sistema aún no ha convergido al estado pseudo-estacionario ya que las variaciones de veloci- dad son demasiado grandes. Lo mismo ocurre con la presión, la cual mantiene niveles razonables tendiendo a la estratificación, pero sin conseguir alcanzar el pseudo-equilibrio esperado. 4. Conclusiones y trabajo futuro El análisis realizado señala que el modelo SPH puede ser adecuado para modelar corrientes turbulentas en lagos y em- balses en apoyo al estudio del movimiento de cianobacterias. Se ha creado un modelo básico y simplificado de una situación real que ha demostrado un buen comportamiento. Sin embargo, existen elementos en el modelo o en su implementación que ne- cesitan revisión a fin de corregir y prevenir eventos que impiden la convergencia del sistema a un estado de fluido incompresible. Además de la revisión de este ejemplo, la lı́nea de investi- gación avanza gradualmente hacia un modelo más complejo y realista, introduciendo la tercera dimensión, entradas y salidas de agua, geometrı́as más complejas, vientos variables, etc. To- do ello con el fin de obtener un modelo sencillo pero completo de la dinámica general de un embalse y poder combinarlo con el ciclo biológico de las cianobacterias. Agradecimientos Este trabajo está siendo realizado gracias al apoyo de los proyectos IA-GES-BLOOM-CM (Y2020/TCS-6420) del pro- grama de Proyectos Sinergicos de la Comunidad Autónoma de Madrid (CAM); SMART-BLOOMS (TED2021-130123B-I00) de los programas MCIN/AEI/10.13039/501100011033 y Euro- pean Union NextGenerationEU/PRTR; y el Programa Investigo de la CAM cofinanciado por el Fondo Social Europeo a través del Programa Operativo de Empleo Juvenil. Referencias Amati, G., Succi, S., Benzi, R., 1997. Turbulent channel flow simulations using a coarse-grained extension of the lattice boltzmann method. Fluid Dynamics Research 19 (5), 289–302. DOI: 10.1016/S0169-5983(96)00026-3 Bhatnagar, P. L., Gross, E. P., Krook, M., May 1954. A model for collision processes in gases. i. small amplitude processes in charged and neutral one- component systems. Phys. Rev. 94, 511–525. DOI: 10.1103/PhysRev.94.511 Chen, S., Doolen, G. D., 1998. Lattice boltzmann method for fluid flows. An- nual Review of Fluid Mechanics 30 (1), 329–364. DOI: 10.1146/annurev.fluid.30.1.329 Chorin, A. J., 1967. A numerical method for solving incompressible viscous flow problems. Journal of Computational Physics 2 (1), 12–26. DOI: 10.1016/0021-9991(67)90037-X Espinoza-Andaluz, M., Moyón, A., Andersson, M., 2019. A comparative study between d2q9 and d2q5 lattice boltzmann scheme for mass transport pheno- mena in porous media. Computers Mathematics with Applications 78 (9), 2886–2896. DOI: 10.1016/j.camwa.2019.02.012 Flekkøy, E., Herrmann, H., 1993. Lattice boltzmann models for complex fluids. Physica A: Statistical Mechanics and its Applications 199 (1), 1–11. DOI: 10.1016/0378-4371(93)90091-H Gingold, R. A., Monaghan, J. J., 12 1977. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society 181 (3), 375–389. DOI: 10.1093/mnras/181.3.375 Jones, W., Launder, B., 1972. The prediction of laminarization with a two- equation model of turbulence. International Journal of Heat and Mass Trans- fer 15 (2), 301–314. DOI: 10.1016/0017-9310(72)90076-2 Kolditz, O., 2002. Finite Volume Method. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 173–190. DOI: 10.1007/978-3-662-04761-38 Lee, E. S., Moulinec, C., Xu, R., Violeau, D., Laurence, D., Stansby, P., sep 2008. Comparisons of weakly compressible and truly incompressible algo- rithms for the sph mesh free particle method 227 (18), 8417–8436. DOI: 10.1016/j.jcp.2008.06.005 Menter, F. R., 1994. Two-equation eddy-viscosity turbulence models for engi- neering applications. AIAA Journal 32 (8), 1598–1605. DOI: 10.2514/3.12149 Mintu, S., Molyneux, D., Colbourne, B., 2021. Full-scale sph simulations of ship-wave impact generated sea spray. Ocean Engineering 241, 110077. DOI: 10.1016/j.oceaneng.2021.110077 Monaghan, J., 2000. Sph without a tensile instability. Journal of Computational Physics 159 (2), 290–311. DOI: 10.1006/jcph.2000.6439 Multiphysics, C., 1998. Introduction to comsol multiphysics®. COMSOL Mul- tiphysics, Burlington, MA, accessed Feb 9, 2018. Reynolds, O., 1895. Iv. on the dynamical theory of incompressible viscous fluids and the determination of the criterion. Philosophical Transactions of the Royal Society of London. (A.) 186, 123–164. DOI: 10.1098/rsta.1895.0004 Roelvink, J. A., Van Banning, G. K. F. M., 1995. Design and development of delft3d and application to coastal morphodynamics. Oceanographic Litera- ture Review 42 (11), 925. Rouzbahani, F., Hejranfar, K., 2017. A truly incompressible smoothed particle hydrodynamics based on artificial compressibility method. Computer Phy- sics Communications 210, 10–28. DOI: 10.1016/j.cpc.2016.09.008 Shadloo, M., Oger, G., Le Touzé, D., 2016. Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges. Computers Fluids 136, 11–34. DOI: 10.1016/j.compfluid.2016.05.029 Shao, S., Lo, E. Y., 2003. Incompressible sph method for simulating newtonian and non-newtonian flows with a free surface. Advances in Water Resources 26 (7), 787–800. DOI: 10.1016/S0309-1708(03)00030-7 Shi, H., Si, P., Dong, P., Yu, X., 2019. A two-phase sph model for massive sedi- ment motion in free surface flows. Advances in water resources 129, 80–98. Weller, H. G., Tabor, G., Jasak, H., Fureby, C., 11 1998. A tensorial approach to computational continuum mechanics using object-oriented techniques. Computer in Physics 12 (6), 620–631. DOI: 10.1063/1.168744 Wilcox, D. C., 2008. Formulation of the k-w turbulence model revisited. AIAA Journal 46 (11), 2823–2838. DOI: 10.2514/1.36541 398 XLIV Jornadas de Automática 2023. Modelado, Simulación y Optimización. Ferrero-Losada, S. et al., pp. 393-398