Universidad de Costa Rica Sistema de Estudios de Posgrado Modelado de flujo bifásico no lineal por hidrodinámica de part́ıculas suavizadas Tesis sometida a la consideración de la Comisión del Programa de Doctorado en Ingenieŕıa para optar por el grado y t́ıtulo de Doctorado Académico en Ingenieŕıa Juan Gabriel Monge Gapper Ciudad Universitaria Rodrigo Facio, Costa Rica 2024 Dedicatoria Para Rebeca y Ana ii Agradecimiento Agradezco toda la asesoŕıa técnica y el genuino apoyo personal que han prestado los miembros del panel asesor, los investigadores de la Universidad Politécnica de Madrid, el personal de apoyo del Centro Nacional de Computación Avanzada y muchos otros a quienes abrumé con consultas, conjeturas y preocupaciones. iii Eṕıgrafe I am an old man now, and when I die and go to heaven there are two matters on which I hope for enlightenment. One is quantum electrodynamics, and the other is the turbulent motion of fluids. And about the former I am rather optimistic. —H. Lamb iv Esta tesis fue aceptada por la Comisión del Programa de Doctorado en Ingenieŕıa de la Universidad de Costa Rica, como requisito parcial para optar por el grado y t́ıtulo de Doctorado Académico en Ingenieŕıa. Dra. Alicia Rojas Araya Representante de la Decana Sistema de Estudios de Posgrado Dr. Alberto Serrano Pacheco Director de tesis Dr. Juan Gabriel Calvo Alṕızar Asesor Dr. Esteban Meneses Rojas Asesor Dr. Esteban Durán Herrera Representante del Director Programa de Doctorado en Ingenieŕıa Juan Gabriel Monge Gapper Sustentante v Índice general Dedicatoria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Agradecimientos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Eṕıgrafe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Hoja de tribunal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Resumen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Lista de cuadros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Lista de figuras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix Simboloǵıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx Lista de abreviaturas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii Licencia de publicación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiv Lista de publicaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxvi 1. Introducción 1 1.1. Motivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2. Estado del arte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3. Objetivos e hipótesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4. Estructura de la tesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2. Marco teórico 13 2.1. Fundamentación matemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.1. Estado de esfuerzos de un fluido . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.2. Continuidad y conservación de la masa . . . . . . . . . . . . . . . . . . . . 16 2.1.3. Conservación de la cantidad de movimiento . . . . . . . . . . . . . . . . . 18 vi 2.1.4. Ecuación de vorticidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2. Método de SPH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2.1. Aproximación por kernel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2.2. Aproximación de derivadas . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2.3. Aproximación por part́ıculas . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.4. Ecuaciones de estado f́ısico . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.2.5. Ecuaciones de continuidad y enerǵıa . . . . . . . . . . . . . . . . . . . . . 36 2.2.6. Funciones de suavizado (kernel) . . . . . . . . . . . . . . . . . . . . . . . 38 2.2.7. Longitud de suavizado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.2.8. Condiciones de contorno sólidas . . . . . . . . . . . . . . . . . . . . . . . . 43 2.2.9. Mecanismos de estabilización . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.2.10. Listas de part́ıculas vecinas . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.3. Modelos reológicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.4. Consideraciones finales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3. Extensión del método numérico 61 3.1. Esquema numérico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.1.1. Enfoque de viscosidad equivalente . . . . . . . . . . . . . . . . . . . . . . 62 3.1.2. Criterio de fluencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.1.3. Tratamiento de interfase . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.1.4. Estabilización por reposicionamiento de part́ıculas . . . . . . . . . . . . . 73 3.1.5. Detección de part́ıculas de interfase . . . . . . . . . . . . . . . . . . . . . 78 3.1.6. Paso temporal variable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.2. Código fuente SPH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.2.1. Propiedades del código GapSPH . . . . . . . . . . . . . . . . . . . . . . . 82 3.2.2. Paquete AQUAgpuSPH . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.2.3. Aceleración por paralezizado . . . . . . . . . . . . . . . . . . . . . . . . . 94 3.3. Estrategia de validación metodológica . . . . . . . . . . . . . . . . . . . . . . . . 96 4. Validación del esquema numérico 97 4.1. Rotura de presa bidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 vii 4.2. Rotura de presa sobre ŕıo bidimensional . . . . . . . . . . . . . . . . . . . . . . . 107 4.3. Rotura de presa con compuerta bidimensional . . . . . . . . . . . . . . . . . . . . 116 4.4. Colapso viscoplástico bidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.5. Rotura de presa sobre cama de fluido no lineal . . . . . . . . . . . . . . . . . . . 130 4.6. Colapso viscoplástico sumergido . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 4.7. Rotura de presa bifásica sobre lecho mojado . . . . . . . . . . . . . . . . . . . . . 143 4.8. Modelos de prueba tridimensionales . . . . . . . . . . . . . . . . . . . . . . . . . 158 4.8.1. Descarga angular de canal en pileta . . . . . . . . . . . . . . . . . . . . . 159 4.8.2. Choque ortogonal de roturas de presa . . . . . . . . . . . . . . . . . . . . 169 4.9. Convergencia y parámetros numéricos . . . . . . . . . . . . . . . . . . . . . . . . 175 4.9.1. Convergencia de modelos multifásicos . . . . . . . . . . . . . . . . . . . . 175 4.9.2. Ajuste de factores numéricos . . . . . . . . . . . . . . . . . . . . . . . . . 181 4.9.3. Estabilidad de campos de datos . . . . . . . . . . . . . . . . . . . . . . . . 188 4.10. Efecto de la paralelización en el código GapSPH . . . . . . . . . . . . . . . . . . 196 4.11. Recapitulación de resultados esenciales . . . . . . . . . . . . . . . . . . . . . . . . 204 5. Conclusiones 206 5.1. Aportes al modelado de flujos multifásicos . . . . . . . . . . . . . . . . . . . . . . 207 5.2. Aspectos sobre el proceso de implementación . . . . . . . . . . . . . . . . . . . . 214 5.3. Definición de ĺıneas de investigación futura . . . . . . . . . . . . . . . . . . . . . 218 Bibliograf́ıa 221 Anexo 238 viii Resumen Se propuso la extensión para un método numérico basado en hidrodinámica de part́ıculas suavizadas (SPH, por sus siglas en inglés) capaz de modelar la interacción de dos fluidos no miscibles con caracteŕısticas f́ısicas diferentes en un contorno que incluye fronteras sólidas y una frontera libre. Se orienta la metodoloǵıa a la consideración de que uno de los fluidos podrá presentar comportamiento viscoplástico, como el de algunos lodos y el apilamiento de materiales granulados en un medio ĺıquido, como seŕıa arena saturada en agua. Para una primera etapa de desarrollo, se implementó en un programa de cómputo original en lenguaje C de un esquema básico SPH y se valoraron cualitativamente aspectos de estabilidad numérica y se compararon los resultados con experimentos conocidos de hidráulica. Posteriormente, se elaboró una extensión metodológica que incorpora una serie de elementos indispensables para la precisión y estabilidad numérica del modelo y se implementó para integrarse a un paquete de solución SPH que es de estructura modular. El esquema de modelado es una versión modificada del método de SPH con compresibilidad artificial (WCSPH, por sus siglas en inglés), para el que se ideó un enfoque de viscosidad equiva- lente para la fase no lineal, pero con la distinción de que la condición de fluencia, al intervenir en la forma de la ecuación de conservación de cantidad de movimiento, lo haga independiente del valor de viscosidad equivalente en la localidad. Esto resulta imprescindible para garantizar estabilidad numérica y obtener precisión en casos de flujo viscoplástico con valores relativamente altos de es- fuerzo de fluencia, tanto de manera cualitativa en cuanto a comportamiento macroscópico de las interacciones entre fases como la cuantitativa al cotejar con datos experimentales. En esencia, el éxito de la integración de los elementos de este modelo numérico para flujos no lineales complejos se debe a que se segregó por fase f́ısica el cálculo de la densidad local, se aplica un algoritmo de redistribución de part́ıculas espećıfico que evita inestabilidades al campo de velocidades de deformación que de otra manera causan fluencia prematura del fluido, y se utiliza un novedoso factor de atenuación al cálculo de la aceleración del fluido según un criterio de fluencia plástica. Los casos que se utilizaron para validar el método incluyen variedad de configuraciones para ix ilustrar su versatilidad ante estos distintos tipos de fenómenos que pueden aparecer en un flujo no lineal. Se muestran resultados de dinámica de rotura de presa de flujo lineal en una sola fase, flujo no lineal en una sola fase, y flujo multifásico con una de las fases de comportamiento no lineal. Todos los casos corresponden a una configuración de tanque abierto a la atmósfera. Según los resultados de estas simulaciones, se analizan aspectos de sensibilidad a parámetros numéricos y a ciertos parámetros f́ısicos, aśı como las implicaciones que tiene en cuanto al ámbito de aplicabilidad del método WCSPH. Los resultados de calibración indican que este esquema es funcional y razonablemente pre- ciso, con la gran ventaja de que no requiere la implementación de un modelo elástico para la masa que no ha sobrepasado el esfuerzo de fluencia. Es evidente que se reproduce con bastante precisión respecto a ciertos experimentos f́ısicos tanto eventos de erosión como de colapso sin que los elementos metodológicos destinados a replicar este comportamiento se anulen entre śı. Adicionalmente, a pesar de las distorsiones que introducen los elementos no lineales de las ecua- ciones utilizadas, con la estratagema propuesta se logra excelente estabilidad en los campos de presiones y de razones de deformación del flujo. Esto no se hab́ıa logrado en estudios previos orientados a flujos de una fase granular interactuando con otra fase ĺıquida y alguna frontera libre sin recurrir a hibridización con modelos elásticos, que no son necesariamente aplicables a la naturaleza de su dinámica, que en realidad está dirigida por fuerzas de fricción intergranular. La ĺınea de trabajo de esta investigación, centrada en un enfoque de viscosidad aparente con una ecuación de conservación modificada para modelar un flujo granular saturado y el conjunto de módulos de cómputo elaborados provee una excelente herramienta para llevar a cabo estudios exhaustivos de caracterización de flujos multifásicos de este tipo. x Abstract An extension was proposed for a numerical method based on smoothed particle hydrody- namics (SPH) that can model the interaction of two immiscible fluids with di↵erent physical characteristics at a boundary that includes solid boundaries and a free boundary. The metho- dology is oriented to the consideration that one of the fluids may present viscoplastic behavior, similar to that of some sludges and the stacking of granulated materials in a liquid medium, such as water-saturated sand. For a first stage of development, a basic SPH scheme was implemented in an original computer program in C language, aspects of numerical stability were qualitati- vely assessed, and the results were compared with known hydraulic experiments. Subsequently, a methodological extension was developed that incorporates a series of essential elements for the precision and numerical stability of the model and it was implemented to be integrated into an SPH solution package that is modular in structure. The modeling scheme is a modified version of the Weakly Compressible SPH (WCSPH) method for which an equivalent viscosity approach was devised for the non-linear phase, but with the distinction that the yield criteria, by intervening in the form of the momentum conservation equation, makes it independent of the locally equivalent value of viscosity. This is essential to guarantee numerical stability and obtain good precision in cases of viscoplastic flow with relatively high values of yield stress, both qualitatively in terms of the macroscopic behavior of the interactions between phases and quantitatively when comparing with experimental data. In essence, the success of the integration of the elements of this numerical model for complex nonlinear flows is due to the fact that the calculation of the local density was segregated by physical phase, a specific particle redistribution algorithm is applied that avoids instabilities in the field. of deformation rates that would otherwise cause premature fluid yielding, and a novel attenuation factor is used to calculate fluid acceleration based on a yield criterion. The cases used to validate the method include a variety of configurations to illustrate its versatility in the face of these di↵erent types of phenomena that can appear in a nonlinear flow. xi Results of dam break dynamics of linear flow in a single phase, nonlinear flow in a single phase, and multiphase flow with one of the phases of nonlinear behavior are shown. All cases correspond to a tank configuration open to the atmosphere. Based on the results of these simulations, aspects of sensitivity to numerical parameters and to certain physical parameters are analyzed, as are the implications that these results have in terms of the field of applicability of the WCSPH method. The calibration results indicate that this scheme is functional and reasonably accurate, with the great advantage that it does not require the implementation of an elastic model for the mass that has not exceeded the yield stress. It is evident that both erosion and collapse events are reproduced with good precision with respect to certain physical experiments without the met- hodological elements intended to replicate this behavior canceling each other out. Additionally, despite the distortions introduced by the non-linear elements of the equations used, excellent stability in the fields of pressures and flow deformation ratios is achieved with the proposed stratagem. This had not been achieved in previous studies oriented to flows of a granular phase interacting with another liquid phase and some free boundary without resorting to hybridization with elastic models, which are not necessarily applicable to the nature of its dynamics, that are in reality driven by intergranular friction forces. The line of work of this research, focused on an apparent viscosity approach with a modified conservation equation to model a saturated granular flow, and the set of computational modules developed provide an excellent tool to carry out exhaustive flow characterization studies for these families of multiphase flow. Ćıtese este trabajo como: Monge-Gapper, Juan Gabriel. (2024) Modelado de flujo bifásico no lineal por hidrodinámica de part́ıculas suavizada Tesis de grado para optar por el t́ıtulo de Doctorado Académico en Ingenieŕıa. Universidad de Costa Rica: San Pedro de Montes de Oca, San José, COSTA RICA. xii Índice de cuadros 4.1. Condiciones generales de casos para validación . . . . . . . . . . . . . . . . . . . 100 4.2. Parámetros del modelo numérico para rotura de presa bidimensional . . . . . . . 102 4.3. Parámetros del modelo numérico para rotura de presa sobre ŕıo . . . . . . . . . . 108 4.4. Parámetros del modelo numérico para presa con compuerta . . . . . . . . . . . . 118 4.5. Parámetros del modelo numérico para rotura de presa viscoplástica (A) . . . . . 125 4.6. Parámetros del modelo numérico para rotura de presa viscoplástica (B) . . . . . 125 4.7. Parámetros del modelo numérico para rotura de presa sobre cama de arena . . . 131 4.8. Parámetros del modelo numérico para rotura de presa viscoplástica sumergida . . 137 4.9. Parámetros del modelo numérico para rotura de presa bifásica (a) . . . . . . . . 145 4.10. Parámetros del modelo numérico para rotura de presa bifásica (b) . . . . . . . . 145 4.11. Parámetros del modelo numérico para descarga angular de canal en pileta . . . . 161 4.12. Parámetros del modelo numérico para choque ortogonal de roturas de presa . . . 170 4.13. Directivas de paralelización en funciones de GapSPH . . . . . . . . . . . . . . . 200 4.14. Tiempos de procesamiento para rotura de presa sobre ŕıo poco profundo . . . . . 202 4.15. Distribución de tiempos de procesamiento por función en GapSPH . . . . . . . . 202 xiii Índice de figuras 1.1. Modelado computacional de interacción entre fases inmiscibles . . . . . . . . . . 1 2.1. Estado de esfuerzos de un fluido . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2. Volumen de control para conservación de masa . . . . . . . . . . . . . . . . . . . 16 2.3. Volumen de control para conservación de cantidad de movimiento . . . . . . . . . 18 2.4. Discretización por malla o euleriana . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.5. Discretización por part́ıculas o lagrangiana . . . . . . . . . . . . . . . . . . . . . 27 2.6. Modelo de part́ıculas de fluido no viscoso en el paquete Blender . . . . . . . . . . 28 2.7. Interpretacion conceptual de la función de suavizado W . . . . . . . . . . . . . . 30 2.8. Concepto de part́ıculas virtuales de frontera . . . . . . . . . . . . . . . . . . . . 44 2.9. Concepto del algoritmo de búsqueda directo . . . . . . . . . . . . . . . . . . . . . 52 2.10. Concepto del algoritmo de búsqueda por listas de correlaciones . . . . . . . . . . 53 2.11. Concepto del algoritmo de búsqueda por zonaje . . . . . . . . . . . . . . . . . . . 53 2.12. Concepto del algoritmo de búsqueda árbol . . . . . . . . . . . . . . . . . . . . . . 54 2.13. Morfoloǵıa general de curvas de variación de viscosidad (adaptado de [2]) . . . . 56 2.14. Variación de viscosidad de una suspensión de bentonita . . . . . . . . . . . . . . 57 3.1. Composición de ecuaciones de viscosidad efectiva . . . . . . . . . . . . . . . . . . 63 3.2. Variación del esfuerzo de fluencia local por apilamiento . . . . . . . . . . . . . . . 66 3.3. Ausencia de part́ıculas cerca de fronteras sólidas o libres . . . . . . . . . . . . . . 69 3.4. Redistribución de part́ıculas en flujo 2D en una cavidad . . . . . . . . . . . . . . 74 3.5. Variables geométricas asociadas a fronteras libres . . . . . . . . . . . . . . . . . . 78 3.6. Diagrama de flujo general del código fuente Gap-SPH v.1.0 . . . . . . . . . . . . 83 xiv 3.7. Diagrama de flujo general del código fuente AQUAgpuSPH v. 3.0 . . . . . . . . . 90 3.8. Diagrama de flujo del módulo Rates.cl para flujo no lineal . . . . . . . . . . . . . 92 4.1. Evolución del flujo en experimento de rotura de presa [77] . . . . . . . . . . . . . 99 4.2. Condiciones de frontera para rotura de presa bidimensional . . . . . . . . . . . . 101 4.3. Campo de velocidades para rotura de presa 2D . . . . . . . . . . . . . . . . . . . 103 4.4. Campo de razón de deformación para rotura de presa 2D . . . . . . . . . . . . . 103 4.5. Campo de presiones para rotura de presa 2D (t = 3.0 s) . . . . . . . . . . . . . . 104 4.6. Comparación con perfiles experimentales para rotura de presa (td < 10) . . . . . 105 4.7. Hidrograma para posición X1 de rotura de presa y soluciones numéricas . . . . . 106 4.8. Hidrograma para posición X2 de rotura de presa y soluciones numéricas . . . . . 106 4.9. Fronteras para rotura de presa sobre ŕıo . . . . . . . . . . . . . . . . . . . . . . . 108 4.10. Evolución del campo de velocidades para rotura de presa y ŕıo con ↵ = 0.1 . . . 109 4.11. Evolución del campo de velocidades para rotura de presa y ŕıo con ↵ = 0.4 . . . 110 4.12. Comparación con fotograf́ıas del experimento f́ısico según [61] . . . . . . . . . . . 111 4.13. Hidrograma para posición P1 de rotura de presa sobre ŕıo (↵ = 0.1) . . . . . . . 112 4.14. Hidrograma para posición P1 de rotura de presa sobre ŕıo (↵ = 0.4) . . . . . . . 112 4.15. Hidrograma para posición P2 de rotura de presa sobre ŕıo (↵ = 0.1) . . . . . . . 112 4.16. Hidrograma para posición P2 de rotura de presa sobre ŕıo (↵ = 0.4) . . . . . . . 113 4.17. Hidrograma para posición P3 de rotura de presa sobre ŕıo (↵ = 0.1) . . . . . . . 113 4.18. Hidrograma para posición P3 de rotura de presa sobre ŕıo (↵ = 0.4) . . . . . . . 113 4.19. Geometŕıa a escala distorsionada de las crestas que viajan de derecha a izquierda para la rotura de presa sobre ŕıo para td = 56,4 (↵ = 0.4) . . . . . . . . . . . . . 114 4.20. Evolución del campo de presiones para rotura de presa sobre ŕıo (↵ = 0.4) . . . . 115 4.21. Evolución del campo de presiones para rotura de presa sobre ŕıo (↵ = 0.1) . . . . 115 4.22. Detalle del campo de razón de deformación angular para rotura de presa sobre ŕıo (↵ = 0,4) a td = 56,4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.23. Fronteras para presa con compuerta . . . . . . . . . . . . . . . . . . . . . . . . . 117 4.24. Detalle a escala distorsionada de transición en rotura de presa larga con 3500 part́ıculas (td = 16) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 xv 4.25. Detalle a escala distorsionada de transición en rotura de presa larga con 14 000 part́ıculas (td = 16) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 4.26. Solución de rotura de presa para td = 15.6 . . . . . . . . . . . . . . . . . . . . . . 119 4.27. Evolución del campo de velocidades para apertura de compuerta . . . . . . . . . 120 4.28. Evolución del campo de presiones para apertura de compuerta . . . . . . . . . . 120 4.29. Solución de presa con compuerta de 20% de altura para td = 15.6 . . . . . . . . 121 4.30. Detalle del campo de razón de deformación cortante para apertura de compuerta a td = 15.0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.31. Condiciones de frontera para rotura de presa viscoplástica . . . . . . . . . . . . . 123 4.32. Perfiles para rotura de presa viscoplástica para t/ p g/Hf  8,7 (t  1.0 s) . . . . 126 4.33. Viscosidades efectivas para rotura de presa viscoplástica . . . . . . . . . . . . . . 127 4.34. Campo de presión para rotura de presa viscoplástica . . . . . . . . . . . . . . . . 127 4.35. Razón de deformación angular para rotura de presa viscoplástica a t = 1.000 s . 128 4.36. Modelo tridimensional de rotura de presa viscoplástica a td = 8.7 (t = 1.0 s) . . 129 4.37. Geometŕıa para el caso de rotura de presa sobre cama de arena . . . . . . . . . . 131 4.38. Perfiles de superficie para arena de ruptura de presa sobre cama no lineal en td = 2.9 (t =0.50 s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 4.39. Perfiles para arena de ruptura de presa sobre cama no lineal en td = 5.7(t = 1.00 s)133 4.40. Perfiles para agua de ruptura de presa sobre cama no lineal en td = 2.9 (t = 0.50 s)133 4.41. Perfiles para agua de ruptura de presa sobre cama no lineal en td = 5.7 (t = 1.00 s)134 4.42. Campo de velocidades para ruptura de presa sobre fondo erosionable . . . . . . . 135 4.43. Viscosidad efectiva para ruptura de presa sobre fondo erosionable . . . . . . . . . 135 4.44. Campo de presiones para ruptura de presa sobre fondo erosionable . . . . . . . . 136 4.45. Campo de viscosidades para ruptura de presa sobre fondo erosionable (t = 6.000 s)136 4.46. Geometŕıa para el caso de ruptura de presa sumergida . . . . . . . . . . . . . . . 138 4.47. Campo de presiones para colapso viscoplástico sumergido (t = 1.5 s) . . . . . . . 138 4.48. Proceso de rotura de presa viscoplástica sumergida en instantes significativos . . 139 4.49. Perfiles de superficie de presa viscoplástica sumeergida para 0.5 s < t < 2.5 s . . 140 4.50. Campo de razones de deformación angular para colapso viscoplástico sumergido (t = 0.500 s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 xvi 4.51. Campo de razones de deformación angular para colapso viscoplástico sumergido (t = 1.500 s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 4.52. Cohesión de origen numérico en un caso de prueba preliminar de colapso viscoplástico143 4.53. Geometŕıa para ruptura de presa bifásica sobre lecho mojado . . . . . . . . . . . 144 4.54. Evolución de perfiles de flujo para rotura de presa bifásica con 80% de arena . . 147 4.55. Evolución de perfiles de flujo para rotura de presa bifásica con 10% de arena . . 147 4.56. Profundidad del agua para presa con 66% de arena en x/L = 0.5 . . . . . . . . . 148 4.57. Profundidad del agua para presa con 66% de arena en x/L = 1.0 . . . . . . . . . 149 4.58. Profundidad del agua para presa con 66% de arena en x/L = 1.7 . . . . . . . . . 149 4.59. Profundidad del agua para presa con 80% de arena en x/L = 0.5 . . . . . . . . . 150 4.60. Profundidad del agua para presa con 80% de arena en x/L = 1.0 . . . . . . . . . 150 4.61. Profundidad del agua para presa con 80% de arena en x/L = 1.7 . . . . . . . . . 151 4.62. Campo de presiones a t = 1.0 s para presa con 80% de arena . . . . . . . . . . . 152 4.63. Campo de presiones a t = 1.0 s para presa con 20% de arena (sin redistribución de part́ıculas) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 4.64. Evolución de niveles de fase no lineal para presa con 80% de arena . . . . . . . . 153 4.65. Evolución de niveles de fase no lineal para presa con 10% de arena . . . . . . . . 153 4.66. Profundidad de fase no lineal para presa con 50% de arena . . . . . . . . . . . . 154 4.67. Profundidad de fase lineal para presa con 50% de arena . . . . . . . . . . . . . . 155 4.68. Evolución del campo de viscosidades para rotura de presa bifásica con 50% de arena156 4.69. Hidrogramas para rotura de presa bifásica con 50% de arena . . . . . . . . . . . 156 4.70. Campo de razón de deformación angular para presa con 50% de arena a t = 1.0 s 157 4.71. Campo de presiones para rotura de presa bifásica con 50% de arena . . . . . . . 157 4.72. Campo de velocidades para rotura de presa bifásica con 50% de arena a t = 1.0 s 157 4.73. Geometŕıa para descarga angular de canal en pileta . . . . . . . . . . . . . . . . . 160 4.74. Velocidad absoluta para descarga angular de canal en pileta (t = 0.5 s y 1.0 s) . 163 4.75. Velocidad absoluta para descarga angular de canal en pileta (t = 1.5 s y 2.5 s) . 164 4.76. Campo de viscosidades para lodo en descarga angular de canal en pileta (t = 2.5 s)165 4.77. Razón de deformación angular para descarga angular de canal en pileta (t = 1.0 s) 165 4.78. Campo de presiones para descarga angular de canal en pileta (t = 1.0 s) . . . . . 166 xvii 4.79. Proceso de mezcla del agua (fase lineal, azul) con el lodo (fase no lineal, rojo) para t = 0.5 s, 1.0 s y 1.5 s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 4.80. Geometŕıa para choque ortogonal de presas . . . . . . . . . . . . . . . . . . . . . 170 4.81. Magnitud de velocidad para choque ortogonal de presas (t = 0.27 s, 0.77 s y 1.33 s)172 4.82. Part́ıculas de fase no lineal separadas de la masa principal a t = 2.0 s . . . . . . 173 4.83. Campo de presiones para choque ortogonal de presas (t = 2.2 s) . . . . . . . . . 173 4.84. Viscosidad equivalente para fase no lineal en choque ortogonal de presas (t = 2.2 s)174 4.85. Razón de deformación angular para choque ortogonal de presas (t = 2.2 s) . . . . 174 4.86. Convergencia de perfiles para rotura de presa viscoplástica a td = 4.3 (t = 0.5 s) 177 4.87. Convergencia de perfiles para rotura de presa viscoplástica a td = 8.7 (t = 1.0 s) 178 4.88. Convergencia de perfiles de profundidad de agua para rotura de presa con 50% de arena (t = 1.0 s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 4.89. Convergencia de perfiles de profundidad de agua para rotura de presa con 50% de arena (t = 2.0 s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 4.90. Convergencia de perfiles de profundidad de arena para rotura de presa con 50% de arena (t = 1.0 s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 4.91. Convergencia de perfiles de profundidad de arena para rotura de presa con 50% de arena (t = 2.0 s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 4.92. Profundidad del agua para presa con 66% de arena en x/L = 0.5 . . . . . . . . . 184 4.93. Profundidad del agua para presa con 66% de arena en x/L = 1.0 . . . . . . . . . 185 4.94. Profundidad del agua para presa con 66% de arena en x/L = 1.7 . . . . . . . . . 185 4.95. Profundidad del agua para presa con 80% de arena en x/L = 0.5 . . . . . . . . . 187 4.96. Profundidad del agua para presa con 80% de arena en x/L = 1.0 . . . . . . . . . 187 4.97. Profundidad del agua para presa con 80% de arena en x/L = 1.7 . . . . . . . . . 188 4.98. Campo de presiones con oscilaciones numéricas en flujo bajo compuerta (t = 6.0 s)190 4.99. Ruido en campo de viscosidad efectiva por inestabilidad de presiones en rotura de presa bifásica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 4.100.Estructuras periódicas en frontera de ruptura de presa no lineal a td = 10 . . . . 191 4.101.Separación de capa superficial de part́ıculas 2D por frontera libre sin renormalización191 4.102.Formación de grietas de origen numérico en la fase no lineal . . . . . . . . . . . . 192 xviii 4.103.Razones de deformación de periodicidad espacial en una ruptura de presa bifásica 193 4.104.Campo de viscosidades con periodicidad espacial en colapso viscoplástico . . . . 193 4.105.Aglomeraciones sin apilamiento en el frente de avance viscoplástico . . . . . . . . 194 4.106.Agrietado artificial y formación de vaćıo circular en fase viscoplástica . . . . . . . 194 4.107.Campo de viscosidades a t = 1.0 s para presa bifásica con 80% de arena . . . . . 194 4.108.Fase no lineal con cohesividad numérica por algoritmo de redistribución . . . . . 195 4.109.Campo de presiones a t = 1.0 s para presa bifásica con 80% de arena . . . . . . 196 4.110.Ciclo de actualización de la densidad con directivas de paralelización . . . . . . . 199 xix Simboloǵıa Las variables que reciben mayor mención en este documento usan la siguiente nomenclatura convencional: hf(x)i Aproximación por part́ıculas SPH de una función f(x) ai Vector aceleración de la part́ıcula i del fluido discretizado B Matriz de renormalización para corrección por redistribución de part́ıculas C Concentración local de part́ıculas cs Velocidad del sonido artificial en el medio � ij Tensor unitario o función delta de Kronecker �x Corrector por redistribución de part́ıculas para la posición �u Corrector por redistribución de part́ıculas para la velocidad D Dt Derivada total respecto de t D Constante difusiva para técnica de redistribución de part́ıculas Di Constante difusiva según técnica � -SPH para corrección de ecuación de conservación �t Espaciamiento temporal del modelo discretizado "ij Deformación de un elemento de fluido en dirección i proyectada sobre la dirección j. "̇ij Razón de deformación angular cortante en el plano ij de un elemento de fluido. fe Vector de fuerza resultante de origen externo que actúa en un elemento de fluido fi Vector de fuerza resultante de origen interno que actúa en un elemento de fluido fij Función de presión artificial para corrección por redistribución de part́ıculas F Función correctora por número de dimensiones para proceso de renormalización �(x) Factor de renormalización de Shepard g Constante de la gravedad h Longitud de suavizado i Número identificador de la part́ıcula de fluido en estudio I Matriz identidad j Número identificador de la part́ıcula vecina a la part́ıcula i xx  Índice de consistencia de un fluido no lineal K Factor de viscosidad equivalente para flujo cuasiestático KR Factor de relajación para ecuación de conservación en flujo cuasiestático � Segunda viscosidad de un fluido o coeficiente de pérdida viscosa por expansión o compresión. L(x) Corrector por frontera para función de suavizado m Razón de espesamiento de un fluido no lineal µeff Viscosidad efectiva local de una part́ıcula i de fluido N Número total de part́ıculas que conforman el fluido discretizado n Índice de flujo de un fluido no lineal rC Gradiente de concentración local de part́ıculas r̂C Aproximación del gradiente de concentración local de part́ıculas rf Gradiente de una función f r⇥ x Divergencia de un vector x r⌦ x Rotacional de un vector x ru↵� i Razón de deformación total de la part́ıcula i en el plano ↵� n Vector normal a un plano definido por un conjunto de part́ıculas vecinas P Presión P0 Presión inicial de referencia Pa Presión atmosférica local r0 Radio de activación de la fuerza de frontera rij Distancia absoluta entre las part́ıculas i y j. ⇢i Densidad local de la part́ıcula i ⇢0 Densidad inicial de referencia de un fluido incompresible s Vector tangencial a un plano definido por un conjunto de part́ıculas vecinas ⌧0 Esfuerzo de fluencia de un fluido no lineal ⌧ij Componente esfuerzo cortante viscoso actuando en un elemento de fluido en dirección que lo deforma en el plano ij ⌧ii Componente esfuerzo axial viscoso actuando en un elemento de fluido en dirección i ui Vector velocidad de la part́ıcula i del fluido discretizado Vi Volumen de la part́ıcula i xxi W Función de suavizado o función kernel (notación compacta) Wii Valor de la función kernel evaluada para la part́ıcula i Wij Valor de la función kernel evaluada para la part́ıcula j respecto de la part́ıcula i W (rij , h) Función de suavizado o función kernel xi Vector posición de la part́ıcula i del fluido discretizado ⌦ Dominio en el que actúan las part́ıculas del fluido discretizado !i Vorticidad vectorial de la part́ıcula i del fluido discretizado xxii Lista de abreviaturas En este documento se utilizan con frecuencia las siguientes abreviaturas de uso convencional en publicaciones en inglés y en español relacionadas con el tema central de esta investigación: CFL Courant-Friedichs-Levy. DNS Simulación Numérica Directa (Direct Numerical Simulation). GPU Unidad de Procesamiento Gráfico (Graphical Processing Unit). HB Herschel-Bulkley. HBP Herschel-Bulkley-Papanastasiou. MDF Método de Diferencias Finitas. MEC Método de Elementos de Contorno. MED Método de Elementos Discretos. MEF Método de Elementos Finitos. MPM Método de Part́ıcula en Movimiento. MPV Método de Part́ıculas Virtuales. MSPM Método Semiimpĺıcito de Part́ıcula en Movimiento. MVF Método de Volúmenes Finitos. PST Técnica de Desplazamiento de Part́ıculas (Particle Shifting Technique). SPH Hidrodinámica de Part́ıculas Suavizadas (Smoothed Particle Hydrodynamics). VOF Volume of Fluid (Volume of Fluid). WCSPH Hidrodinámica de Part́ıculas Suavizadas con Compresibilidad Artificial (Weakly Com- pressible Smoothed Particle Hydrodynamics). xxiii xxv Lista de publicaciones 1. Monge-Gapper, Juan Gabriel, Serrano-Pacheco, Alberto y Calderón-Sánchez, Javier. (2024). Multiphase simulations of non-linear fluids with SPH. Computational Particle Mechanics. DOI: https://doi.org/10.1007/s40571-024-00712-3 2. Monge-Gapper, Juan Gabriel y Serrano-Pacheco, Alberto. (2023) An implementation of the SPH method and its application to two-dimensional dam break cases. Ingenieŕıa, Investigación y Tecnoloǵıa. Volumen 24, número 1. Ciudad de México. DOI: https://doi.org/10.22201/fi. 25940732e.2023.24.1.005 3. Monge-Gapper, Juan Gabriel y Serrano-Pacheco, Alberto (2021). Flujo en compuerta infe- rior de represa bidimensional larga por el método de WSPH. Revista Ingenieŕıa. Volumen 31, número 2. Costa Rica. DOI: https://doi.org/10.15517/RI.V31I2.45850 4. Monge-Gapper, Juan Gabriel; Serrano-Pacheco, Alberto; Duque, Daniel; Calderón-Sánchez, Javier. (2022). Numerical scheme for multifluid open-flow with discontinuous nonlinear viscosity. Bilotta, G. (ed.) Proceedings of the 16th SPHERIC International Workshop. Catania, Italia, 6–9 Junio de 2022. https://spheric2022.it 5. Monge-Gapper, Juan Gabriel, Serrano-Pacheco, Alberto. (2021). Modelado de rotura de presa viscoplástica por método WCSPH. Arrieta-Orozco, O., Schmidt-Dı́az, V. y Castro-Arce, K. (ed.) Memoria de las III Jornadas de Investigación de la Facultad de Ingenieŕıa. San José, Costa Rica, 27–29 octubre de 2021. DOI: https://doi.org/10.15517/ri.v32iNE2.50666. xxvi https://doi.org/10.1007/s40571-024-00712-3 https://doi.org/10.22201/fi.25940732e.2023.24.1.005 https://doi.org/10.22201/fi.25940732e.2023.24.1.005 https://doi.org/10.15517/RI.V31I2.45850 https://spheric2022.it https://doi.org/10.15517/ri.v32iNE2.50666 Caṕıtulo 1 Introducción El modelado matemático del flujo de fluidos es un campo de estudio en que las ecuaciones generales no pueden aplicarse directamente a un caso de interés particular. Se requieren sustan- ciales simplificaciones o bien elaborados métodos numéricos para resolver sistemas de ecuaciones para un contorno f́ısico de particular interés. Los mejores modelos, correlaciones experimentales o modelos determińısticos, se aplican a un fenómeno de dinámica de fluidos muy espećıfico, como el choque en un espacio bidimensional de dos fluidos inmiscibles que se muestra en la Figura 1.1. Los métodos discretizados, entre ellos los Métodos de Elementos Finitos (MEF), los Métodos de Volúmenes Finitos (MVF), y los Métodos de Hidrodinámica de Part́ıculas Suavizadas (SPH, por sus siglas en inglés) han recibido más atención en los últimos dos decenios [53], [107], porque han arrojado buenos resultados sin que el costo computacional resulte tan oneroso que sea poco práctico para la solución de problemas de ingenieŕıa. Se ha probado que los métodos SPH son especialmente versátiles para el modelado de pro- blemas de dinámica de fluidos con alguna frontera abierta como flujo en canales o acciones de vertido. Tienen como ventaja su adaptabilidad relativamente directa a algoritmos de solución que Figura 1.1: Modelado computacional de interacción entre fases inmiscibles 1 2 aprovechan plataformas de cómputo que pueden llevar a cabo procesos de cálculo interdependien- tes pero no necesariamente secuenciales [19]; lo que se conoce como la noción de procesamiento paralelo. Con el advenimiento de procesadores auxiliares tipo GPU que ejecutan cientos o miles de cálculos simultáneos, de bajo costo, compatibles con la arquitectura de computadores per- sonales, tales métodos se han convertido en el sistema de preferencia de modelado de dinámica de fluidos, al reducir los tiempos de convergencia de los métodos numéricos. Además, se facilita la elaboración de código fuente de programación más versátil en caso de que el usuario necesite migrar de plataforma de cómputo [113]. Con el presente estudio se procura aportar a esta ĺınea de investigación al simular mediante el método de hidrodinámica de part́ıculas suavizadas la interacción entre dos fluidos de densidad diferente, en donde uno de ellos es de viscosidad variable y el otro es de viscosidad constante. Esta familia de casos es una a la que el método SPH se ha aplicado para ambas fases con éxito limitado, porque considera únicamente flujo en dos dimensiones y limita la comprobación expe- rimental a casos en que la resistencia a la fluencia es relativamente baja. Además, cuando otros investigadores del método de SPH han tratado casos viscoplásticos con resistencia a la fluencia alta, se han inclinado por modelos h́ıbridos para la fase no lineal; aśı, se tienen combinaciones del método de SPH con MVF [138] o MEF [54] o bien modelos de cuerpos plásticos para la fase sólida [153] (en realidad arena o arcilla con fluido intersticial). En todos los casos se busca una alternativa más eficiente que el uso exclusivo de métodos de malla como MVF o MEF (de uso largamente consolidado en paquetes computacionales comerciales), pues requieren aplicar complejos algoritmos de actualización de malla cuando hay fronteras libres entre fases. La familia de casos que se propone resolver en el presente análisis como objeto de estudio corresponde principalmente al modelado de la interacción de masas de agua con lodos o materia granulada que tienen un comportamiento esencialmente viscoplástico. Esto contribuirá a otras aplicaciones que busquen simular fenómenos de arrastre de sedimentos en canales, transporte de materia granulada, flujo de lodos, y dinámica de la erosión en zonas de derrumbe, entre muchos otros. 3 1.1. Motivación Los métodos de modelado del continuo que se han consolidado en aplicaciones comerciales parten de formulaciones discretizadas asociadas a algún proceso de mallado [107]; por ello un método sin malla como SPH no suele asociarse con paquetes de cómputo comerciales. Sin em- bargo, ofrece ventajas considerables para el tratamiento de fronteras libres y la interacción entre múltiples fases de fluido, sean gaseosas o ĺıquidas, por lo que esta es un área de investigación que ha recibido mucha atención por parte de desarrolladores e investigadores en métodos numéricos aplicados a fluidos. En el último decenio ha sido habitual la alusión a implementaciones novedosas de variantes y modificaciones de SPH para modelado preciso de diferentes fenómenos de flujo, entre las que se destacan el modelado de fenómenos espećıficos de carácter biológico [51], de procesos de manu- factura [114], fenómenos asociados a derrumbes e inundaciones [136] o la interacción de fluidos libres con elementos estructurales [149]. La amplia variedad de parámetros f́ısicos y numéricos implicados en estos asuntos implica que las soluciones deben concentrarse por grupo temático para garantizar su éxito. En el caso del presente estudio, se optó por la interacción de masas de agua con lodos arcillosos o materiales granulados, fenómeno que sirve para analizar procesos de arrastre y erosión de sedimentos, el movimiento de lahares, o el arrastre la carga de fondo en ŕıos y canales. Como queda impĺıcita la presencia de dos fases distintas, necesariamente sus propiedades podrán ser distintas, empezando por la densidad y extendiéndose al otro origen de fuerzas internas: la viscosidad. Es al respecto que surge la consideración indispensable de que la viscosidad de al menos una de las fases será variable según la condición dinámica del fluido. Un lodo arcilloso se resistirá elásticamente al movimiento previo a su punto de fluencia; de igual manera, un material granulado en reposo requerirá una carga mı́nima antes de que inicie su movimiento relativo a fronteras sólidas o material que tampoco ha entrado a una condición dinámica. En suma, la viscosidad de una de las fases se considerará no lineal y la naturaleza de esta reoloǵıa es un factor central al desarrollo de este trabajo. Hay muchos fenómenos de interés cient́ıfico que podrán estudiarse con esta metodoloǵıa, lo que es esencial para comprender fenómenos macroscópicos del flujo de fluidos, pero la meta de 4 este estudio es valorar el potencial que tiene para aplicaciones en ingenieŕıa, algunas de las que se señalan a continuación: Interacción de masas de concentración salina distinta en desembocaduras. Procesos de sedimentación y erosión en desembocaduras. Fenómenos de inundación en terrenos erosionables. Cáıda de deslizamientos de lodos en cuerpos de agua. Procesos de movimiento de tierras en presencia de agua. Arrastre de capas superficiales en suelos irrigados. Comportamiento de trampas de grasa y cajas de registro. Procesos de ensuciamiento de cañeŕıas y desagües. Sin embargo, el modelo numérico no deja de ser una herramienta que se vale de los datos de cualquier evento general en que interactúan dos masas de fluidos con propiedades diferentes, en contacto convectivo y condiciones de contorno cerradas. Fuera de eso, depende casi únicamente de las dificultades de construir la geometŕıa de las fronteras o describir las condiciones iniciales, para convertirse en un problema técnico, a veces muy complejo, pero que no depende de la formulación del modelo numérico del flujo ni de las condiciones de contorno. Como cualquier técnica en su etapa de desarrollo, este trabajo lo propició una serie de retos por afrontar antes de afirmar que una herramienta de este tipo es bastante completa en términos de aplicabilidad a una cierta categoŕıa de fenómenos que pueda ser de interés simular con un método numérico. En el caso de fluidos no lineales, sea en una sola fase o en varias, no solo está la necesidad de implementarlo y cuantificar la precisión esperada respecto a observaciones expe- rimentales, sino detectar aquellas situaciones en que surjan inestabilidades de origen numérico o computacional. Aśı, se pueden agregar algoritmos y estrategias de construcción de modelos que eviten algunos tipos de imprecisión y garanticen, en condiciones claramente delimitadas, que habrá convergencia estable y eficiente al resultado. 5 1.2. Estado del arte En el área de ingenieŕıa, la evolución de un método numérico suele responder al interés par- ticular de un grupo de investigación de usarlo para resolver problemas espećıficos u obtener información adicional para complementar observaciones de campo. Por ello, las aplicaciones de SPH y las formulaciones especiales resultantes evolucionaran de manera distinta, y si bien algu- nas estrategias pueden heredarse de un campo a otro, no es práctico buscar un enfoque universal cuando se trata de un campo de trabajo tan amplio como la dinámica de fluidos. Por ello, esta sección menciona únicamente los estudios en SPH que procuran modelar casos de flujos subsónico con reoloǵıa no lineal; tal es la intención de investigación de este estudio: ofrecer una plataforma de modelado de eventos de fluencia viscoplástica de lodos y granulado cuando interactúan con agua. Desde luego, hay una amplia variedad de casos de estudio y métodos SPH a la medida para flujo potencial, procesos de formación y dinámica de espumas, fenómenos de lubricación intersticial, flujo de poĺımeros en moldes, entre muchos otros, que convencionalmente se conside- ran casos de flujo no lineal. Esta amplia variedad implica comportamientos muy diferentes entre śı, que requieren consideraciones de modelado y de metodoloǵıa particulares para cada cual, lo cual rebasaŕıa los objetivos de este trabajo. La dinámica de fluidos no lineales es una aplicación idónea para los métodos numéricos. Si se considera la formulación propuesta por [89] como la primera formalización generalista de SPH aplicada al problema general de la dinámica de fluidos, tan solo unos pocos años después en [31] se intentó incorporar aspectos no lineales del comportamiento a nivel molecular de los fluidos. Mas adelante, y en especial en el decenio de 2005 a 2015, la cantidad de publicaciones cient́ıficas que se refieren a fenómenos no newtonianos creció significativamente [92], [70]. Algunas revisiones recientes de literatura apuntan a que desde 2018 ha surgido un renovado interés por aplicar SPH a fenómenos de flujo de comportamiento notoriamente no lineal [18], aśı como una atención constante a la interacción entre fluidos y elementos sólidos elásticos [130]. Este es un campo de estudio en que el carácter no lineal de los procesos numéricos lo determinan las fronteras elásticas, no las propiedades f́ısicas del fluido. En todo caso, el tratamiento de flujos en que el fluido es no lineal, por su compresibilidad o por la variación de su viscosidad se le sigue reconociendo como uno de los grandes retos vigentes para los investigadores en SPH [128] por su efecto en la 6 estabilidad y precisión del método. Una particularidad suya es que las publicaciones desde 2007 pueden integrarse entre śı en relación al tópico de modelado de fluidos plásticos usando SPH [115], [68], [108], [105], [51], [40], [145], [142], [96], [66], [50], [118]. La aplicación del método es expĺıcita1, y se proponen muchas variantes y novedosas técnicas de ajuste de parámetros propios del método numérico; también muestran validación con experimentos diseñados especialmente para el efecto. Dentro de este tipo de problemas, entre los primeros intentos por usar el modelo reológico de Herschel-Bulkley (HB) aplicado a SPH para dos fases se encuentra el estudio de [116], pero los resultados carecen de suficiente precisión debido a los modelos de baja resolución que eran computacionalmente factibles en ese momento. Otra aplicación de SPH pero con la aproximación de compresibilidad artificial (Weakly Compressible SPH, frecuentemente denominado por sus siglas WCSPH) lo aplica a un caso particular de sedimentos granulares [110] y emplea el modelo reológico de Bingham y Cross para aproximar una viscosidad local, si bien contrasta los resultados con fotograf́ıas de experimentos en que hay mucha mezcla del agua con la fase sedimentaria, lo que dificulta sacar conclusiones relevantes a valorar la precisión del modelo. Habiendo reconocido las desventajas de usar los modelos reológicos de Bingham y de Mohr- Coulomb para modelar la ı́ndole no lineal del flujo de sedimentos, [40] efectuó una aplicación con SPH usando el modelo reológico de Herschel-Bulkley-Papanastasiou (HBP) y lo aplicó a un caso bidimensional bifásico agua-arena. La meta consistió en modelar procesos de erosión, suspensión y sedimentación de granulados causados por la súbita cáıda de una masa de agua sobre un lecho de arena. Otro estudio derivado para modelar el mismo tipo de fenómeno [153] aplicó un enfoque SPH viscoplástico a algunos casos de flujo bidimensional bifásico en que el interés central era pronosticar procesos de transporte por convección de la fase granular. En tal estudio, el modelo de Bingham se aplica a la fase granular, con la adición de un ĺımite elástico variable para esa fase como función de la presión local. Estos casos bidimensionales muestran buena concordancia con los datos experimentales de perfiles de erosión y asentamiento, pero no resuelven ciertas dificultades de estabilidad numérica. Otros estudios, como [62] y [80] exploran variantes para la aplicación de un criterio de fluencia 1Al encontrarse fuera de los alcances de este trabajo, no se mencionan publicaciones que implican hibridización con otros métodos discretizados, dado que esa técnica se aplica para lidiar con la deformación o transporte de objetos sólidos no particulados, aśı como condiciones de flujo laminar o supersónico. 7 para la fase no lineal de casos similares a los anteriores en que son bidimensionales y bifásicos. Finalmente, con un propósito y una configuración de casos similares, en [42] se agrega la solución de una ecuación de Laplace para calcular un esfuerzo equivalente en la fase no lineal, con el fin de combinar un modelo elástico con una formulación SPH, solo aplicada a la parte de las fases que se comporta como un fluido, sea la granular después de haber iniciado fluencia o el agua con la que interactúa. Las aplicaciones posteriores de esta combinación de técnicas a casos bidimensionales espećıficos de ruptura de presas en dos fases, muestran diferentes aproximaciones sobre la dinámica de flujo exhibida por los modelos numéricos, como la ruptura de una presa de lodo bajo el agua [118], o la ruptura de una presa con fondo sedimentado [135] con alguna comparación con datos experimentales o incluso la presentación de casos tridimensionales con resolución espacial limitada. En estos estudios, que se ocupan de la arena y otras fases granulares, el uso del modelo reológico de Bingham es más común. Para acoplarlos a los métodos SPH bidimensionales, los modelos reológicos de HB o HBP se han usado para estudiar casos monofásicos, como en [47]. En publicaciones que tratan el tema de flujo bifásico pero con los métodos tradicionales de métodos de malla en lugar de SPH, es habitual encontrar estudios recientes que igualmente acuden a estos modelos reológicos; por ejemplo, [56] usa el método de volumen de fluido (VOF) junto con la ecuación de HBP, con éxito moderado. Otro ejemplo, en una sola fase de un lodo arcilloso, es el de [67] que aplica el MVF al modelado del flujo. En tales estudios, como en otros que emplean métodos de malla a la fase fluida, se concentran en la dinámica de la interacción entre el agua como fase lineal y la arena húmeda saturada o el lodo arcilloso como fase no lineal. Otro enfoque que afronta el modelado de flujo granular en suspensiones usando SPH conside- ra que la fase fluida se mezcla completamente con el sólido granulado. Según esta consideración, ambas fases tienen sus propias ecuaciones constitutivas, de continuidad y de conservación afec- tadas por una proporción de presencia de fase i en la mezcla cuando corresponda. Cualquier transferencia de enerǵıa que ocurra entre fases se da mediante un coeficiente de arrastre efectivo entre part́ıculas sólidas y fluidas. Este concepto es mejor conocido en [127] y se ha aplicado en otros desarrollos recientes como lo describe [147], siempre para casos bidimensionles. Las aplicaciones en tres dimensiones de un modelo numérico con una fase no lineal se men- cionan con frecuencia, especialmente porque la mayoŕıa de las formulaciones SPH se pueden 8 implementar fácilmente para tal situación, aunque sea limitadas por el costo computacional u otros problemas prácticos referidos a la estabilidad numérica. Algunos autores prefieren mostrar resultados con casos simples de resolución limitada, como el modelo de flujo de detritos monofási- co de [48]; no obstante, debido a las dificultades en la configuración y medición de experimentos f́ısicos, la caracterización cuantitativa de estos modelos parece ser muy escasa, sea dentro del aplicación de SPH o con otros métodos numéricos. Con alguna frecuencia, como en el estudio de [126] o el de [23], aunque se presenta un modelo tridimensional, la geometŕıa del flujo es predominantemente bidimensional para hacerla coincidente con la configuración del experimento f́ısico. Otros estudios, centrados en la generación y evolución de ondas superficiales en el fluido, son bifásicos pero con la consideración de que SPH sólo procede para la fase lineal, mientras que en la fase no lineal se aplica un modelo totalmente diferente, como el trabajo de [68], que usa un modelo de Boussineq para evolucionar la dinámica de la onda en movimiento. cuerpo de suelo submarino. Otros casos que combinan de manera similar SPH lineal con aproximaciones no lineales para el comportamiento de una fase elástica o perfectamente plástica, parten del Método de Elementos Discretos (MED) [131], MVF [11] o el de Sistema de Part́ıculas en Movimiento (SPM) [57] para desarrollar las ecuaciones constitutivas y la dinámica correspondiente. Aśı, en este tipo de solución en rigor no se está aplicando SPH a la fase no lineal, de manera que sus resultados son de interés para hacer cotejos cuantitativos, pero no proponen novedad en la utilización de SPH para la fase no lineal. Dado su impacto en la viscosidad equivalente de la fase no lineal, otra ĺınea de investigación con implicaciones para el presente trabajo tiene que ver con la modernización del método me- diante propuestas que buscan eliminar el término de viscosidad artificial que convencionalmente se agrega a la formulación SPH para garantizar estabilidad [4]. En muchos casos en que las fuerzas viscosas reales son pequeñas en comparación con las cargas inerciales, no le disminuye precisión al modelo; empero, cuando las fuerzas viscosas son grandes o bien su valor adelantará o atrasará fenómenos macro del flujo, pueden verse afectados negativamente los resultados. En si- tuación similar se cuentan las aproximaciones que proponen técnicas para mejorar la estabilidad del método de SPH con fenómenos de impacto o fronteras con discontinuidades muy marcadas [150], como ocurre en el caso de aguas poco profundas o en la frontera entre fluidos de fases diferentes [112],[147]. 9 Son abundantes las publicaciones que se pueden clasificar, en términos generales, en la cate- goŕıa de hibridización de SPH con otros métodos, sean de malla o de part́ıculas. Al considerar una muestra arbitraria de art́ıculos cient́ıficos fechados entre 2005 y 2015, que tratan directa- mente el tema de SPH2, predominan los tópicos de aplicaciones especiales y su validación, la proposición de variantes especiales de SPH3, y la hibridización del método con otras técnicas de modelado numérico como MVF, MEF, MDF y el Método de Elementos de Contorno (MEC), entre otras [114]. En este trabajo no se proponen esquemas h́ıbridos, por lo que no se proveerá mayor detalle, pero śı es una vertiente de acción muy rica en la actualidad. En concordancia con el propósito del presente estudio y con los resultados a los que se llegó al término de la investigación, esta sección sobre el estado del arte se elaboró orientada a los modelos WCSPH en que se da preferencia a un enfoque de variación de la viscosidad equivalente como aproximación al comportamiento de la fase no lineal. Por supuesto, muchos de los estudios mencionados en esta sección de casos bidimensionales abren interrogantes sobre la influencia de los efectos de pared, las interacciones en el ĺımite entre las fases, la aparición de inestabilidades de presión y el error en la forma de la interfaz debido a artefactos numéricos. Dado que tales preguntas se pueden aplicar a una amplia variedad de problemas, este trabajo se centró en ofrecer soluciones a los retos abiertos para aplicaciones práctica que tiene este esquema general de modelado. 1.3. Objetivos e hipótesis El objetivo general de esta investigación consiste en formular un modelo tridimensional de dinámica de fluidos mediante un método de hidrodinámica de part́ıculas suavizadas (SPH), que simule la interacción entre dos fluidos inmiscibles de densidad distinta, uno de ellos de compor- tamiento no lineal. Para alcanzarlo, se ha formulado una metodoloǵıa orientada por los siguientes objetivos espećıficos: 2Dado el enfoque de ciencia aplicada a la ingenieŕıa de este trabajo y para mantener esfuerzos al aspecto de novedad requerido para el estudio, no se explorarán publicaciones que analizan el método desde un punto de vista estrictamente matemático. 3Las que más mención reciben son SPH con factor desconocido (XSPH) [88], �-SPH [7] y SPH-✏ [91], entre otros. 10 Fundamentar el modelo matemático mediante la formulación de las ecuaciones de estado para las respectivas fases de fluido. Transformar las ecuaciones de continuidad y de estado, para implementarlas conforme a las estrategias de modelado del método de part́ıculas hidrodinámicas libres de malla (SPH). Formular un algoritmo de solución de los sistemas de ecuaciones resultantes compatible con una estrategia de solución computacionalmente paralelizable. Validar los resultados de las simulaciones con casos de estudio particulares previamente publicados y con datos experimentales de carácter espećıfico. Estos objetivos se plantean partiendo de que a lo largo de la investigación seŕıan válidas las siguientes hipótesis de trabajo, algunas de las que, de comprobarse con prueba sustancial, podrán conformar el aporte de novedad de este estudio: El método de SPH seŕıa aplicable a la interacción tridimensional no miscible de dos fluidos distintos sin que resulten sistemas de ecuaciones computacionalmente inestables. Es factible formular un algoritmo de solución en paralelo al sistema de ecuaciones que resulte de la aplicación del método de SPH a fluidos con comportamiento no lineal. El sistema de ecuaciones seleccionadas constituye un modelo congruente con la condición dinámica macroscópica de la interacción tridimensional de dos fluidos inmiscibles con ca- racteŕısticas f́ısicas determinadas dentro del margen de error aceptable para aplicaciones en ingenieŕıa. La arquitectura de computadores personales actuales tiene capacidad de procesamiento de instrucciones suficiente para llevar a buen término una implementación del algoritmo de solución para un caso de estudio de interés particular. Para circunscribir la investigación a los ĺımites pertinentes a un trabajo de ı́ndole doctoral, y en concordancia con la naturaleza aplicada de un estudio en el área de ingenieŕıa, estos objetivos e hipótesis no se ocupan de considerar aspectos de completitud, generalidad, ni consistencia numérica, sino que se centra en comprobar su utilidad para una cierta familia de fenómenos de 11 flujo, claramente identificada a la hora de seleccionar los casos muestra. Ademas, los estudios sobre la sensibilidad a parámetros numéricos o f́ısicos presentes en el modelo se harán dentro de la variabilidad aplicable a los fenómenos estudiados en esta investigación, que podrán ser extensibles a otras familias de casos pero para los que no hay garant́ıa de que necesariamente producirán igual rendimiento en términos de estabilidad o precisión. 1.4. Estructura de la tesis Este documento es un informe final de investigación para el que se ha elegido una estructura de tesis doctoral en el área de métodos numéricos aplicados a simular problemas abiertos de ingenieŕıa. En medida de ello, lo conforman las secciones que se describen a continuación. El primer caṕıtulo ofrece una sinopsis general del propósito y contexto general del estudio, aśı como un estado del arte en el área de estudio espećıfico al objetivo general, con especial énfasis en trabajos previos en modelado por WCSPH de flujo bifásico no lineal que son el antecedente inmediato al presente trabajo. Se afirman de manera expĺıcita los objetivos de investigación y las hipótesis de trabajo que guiaron todo el proceso. Seguido está el segundo caṕıtulo, que conforma el marco teórico de este trabajo, se presenta el origen de las ecuaciones fundamentales de continuidad y de conservación de la cantidad de movimiento de la Mecánica de Fluidos, se describe la formulación básica del método de SPH con compresibilidad artificial, y se muestran los modelos reológicos aplicables a algunos tipos de comportamiento no lineal de los fluidos. Con esto se da paso al tercer caṕıtulo, que corresponde a la metodoloǵıa aplicada para cumplir con los objetivos de esta investigación, lo que incluye el proceso de implementación de un progra- ma que aplica el método WCSPH en lenguaje C y la descripción de los algoritmos y extensiones para lidiar con flujos no lineales. En el cuarto caṕıtulo se describen los casos de estudio utilizados para generar una serie de resultados numéricos para la validación del método, y se analizan centrados en aspectos parti- culares de la metodoloǵıa que se propone. Estos experimentos incluyen casos convencionales de hidráulica de una sola fase lineal, otros con una fase no lineal y multifásico con una fase no lineal; la mayor parte son flujos bidimensionales, pero se presentan algunos casos tridimensionales. Para 12 casi todos se dispońıa de datos provenientes de experimentos f́ısicos, pero de especial valor es la publicación de la base de datos de rotura de presa bifásica de [134]. El análisis de caso por caso se cierra con una sección dedicada a comentar aspectos comunes del proceso de convergencia, ajuste de parámetros numéricos, y del rendimiento del modelo propuesto. En este caṕıtulo también se tratan temas derivados del proceso de implementación del método de cálculo. El cuerpo del documento se cierra con el quinto caṕıtulo, correspondiente a las conclusiones, en que se hace una valoración final de los aportes del trabajo y se proponen ĺıneas de investigación futuras con las que se puede contribuir directamente al desarrollo de la técnica de SPH a nivel internacional. En la sección de anexos se incluyen los manuscritos de las publicaciones derivadas de esta investigación y el código fuente original del autor en lenguaje de programación C. Caṕıtulo 2 Marco teórico A continuación se ofrece un desarrollo de las ecuaciones de conservación de masa y de conser- vación de la cantidad de movimiento, en coordenadas cartesianas tridimensionales para fluidos compresibles viscosos sin cambio de fase; con ello, se establece el marco de trabajo en términos de fundamentación f́ısica del método. Luego, se presenta la formulación SPH con compresibili- dad artificial (WCSPH) aplicada a esas ecuaciones de conservación junto con las consideraciones metodológicas convencionales que requiere el método para su implementación práctica. Final- mente, se muestran los modelos reológicos aplicables al tipo de flujo no lineal de interés para este trabajo. 2.1. Fundamentación matemática de la mecánica de fluidos La formulación de modelos matemáticos para fenómenos de dinámica de fluidos se puede orientar a casos muy espećıficos para obtener información anaĺıtica. Esto es útil, aunque con el inconveniente de restringirse a aplicaciones muy sencillas o delimitada por contornos ideales. En cursos básicos de mecánica de fluidos [64], estos modelos se orientan a crear estrategias de pensamiento y criterio anaĺıtico, y en términos generales se incluyen los casos de placas paralelas, tubeŕıa circular, tambores rotatorios y discos rotatorios, entre otros. Las ecuaciones resultantes de los análisis son demostrables anaĺıticamente mediante un razonamiento espećıfico, aunque tam- bién se pueden derivar como simplificaciones de las ecuaciones generales de dinámica de fluidos. Estas ecuaciones generales, que rigen para una masa de fluido continuo, se formulan mediante un 13 14 análisis de una masa afectada por los fenómenos f́ısicos correspondientes a su comportamiento observado. Una presentación convencional del desarrollo de las ecuaciones generales es la aproximación diferencial [8], [45], que deriva para una part́ıcula diferencial las correlaciones producto de la aplicación de las leyes de conservación de la masa y de conservación de la cantidad de movimiento. La aproximación integral es igualmente válida y su desarrollo por esa forma produce las mismas ecuaciones diferenciales de flujo de fluidos [8], [64]. En esta sección se muestra de forma sintética el desarrollo de las ecuaciones generalizadas de la mecánica de fluidos compresibles1, conocidas como ecuaciones de Navier-Stokes para flujo compresible. Esto en parte sirve de fundamento al método de SPH en que aplican como parte del método, pero permite su trazabilidad al origen de cada uno de los términos, con lo que podrá estimarse su impacto en la estabilidad del método numérico y en la exactitud de las soluciones numéricas. 2.1.1. Estado de esfuerzos de un fluido Si se toma como punto de partida el análisis de las fuerzas que actúan en un elemento infinitesimal de fluido tridimensional ubicada en el espacio, es posible desarrollar un conjunto de ecuaciones fundamentales, que en este caso se elaborará para un sistema de coordenadas cartesianas. Este elemento, de dimensiones dx, dy, y dz está localizada en la posición x = [x, y, z] y se mueve a una velocidad u = [u, v, w] con aceleración Du Dt (derivada total de la velocidad en el tiempo) por la acción de la suma vectorial de fuerzas internas fi y las fuerzas externas fe. Esta condición cinética se completa con escalares de interés asociados a la condición termodinámica de la part́ıcula, entre los que se encuentra su viscosidad µ, densidad ⇢, temperatura �, enerǵıa interna e y entroṕıa h. Las fuerzas internas corresponden a la acción resultante de los aportes de las presiones y de la fricción viscosa en las caras. La presión, que actúa en dirección perpendicular a las caras del elemento genera los esfuerzos axiales netos �xx, �yy, �zz. Por otra parte, la fricción viscosa, que actúa en dirección paralela a las caras del elemento, genera en cada plano esfuerzos cortantes 1Si bien el caso general de interés de este trabajo está en el ámbito subsónico y además no es un gas, puede suponerse flujo incompresible; sin embargo, es importante conocer qué términos de las ecuaciones generalizadas podŕıan verse afectados por incorporar una compresibilidad artificial, que en WCSPH es la Ecuación de Estado (denominada frecuentemente como EOS por sus siglas en inglés) y forma parte indispensable de la solución como se verá en la sección 2.2.9. 15 Figura 2.1: Estado de esfuerzos de un fluido netos con componentes ortogonales: ⌧yz y ⌧zy en el plano yz; ⌧yx y ⌧xy en el plano xy; y los esfuerzos cortantes netos ⌧xz y ⌧zx en el plano xz. Cada uno de estos esfuerzos se indican relativo a la geometŕıa del elemento cuadrangular en la Figura 2.1. En concordancia con tales definiciones, el estado de esfuerzos que afecta el elemento, que incluye el término de presión p, los esfuerzos viscosos axiales ⌧ii y los esfuerzos viscosos cortantes ⌧ij podŕıa expresarse de manera compacta como � = pI+ ⌧ . Si se suman la presión y el esfuerzo viscoso axial en un solo término �ii = p+ ⌧ii, el estado de esfuerzos se podŕıa describir con una matriz como la que aparece en la ec. (2.1). � = 2 66664 �xx ⌧xy ⌧xz ⌧yx �yy ⌧yz ⌧zx ⌧zy �zz 3 77775 (2.1) Por otra parte, las fuerzas externas corresponden a la acción de la gravedad u otras acelera- ciones2, a reacciones electromagnéticas y otros efectos proporcionales a la masa de la part́ıcula que no sean los ya mencionados hasta aqúı. 2Algunos autores se refieren a estas acciones como cargas inerciales, lo que no es conceptualmente correcto por cuanto la inercia de un cuerpo no es una fuerza, pero es conveniente tratar el término como una carga para efectos de cálculo. En este trabajo se evitará esta denominación. 16 Figura 2.2: Volumen de control para conservación de masa 2.1.2. Continuidad y conservación de la masa La condición de carga de la part́ıcula n de la Figura 2.1 se presta a un desarrollo convencional de las ecuaciones generales de la mecánica de fluidos. En el presente estudio se utiliza la aproxi- mación diferencial aplicando las leyes de conservación de la masa y de conservación de la cantidad de movimiento [8], [45] para obtener las ecuaciones que describiŕıan la dinámica del flujo. En esta sección se desarrolla la ecuación de conservación de masa o ecuación de continuidad. Según se muestra en la Figura 2.2, dado que el volumen de control rectangular de la part́ıcula de fluido de volumen dV = dx dy dz y masa dm = ⇢ dV , se puede obtener la cantidad de masa que atraviesa cada una de las caras, multiplicando la densidad del fluido por el caudal, que seŕıa a su vez la velocidad del flujo multiplicada por el área de sección transversal. Para las caras perpendiculares al eje x, la masa de entrada seŕıa ⇢u dy dz, mientras que en la salida la densidad habrá variado a ⇣ ⇢+ @⇢ @xdx ⌘ y la velocidad será � u+ @u @xdx � . Este razonamiento es aplicable a las caras en los otros ejes, de manera que el flujo de salida correspondiente a los ejes x, y y z se expresaŕıan con las ecs. (2.2), (2.3) y (2.4) respectivamente. 17 ✓ ⇢+ @⇢ @x dx ◆✓ u+ @u @x dx ◆ dy dz ⇡ ✓ ⇢u+ @ @x (⇢u) dx ◆ dy dz (2.2) ✓ ⇢+ @⇢ @y dy ◆✓ v + @v @y dy ◆ dx dz ⇡ ✓ ⇢v + @ @y (⇢v) dy ◆ dx dz (2.3) ✓ ⇢+ @⇢ @z dz ◆✓ w + @w @z dz ◆ dx dy ⇡ ✓ ⇢w + @ @z (⇢w) dz ◆ dx dy (2.4) El lado derecho de las ecuaciones anteriores es una simplificación del lado izquierdo, para lo que se han despreciado los términos diferenciales de orden superior. Por otra parte, la razón de acumulación de masa en el volumen de control se puede escribir en función de la densidad según la ec. (2.5). @m @t = @ @t ⇢ (dx dy dz) = @⇢ @t dV (2.5) Con estos componentes, se puede escribir la ecuación de balance de masa del volumen de control; en suma, que el flujo másico de entrada equivale a la razón de acumulación de masa sumado al flujo másico de salida, con lo que queda la ec. (2.6). ⇢u dy dz + ⇢ v dx dz + ⇢w dx dy = @⇢ @t dV + ✓ ⇢+ @⇢ @x dx ◆✓ u+ @u @x dx ◆ dy dz + ... ...+ ✓ ⇢+ @⇢ @y dy ◆✓ v + @v @y dy ◆ dx dz + ✓ ⇢+ @⇢ @z dz ◆✓ w + @w @z dz ◆ dx dy (2.6) Si todos los términos se transfieren al lado derecho y se desarrollan los productos despreciando los términos de segundo orden, queda la ec. (2.7), que a su vez se puede expresar de una forma más compacta como aparece en la ec. (2.8) a sabiendas de que el grupo dx dy dz = dV , que seŕıa un diferencial de volumen. 0 = ✓ ⇢ @u @x + u @⇢ @x ◆ + ✓ ⇢ @v @y + v @⇢ @y ◆ + ✓ ⇢ @w @z + w @⇢ @z ◆� dx dy dz + @⇢ @t dx dy dz (2.7) 18 Figura 2.3: Volumen de control para conservación de cantidad de movimiento  @⇢ @t + @ @x (⇢u) + @ @y (⇢ v) @ @z (⇢w) � dV = 0 (2.8) Esta es la ecuación de continuidad para el caso de un fluido compresible en estado transitorio. El término dV desaparece y con ello se puede expresar en notación vectorial como en la ec. (2.9). @⇢ @t +r · (⇢u) = 0 (2.9) 2.1.3. Conservación de la cantidad de movimiento A partir del estado de esfuerzos descrito en la Figura 2.3 se elabora una ecuación para agregar la conservación de la cantidad de movimiento [8], [64], [45] si se le aplica la segunda Ley de Newton según se expresa en la ec. (2.10) en donde la notación D Dt se refiere a la derivada total. Las fuerzas fe y fi corresponden a las cargas netas volumétricas (reacción a aceleraciones) y superficiales (presión y viscosidad), respectivamente. ⇢ dV Du Dt = fe + fi (2.10) Para expresar las cargas triaxiales en cada una de las seis caras del volumen de control, 19 conviene emplear la notación tensorial, con lo que es posible definir tensores de esfuerzos axiales �ij (que incluye términos de presión y de esfuerzos viscosos), el de viscosidad cortante ⌧ij y el de razón de deformación3 " como se enuncia en las ecs. (2.11), (2.12) y (2.13). �ij = 2 66664 �xx �xy �xz �yx �yy �yz �zx �zy �zz 3 77775 (2.11) ⌧ij = 2 66664 ⌧xx ⌧xy ⌧xz ⌧yx ⌧yy ⌧yz ⌧zx ⌧zy ⌧zz 3 77775 (2.12) ✏ij = 2 66664 ✏xx ✏xy ✏xz ✏yx ✏yy ✏yz ✏zx ✏zy ✏zz 3 77775 = 2 66664 @u @x 1 2 ⇣ @u @y + @v @x ⌘ 1 2 � @u @z + @w @x � 1 2 ⇣ @v @x + @u @y ⌘ @v @y 1 2 ⇣ @v @z + @w @y ⌘ 1 2 � @w @x + @u @z � 1 2 ⇣ @w @y + @v @z ⌘ @w @z 3 77775 (2.13) Las cargas que actúan en cada uno de los ejes coordenados se pueden analizar por separado. Inspeccionando en primera instancia las fuerzas de superficie que tienen componente en el eje x, se puede ver que hay una carga FpA por acción de la presión, una carga FpB por acción axial de la viscosidad, y cuatro cargas Fma, Fmb, Fmc y Fmd que se deben a acción cortante de la viscosidad. Cada una se calcula como el producto del área del elemento por el esfuerzo efectivo en esa cara; el término de esfuerzo, como podrá variar en el continuo, se obtiene como una aproximación de Taylor de primer orden centrada en los elementos de �ij que están presentes. El resultado de este razonamiento queda descrito en las ecs. (2.14), (2.15), (2.16), (2.17), (2.18) y (2.19). FpA = � ✓ �xx � dx 2 @�xx @x ◆ dy dz (2.14) FpB = ✓ �xx + dx 2 @�xx @x ◆ dy dz (2.15) 3En este desarrollo es importante tener en mente que este término contiene razones de deformación axial "ii y razones de deformación angular cortante "ij para i 6= j. 20 Fma = � ✓ �yx � dy 2 @�yx @y ◆ dx dz (2.16) Fmb = ✓ �yx + dy 2 @�yx @y ◆ dx dz (2.17) Fmc = � ✓ �zx � dz 2 @�zx @z ◆ dx dy (2.18) Fmd = ✓ �zx + dz 2 @�zx @z ◆ dx dy (2.19) Si se agrega la componente en dirección x de la carga volumétrica que afecta al fluido, se pueden ingresar las equivalencias de las ecs. (2.14) a (2.19) en la ec. (2.10) para obtener la componente en x de la ecuación de equilibrio dinámico. Con fines ilustrativos, se partiŕıa de que la componente de la gravedad en dirección x ejerce la única carga volumétrica, que equivaldŕıa a ⇢ dV gx, con lo que queda preliminarmente la ec. (2.20); esta expresión entonces se transforma en la ec. (2.21). ⇢ dV Du Dt = ⇢ dV gx + (FpA + FpB + Fma + Fmb + Fmc + Fmd); (2.20) ⇢ dx dy dz Du Dt = ⇢ dx dy dz gx + ✓ @�xx @x + @�yx @y + @�zx @z ◆ dx dy dz (2.21) Al eliminar el grupo dx dy dz y desarrollar la derivada total Du Dt , queda la ec. (2.22), que será la componente x de la ecuación de conservación de cantidad de movimiento para un fluido en coordenadas cartesianas tridimensionales. Las componentes y y z se pueden escribir por analoǵıa para completar el conjunto con las ecs. (2.23) y (2.24) correspondientes. Estas tres ecuaciones se conocen también como la ecuación de Cauchy. ⇢ ✓ @u @t + u @u @x + v @u @y + w @u @z ◆ = ⇢ gx + ✓ @�xx @x + @�yx @y + @�zx @z ◆ (2.22) 21 ⇢ ✓ @v @t + u @v @x + v @v @y + w @v @z ◆ = ⇢ gy + ✓ @�xy @x + @�yy @y + @�zy @z ◆ (2.23) ⇢ ✓ @w @t + u @w @x + v @w @y + w @w @z ◆ = ⇢ gz + ✓ @�xz @x + @�yz @y + @�zz @z ◆ (2.24) Teniendo la previsión de separar las componentes de presión de las componentes de viscosidad previo al ensamble de la ec. (2.22), quedaŕıa, para las fuerzas de superficie en dirección X las ecs. (2.25), (2.26) y (2.27). @�xx @x = �@P @x + @⌧xx @x (2.25) @�yx @y = @⌧yx @y (2.26) @�zx @z = @⌧zx @z (2.27) Las componentes en los ejes Y y Z reciben tratamiento análogo, con lo que las ecs. (2.22), (2.23) y (2.24) quedan transformadas al conjunto de ecs. (2.28), (2.29) y (2.30). ⇢ ✓ @u @t + u @u @x + v @u @y + w @u @z ◆ = ⇢ gx � @P @x + ✓ @⌧xx @x + @⌧yx @y + @⌧zx @z ◆ (2.28) ⇢ ✓ @v @t + u @v @x + v @v @y + w @v @z ◆ = ⇢ gy � @P @y + ✓ @⌧xy @x + @⌧yy @y + @⌧zy @z ◆ (2.29) ⇢ ✓ @w @t + u @w @x + v @w @y + w @w @z ◆ = ⇢ gz � @P @x + ✓ @⌧xz @x + @⌧yz @y + @⌧zz @z ◆ (2.30) Si se desarrolla el tensor ⌧ij , asumiendo que la viscosidad dinámica µ es la constante de proporcionalidad entre entre la resistencia al corte y la razón de deformación (en la ec. (2.13) esto se hab́ıa definido de previo), se podrán sustituir sus componentes en las ecs. (2.28), (2.29) y (2.30). Conviene señalar que en este desarrollo se ha decidido conservar el término � que es la 22 viscosidad volumétrica4 con el propósito de que posteriormente se pueda valorar su significancia en comparación con la viscosidad artificial que implica el uso del método de WCSPH. ⌧ = 2 66664 2µ@u @x + � ⇣ @u @x + @v @y + @w @z ⌘ µ ⇣ @u @y + @v @x ⌘ µ � @u @z + @w @x � µ ⇣ @u @y + @v @x ⌘ 2µ@v @y + � ⇣ @u @x + @v @y + @w @z ⌘ µ ⇣ @v @z + @w @y ⌘ µ � @u @z + @w @x � µ ⇣ @v @z + @w @y ⌘ 2µ@w @z + � ⇣ @u @x + @v @y + @w @z ⌘ 3 77775 (2.31) De esta manera, quedan las ecuaciones de Navier-Stokes para flujo compresible, variante bastante general de la ecuación de conservación de la cantidad de movimiento de la dinámica de fluidos. ⇢ ✓ @u @t + u @u @x + v @u @y + w @u @z ◆ = ⇢gx � @P @x + @ @x ✓ 2µ @u @x + � ✓ @u @x + @v @y + @w @z ◆◆ + @ @y  µ ✓ @u @y + @v @x ◆� + @ @z  µ ✓ @u @z + @w @x ◆� (2.32) ⇢ ✓ @v @t + u @v @x + v @v @y + w @v @z ◆ = ⇢gy � @P @y + @ @x  µ ✓ @u @y + @v @x ◆� + @ @y ✓ 2µ @v @y + � ✓ @u @x + @v @y + @w @z ◆◆ + @ @z  µ ✓ @v @z + @w @y ◆� (2.33) ⇢ ✓ @w @t + u @w @x + v @w @y + w @w @z ◆ = ⇢gz � @P @z + @ @x  µ ✓ @u @z + @w @x ◆� + @ @y  µ ✓ @v @z + @w @y ◆� + @ @z ✓ 2µ @w @z + � ✓ @u @x + @v @y + @w @z ◆◆ (2.34) 4La viscosidad volumétrica o segunda viscosidad se refiere a las fuerzas disipativas que dependen de la razón de expansión o compresión del fluido. En fluidos incompresibles, este término se puede despreciar, pero debe tenerse en cuenta para modelar fenómenos de atenuación de ondas en gases o para el modelado de flujos bifásicos ĺıquido-gas donde las velocidades son cercanas o superiores a las del sonido en el medio[64]. 23 Para más claridad, conviene señalar que estas ecuaciones se pueden escribir en una notación más compacta usando las definiciones de derivada total D Dt y acudiendo al operador r. Para los términos multiplicados por la densidad, si se considera que las componentes de la velocidad u son la derivada en el tiempo de la posición en cada eje correspondiente se tiene que u = @x @t ; v = @y @t ; y w = @z @t ; aśı, es válida la equivalencia descrita por la ec. (2.35), ⇢ 0 BBBB@ @u @t + @x @t @u @x + @y @t @u @y + @z @t @u @z @v @t + @x @t @v @x + @y @t @v @y + @z @t @v @z @w @t + @x @t @w @x + dy dt @w @y + @z @t @w @z 1 CCCCA = ⇢ 0 BBBB@ Du Dt Dv Dt Dw Dt 1 CCCCA = ⇢ Du Dt (2.35) El vector f e, definido en la ec. (2.36) contiene todas las fuerzas de origen externo, que en este caso únicamente se refiere a la carga gravitacional en las componentes que pueda tener en el sistema de coordenadas elegido. fe = ⇢ 0 BBBB@ gx gy gz 1 CCCCA (2.36) Los términos que contienen la viscosidad µ requieren un tratamiento distinto. Si se desarrollan las derivadas como en la ec. (2.37), se pueden separar en dos componentes para cada uno de los ejes coordenados, como aparece en la ec. (2.38). µ 0 BBBB@ @ @x � 2@u @x � + @ @y ⇣ @u @y + @v @x ⌘ + @ @z � @u @z + @w @x � @ @x ⇣ @u @y + @v @x ⌘ + @ @y ⇣ 2@v @y ⌘ + @ @z ⇣ @v @z + @w @y ⌘ @ @x � @u @z + @w @x � + @ @y ⇣ @v @z + @w @y ⌘ + @ @z � 2@w @z � 1 CCCCA = µ 0 BBBB@ 2@2u @x2 + @2u @y2 + @2v @y@x + @2u @z2 + @2w @x@z @2v @x2 + @2u @y@x + 2@2v @y2 + @2v @z2 + @2w @y@z @2w @x2 + @2w @y2 + @2u @x@z + @2v @y@z + 2@2w @z2 1 CCCCA (2.37) µ 0 BBBB@ @2u @x2 + @2u @y2 + @2u @z2 + @2u @x2 + @2v @y@x + @2w @x@z @2v @x2 + @2v @y2 + @2v @z2 + @2u @y@x + @2v @y2 + @2w @y@z @2w @x2 + @2w @y2 + @2w @z2 + @2u @x@z + @2v @y@z + @2w @z2 1 CCCCA = µ 0 BBBB@ r2 u+ @ @x ⇣ @u @x + @v @y + @w @z ⌘ r2 v + @ @y ⇣ @u @x + @v @y + @w @z ⌘ r2 w + @ @z ⇣ @u @x + @v @y + @w @z ⌘ 1 CCCCA (2.38) 24 Como la derivada del escalar que resulta de la divergencia es cero, quedan únicamente los términos que aparecen en la ec. (2.39). µ 0 BBBB@ r2 u+ @ @x (r · u) r2 v + @ @y (r · u) r2 w + @ @z (r · u) 1 CCCCA = µ 0 BBBB@ r2 u r2 v r2 w 1 CCCCA = µr2u (2.39) Aplicando estas equivalencias, se puede escribir la ec. (2.40), en la que se pueden identificar todos los términos presentes en las ecs. (2.32), (2.33) y (2.34). ⇢ Du Dt = f �rP + �(r · u)13 + µr2u (2.40) Por tener un carácter general, estos modelos matemáticos tienen la versatilidad conceptual para modelar cualquier fenómeno conocido de dinámica de fluidos dentro de las definiciones y las suposiciones que se hayan expuesto. Existe el inconveniente de que no hay solución anaĺıtica conocida, ni siquiera para las familias de problemas que ocurren con más frecuencia en aplicacio- nes de ingenieŕıa. Esto se hace aun más complejo al agregar condiciones de frontera reales, que con frecuencia implican discontinuidades matemáticas tanto en las fronteras como en el medio. 2.1.4. Ecuación de vorticidad La ecuación de cantidad de movimiento, según se muestra en este trabajo, implica que la cantidad de movimiento angular se conserva, pero sólo porque se ha usado la aproximación de considerar el flujo como no rotacional (cada elemento diferencial con velocidad angular ✓̇ = 0). Aunque se pueden formular las ecuaciones adicionales para considerar intercambio de enerǵıa rotacional, los métodos generales de SPH utilizados en este trabajo no lo contemplan, además de que tal simplificación no afecta en demaśıa el comportamiento macroscópico de un fluido [64] en las escalas de tamaño y velocidad delimitados por los objetivos de esta investigación5. No es indispensable considerar estos grados de libertad adicionales (y la consecuente complejidad adicional de solución, cualesquiera que sea su naturaleza), pero esto no anula la posibilidad de 5La formación de vórtices a microescala en este momento es un problema abierto en dinámica de fluidos, dado que es un fenómeno en donde las fuerzas electromagnéticas intermoleculares son de magnitud comparable a las fuerzas viscosas e inerciales. 25 aplicar modelos de turbulencia6 o cuantificar el tamaño y velocidad de vórtices macroscópicos en el fluido. Por eso, aun para el caso de flujo no rotacional, es conveniente calcular a partir del campo vectorial de velocidades un campo de potencial de rotación conocido como la vorticidad. La ecuación de vorticidad [8], se puede derivar de la ec. (2.40), de conservación de la cantidad de movimiento lineal, sabiendo que la definición de vorticidad ! es el rotacional7 de la velocidad u según se describe en la ec. (2.41). ! = r⇥ u (2.41) Si se calcula el rotacional de la ecuación de conservación de la cantidad de movimiento para flujo compresible viscoso, se obtiene [45] la ec. (2.42), recordando que ⌧ es el tensor de esfuerzos cortantes por acción de la viscosidad µ y de la segunda viscosidad �. @! @t +r⇥ [(u ·r)u] = (! ·r)� !(r · u) + 1 ⇢2 r⇢⇥rP +r⇥ ✓ r · ⌧ ⇢ ◆ +r⇥ ✓ f ⇢ ◆ (2.42) Si se considera el caso de flujo viscoso incompresible (densidad ⇢ constante y segunda vis- cosidad � = 0), el campo de presiones P y el campo de fuerzas constantes f no aportan a la condición de rotación, los términos asociados a esas magnitudes son cero, con lo que la anterior expresión se simplifica a la ec. (2.43) considerando la expansión de r⇥ [(u ·r)u], @! @t +r⇥ h r ⇣u · u 2 ⌘ + ! ⇥ u i = µ ⇢ r2 ! (2.43) @! @t +r⇥ (! ⇥ u) = µ ⇢ r2 ! (2.44) Utilizando identidades para desarrollar r⇥(!⇥u) y eliminando los términos nulos, se obtiene la ec. (2.45), 6Los modelos de turbulencia que se aplican en algunos métodos numéricos para flujos son una forma emṕırica de compensar, entre otras cosas, el efecto de la ausencia de las ecuaciones de conservación de cantidad de momento angular que, de todas maneras, no son suficientes para explicar el origen aparentemente estocástico de los vórtices a nivel microscópico. 7El rotacional de un campo vectorial muestra la intensidad con que el campo tiende a inducir rotación. Es un ope- rador que, aplicado a un campo vectorial arbitrarioB, se denota conr⇥B y está definido como ĺım�S!0 1 �S H Bdr. 26 @! @t + (u ·r)! = (! ·r)u+ µ ⇢ r2 ! (2.45) No debe confundirse la ecuación de vorticidad con los modelos de turbulencia. La solución de las ecuaciones de Navier-Stokes podŕıa dar como resultado un campo de velocidades determińısti- co, del que podŕıa calcularse el rotacional para obtener un indicador del grado de turbulencia del fluido. Aśı, el concepto de vorticidad se puede aplicar como un posproceso a los resultados de un modelo numérico discretizado, se haya usado o no un modelo de turbulencia. Tampoco es un sustituto para las ecuaciones de conservación de cantidad de momento angular8, que no se presentan en este trabajo dada la naturaleza irrotacional del método de SPH9. Esta sección, dedicada a presentar las ecuaciones generales de la dinámica de fluidos, da el fundamento conceptual para el pleno conocimiento de las alteraciones que implica el uso de métodos numéricos para resolverlas para determinadas condiciones y fronteras. A esto hay que agregar que como los fluidos cambian de condición según su estado energético, las ecuaciones diferenciales parciales que describen su movimiento se deben acoplar a las ecuaciones de estado de la materia que sean pertinentes. Las variantes que esto implica inciden en la exactitud e incluso en la factibilidad de los métodos numéricos que se apliquen, lo que impulsa constantemente las técnicas novedosas en métodos numéricos y el aprovechamiento de los recursos computacionales accesibles a los desarrolladores. 8El desarrollo de las ecuaciones generalizadas de la mecánica de fluidos que se presenta en este trabajo única- mente considera las cargas lineales con componentes en los ejes ortogonales (tres grados de libertad) y la reacción inercial correspondiente, pero no considera los grados de libertad rotacionales ni la inercia rotacional. Aparte del hecho de que agregarlas vuelve excesivamente complejo el sistema de ecuaciones diferenciales resultante, el método de SPH en la formulación que se usa en este trabajo por definición no incorpora su consideración. Un desarrollo conciso de las ecuaciones de conservación de momento angular para flujo rotacional se puede ver en notación vectorial en [64]. 9Sin prejuicio de lo anterior, hay modificaciones del método de SPH que incorporan entre sus variables de- pendientes un escalar de vorticidad ! que puede usarse como aproximación emṕırica del efecto disipativo de la rotación [49]. 27 Figura 2.4: Discretización por malla o euleriana Figura 2.5: Discretización por part́ıculas o lagrangiana 2.2. Método de hidrodinámica de part́ıculas suavizada (SPH) El método de hidrodinámica de part́ıculas suavizada10 (SPH) consiste en una técnica en la que el volumen de fluido se divide en una cantidad finita de part́ıculas cuya posición no es fija11, como en el esquema de la Figura 2.5. Esto contrasta con los métodos de malla12 en que los nodos de cálculo no cambian de posición, tal como se esboza en la Figura 2.4. El movimiento de estas part́ıculas de tamaño discreto y finito, como las de la Figura 2.6, se modela con las ecuaciones de dinámica de fluidos y de estado f́ısico que se elijan considerar, como cambios en densidad, viscosidad o estado energético. La interacción con otras part́ıculas no es uńıvoca; las variaciones en su estado toman en cuenta el estado de las part́ıculas cercanas o vecinas que se encuentren dentro de una zona de influencia de tamaño arbitrario, y esta ponderación se hace utilizando una función de atenuación W [89], llamada kernel de suavizado especialmente definida para ello. 10Esta es la traducción más leal al concepto contenido en el nombre original en inglés smoothed particle hydrody- namics, denominado en publicaciones y monograf́ıas por sus siglas SPH. 11En esta familia de métodos numéricos el medio lo componen part́ıculas discretas que tienen grados de libertad definidos y cuyo estado y posición es variante en el tiempo. En esta categoŕıa, conocida como los de formulación lagrangiana, se encuentran los métodos de Hidrodinámica de Part́ıculas Suavizada (SPH), los de Petrov-Galerkin (P-G), algunos métodos de elementos de contorno y los métodos de turbulencia intracelular (IP), entre otros. 12En los métodos que pertenecen a esta clasificación, la geometŕıa del medio se ha dividido en una serie de elementos discretos de localización fija para los que el sistema de ecuaciones describe el estado de cada punto como función de su posición y del tiempo. Entre estos métodos, conocidos como los de formulación euleriana, se encuentran los métodos de Diferencias Finitas (MDF), Elementos Finitos (MEF), Volúmenes Finitos (MVF), y la Simulación Numérica Directa (DNS), entre otros. 28 Figura 2.6: Modelo de part́ıculas de fluido no viscoso en el paquete Blender El propósito de esta técnica, que equivale conceptualmente a calcular una integral de volumen delimitada por la zona de interacción para determinar el valor medio de una variable activa (como la velocidad o temperatura de una part́ıcula), consiste en que el resultado no dependa solo del estado de las part́ıculas inmediatamente adyacentes, sino de lo que ocurre en toda una zona cercana a ella [89]. Esta estrategia evita ciertos tipos de inestabilidad numérica y se podŕıa decir que suaviza la trayectoria de las part́ıculas, lo que en concepto evita que en la solución haya fluctuaciones o discontinuidades de origen numérico y que el movimiento simulado de las part́ıculas también será más congruente al de un fluido real en condición de flujo subsónico. A esta función de atenuación se le conoce como el kernel13 del modelo de SPH. Una vez discretizado el volumen de fluido delimitado por su contorno (en adelante el dominio), se incorporan al modelo los efectos de interés. Esto significa que deben incorporarse las variables de campo propias de la mecánica de fluidos como densidad, enerǵıa y otras que estarán descritas con ecuaciones diferenciales parciales. Para ello se elabora una versión discretizada del efecto de estas ecuaciones diferenciales, de manera que se genere un conjunto de ecuaciones diferenciales ordinarias discretizadas, con el tiempo como única variable independiente. Esto se logra al aplicar una función de aproximación, para obtener tanto el valor aproximado de la variable de campo como el de su primera derivada [72]. El sistema de ecuaciones diferenciales ordinarias discretizadas 13Resulta especialmente confuso este nombre dado que es el mismo término que se usa en informática para ciertos componentes del sistema operativo. Para sumar a la ambigüedad, también se denomina como kernel al segmento de código fuente del que se ejecuta en múltiples copias en sistemas de procesamiento paralelo como los que se utilizan con frecuencia en implementaciones de métodos SPH. 29 resultante se puede resolver con esquemas de integración numérica convencionales [19]. La discretización del dominio se logra representándolo con una distribución arbitraria de part́ıculas que pueden entrar en contacto, pero no tienen posición fija. Posteriormente, como la función de aproximación se puede aplicar a una zona más allá de la frontera de una sola part́ıcula, si se elige el tamaño apropiado de esta zona, esta etapa tiende a estabilizar el método de cálculo al tratarse de una integral de volumen que suaviza las variaciones numéricas que haya en su interior. Este enfoque se conoce como aproximación por kernel [72], y se ejecuta numéricamente en forma sucesiva con una nueva aproximación, en que la integral y las derivadas de las variables de interés se sustituyen por aproximaciones numéricas convencionales. A este paso subsecuente se le conoce como aproximación por part́ıculas [72]. Este proceso ocurre para cada paso temporal �t en el que se haya discretizado el dominio de tiempo, de manera que se reevalúa el valor de cada una de las variables en la nueva distribución de part́ıculas, similar en naturaleza (pero sin la necesidad de algoritmos de mallado) a las mallas que se adaptan a cambios en el volumen de control. Dado que se usa un algoritmo de integra- ción expĺıcito para el conjunto de ecuaciones diferenciales ordinarias discretizadas resultante, la selección del tamaño de paso temporal no queda restringida a valores tan pequeños que resulten poco prácticos14. 2.2.1. Aproximación por kernel La formulación de esta fase de aproximación es uno de los pilares conceptuales del método de SPH. Consiste en la integral del producto de una función f(x) (que podŕıa representar una determinada variable dependiente de la posición x, como velocidad o temperatura) y una función de atenuación, denominada formalmente como función kernel de suavizado [72]. En una etapa posterior, la aproximación por part́ıculas, la integral se aproxima numéricamente utilizando una sumatoria de los valores de la variable de interés f(x) correspondientes a las part́ıculas más cercanas a aquella para la que se busca el valor aproximado. El ĺımite de integración es un subsector espacial ⌦ del dominio completo D. Para definir la función de aproximación, se parte de la identidad descrita por la ec.(2.46) [72], 14En algunos métodos el tamaño máximo de �t para garantizar estabilidad numérica es tan pequeño que un computador personal convencional duraŕıa semanas o meses en llevar a cabo los cálculos para simular casos relativamente triviales. 30 Figura 2.7: Interpretacion conceptual de la función de suavizado W con la función delta Dirac �d centrada en x = x0. Esta equivalencia es exacta siempre que f(x) sea continua y esté definida para todo el subsector ⌦. f(x) = Z ⌦ f(x0)�d( ��x� x0 ��)dx0 (2.46) Se sustituye la función �d(x) por una función arbitraria W (r, h) (denominada en adelante como función kernel en este documento) para la que r = kx� x0k y cuya entrada adicional h es la longitud de suavizado, para definir junto con el atenuador  el tamaño de la zona de influencia15 (una zona circular o esférica de radio h) de la función [72], tal como se muestra en la Figura 2.7. Como es una aproximación, en la bibliograf́ıa acerca de métodos SPH suele representarse con el operador h·i que se aplica a la función f(x) como aparece en la ec. (2.47). hf(x)i = Z ⌦ f(x0)W ( ��x� x0 �� , h)dx0 (2.47) La selección de la función kernel W (r, h) depende tanto de las propiedades matemáticas que debe tener para no introducir alteraciones al fenómeno f́ısico que se está modelando como de la influencia que podrá tener en la est