Técnicas para la mejora de la estabilidad en el Método SPH para la simulación de fluidos
Abstract
En el campo de la Informática Gráfica existen diversos problemas que siguen investigándose por no haber encontrado aún una solución suficientemente apropiada. En el caso de la Simulación de Fluidos, aunque existen herramientas y técnicas para simular el comportamiento dinámico de un fluido, los resultados obtenidos deben ser posprocesados para eliminar las perturbaciones que aparecen y que disminuyen el realismo visual. Simular un fluido, o más técnicamente el flujo de un fluido, requiere de un procesamiento en el que se resuelven el conjunto de ecuaciones en derivadas parciales con las que se describen, cuantitativamente, el flujo. En general, no puede simularse el fluido a partir de los resultados obtenidos mediante técnicas analíticas, en parte debido a que sólo es posible obtener estas soluciones bajo condiciones determinadas muy restrictivas. Además, esto supondría por una parte un esfuerzo computacional enorme y por otra parte unas simulaciones poco versátiles, ya que no sería posible introducir elementos dinámicos en el escenario de simulación. Para sobrepasar estas limitaciones es necesario utilizar técnicas discretas basadas en cálculo numérico que nos ofrecen mayor capacidad de intervención en la generación de resultados en un tiempo de cálculo plausible. Una de las técnicas más destacables es Smoothed Particles Hydrodynamic (SPH). Descriptivamente, esta técnica numérica opera sobre discretizaciones de todo el volumen del fluido mediante un conjunto ligado de puntos de control, que en esta técnica son partículas. La interacción mutua entre cada par de partículas es ponderada mediante una función continua y decreciente, conocida como función kernel. Algunas ventajas del método SPH son su adaptabilidad a contornos complejos, su capacidad para simular superficies libres, el posible tratamiento de las salpicaduras o gotas, la simplificación de las ecuaciones a resolver o la conservación de las magnitudes dinámicas que definen el fluido. A pesar de estas ventajas, el método SPH presenta algunas deficiencias, como es el caso de la incompresibilidad, que al no ser modelada correctamente induce unos resultados visuales poco realistas en los que se muestran sobreamortiguamiento de las partículas. Además, el método SPH presenta problemas de inestabilidad e imprecisión que no han sido resueltos de manera satisfactoria. Para resolver estos problemas proponemos tres mecanismos que mitigan las deficiencias enunciadas anteriormente. Estas son: 1. Formulamos una nueva función de interacción entre partículas (conocida como función kernel) que mejora la estabilidad incluso en distribuciones de partículas no homogéneas. 2. Modelamos la fuerza mutua entre partículas de modo que es posible imponer la incompresibilidad. 3. Desarrollamos un nuevo mecanismo de modelado de las condiciones de contorno, que amplía la versatilidad de SPH al permitirnos introducir en la escena de simulación contornos no suaves. Para demostrar nuestras mejoras desarrollaremos un conjunto de pruebas tanto cuantitativas, como cualitativas, en las que compararemos nuestra propuesta con otras técnicas también basadas en SPH.
Full text
Técnicas para la Mejora de la Estabilidad en el Método SPH para la Simulación de Fluidos Juan José Perea Rodríguez Director: Prof. Dr. Juan M. Cordero Tutora: Profª. Dra. Isabel Ramos Departamento de Lenguajes y Sistemas Informáticos Universidad de Sevilla Esta tesis se ha desarrollado para la obtención del grado de Doctor en Ingeniería Informática Universidad de Sevilla Octubre 2019
Agradecimientos Mi más sincero agradecimiento a mi familia, por todo lo que me ha soportado, por su comprensión, por su apoyo y su abnegación día tras día. A mi tutora Isabel Ramos por todas sus recomendaciones y apoyo. Especialmente quiero dar las gracias a mi director Juan M. Cordero por toda su dedicación, enseñanza, consejos y tanta ayuda que me ha prestado a lo largo de todos estos años. Sin duda, nada de esto se hubiese podido cumplir sin él. Gracias.
iv
Abstract In Computer Graphics scope there are several problems that are still being researched because a sufficiently suitable solution has not yet been found. In fluid flow simulation, although there are tools and techniques to simulate the fluid dynamic behaviour, the obtained results must be postprocessed to clear away the disturbances that appear and whom are responsible for reducing visual realism. The fluid simulation, or more specifically the fluid flow, require a processing to which a partial differential equation set, that quantitative describe the flow, must be solved. In general, it is not possible to simulate the fluid from the obtained results by analytic techniques, owing in partly that is only possible to obtain these solutions type under specifically conditions too restrictive. In addition, this would involve: on the one hand, a huge computational effort and not very versatile simulations on the other hand, as it would not be possible to introduce dynamic elements into the simulation environment. To overcome these limitations, it is necessary to use discrete techniques based on numerical calculation, techniques that offer us greater ability to intervene in the results generation in a plausible processing time. One of the most highlight techniques is Smoothed Particles Hydrodynamic (SPH). Descriptively, this numerical technique operates on discretization of whole fluid volume by linked set of control points, the particles. The mutual interaction between each particles pair is weighed by a continuous and decreasing function, known as kernel function. Some advantages of SPH method are its adaptability to complex boundaries, its ability to simulate free surfaces, the possible treatment of splashes or droplets, the equations simplification that must be solved or the dynamic magnitudes conservation that define the fluid. In spite of these advantages, SPH method shows some deficiencies, such as the incompressibility case which, if is not properly modelled , induces unrealistic visual results in which overdamping of the particles is shown. All this leads, SPH method also shows instability and imprecision problems that have not satisfactorily been solved. In order to solve these problems, we propose three mechanisms that mitigate the deficiencies listed above, these ones are:
vi 1. We formulate a new particle interaction function (known as kernel function) that improves stability even in inhomogeneous particle distributions. 2. We model the mutual force between particles so that it is possible to impose incompressibility. 3. We develop a new technique for modelling the boundary constraints, which extends the SPH versatility by allowing us to introduce non–soft boundaries into simulation environment. In order to show our improvements, we will develop a set of both the quantitative and the qualitative tests, where we will compare our proposal versus other techniques also based on SPH.
vii
Resumen En el campo de la Informática Gráfica existen diversos problemas que siguen investigándose por no haber encontrado aún una solución suficientemente apropiada. En el caso de la Simulación de Fluidos, aunque existen herramientas y técnicas para simular el comportamiento dinámico de un fluido, los resultados obtenidos deben ser posprocesados para eliminar las perturbaciones que aparecen y que disminuyen el realismo visual. Simular un fluido, o más técnicamente el flujo de un fl uido, requiere de un procesamiento en el que se resuelven el conjunto de ecuaciones en derivadas parciales con las que se describen, cuantitativamente, el flujo. En general, no puede simularse el fluido a partir de los resultados obtenidos mediante técnicas analíticas, en parte debido a que sólo es posible obtener estas soluciones bajo condiciones determinadas muy restrictivas. Además, esto supondría por una parte un esfuerzo computacional enorme y por otra parte unas simulaciones poco versátiles, ya que no sería posible introducir elementos dinámicos en el escenario de simulación. Para sobrepasar estas limitaciones es necesario utilizar técnicas discretas basadas en cálculo numérico que nos ofrecen mayor capacidad de intervención en la generación de resultados en un tiempo de cálculo plausible. Una de las técnicas más destacables es Smoothed Particles Hydrodynamic (SPH). Descriptivamente, esta técnica numérica opera sobre discretizaciones de todo el volumen del fluido mediante un conjunto ligado de puntos de control que en esta técnica son partículas. La interacción mutua entre cada par de partículas es ponderada mediante una función continua y decreciente, conocida como función kernel. Algunas ventajas del método SPH son su adaptabilidad a contornos complejos, su capacidad para simular superficies libres, el posible tratamiento de las salpicaduras o gotas, la simplificación de las ecuaciones a resolver o la conservación de las magnitudes dinámicas que definen el fluido. A pesar de estas ventajas, el método SPH presenta algunas deficiencias, como es el caso de la incompresibilidad, que al no ser modelada correctamente induce unos resultados visuales poco realistas en los que se muestran sobreamortiguamiento de las partículas. Además, el método SPH presenta problemas de inestabilidad e imprecisión que no han sido resueltos de manera satisfactoria.
ix Para resolver estos problemas proponemos tres mecanismos que mitigan las deficiencias enunciadas anteriormente. Estas son: 1. Formulamos una nueva función de interacción entre partículas (conocida como función kernel) que mejora la estabilidad incluso en distribuciones de partículas no homogéneas. 2. Modelamos la fuerza mutua entre partículas de modo que es posible imponer la incompresibilidad. 3. Desarrollamos un nuevo mecanismo de modelado de las condiciones de contorno, que amplía la versatilidad de SPH al permitirnos introducir en la escena de simulación contornos no suaves. Para demostrar nuestras mejoras desarrollaremos un conjunto de pruebas tanto cuantitativas, como cualitativas, en las que compararemos nuestra propuesta con otras técnicas basadas también en SPH.
Lista de Tablas 6.1 Valores del parámetro de orden, ς , asociado a la prueba de sedimentación de una distribución aleatoria de partículas. . . . . . . . . . . . . . . . . . . . 76 6.2 Valores relevantes obtenidos de las gráficas de la densidad para la prueba Sod’s shock–tube, figura 6.2a. . . . . . . . . . . . . . . . . . . . . . . . . 79 6.3 Valores relevantes obtenidos de las gráficas de la energía térmica para la prueba Sod’s shock–tube, figura 6.2b. . . . . . . . . . . . . . . . . . . . . 79 6.4 Valores relevantes obtenidos de las gráficas de la velocidad para la prueba Sod’s shock–tube, figura 6.2c. . . . . . . . . . . . . . . . . . . . . . . . . 80 6.5 Valores relevantes obtenidos de las gráficas de la densidad para la prueba de onda expansiva, figura 6.3a. . . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.6 Valores relevantes obtenidos de las gráficas de la energía térmica para la prueba de onda expansiva, figura 6.3b. . . . . . . . . . . . . . . . . . . . . 82 6.7 Valores relevantes obtenidos de las gráficas de la velocidad para la prueba de onda expansiva, figura 6.3c. . . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.8 Valores relevantes obtenidos de las gráficas de la densidad para la prueba de impacto sobre muro, figura 6.4a. . . . . . . . . . . . . . . . . . . . . . . . 85 6.9 Valores relevantes obtenidos de las gráficas de la energía térmica para la prueba de impacto sobre muro, figura 6.4b. . . . . . . . . . . . . . . . . . 85 6.10 Valores relevantes obtenidos de las gráficas del momento para la prueba de impacto sobre muro, figura 6.4c. . . . . . . . . . . . . . . . . . . . . . . . 85 6.11 Valores relevantes obtenidos de las gráficas de la densidad para el test Sjögreen,figura6.5a................................ 87 6.12 Valores relevantes obtenidos de las gráficas de la energía térmica para el test Sjögreen,figura6.5b.............................. 88 6.13 Valores relevantes obtenidos de las gráficas de la velocidad para el test Sjögreen,figura6.5c.............................. 88
Lista de Tablas xvii 6.14 Rangos de valores de estabilidad. Las magnitudes se expresan en unidades del sistema internacional. . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Capítulo 1 Introducción 1.1 Introducción La simulación de fluidos juega un papel fundamental en el ámbito de la ingeniería informática. Por ello, a lo largo de las últimas décadas se han desarrollado, y mejorado, distintas técnicas computacionales con las que obtener simulaciones realistas, esto es, simulaciones estables y precisas acordes con el comportamiento esperado del fluido. En términos cuantitativos un fluido se conforma como un medio continuo en el que se definen variables de estado, como son la densidad de masa, la presión, la velocidad, energía térmica, etc. ypropiedades constitutivas, tales como la viscosidad y la incompresibilidad entre otras. Por otro lado, la dinámica del flujo del fluido, se describe mediante un sistema de ecuaciones en derivadas parciales, donde se establecen relaciones entre las variables de estado y las propiedades constitutivas. A través de estas ecuaciones, basadas en los principios de conservación, se cuantifican y relacionan la variación temporal y espacial de las variables de estado, como se describe en [ 91 ]. En términos operativos, simular el flujo de un fluido consiste en conocer el valor de las variables de estado en todos los puntos del fluido y en cualquier instante de tiempo. Esto se consigue resolviendo las ecuaciones en derivadas parciales que describen el flujo. Sin embargo, sólo es posible obtener soluciones analíticas, de dichas variables de estado, para problemas muy concretos con alta simetría espacial. Para casos más generales, es necesario el uso de métodos numéricos [91]. En síntesis, estos métodos numéricos operan sobre discretizaciones del fluido continuo. Estas discretizaciones consisten en un conjunto ligado de puntos de muestreo distribuidos a lo largo de todo el volumen del fluido. La ventaja operativa de los métodos numéricos se basa en el hecho de que las variables de estado se particularizan sobre los puntos de muestreo. De este modo, las ecuaciones en derivadas parciales que describen el flujo, se transforman en un conjunto de ecuaciones algebraicas resolubles mediante procesamiento computacional.
1.1 Introducción 2 Dependiendo de los objetivos que se establezcan para la simulación, se puede prescindir de la implementación de algunas de las ecuaciones en derivadas parciales que describen el flujo del fluido. En el caso particular de la informática gráfica (que designaremos por sus siglas en inglés CG), donde se busca un equilibrio entre el realismo y la eficiencia computacional, como se describe en [ 35 , 62 , 80 ], suele ser habitual resolver sólo la ecuación de Navier–Stokes (NS). Para simular otros fenómenos que aparecen en el fluido, tales como la incompresibilidad o el balance de energía térmica, es necesario el modelado de restricciones anexas sobre la ecuación de NS, tal como se propone en [ 3 , 8 , 17 , 68 ], ecuación que en su formulación más compacta se define a través de la expresión 1.1: ρDv Dt =−∇P+µ∇2v+ρg,(1.1) donde ρ es la densidad del fluido, D/Dt hace referencia a la derivada sustancial, v es la velocidad del fluido, P es la presión del fluido, g es la gravedad y µ es la viscosidad del fluido. En el ámbito de la simulación de fluidos, el término −∇P suele ser designado como fuerza de presión y representa la tasa de cambio de presión que aparecen entre regiones contiguas del fluido; µ∇2v se denomina fuerza de viscosidad que es de naturaleza disipativa y se relaciona con las propiedades constitutivas del fluido a través de la viscosidad µ ; y el término ρges el peso. Existen diferentes métodos con los que resolver numéricamente la ecuación 1.1 y simular el fluido, si bien suelen agruparse en dos categorías: euleriana ylangrangiana. En una descripción breve de sus características podemos afirmar que los métodos eulerianos operan sobre discretizaciones del fluido a través de mallas conformadas mediante triángulos o cuadriláteros. En cada uno de estos polígonos se ubican, generalmente en sus vértices y aristas, un conjunto de puntos de muestreo, sobre los que se particularizan las variables de estado [ 47 , 91 ]. Por otro lado, en los métodos lagrangianos, el fluido se modela mediante un conjunto de partículas sobre las que se definen las variables de estado. La interacción mutua entre estas partículas, es ponderada a través de una función que depende de la distancia [13, 71]. En la figura 1.1 se ilustran las dos metodologías que acabamos de describir. Simular un fluido utilizando métodos eulerianos, transforma las ecuaciones en derivadas parciales, que describen el flujo del fluido, en un sistema de ecuaciones algebraicas. Esta formulación presenta algunas ventajas operativas como es el hecho de que permite acoplar directamente presión y velocidad, con lo que es relativamente sencillo imponer la rectricción de incompresibilidad [ 91 ]. Además, fenómenos como la advección o la ruptura de flujo laminar son fácilmente modelados y simulados. A pesar de estas ventajas, los métodos eulerianos tienen algunos inconvenientes, los más destacables son: (1) la distorsión en los resultados obtenidos, (2) la reestructuración reiterada de la malla de polígonos para
1.1 Introducción 3 Figura 1.1 Representación de los modelados del fluido para la implementación de técnicas numéricas. La imagen superior representa el fluido a simular, mientras que en las inferiores se muestra del dominio discretizado. La imagen inferior izquierda muestra la discretización asociadas a la metodología lagrangiana, en la inferior derecha se muestra la euleriana.
1.1 Introducción 4 poder simular el flujo, (3) la complejidad del sistema de ecuaciones algebraicas modelado, según se describe en [ 91 ]. En referencia (1), la distorsión se relaciona con la alta tasa de cambio que puede aparecer en las variables de estado, especialmente en la velocidad y la presión. La consecuencia de la distorsión, es que los valores de estas magnitudes acaban siendo multivaluados, lo que induce alta inestabilidad en resultados [ 91 ]. Para resolver este problema se requiere reestructurar la malla de polígonos y reformular las magnitudes dinámicas. Respecto a la limitación (2), la malla requiere ser reestructurada por dos causas, la primera es para evitar la distorsión, que hemos descrito previamente, la segunda es para que se adapte a los cambios del contorno, especialmente en la superficie libre, cuando exista interfase en el fluido. En la reestructuración se induce una pérdida de eficiencia computacional, ya que cada vez que es necesario llevarla a cabo, el cálculo de las variables de estado se paraliza. Así pues, la eficiencia se reduce drásticamente. A todo ello hay que añadir la limitación (3), esto es, la complejidad del sistema de ecuaciones algebraicas ya que del proceso de reestructuración induce el incremento del número de ecuaciones e incógnitas en el sistema de ecuaciones, lo que dificulta su resolución. Es por ello por lo que, aunque estos inconvenientes no perjudican los resultados de la simulación, no son idóneos para uso en el ámbito de CG, especialmente por su excesivo coste computacional, tal como se advierte en [63, 80]. Por otro lado, utilizar métodos lagrangianos para modelar la dinámica del fluido, transforma las ecuaciones en derivadas parciales del fluido en un sistema de ecuaciones algebraicas, en el que sus ecuaciones habitualmente están desacopladas, como se indica en [ 4 , 71 ], por lo cual, en comparación con los métodos eulerianos, los lagrangianos son mucho más eficientes. Además, otras ventajas que presentan estos métodos son: (1) su adaptabilidad y versatilidad para simular entornos complejos, según [ 62 ], (2) la conservación de masa, según se justifica en [ 54 , 62 , 71 ], (3) la simplificación de la formulación para resolver las ecuaciones dinámicas del fluido [ 45 , 54 ], (4) la obtención directa de la superficie libre, ya que al estar modelado mediante partículas, la interfase se simula directamente, según la descripción realizada en [ 54 , 62 ]. Por estas características, que representan una ventaja operativa frente a los métodos eulerianos, los métodos lagrangianos han sido intensivamente utilizados en el ámbito de CG, siendo el método Smoothed Particle Hydrodynamics (SPH) uno de los más implementados. Sin embargo, a pesar de las ventajas previamente descritas, es habitual que en las simulaciones mediante SPH aparezcan comportamientos erráticos, tales como agrupaciones irreales de partículas, dispersiones espurias en los contornos o sobreamortiguamiento. La consecuencia son simulaciones en las que se muestra un comportamiento errático del fluido, que no es acorde con su comportamiento real.
1.2 Limitaciones del Método SPH 5 El origen de este comportamiento errático suele relacionarse con problemas de inestabilidad e imprecisión del método SPH, como se destaca en [21, 74, 90]. Con el objetivo de resolver los problemas de inestabilidad e imprecisión que aparecen en las simulaciones de fluidos mediante SPH, y obtener así simulaciones realistas, hemos desarrollado un conjunto de mecanismos que serán descritos a lo largo de esta tesis. Comenzaremos analizando las causas que generan esos problemas y justificaremos cómo nuestras propuestas las resuelven. En la siguiente sección vamos a describir el origen de esos problemas de inestabilidad e imprecisión y las consecuencias que tienen en la simulación. 1.2 Limitaciones del Método SPH Como ya hemos mencionado con anterioridad, existe un consenso en relacionar las simulaciones erráticas de SPH, con los problemas de imprecisión e inestabilidad del método. En esta sección vamos a describir, brevemente, estos problemas inestabilidad e imprecisión y su repercusión en la simulación. Nos centraremos en las diversas metodologías propuestas para resolverlos, así como las limitaciones que dichas propuestas presentan. El objetivo es hacer una introducción sobre los problemas que el método SPH presenta y cómo su resolución es un campo de investigación aún abierto. Diversos estudios [ 36 , 57 , 71 , 74 , 90 ] concluyen que el método SPH presenta problemas de inestabilidad e imprecisión que acaba reduciendo el realismo de la simulación. Dichos problemas suelen relacionarse con: 1. El número y distribución de partículas vecinas. En las simulaciones mediante SPH suelen aparecer dos problemas que, a priori, podrían parecer contradictorios, estos son: (1) dispersiones irreales de partículas, (2) agrupamiento erróneo de partículas. El origen suele relacionarse con el número y la distribución de las partículas con las que cada una de ellas interactúa (las partículas vecinas), según los estudios de [ 71 , 74 ]. Para resolver estos problemas se han desarrollado técnicas centradas en incrementar el número de partículas vecinas, según la recomendación dada por [ 21 , 74 ]. Sin embargo, existen regiones en el fluido donde este incremento no puede llevarse a cabo, como ocurre en las interfases y en las regiones de interacción con el contorno, como se destaca en [ 80 ]. Además, este incremento no resuelve los problemas relacionados con el agrupamiento erróneo de las partículas, según el estudio desarrollado por [ 14 , 21 ]. Es por ello, por lo que se han desarrollado técnicas que afrontan el problema desde un punto de vista dinámico y no de distribución de partículas. En síntesis, estas técnicas resuelven la agrupación irreal o la dispersión introduciendo fuerzas de control
1.2 Limitaciones del Método SPH 6 ficticias que inducen una repulsión artificiosa de las partículas, cuando aparece la agrupación irreal, o amortigüen en el caso de que se manifieste la dispersión errónea. Con estas metodologías se puede paliar, en parte, los problemas de agrupamiento y dispersión. Sin embargo, estas fuerzas ficticias suelen acabar desestabilizando el fluido, ya que habitualmente violan la conservación del momento, lo que favorece el sobreamortiguamiento, tal como se indica en [21]. Así pues, es necesario desarrollar una técnica que sea inmune a la distribución y el número de partículas vecinas. 2. La ponderación de la dinámica entre las partículas. Desde el punto de vista de la interacción, en SPH, el pilar fundamental es que la interacción entre partículas decae continuamente con la distancia. Dicha interacción es modelada mediante una función continua y monótona decreciente con la distancia, (la función kernel). Diversos estudios se han centrado en formular funciones kernel que garanticen estabilidad y precisión, como los llevados a cabo por [ 21 , 88 ]. En estos estudios se destaca la relación entre la función kernel, el agrupamiento irreal de las partículas y el sobreamortiguamiento, por lo que se proponen distintas funciones para solventar estos problemas. Los resultados que se obtienen son satisfactorios en la mayoría de los casos. Sin embargo, requiere que el número de partículas vecinas sea relativamente alto, según se indica en [ 21 ], restricción que no se puede satisfacer ni en las interfases, ni en el contorno, tal como se advierte en [ 74 , 80 ]. Es por ello que se hace necesario, formular nuevas funciones kernel que garanticen estabilidad aun cuando el número de partículas vecinas sea bajo. 3. La imposición de la incompresibilidad. Una de las restricciones más destacables que tienen los líquidos es la incompresibilidad. Cuantitativamente, la incompresibilidad se describe mediante una ecuación de continuidad, donde dos varaibles de estado como son la velocidad y la tasa de cambio de la densidad, no pueden variar libremente [ 9 ], es decir, que la densidad de masa se mantiene constante en todo el flujo del fluido. Cualitativamente, no imponer la incompresibilidad da como resultado simulaciones poco realistas, donde el efecto más destacable es el sobreamortiguamiento, según [ 57 ]. Para evitar este problema se han desarrollado diversas estrategias que dependiendo del proceso de resolución se suelen clasificar en implícitas yexplícitas. En las técnicas implícitas, se suele modelar la condición de incompresibilidad que se acopla con la de ecuación de NS, como se describe en [ 36 ], el resultado es un sistema de ecuaciones que satisface la condición de incompresibilidad, a partir del cual, se calcula el campo de presiones o de densidades. Por otro lado, en los métodos explícitos, o bien se establece una relación entre la densidad y la presión, mediante la ecuación de estado
1.2 Limitaciones del Método SPH 7 [ 9 ], y se fijan fuerzas que impidan el sobreamortiguamiento [ 57 ]; o bien se relaciona la densidad y la presión, mediante un método predictor–corrector, de modo que se satisfaga la incompresibilidad, según [ 80 ]. Si bien ambas estrategias, implícitas y explícitas, pueden dar resultados satisfactorios desde el punto de vista del realismo de la simulación, cada una de ellas presentan limitaciones específicas. Por ejemplo, en las técnicas implícitas, la resolución del sistema de ecuaciones formado por la ecuación de continuidad y la de NS, requiere de una simplificación de la ecuación de continuidad, lo que puede reducir el realismo, tal como se indica en [ 8 , 36 ]. Además, resolver un sistema de ecuaciones en derivadas parciales incrementa el coste computacional. Por otro lado, las técnicas explícitas, aunque son eficientes desde el punto de vista del coste computacional, no consiguen controlar la incompresibilidad por debajo de un 2% de cambio en la densidad, según se advierte en el estudio de [ 80 ]. Además, con estas técnicas explícitas, no se solventa totalmente los problemas de estabilidad en las regiones del contorno del fluido, donde el número de partículas vecinas no puede controlarse. Así pues, se hace necesario desarrollar una técnica que aproveche la eficiencia de los métodos explícitos y con la que se obtengan resultados estables en todas las áreas del fluido, especialmente en las regiones límites del contorno. 4. El balance de energía térmica en el contorno. En el flujo del fluido se producen balances de energía térmica que modifican su estado interno [ 74 ]. Cuantitativamente, el cambio de energía térmica se describe mediante la ecuación diferencial del calor [ 9 , 74 ], donde los cambios de energía térmica se compensan con variaciones de la presión y la densidad [ 74 , 91 ]. La complejidad de resolver esta ecuación junto con la de NS lo hace poco atractivo para su uso en el ámbito de CG, por lo que el balance de energía térmica en el fluido no suele ser simulado. Sin embargo, es necesario compensar los cambios de energía térmica para obtener simulaciones realistas. El modo habitual de hacerlo es implementar modelos viscosos en los que la viscosidad amortigua el exceso de energía térmica. Sin embargo, estos modelos suelen inducir sobreamortiguamiento [ 57 ]. Además, existe un problema entre la incompresibilidad y el balance de energía térmica en el contorno. Sucintamente, lo que ocurre es que la incompresibilidad favorece un trasvase de energía térmica en el contorno del fluido. Dicho balance de energía térmica, de no ser disipado, genera dispersiones erróneas de las partículas, lo que afecta a la estabilidad y precisión de la simulación, como se destaca en [ 74 ]. La manera habitual de tratar estas dispersiones de partículas, es introduciendo un coeficiente de disipación, que actúa sobre la velocidad de las partículas del fluido que interactúan con el contorno. El problema es que el balance de energía térmica en el contorno, no sólo se debe al cambio de velocidad, sino más
1.2 Limitaciones del Método SPH 8 bien a los cambios de presión al “impactar” el fluido sobre el contorno. Es por ello que se hace necesario desarrollar una técnica que relacione la incompresibilidad y el balance de energía térmica, especialmente en el contorno. 5. La interacción fluido contorno. En la simulación del fluido tiene especial interés la interacción con el contorno, según [ 64 ]. Dos son los principales retos a la hora de simular la interacción fluido–contorno: el primero es la detección y el segundo es el modelado de las condiciones de contorno. Para la detección se han desarrollado diversas técnicas, geométricas [ 18 ], discretas [ 51 ] o analíticas [ 90 ]. El principal objetivo de la detección es determinar las partículas del fluido que contactan con el contorno. Con el modelado de las condiciones de contorno, el objetivo es simular los diferentes efectos, que en las variables de estado, tiene el contorno. Estos efectos suelen estar descritos mediante ecuaciones en derivadas parciales. El problema que plantea el modelado de las condiciones del contorno es que, los usados habitualmente en CG, no pueden modelar condiciones de contorno complejas, como se advierte en [ 84 ], lo cual limita el espectro de simulaciones. Para poder simular condiciones de contorno complejas, aunque existen estudios que afrontan este problema, las soluciones que proponen son poco adecuadas para su uso en CG, por la complejidad numérica que presentan. Complejidad que induce un incremento de coste computacional que no compensa los efectos simulados.Por ello, es necesario desarrollar una técnica que permita una correcta gestión tanto de la detección, como de la imposición de condiciones de contorno complejas. Técnica que debe tener dos objetivos: el primero es conseguir estabilidad y precisión y el segundo es que sea útil para su uso en el ámbito de CG. 6. Los errores en las magnitudes dinámicas en el contorno. Uno de los problemas recurrentes en SPH se relaciona con el bajo número de partículas vecinas. Cuando el número de partículas vecinas es bajo, magnitudes como la velocidad, densidad y presión acumulan errores numéricos, lo que afecta a la estabilidad de la simulación, que se deduce de los estudios [ 8 , 80 ]. Este problema se ve agudizado en el contorno, donde además de los problemas relacionados con el número de partículas vecinas se le añade el modelado de magnitudes como la densidad, presión y velocidad, al interactuar con el contorno [ 84 ]. Así pues, se debe desarrollar una metodología acorde con la imposición de las condiciones de contorno, y que permita calcular magnitudes como la densidad o la presión, aunque no exista una condición de contorno específica para estas magnitudes.
2.2 Características del Método SPH 15 Figura 2.1 Representación esquemática de los elementos relevantes en la formulación SPH. W(r−r′,h) representa a la función kernel, h hace referencia al radio del dominio soportado. Las partículas cuya distancia cumple r−r′<hson partículas vecinas. simular su flujo. Así mismo, esbozaremos algunas de las limitaciones que este método presenta. Para resolver estas limitaciones, hemos desarrollado un conjunto de mecanismos que serán descritas en los capítulos 3, 4 y 5. 2.2 Características del Método SPH Cualitativamente, el hecho de que el método SPH se base en el modelado del fluido mediante partículas, implica que los desplazamientos de éstas representen al flujo del fluido. Por otro lado, desde un punto de vista cuantitativo, sobre cada una de las partículas se particularizan las propiedades constitutivas y las variables de estado, tales como: la presión, la densidad de masa, la energía térmica, la velocidad, la posición, etc. Así pues, las partículas juegan el papel de los puntos de muestreo que “fluyen” con el movimiento del fluido. Esta es la causa por la que el método SPH se basa en la formulación lagrangiana de la ecuación de NS, según se describe en [55, 70], ecuación 1.1. Pasar de la formulación continua de las variables de estado, a un modelado discreto con el que opera SPH, implica una pérdida de información dinámica ya que, a priori, no se modela la interacción definida a través de la ecuación de NS. Para evitar este problema, en el método SPH, se formula una función monótona decreciente con la distancia, denominada función kernel, que constituye un pilar fundamental de este método. Analíticamente, la función kernel debe satisfacer un conjunto de características, las cuales serán descritas en profundidad en el capítulo 3, entre la que podemos destacar que debe tener un dominio de definición limitado. Esta condición implica restringir el rango de interacción entre las
2.2 Características del Método SPH 16 partículas a una distancia máxima, por encima de la cual, dos partículas no interactúan. Esta distancia máxima se fija en la función kernel mediante el radio del dominio soportado, que ya ha sido introducido en la sección 2.1. Así pues, a través del radio del dominio soportado se controlan las partículas que interactúan mutuamente, esto es, con el radio del dominio soportado se establece la vecindad de cada partícula. Además de las características previamente descritas, en el método SPH existen otras que han sido destacadas en diversos estudios [ 7 , 21 , 54 , 55 , 71 ] por representar una ventaja tanto cuantitativa, como cualitativa en la simulación de fluidos. Estas características relevantes son: 1. Satisface directamente la conservación de la masa, lo que representa una ventaja operativa y simplifica los cálculos. Ello se debe a que el método SPH opera sobre discretizaciones del fluido mediante partículas, en las que la masa de cada una de ellas, representa una fracción discreta del volumen del fluido. Por lo tanto, la suma de las masas de todas las partículas, conformará la masa total del volumen del fluido. 2. Por su formulación lagrangiana, el método satisface la invariancia galileana. Debido a que los desplazamientos de las partículas representan el flujo del fluido, la formulación con la que opera el método es la formulación lagrangiana. Desde un punto de vista dinámico, este hecho implica la invariancia galileana de las magnitudes vectoriales (aceleración, velocidad, posición, gradiente de presión...) lo que simplifica las operaciones y el cálculo de dichas magnitudes vectoriales. 3. El modelado de las ecuaciones dinámicas mediante SPH es conceptualmente sencillo y transforma complejas ecuaciones en derivadas parciales en ecuaciones algebraicas independientes. De este modo, es relativamente fácil su implementación y su resolución. Esta ventaja es especialmente destacable cuando se simulan entornos tridimensionales. 4. Modela directamente flujos con superficie libre e interfases. Al operar sobre sistemas de partículas, la formación de la superficie libre o de las interfases se obtiene de manera directa. 5. Permite calcular la advección y los fenómenos de transporte en el fluido. A pesar de que el término de advección del fluido no aparece explícitamente en la formulación lagrangiana, el método permite “reconstruirlos” a partir de la tasa de cambio de la velocidad y de la densidad de masa durante el estado transitorio del fluido. 6. Simplifica la simulación de fenómenos tanto de muy alta como de muy baja energía térmica, como mostraremos en el capítulo 6. Si bien, existen zonas en estas simulaciones donde los resultados obtenidos son imprecisos.
2.3 Formulación del Método SPH 17 A pesar de estas ventajas, el método SPH presenta limitaciones de estabilidad y precisión. Limitaciones que fueron esbozadas en el capítulo 1 y que serán descritas, con mayor profundidad, en los capítulos 3, 4 y 5, donde explicaremos nuestras propuestas para superarlas. 2.3 Formulación del Método SPH La ecuación fundamental de SPH, se basa en la formulación integral de la función delta de Dirac [ 54 ] y queda abalada por la teoría de interpolación de funciones [ 7 ]. Operativamente, la obtención de la formulación de SPH, puede describirse en dos pasos: 1. Se generaliza la formulación integral de la función δ de Dirac, para “aproximar” cualquier magnitud en un dominio continuo. En este proceso, la función δ es sustituida por una función con extensión espacial, a esta función se la conoce cómo función kernel. Con este primer paso se obtiene la formulación integral de SPH. 2. Se modela, o sustituye, la formulación integral mediante una sumatoria que se extiende a todo el dominio soportado de la función kernel. Esta nueva formulación discretizada se concreta sobre cada partícula, de modo que el valor final de cualquier magnitud, se obtiene como resultado de la suma ponderada del valor que, dicha magnitud, tiene en cada partícula vecina. Con esta operación se obtiene la formulación discreta de SPH. En las secciones 2.3.1 y 2.3.2 se describirán más detalladamente cada uno de los dos puntos anteriores. Así mismo, se obtendrán las expresiones con las que se modelan las derivadas de las funciones y se destacará el problema de conservación del momento que presentan las ecuaciones obtenidas. 2.3.1 Formulación integral de SPH Según la formulación integral de la función delta Dirac, cualquier función f puede aproximarse cómo: f(r) = Zf(r′)δ(r−r′)dr′,(2.1) donde r hace referencia a la posición, r′ es la posición donde está definida f(r′) , f representa a cualquier función escalar, que en el caso del fluido será cualquier variable de estado, también designadas como magnitudes dinámicas, (velocidad, presión, densidad...), y δ hace referencia a la función delta Dirac.
2.3 Formulación del Método SPH 18 En base a la teoría de interpolación [ 7 , 54 ], la función δ puede ser modelada mediante cualquier función continua con extensión espacial. En base a ello, la ecuación 2.1 puede ser reformulada mediante la expresión: ⟨f(r)⟩=ZΩ f(r′)W(r−r′,h)dr′(2.2) donde ⟨⟩ indica que el valor obtenido es una aproximación de f en la posición r , W es la función kernel, Ω hace referencia al volumen del dominio de definición de W y h es el radio del dominio soportado, a partir del cual, se define la distancia máxima de interacción entre partículas. Sucintamente, desde un punto de vista analítico, la función kernel debe cumplir, al menos, las siguientes características: 1. Debe estar normalizada y tender a una pulsiana, δ , cuando h tiende a cero, según [ 54 ], esto es: lim h→0W(r−r′,h) = δ(r−r′)yZΩ Wdr′=1.(2.3) 2. Debe ser defina como una función positiva, par y con un decrecimiento continuo con la distancia (r−r′) . Además, tanto la función kernel como sus derivadas deben presentar gráficas con perfiles suaves, como se indica en [54, 71]. 3. Debe ser definida como función simétrica, y compacta, respecto a (r−r′) , según la conclusión de [45, 54] Pasar de la formulación integral de la función delta de Dirac, ecuación 2.1, a la formulación basada en la teoría de interpolación, ecuación 2.2; esto es, sustituir la función δ por la función W , genera errores de convergencia que afectan a la estabilidad y a la precisión de los resultados. Diversos estudios se centran en analizar y cuantificar dicho error de convergencia [ 21 , 71 , 88 ]. El modo habitual de demostrar el orden de dicho error, es mediante el desarrollo en serie de Taylor de la ecuación de interpolación 2.2. Para concretar el orden mínimo de error, se va a calcular el desarrollo en serie de Taylor de la ecuación 2.2 para el caso monodimensional. Este desarrollo no implica una pérdida de generalidad, ya que el mismo proceso puede extenderse para los casos bidimensional y tridimensional. Así pues, la aproximación de una magnitud arbitraria f(x) , en el ámbito de formulación de integral de SPH, se puede escribir como: ⟨f(x)⟩=ZΩ f(x′)W(x−x′,h)dx′=f(x′)ZΩ W(x−x′,h)dx′+f′(x′)ZΩ(x−x′)W(x−x′,h)dx′
2.3 Formulación del Método SPH 19 +1 2f′′(x′)ZΩ(x−x′)2W(x−x′,h)dx′+O(h3)(2.4) donde Ohace referencia al resto del desarrollo en serie de Taylor. Analizando la expresión 2.4, ésta se puede simplificar teniendo en cuenta dos de las tres características, previamente descritas, que la función kernel tiene que satisfacer: la primera es la condición de normalización, ecuación 2.3, con lo que el valor de la integral del primer sumando de la ecuación 2.4 vale 1 ; la segunda es la condición de paridad que se le exige a la función kernel, es decir RxW(x,h)dx =0 , con lo que el segundo sumando de la ecuación 2.4 se anula. Teniendo en cuenta estas dos características, se puede concluir que: ∆f(x)∼Z(x−x′)2W(x−x′,h)dx′(2.5) donde ∆ hace referencia al error que se comete al calcular la función f(x) y con ∼ indicamos el orden de dicho error. Generalizando la ecuación 2.5 para el caso tridimensional, se obtiene la expresión 2.6. ∆f(r)∼Zr2W(r,h)dr.(2.6) Es habitual relacionar la ecuación 2.6 con el momento de segundo orden de la función kernel, lo que permite cuantificar el error cometido al implementar una determina función kernel. Entonces, en base al resultado obtenido en la expresión 2.6, se puede concluir que la formulación fundamental de SPH garantiza, en el caso óptimo, una precisión de hasta el segundo orden. Si bien, cabe destacar que este orden de error se da en las mejores condiciones de simulación. Sin embargo, en situaciones como la interacción con el contorno, las interfases o las regiones donde aparece discontinuidad para la densidad de masa, la energía térmica o la presión, el error es muy superior y requiere seleccionar la mejor función kernel que ofrezca el mínimo error (como describiremos en el capítulo 3). Así mismo, a partir de la ecuación 2.6, es posible cuantificar el error esperado al utilizar una determinada función kernel. Para funciones kernel como la Gaussiana o el Quartic B–spline, muy utilizados en SPH, el error mínimo que se puede cometer es 0.6h2 y 3h2 respectivamente. Especial mención merece el caso de funciones kernel formuladas mediante la composición de funciones coseno, en este caso el rango de error oscila entre 0.5h2−0.65h2 [ 88 ]. Así pues, se puede afirmar que una función kernel basada en la composición de funciones coseno, es una buena opción. Esta va a ser una de nuestras premisas para obtener la función kernel propuesta, tal como describiremos en el capítulo 3.
2.3 Formulación del Método SPH 20 Formulación integral de las derivadas en SPH Una vez formulada la ecuación general con la que es posible aproximar cualquier magnitud mediante el método SPH, ecuación 2.2, el siguiente paso es calcular las expresiones con las que aproximar las derivadas de las funciones. Esto se consigue geneneralizando la ecuación 2.2 para ∇·f. La expresión que se obtiene se expresa en términos de la ecuación 2.7: ⟨∇·f(r)⟩=ZΩ∇·f(r′)W(r−r′,h)dr′,(2.7) donde la integral de [∇·f(r′)] se calcula sobre la posición r′. Aplicando la regla de la cadena sobre el primer multiplicador del argumento de la integral 2.7, se obtiene la ecuación 2.8: ∇·f(r′)W(r−r′,h) = ∇·f(r′)W(r−r′,h)−f(r′)·∇W(r−r′,h),(2.8) Sustituyendo la ecuación 2.8, en la ecuación 2.7, es posible reformularla en términos de ecuación 2.9: ⟨∇·f(r)⟩=ZΩ ∇·f(r′)W(r−r′,h)dr′−ZΩ f(r′)·∇W(r−r′,h)dr′.(2.9) Para simplificar la ecuación 2.9, se puede utilizar el Teorema de la Divergencia sobre la primera integral, con lo que se obtiene: ⟨∇·f(r)⟩=ZS f(r′)W(r−r′,h)·ndS −ZΩ f(r′)·∇W(r−r′,h)dr′.(2.10) dondenhace referencia al vector normal a la superficie que recubre el dominio Ω. La importancia de la ecuación 2.10, radica en el hecho de que su primer término, esto es, el término de la integral de superficie, es nulo siempre que se cumpla la condición de compacidad de la función kernel (tercera característica que se definió en la 2.3.1). En el caso de que se cumpla dicha condición, la ecuación 2.10 se simplifica cómo: ⟨∇·f(r)⟩=−ZΩ f(r′)·∇W(r−r′,h)dr′.(2.11) Sin embargo, existen situaciones en las que esa simetría no se satisface, como sucede cuando la partícula se aproxima al contorno, según se destaca en [ 7 , 54 ]. En esta situación,
2.3 Formulación del Método SPH 21 el contorno es el responsable de romper dicha simetría. En estas circunstancias, la ecuación 2.11 no es una buena aproximación del valor de ⟨∇·f(r)⟩. Analizando la ecuación 2.11 se puede afirmar que, en el método SPH, el gradiente de cualquier magnitud f(r) se puede obtener como el producto de dicha magnitud, por el gradiente de la función kernel extendida sobre el dominio de interacción (el dominio soportado). Siempre que se satisfaga la condición de compacidad, a través de la cual, la integral de superficie de la ecuación 2.10 se anula. Esta simplificación de la derivada, mostrada en la ecuación 2.11, representa una ventaja operativa del método SPH frente a otras técnicas numéricas, como es destacado en [54]. Generalizando las operaciones desarrolladas entre las ecuaciones 2.7–2.11, se pueden formular las derivadas de orden superior. Así pues, la ecuación general con la que se calculan las derivadas de cualquier orden en SPH, se puede expresar en términos de la ecuación 2.12: ⟨∇lf(r)⟩=−ZΩ f(r′)·∇lW(r−r′,h)dr′.(2.12) donde lhace referencia la orden diferencial. Sin embargo, aunque la ecuación 2.12, establece una formulación general para calcular las derivadas de cualquier función mediante SPH, dicha formulación no es adecuada cuando se utiliza en el cálculo de magnitudes dinámicas y suele generar inestabilidades, como se destaca en [ 54 ]. Dichas inestabilidades inducen un comportamiento errático de las partículas, como se advierte en [ 7 ], lo que afecta al realismo de la simulación. El origen de este problema, en el caso de la simulación de fluidos, está relacionado con el hecho de que la ecuación 2.12, viola la conservación del momento, según el estudio de [ 54 ]. Para solventar este problema, se ha desarrollado una metodología que, basándose en la formulación de SPH, permite obtener expresiones simétricas. De este modo, se podrá satisfacer la conservación del momento, ya que las magnitudes dinámicas aparecen por pares [ 54 ]. El proceso seguido para obtener dichas ecuaciones simétricas, así como sus fundamentos, se describirá en la sección 2.3.2, todo ello según la investigación desarrollada en [54]. 2.3.2 Formulación discreta de SPH Cuantitativamente, modelar un medio continuo mediante un sistema de partículas, implica transformar las integrales descritas en las secciones anteriores, ecuaciones 2.2–2.11, en un sumatorio, donde la masa de cada partícula representa un elemento de volumen cuyo valor es ρdV . Así pues, la formulación discreta de la ecuación 2.2 queda expresada en términos de la ecuación 2.13:
2.3 Formulación del Método SPH 22 ⟨f(ri)⟩=∑ j∈N(i) mj ρj f(xj)W(ri−rj,h),(2.13) donde N(i) hace referencia a la vecindad de cada partícula i . Vecindad que queda restringida por el radio del dominio soportado. La expresión de la suma extendida sobre las partículas vecinas, establecida cuantitativamente a través de la ecuación 2.13, permite obtener el valor de cualquier magnitud f definida sobre la partícula i , como un promedio ponderado por la función kernel, de los valores que esa misma magnitud ftiene en cada unas de sus jpartículas vecinas. Conocida la formulación discreta con la que se puede calcular cualquier función, o magnitud, f el siguiente paso es obtener la formulación discreta de sus derivadas. Para ello se va a reformular la ecuación 2.7 en términos de derivadas parciales. La expresión obtenida es la ecuación 2.14. ⟨∇·f(r)⟩=ZΩ∇·f(r′)W(r−r′,h)dr′=∂ ∂rZf(r′) ρ(r′)W(r−r′,h)ρ(r′)dr′(2.14) Considerando la condición de compacidad, la ecuación 2.14 se simplifica y permite aproximar el gradiente de cualquier magnitud escalar f(r) , en la forma discreta expresada a través de la ecuación 2.15: ⟨∇f(r)⟩=∂ ∂rZf(r′) ρ(r′)W(r−r′,h)ρ(r′)dr′≈∑ j∈N(i) mj fj ρj ∇W(rj−ri,h).(2.15) Reordenando la ecuación 2.15, es posible formular una expresión similar a la ecuación 2.11. Tras este proceso la ecuación obtenida es 2.16: ⟨∇·f(ri)⟩=−∑ j∈N(i) mj ρj f(rj)·∇W(ri−rj,h),(2.16) donde cabe destacar que ∇Wse particulariza para la partícula i. La ecuación 2.16 establece que el valor del gradiente de la función f , en la posición de la partícula i , se obtiene a partir de la suma ponderada, extendida a todas sus partículas vecinas j , del producto del valor que dicha función f tiene en la posición de la partícula i , por el gradiente de la función kernel. Siempre que se satisfaga la condición de compacidad, al igual que ocurría en la formulación integral (ecuación 2.11).
2.3 Formulación del Método SPH 23 Generalizando la ecuación 2.16, es posible formular cualquier derivada de orden superior. La ecuación general que se obtiene, se expresa en términos de la ecuación 2.17: ⟨∇l·f(r)⟩=−∑ j∈N(i) mj ρj f(rj)·∇lW(r−rj,h),(2.17) Cabe destacar que la ecuación 2.17, al igual que la ecuación 2.16, sólo ofrece resultados satisfactorios, desde el punto de vista de minimizar el error, si se cumple la condición de compacidad. Condición que no se cumple para las partículas cercanas al contorno, como se advierte en [ 71 ]. Otra limitación asociada a la ecuación 2.17, es que no permite que se cumpla la conservación del momento, según se destaca en [ 54 ], lo cual no la hace adecuada, con esta formulación, para su uso en simulaciones de fluidos. La consecuencia inmediata es que, aunque la ecuación 2.17 cumple las “premisas” establecidas en SPH, la precisión de los resultados obtenidos suele verse reducida. Para solventar este problema, se hace necesario modificar dicha ecuación, de manera que se obtenga una expresión simétrica, que satisfaga la conservación del momento y, por lo tanto, pueda ser utilizada en la simulación del flujo del fluido. El proceso para obtener dicha ecuación se describirá a continuación. Formulación de derivadas simétricas en SPH En el método SPH existen diferentes técnicas para formular expresiones simétricas de la ecuación 2.17, de modo que se satisfaga la conservación del momento. En este estudio, nos vamos a centrar en la propuesta desarrollada por Monaghan [ 54 ], cuyo uso está muy extendido en la simulación de fluidos dentro del ámbito de CG. Según se destaca en [ 54 ], el punto de partida son dos identidades funcionales basadas en la regla de la cadena. Con estas identidades, es posible introducir la densidad, con la que las ecuaciones obtenidas están optimizadas para su uso en la simulación de fluidos. Con esa formulación es posible obtener dos expresiones válidas, las ecuaciones 2.18 y 2.19. ∇lf(r) = 1 ρh∇l(ρf(r))−f(r)∇lρi,(2.18) ∇lf(r) = ρ∇lf(r) ρ+f(r) ρ2∇lρ.(2.19) Siguiendo el mismo procedimiento de aproximación, utilizado para obtener la ecuación 2.16, sobre las ecuaciones 2.18 y 2.19 se consigue la formulación de las ecuaciones 2.20 y 2.21.
2.3 Formulación del Método SPH 24 ⟨∇lf(ri)⟩=1 ρi"N ∑ j=1 mjf(rj)−f(ri)∇l iWi j#(2.20) ⟨∇lf(ri)⟩=ρi"N ∑ j=1 mj" f(rj) ρ2 j!+ f(ri) ρ2 j!#·∇l iWi j#,(2.21) Ambas ecuaciones son válidas, si bien la ecuación 2.21 es la que habitualmente se usa en el ámbito de CG. Como puede comprobarse, tanto la ecuación 2.20, como la ecuación 2.21 se formulan con magnitudes que aparecen por pares, tanto simétrica (referida a la suma de las magnitudes), como asimétricamente (asociada a la resta). De este modo, se solventa el problema de conservación del momento y, en consecuencia, se favorece la estabilidad de los resultados obtenidos, como se justifica en [54].
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 76 A partir de los resultados mostrados en la figura 6.1, podemos afirmar que las cuatro funciones kernel favorecen la obtención de resultados estables y precisos. Sin embargo, los mejores resultados se obtienen utilizando con la función Ci j –Cosine propuesta, ya que es la que permite la obtención de resultados con el mayor grado de reticularidad. Para cuantificar la reticularidad obtenida en el instante final, hemos formulado un parámetro de orden, ς , que satisface: ς=∑s k=1(|dl−dk i j|) s,(6.2) donde s= (m−1)(n−1) , y m y n son el número de partículas en cada dirección del espacio, dl es la longitud óptima de regularidad y dk i j es la distancia entre cada par de partículas vecinas. A partir de la formulación del parámetro ς , podemos afirmar que cuanto menor sea el valor, mayor será el grado de reticularidad alcanzado y, en consecuencia, mayor es la estabilidad y la precisión obtenida. Los valores obtenidos de este parámetro se muestran en la tabla 6.1. Quartic B–Spline Wendland6 Gaussiana C24–Cosine ς35.6 20.1 16.2 3.9 Tabla 6.1 Valores del parámetro de orden, ς , asociado a la prueba de sedimentación de una distribución aleatoria de partículas. Finalmente, podemos destacar que, de las distintas funciones kernel que podemos obtener dando valores a sus parámetros i j a la función propuesta, para esta prueba hemos utilizado la C24–Cosine ya que ha demostrado dar los mejores resultados. 6.2.2 Prueba Sod’s shock–tube La prueba Sod’s shock–tube nos permite analizar los resultados de la simulación de un fluido cuando, debido a diferencias de presión y densidad, aparecen ondas de choque y discontinuidades. Del test Sod’s shock–tube se conoce su solución exacta, por lo que es posible cuantificar la estabilidad y precisión, de cualquier técnica numérica, comparando resultado exacto y numérico, tal como se describe en los estudios [46, 56, 77]. Descriptivamente, el test Sod’s shock–tube se implementa sobre un tubo en el que se introduce una membrana no porosa. Esta membrana define dos secciones que están rellenadas de un fluido, un gas. En una sección, que se denomina sección de alta, el gas es alrededor de 10 más denso y más confinado que en la otra sección que se designa como sección de
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 77 baja. Con esta distribución, en un determinado instante considerado como instante inicial, la membrana se retira instantáneamente. A partir de ese momento se genera, al menos, tres regiones diferentes: (1) una onda de choque, como consecuencia de las diferencias de presión y densidad entre las secciones de alta y de baja. Esta onda de choque evoluciona hacia la sección de baja; (2) una onda de rarefacción que aparece en la sección de alta y al mismo tiempo que la onda de choque (3) una zona de discontinuidad que surge entre lo que fueron las secciones alta y baja iniciales, que suele denominarse como discontinuidad de contacto. A partir de la precisión obtenida en la simulación de estas tres regiones, es posible cuantificar la idoneidad de la función kernel utilizada, tal como es descrito por Sod [79]. Para llevar a cabo esta prueba, hemos considerado un dominio de dimensiones [-0.5,0.5] . La sección de alta está limitada entre [-0.5,0.0[ y la sección de baja está limitada entre ]0.0,0.5] . Inicialmente el número de partículas es nh=1200 y nl=300 , donde los subíndices hylse refieren a las secciones de alta y de baja, respectivamente. La distancia inicial entre partículas es ∆xh=5.8610−3 , para la sección de alta, y ∆xl=23.4410−3 , para la sección de baja y el radio inicial del dominio soportado es hh=hl=1.75∆xr . Para controlar el tamaño del dominio soportado del Ci j –Cosine, hemos asignado los valores de η=1.28 y K=1.85 . Los valores iniciales de la densidad de masa, presión y velocidad de cada sección son (ρh,ph,vh) = (1.0,1.0,0.0.0) y (ρl,pl,vl) = (0.125,0.1,0.0) . Para establecer el valor de la masa de las partículas, hemos considerado la sugerencia dada por [ 74 ]; por lo tanto, m=0.5(ρh∆xh+ρl∆xl) = 4.3910−3. Por otro lado, los valores de los parámetros de control de las ecuaciones de Euler (ecuaciones A.6 y A.7, véase Anexo A), son α=1.08 , αU=0.145 y β=1.85697 . Estos valores se han considerado, teniendo en cuenta las sugerencias dadas por [ 5 , 77 ]. Finalmente, hemos utilizado un paso de tiempo de ∆t=510−4. A partir de estos valores iniciales, hemos obtenido los resultados mostrados en la figura 6.2 en el instante t=0.15 para las magnitudes densidad, energía térmica y velocidad como se sugiere en los estudios [46, 74, 77]. En la figura 6.2 mostramos una comparativa entre la solución exacta, línea sólida, y la solución numérica, puntos, usando las cuatro funciones kernel. En cada gráfica se pueden localizar las tres regiones que aparecen cuando se retira la membrana. La discontinuidad del contacto aparece alrededor de xcd =1.004 , la onda de rarefacción aparece desde x<xcd y la onda de choque desde x>xcd . Aunque la posición de la discontinuidad de contacto es la misma en todas las gráficas, tanto en el ancho espacial de la discontinuidad de contacto, como en la existencia de una ondulación, existen variaciones entre los distintos resultados obtenidos para cada función kernel. Varias investigaciones analizan la causa de este error [ 46 , 56 , 71 , 77 ], de sus conclusiones se pueden deducir dos causas principales. La primera,
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 78 Energía Térmica (a) Densidad (b) Energía Térmica (c) Velocidad Figura 6.2 Comparativa entre la solución exacta y las soluciones numéricas, de los resultados obtenidos para la prueba Sod’s shock–tube. El valor de los parámetros de la función Ci j – Cosine son i=3 y j=5 para la densidad y la velocidad y i=2 y j=4 para la energía térmica, con los que se obtienen los mejores resultados. es el calentamiento en la discontinuidad del contacto, asociado al balance de energía térmica que no se simula completamente por la formulación variacional de SPH. La segunda es la imprecisión al calcular la presión a través de la discontinuidad del contacto, lo que induce un aumento de la velocidad y la energía térmica en esa posición. Para evitar esta imprecisión se propone un aumento de resolución, según la sugerencia de [ 46 ], pero este aumento debe realizarse sin incrementar el radio del dominio soportado. Este requisito se satisface con la función Ci j –Cosine propuesta, que nos permite obtener el menor ancho y amplitud de la onda en la discontinuidad del contacto, como se muestra en la figura 6.2. Asimismo, la cola de la onda de rarefacción, ubicada alrededor de x=−0,25 , muestra un “desbordamiento” que no debería aparecer. Sin embargo, tanto la amplitud como el ancho espacial de los resultados obtenidos disminuyen, desde los resultados obtenidos utilizando el Quartic B–Spline hasta
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 79 los resultados de nuestra kernel Ci j –Cosine. Este mismo decrecimiento se obtiene para la onda de choque antes de la discontinuidad del contacto. Para cuantificar la precisión que obtenemos de cada función kernel, hemos seleccionado los valores numéricos de cada parte del flujo, tales como onda de rarefacción, discontinuidad de contacto y onda de choque, tanto la amplitud máxima, como su anchura. Además, para cuantificar la precisión global, calcularemos el coeficiente de determinación, r2 . Los resultados obtenidos para la densidad, la energía térmica y la velocidad se muestran en las tablas 6.2, 6.3 y 6.4, respectivamente. En estas tablas, Art es la amplitud máxima de la cola de rarefacción, Wrt es el ancho de la cola de rarefacción, Acd es la amplitud máxima de la discontinuidad del contacto, Wcd es el ancho de la discontinuidad del contacto, xsw es la posición de la onda de choque y Wsw es la anchura de la onda de choque. r2Art Wrt Acd Wcd xsw Wsw Quartic B–Spline 0.9483 0.0576 0.131 0.021 0.098 0.3803 0.217 Wendland6 0.9701 0.0412 0.124 0.012 0.082 0.3721 0.124 Gaussiana 0.9728 0.0314 0.093 0.009 0.068 0.3689 0.103 Ci j–Cosine 0.9899 0.0091 0.0162 0.003 0.029 0.3618 0.059 Tabla 6.2 Valores relevantes obtenidos de las gráficas de la densidad para la prueba Sod’s shock–tube, figura 6.2a. r2Art Wrt Acd Wcd xsw Wsw Quartic B–Spline 0.9652 0.0374 0.131 0.0542 0.156 0.3792 0.186 Wendland6 0.9784 0.0104 0.124 0.042 0.097 0.3718 0.114 Gaussiana 0.9875 0.0169 0.092 0.031 0.056 0.3685 0.084 Ci j–Cosine 0.9941 0.0098 0.016 0.009 0.025 0.3623 0.019 Tabla 6.3 Valores relevantes obtenidos de las gráficas de la energía térmica para la prueba Sod’s shock–tube, figura 6.2b. 6.2.3 Prueba de onda expansiva La prueba onda expansiva (blast–wave test) es una versión más extrema del test Sod’s shock– tube y se formula para un muy alto número de Mach, alrededor M=200 . Fue desarrollada por Woodward y Colella [ 85 ] para analizar la precisión de un determinado método numérico a la hora de simular el comportamiento del fluido en una situaciones con alta discontinuidad de contacto. Esta prueba esta diseñada, especialmente para analizar la precisión de la velocidad,
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 80 r2Art Wrt Acd Wcd xsw Wsw Quartic B–Spline 0.9214 0.0312 0.158 0.048 0.206 0.3803 0.146 Wendland6 0.9431 0.0223 0.136 0.027 0.101 0.3721 0.089 Gaussiana 0.9648 0.015 0.098 0.017 0.072 0.3694 0.065 Ci j–Cosine 0.9895 0.0097 0.041 0.004 0.056 0.3618 0.052 Tabla 6.4 Valores relevantes obtenidos de las gráficas de la velocidad para la prueba Sod’s shock–tube, figura 6.2c. el comportamiento de la presión y la densidad a altos valores. En el test de onda expansiva, la evolución del flujo muestra las mismas características que la prueba de Sod’s shock–tube, es decir, una discontinuidad de contacto aparecerá insertada entre la onda de choque y la onda de rarefacción. Debido a las similitudes con el test Sod’s shock–tube, en esta implementación sólo hemos cambiado los valores iniciales de algunos parámetros de simulación. Los nuevos valores son para la densidad de masa, la presión y la energía térmica, estos son: (ρh, ph,uh) = (1.0,1103,0.0) y (ρl, pl,ul) = (1.0,0.01,0.0) , según las recomendaciones dadas por [ 77 , 85 ]. Así mismo, reduciremos el paso de tiempo, debido al alto valor de la velocidad inicial, en este caso hemos utilizado ∆t = 810−5 . Usando estos valores, hemos obtenido los resultados que se muestran en la figura 6.3 en el instante de tiempo t = 0.0075 . En dichas gráficas se muestran las soluciones exactas, representada por líneas continuas, mientras que los resultados numéricos se indican mediante puntos. Análogamente a la fenomenología producida en el test de Sod’s shock–tube, en esta simulación, también aparecen tres regiones que aparecen a partir del instante en el que se retira la membrana, mostradas en las gráficas de la figura 6.3, estas son: la discontinuidad de contacto, que se ubica xcd = 0.1081 , la onda de rarefacción x < xcd y la onda de choque x > xcd . Así mismo, en cada una de las gráficas mostradas en la figura 6.3, el ancho de discontinuidad de contacto cambia para cada función kernel, siendo el menor el que se obtiene con la función Ci j –Cosine propuesta, por lo que es la que ofrece los mejores resultados, ya que idealmente este debería ser cero, como se deduce de los resultados exactos. Cabe destacar, en las gráficas de la energía térmica y de la velocidad de la figura 6.3, el valor sobre–estimado que aparece en la energía térmica y la oscilación de la velocidad en la posición de la discontinuidad de contacto. Esta deficiencia se relaciona con la imprecisión en el modelado de las ecuaciones de Euler, que no permiten determinar, con precisión, el punto donde se produce el calentamiento, como se describe en [ 74 ]. Así mismo, estas perturbaciones también están relacionadas con los valores de la presión entorno a la discontinuidad, que no son precisos en simulaciones SPH, lo que afecta a la energía térmica. Por otro lado, la cola
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 81 Energía Térmica (a) Densidad (b) Energía Térmica (c) Velocidad Figura 6.3 Comparativa entre la solución exacta y las soluciones numéricas, de los resultados obtenidos para la prueba de onda expansiva. El valor de los parámetros de la función Ci j – Cosine son i=3 y j=5 para la densidad y la velocidad y i=2 y j=4 para la energía térmica, con los que se obtienen los mejores resultados. de rarefacción muestra un desbordamiento en las gráficas de velocidad; desbordamiento que aparece, en menor medida, en las gráficas de presión. Este desbordamiento tiene su origen en los problemas de precisión asociados a la energía térmica y, más concretamente, con el calentamiento no disipado. Aunque estas perturbaciones están relacionadas principalmente con la formulación de la ecuación de Euler mediante SPH, como se advierte en [ 74 ], la precisión que se obtiene depende de la función kernel utilizada. Este hecho queda confirmado por los resultados obtenidos y que se muestran en la figura 6.3, donde la mejor precisión se obtiene con la función Ci j–Cosine propuesta. Para demostrar la precisión que se obtiene en la simulación, con cada función del kernel, hemos seleccionado los valores más relevantes asociados a las gráficas de la figura 6.3. Estos
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 82 valores se muestran en la tabla 6.5 para la densidad de masa, en la tabla 6.6 para la energía térmica y en la tabla 6.7 para la velocidad. r2Wrt Acd Wcd xsw Wsw Quartic B–Spline 0.9304 0.141 0.613 0.065 0.193 0.142 Wendland6 0.9518 0.104 0.491 0.039 0.190 0.094 Gaussiana 0.9681 0.081 0.285 0.021 0.194 0.079 Ci j–Cosine 0.9887 0.056 0.146 0.009 0.196 0.038 Tabla 6.5 Valores relevantes obtenidos de las gráficas de la densidad para la prueba de onda expansiva, figura 6.3a. r2Wrt Acd Wcd xsw Wsw Quartic B–Spline 0.9187 0.079 218.063 0.0791 0.193 0.087 Wendland6 0.9498 0.038 114.277 0.0639 0.189 0.064 Gaussiana 0.9636 0.028 72.843 0.0453 0.191 0.051 Ci j–Cosine 0.9873 0.011 20.442 0.0206 0.196 0.021 Tabla 6.6 Valores relevantes obtenidos de las gráficas de la energía térmica para la prueba de onda expansiva, figura 6.3b. r2Art Wrt Acd Wcd xsw Wsw Quartic B–Spline 0.9098 0.4986 0.1066 0.491 0.1116 0.193 0.159 Wendland6 0.9321 0.3469 0.0945 0.398 0.0826 0.189 0.087 Gaussiana 0.9658 0.2877 0.0813 0.294 0.0582 0.191 0.026 Ci j–Cosine 0.9891 0.0934 0.0585 0.169 0.0348 0.196 0.009 Tabla 6.7 Valores relevantes obtenidos de las gráficas de la velocidad para la prueba de onda expansiva, figura 6.3c. 6.2.4 Prueba de impacto sobre muro En esta prueba simularemos la colisión con una pared rígida de un fluido frío uniformemente distribuido. Esencialmente, la prueba puede ser modelada mediante dos flujos idénticos que viajan en direcciones opuestas y que colisionan en la ubicación donde se situaría la pared rígida, según la propuesta de [ 56 ]. Como consecuencia de la interacción entre las dos corrientes, se generan dos ondas de choque simétricas que surgen desde la posición del punto
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 83 de colisión. Desde esta posición, hasta los límites de la simulación, aparecen un conjunto de fenómenos disipativos. Estos fenómenos son más evidentes alrededor de la región de colisión y en la cola de la onda de rarefacción. La fidelidad con la que se simulan estos fenómenos permite cuantificar la exactitud de los resultados obtenidos. Para implementar esta prueba hemos utilizado un dominio restringido entre [−0.4,0.4] , según la recomendación dada por [ 74 ], la ubicación del punto de colisión (la posición de la pared) será en x=0 . Así, se definen dos secciones: una sección izquierda y otra sección derecha que se extienden [−0.4,0.0[ y ]0.0,0.4] respectivamente. Las magnitudes dinámicas iniciales son (ρL,vL,UL) = (1.0,1.0,10−6) y (ρR,vR,UR) = (1.0,−1.0,10−6) , según la recomendación dada por [ 56 ], donde los subíndices L y R se refieren a la sección izquierda y derecha, respectivamente. El número de partículas utilizado es n=1500 , estas partículas están distribuidas equitativamente, a saber, nL=750 y nR=750 y el valor de la masa de partículas es mi=4.2910−3 . La distancia inicial entre partículas es ∆x=4.3910−3 y el radio inicial del dominio soportado es hL=hR=1.5∆x . Hemos utilizado un paso de tiempo ∆t=510−5 para la integración de tiempo y los valores de los parámetros disipativos, asociados a las ecuaciones Euler (ecuaciones A.6 y A.7 del Anexo ) son α=1.1 , αU=0.125 yβ=1.85. A partir de estas condiciones iniciales, hemos simulado la densidad, la energía térmica y el momento, según se recomienda en [74]. Los resultados se muestran en la figura 6.4 en el instante t=0.45 , donde las soluciones exactas y numéricas están representadas por líneas sólidas y puntos, respectivamente. Fenomenológicamente, la región más relevante se sitúa en torno a x=0, donde aparece el calor asociado a la pared, calor que no se disipa. Este calentamiento induce el error que muestra como un “pico” en los gráficos de la figura 6.4. Tanto el ancho como la altura del pico permiten establecer una métrica de la precisión obtenida a partir de cada función kernel, como se describe en [ 74 ]. Por otro lado, alrededor de cada cola de rarefacción, x=−0,16 y x=0,16 , aparece una caída en el valor de las gráficas del momento asociado al cambio brusco de esta magnitud. Analizando la figura 6.4, podemos afirmar que los resultados numéricos se aproximan a las soluciones exactas, excepto en el punto donde se ubicaría la pared (el punto en el que colisionan las dos corrientes) y las colas de rarefacción. Sin embargo, los resultados numéricos obtenidos a partir de la función kernel Ci j –Cosine propuesta, son más precisos que los obtenidos con las otras tres funciones kernel. Para resaltar la precisión obtenida utilizando cada función kernel, hemos seleccionado, en cada gráfica, los valores asociados a las partes más relevantes de flujo.
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 84 Energía Térmica Momento Energía Térmica Momento (a) Densidad Energía Térmica (b) Energía Térmica Momento (c) Momento Figura 6.4 Comparativa entre la solución exacta y las soluciones numéricas, de los resultados obtenidos para la prueba impacto sobre muro. El valor de los parámetros de la función Ci j –Cosine son i=2 y j=6 para la presión y la velocidad y i=2 y j=5 para la energía térmica, con los que se obtienen los mejores resultados. Para las gráficas de la densidad de masa y de la energía térmica, figuras 6.4a y 6.4b, los valores seleccionados son: la amplitud, la anchura y la posición del máximo asociado a la onda de choque, que designamos cómo Asw , Wsw y xsw , respectivamente. Los valores obtenidos se muestran en las tablas 6.8 y 6.9. Por otro lado, a partir de las gráficas del momento, figura 6.4c, podemos afirmar que las zonas más relevantes del flujo se encuentran en la cola de rarefacción. Por lo tanto, los valores seleccionados son la amplitud máxima, el ancho y la posición de cada cola de rarefacción. Estos valores los identificamos como AL rt , WL rt ,xL rt,AR rt,WR rt yxR rt en la tabla 6.10.
6.2 Pruebas de Estabilidad y Precisión para la Función Kernel 85 r2Asw Wsw xsw Quartic B–Spline 0.9614 0.682 0.059 2.4506 10−5 Wendland6 0.9761 0.529 0.041 1.8314 10−5 Gaussiana 0.9814 0.421 0.031 9.265 10−6 Ci j–Cosine 0.9928 0.187 0.014 2.05 10−7 Tabla 6.8 Valores relevantes obtenidos de las gráficas de la densidad para la prueba de impacto sobre muro, figura 6.4a. r2Asw Wsw xsw Quartic B–Spline 0.9671 0.708 0.0471 1.06024 10−6 Wendland6 0.9781 0.518 0.0317 0.05020 10−6 Gaussiana 0.9831 0.431 0.0271 0.7509 10−7 Ci j–Cosine 0.9903 0.207 0.0178 2.05 10−9 Tabla 6.9 Valores relevantes obtenidos de las gráficas de la energía térmica para la prueba de impacto sobre muro, figura 6.4b. r2AL rt WL rt xL rt AR rt WR rt xR rt Quartic B–Spline 0.9698 0.171 0.0916 -0.1911 0.146 0.0843 0.1906 Wendland6 0.9756 0.121 0.039 -0.1805 0.103 0.041 0.1865 Gaussiana 0.9884 0.082 0.026 -0.1830 0.091 0.023 0.1811 Ci j–Cosine 0.9926 0.055 0.016 -0.1919 0.057 0.010 0.1936 Tabla 6.10 Valores relevantes obtenidos de las gráficas del momento para la prueba de impacto sobre muro, figura 6.4c. 6.2.5 Prueba de Sjögreen La prueba de Sjögreen comparte ciertas similitudes con el test de impacto sobre muro, aunque este test fue diseñado para determinar la precisión y la estabilidad, de un determinado método numérico, en simulaciones a muy baja densidad. Descrito por Einfeldt et al. [ 25 ] su configuración inicial implica dos estados con baja densidad y presión. Configuración en la que, durante la evolución de la prueba, aparecen dos ondas de rarefacción que viajan en direcciones opuestas desde la posición de la discontinuidad. Durante el estado transitorio, se desarrolla una expansión mutua en la que no se genera un vacío en la posición inicial de la discontinuidad. Así pues la discontinuidad no avanza, como ocurría en las pruebas anteriores, esto genera una discontinuidad de contacto estacionario.
6.4 Resultados del Modelado de Interacción Fluido–Contorno 92 Sin embargo, a pesar de las semejanzas con estas metodologías, existen diferencias que deben ser implementadas “ad–hoc”. Este es el caso de las densidad asociada a las partículas de contorno, que no es necesario utilizar en nuestra propuesta, el número de iteraciones máxima en el proceso predictor–corrector de la técnica propuesta por [ 11 ] o la función kernel utilizada por [ 3 ] para calcular la densidad de las partículas de contorno. Para llevar a cabo la implementación de cada una de ellas, hemos utilizado los valores y recomendaciones dada por los autores de cada técnica. Para cada simulación hemos utilizado 5.0104 partículas fluidas. Con estas condiciones hemos obtenido los resultados que mostramos en la figura 6.9, todos ellos correspondiente al mismo instante de simulación. Las posibles diferencias son debidas a las distintas características de cada técnica. Analizando estas imágenes, se puede deducir que la técnica propuesta, figura 6.9d, ofrece unos buenos resultados desde el punto de vista de la estabilidad y de la adaptabilidad al contorno. Analizando cada imagen, se puede afirmar que el resultado obtenido a partir de la técnica propuesta por [ 66 ], figura 6.9a, si bien evita el comportamiento errático de las partículas, que es uno de los problemas más destacable en la interacción fluido–contorno en simulaciones mediante SPH, presenta dos importantes limitaciones: la primera es que muestra un sobreamortiguamiento artificial, la segunda es que se produce una ligera penetración en el contorno. Ello se debe a que el modelo de [ 66 ] necesita del uso de fuerzas de fricción no simétricas, cuya elección condiciona el realismo de la simulación. Por otro lado, en el resultado obtenido a partir de la propuesta de Becker et al. [ 11 ], figura 6.9b, se solventa este problema, es decir, no se muestra sobreamortiguamiento, ni penetración en el contorno. Sin embargo, muestra dispersión de las partículas y cierta separación del contorno. Este comportamiento se debe a que, aunque controla la fuerza de presión, algunas fuerzas dinámicas son despreciadas, como por ejemplo la viscosidad, lo que afecta a la estabilidad y favorece la dispersión. Con respecto al resultado obtenido a partir de la técnica propuesta por Akinci et al. [ 3 ], figura 6.9c, se puede afirmar que mejora los resultados de la técnica anterior. No obstante, presenta el mismo problema de la separación del contorno, limitación que es solventada en el resultado que se ha obtenido con nuestra propuesta, como se puede ver en la figura 6.9d. 6.4.2 Comparativa entre diferentes condiciones de contorno Una vez que hemos constatado las ventajas que ofrece nuestra propuesta de modelado de la interacción fluido–contorno, pasamos a mostrar las capacidades que tiene a la hora de implementar condiciones de contorno complejas. Condiciones de contorno que no pueden ser implementadas con las metodologías habitualmente empleadas en SPH, dentro del ámbito
6.4 Resultados del Modelado de Interacción Fluido–Contorno 93 de CG, y que afecta a la posibilidad de simular la interacción del fluido con diferentes materiales. En nuestro caso, vamos a simular la interacción con dos materiales estas son: mármol y madera. Cada uno de estos materiales van a tener diferentes coeficientes de rozamiento. Para el mármol hemos considerado una variación proporcional a la velocidad tanto para la dirección normal como para la tangencial –una pérdida del 50% para la normal y un 20% para tangencial (que designaremos como tipo 1). Para la madera hemos implementado la condición de Robin para la dirección tangencial, y de Dirichlet para la normal; para la tangencial los parámetros que hemos utilizado son: γ=0.35 , σ=0.8 y ζδΩ=3.5 (que designaremos como tipo 2). Esta elección es un modo de modelar la absorción que presenta la madera, lo que reducirá el flujo laminar sobre el elipsoide. Los resultados obtenidos se muestra en las figura 5.1. Analizando comparativamente las imágenes mostradas en la figura 6.10, podemos ver las diferencias en el comportamiento del fluido. Esto se debe a la diferencia en las condiciones de contorno, aunque el contorno de interacción sea el mismo. Con nuestra propuesta para gestionar la interacción fluido–contorno, hemos podido implementar condiciones como la tipo 2, modelada a partir de formulación diferencial. Tal como se modela dicha condición, se impide la formación de una capa límite, como se puede ver en las imágenes de 6.10b, a diferencia de la condición tipo 1. La importancia de simular la tipo 2 radica en que, con las metodologías utilizadas en SPH en el ámbito de CG, requiere de una simplificación para poder ser simulada. Simplificación que suele reducir el realismo.
6.4 Resultados del Modelado de Interacción Fluido–Contorno 94
6.4 Resultados del Modelado de Interacción Fluido–Contorno 95 (a) (b) (c) (d) (e) (f) Figura 6.6 Imágenes de simulaciones de fluidos obtenidas mediante diferentes técnicas. Las imágenes 6.6a y 6.6b se han obtenido a partir de la implementación del modelo WCSPH; las imágenes 6.6c y 6.6d se han obtenido mediante PCISPH; las imágenes 6.6e y 6.6f son el resultado de nuestro método híbrido Potencial–SPH.
6.4 Resultados del Modelado de Interacción Fluido–Contorno 96 Figura 6.7 Imágenes de la simulación de fluidos en dos momentos diferentes, activando y desactivando la disipación de calor. Las imágenes superiores se obtienen al activar la disipación de calor, en las imágenes inferiores se desactiva la disipación de calor.
6.4 Resultados del Modelado de Interacción Fluido–Contorno 97 (a) (b) (c) (d) Figura 6.8 Imágenes de la simulación del fluido en el instante t=1.25s obtenido para diferentes instantes de tiempo: figura 6.8a ∆t=0.001s , figura 6.8b ∆t=0.008s , figura 6.8c ∆t=0.02sy figura 6.8d ∆t=0.07s.
6.4 Resultados del Modelado de Interacción Fluido–Contorno 98 (a) (b) (c) (d) Figura 6.9 Comparativa de los resultados obtenidos con diferentes técnicas basadas en partículas de contorno. En la figura 6.9a, se muestra el resultado de la técnica propuesta por Becker et al. [ 11 ]. En la figura 6.9b, el obtenido a partir del modelo de por Oger et al [ 66 ]. En la figura 6.9c, el de la propuesta de Akinci et al. [ 3 ]. Finalmente, en la figura 6.9d, se muestra el resultado obtenido con nuestra propuesta. En cada imagen, se ha destacado la región donde mejor se muestra la capacidad de cada técnica.
6.4 Resultados del Modelado de Interacción Fluido–Contorno 99 (a) Condición de contorno tipo 1(mármol) (b) Condición de contorno tipo 2 (madera) Figura 6.10 Resultados de interacción del fluido con un elipsoide, cada dupla (fila de imágenes) se obtiene para el mismo instante. Comparativamente se pueden ver las diferencias.
Capítulo 7 Conclusiones 7.1 Conclusiones En esta tesis, hemos estudiado los problemas de inestabilidad e imprecisión que aparecen en las simulaciones de fluidos mediante SPH. Hemos analizado el origen de estos problemas y hemos deducido que si el número de partículas vecinas es bajo, del orden de algunas decenas, o su distribución en el dominio soportado es inhomogenea, la simulación se vuelve inestable y los resultados imprecisos. También hemos analizado cómo la distribución y el número de partículas vecinas bajo reduce la precisión en el gradiente de presión y en la energía térmica; imprecisión que se ve incrementada al imponer la condición de incompresibilidad. Así mismo, hemos descrito las distintas propuesta desarrolladas para evitar estos problemas, cuyo objetivo es incrementar el número y mejorar la distribución de partículas vecinas. Sin embargo, existen regiones límites, como son el contorno o las interfases, donde no es posible incrementar el número de partículas vecinas, por lo que los problemas de imprecisión e inestabilidad no están resueltos en estas regiones. Además, hemos detectado que, en el ámbito de CG, el modelado de la incompresibilidad no está completamente resuelto al simular líquidos y los modelos existentes pueden inducir inestabilidad e imprecisión, ya que no consideran los balances de energía fluido–contorno. La consecuencia es la dispersión errónea de las partículas en estas regiones. Con el objetivo de solventar estas limitaciones, hemos desarrollado tres propuestas, estas han sido: (1) una nueva función kernel, (2) un nuevo modelado del gradiente de presión para imponer la incompresibilidad y (3) una nueva gestión de la interacción fluido–contorno . En cada una de estas propuestas: 1. Hemos identificados las características analíticas de la función kernel que favorecen la estabilidad y la precisión. A partir de estas características hemos formulado una
7.1 Conclusiones 101 nueva función kernel paramétrica desarrollada como una composición de funciones armónicas. Por su formulación las hemos denominado como Ci j–Cosine. Para demostrar las mejoras en precisión y estabilidad que se obtienen con nuestra función kernel, hemos implementado un conjunto de pruebas estandarizadas tal como describimos en el capítulo 6. El objetivo de estas pruebas ha sido desarrollar un entorno comparativo entre la función propuesta y otras funciones kernel muy utilizadas en las simulaciones mediante SPH. A partir de los resultados obtenidos en la sección 6.2 podemos concluir que: • El kernel Ci j –Cosine nos permite obtener una alta estabilidad. Esta conclusión puede deducirse de los resultados obtenidos en el test de sedimentación, figura 6.1, donde los resultados de Ci j –Cosine muestran un mayor grado de reticularidad en comparación con los obtenidos a partir las otras tres funciones kernel. Desde el punto de vista cuantitativo, esta mejora se refleja en los resultados del parámetro ς , tabla 6.1, donde los valores obtenidos por las otras tres funciones kernel son de un valor mayor que los ofrecidos por nuestra Ci j –Cosine. Lo que significa, según las características del parámetro ς , que Ci j –Cosine es la que mayor estabilidad garantiza. • El kernel Ci j –Cosine nos permiten obtener resultados altamente precisos en simulaciones con condiciones extremas, en las que aparecen ondas de choque y la discontinuidades en las magnitudes que describen al fluido. Podemos hacer esta afirmación a partir de los resultados obtenidos en la prueba de Sod (Sod’s shock–tube), sección 6.2.2. En las gráficas de la figura 6.2 se muestran los resultados obtenidos para la densidad, la energía térmica y la velocidad. A partir de estos resultados, podemos inferir que con nuestra propuesta se obtienen los resultados más precisos. Esta precisión es más evidente en la cola de rarefacción, la discontinuidad del contacto y la onda de choque. Para reforzar este hecho, hemos seleccionado los valores asociados a estas regiones que se muestran en las tablas 6.2, 6.3 y 6.4. Todos los valores obtenidos evidencian que la mejor precisión se obtiene a partir de Ci j–Cosine. • Con la función kernel Ci j –Cosine propuesta conseguimos precisión incluso con altos valores de velocidad que induce la aparición de regiones con alta discontinuidad, como se ha demostrado en el test de onda expansiva (blast–wave test), que representa una versión más extrema del test de Sod. Al igual que en los resultados obtenidos en esa prueba, la precisión más alta se obtiene con la función
8.1 Conclusions 108 that the smaller these two magnitudes of dip, the better the obtained accuracy, since ideally, this disturbance should not appear. Based on this fact and analysing the graphs shown in figure 6.4, we can affirm the proposed kernel function, Ci j – Cosine, allows us to obtain the most accurate results. Numerically, this fact is also confirmed by the values shown in tables 6.8, 6.9 and 6.10 where it is shown the greatest accuracy is obtained from the Ci j–Cosine function. • The proposed kernel function also ensures accuracy when density and pressure are low, conditions that are simulated in Sjögreen’s test. The obtained results by our Ci j –Cosine are the most accurate, as can be deduced from the graphs shown in the figure 6.5. Furthermore, this statement is confirmed by the results shown in tables 6.11, 6.12 and 6.13, from which it can be deduced that, again, the highest precision is obtained by our kernel function. On the basis of these facts, we can state that the proposed kernel function, Ci j –Cosine, offers the most accurate results in all tests we have submitted, even under extreme simulation constraints. 2. We have identified the factors that relate instability to the pressure gradient. We have concluded the incompressibility enforcing in liquids simulation, exacerbates the instability problems, especially when the number of neighbouring particles is low. The consequence in fluid simulation is that instability results in erratic particle dispersions, especially on fluid boundaries. In addition, we have found the energy balance in fluid–boundary interaction must be considered to avoid instability. To solve these problems related to incompressibility enforcing, we have developed a hybrid Potential—SPH technique, which guarantees stability even to a very low number of particles. Particularly we have achieved to: • Improve stability in pressure force, we have formulated a pressure force independent of the density, that is defined on each pair of particles, so that this stability is maintained even to situations of very low number of neighbouring particles. Thus, we have obtained stable results in all fluid regions, especially in boundaries. This statement is confirmed by the results shown in figure 6.6, where it can be seen the stability achieved especially in boundaries. • Enforce the incompressibility an efficient way. Efficiency is since we attain to enforce explicitly incompressibility and in a single step, since our proposal is not a predictive–corrective process. The obtained improvement can be deduced from the obtained results in the simulation shown in figure 6.6.
8.1 Conclusions 109 • Detect the heating into the fluid and to dissipate it in an efficient and accurate way. This constitutes a highlight novelty since the existing methods in Computer Graphics, based on SPH, do not consider heating. The importance of heating, especially in boundaries, is shown in the simulations shown in figure 6.7 where deactivating heat dissipation leads to instabilities in boundaries. • Improve the results obtained by other techniques of explicit incompressibility enforcement. This improvement is shown in figure 6.6 where we have performed a comparison to other explicit incompressibility enforcement techniques. • Develop a stable and efficient model for a wide values range. As it is shown in figure 6.8 where the result of obtained simulation appears, for the same moment time, but with different time steps. In addition, in table 6.14 we show the result of a comparative study of the stability ranges for different methods. Particularly important is the number range of neighbouring particles in which we can operate and the time steps in which the system is stable. Thus, we can conclude that our proposal to explicitly incompressibility enforcement improves the stability in fluid simulations. In addition, we have developed an efficient model that allows us to manage the energy balances, especially in boundaries. 3. We have identified the stability problems that are related to fluid–boundary interaction. Likewise, we have analysed the different existing techniques for modelling the boundary constraints that govern the fluid behaviour in these regions. Based on these studies, we have developed a new technique that guarantees stability and allows the simulation of more complex boundary constraints than can be simulated by existing Computer Graphics techniques. Particularly, the obtained achievements by our technique of fluid–boundary interaction are: • An optimized method for particle–modelled boundaries in which only particle interaction of the fluid–particle boundary is managed. In this way, it is possible to manage the interaction without the need to generate artificial magnitudes, such as: density or pressure, on the boundary particles. The obtained advantage is not only operational, but also avoids the particles density in fluid, which interact on boundary, from acquiring erroneous values that could introduce instability and inaccuracy. • We have managed to optimize the adaptability fluid to boundaries. Due to the structure with which we model our boundary and the developed formulation, the fluid particles stick to the boundary without the need to implement artificial
8.1 Conclusions 110 operations. This fact can be deduced from the results shown in figure 6.9, where the obtained results by implementing some of the methodologies usually used in SPH are illustrated. • We have developed a formulation that allows us to model constraints in Computer Graphics scope which until now could only be simulated in the Computational Fluid Dynamic scope. By this formulation, we have expanded the spectrum of fluid simulations. In this way, it is possible to simulate complex fluid behaviours, as it is shown in figure 6.10 where it can be seen the effect of modifying boundary constraint with directional gradient. Simulation that it is not possible to carry out by the existing Computer Graphics methodologies. Therefore, we can conclude that we have developed a series of techniques, with which stability and precision in fluid simulations are improved by means of SPH and that, due to their simplicity, they can be used in simulations into Computer Graphics scope, without restricting their use to more accuracy implementations.
8.1 Conclusions 111
Referencias [1] Abel, T. (2011). rpsph: A novel smoothed particle hydrodynamics algorithm. Monthly Notices of the Royal Astronomical Society, 413(1):271–285. [2] Afanasiev, K., Makarchuk, R., and Popov, A. (2008). Comparative analysis of the sph and isph methods. Notes on Numerical Fluid Mechanics and Multidisciplinary Design, 101:206–223. [3] Akinci, N., Ihmsen, M., Akinci, G., Solenthaler, B., and Teschner, M. (2012). Versatile rigid–fluid coupling for incompressible sph. ACM Trans. Graph., 31(4):62:1–62:8. [4] Alduán, I., Tena, A., and Otaduy, M. A. (2015). Efficient and robust position-based fluids for vfx. In Proc. of Congreso Español de Informática Gráfica. [5] Amini, Y., Emdad, H., and Farid, M. (2011). A new model to solve fluid-hypo-elastic solid interaction using the smoothed particle hydrodynamics (sph) method. European Journal of Mechanics - B/Fluids, 30(2):184–194. [6] Ata, R. and Soulaïmani, A. (2005). A stabilized sph method for inviscid shallow water flows. International Journal for Numerical Methods in Fluids, 47:139–159. [7] B. Liu, M. and Liu, G. (2010). Smoothed particle hydrodynamics (sph): an overview and recent developments. Archives of Computational Methods in Engineering, 17(1):25–76. [8] Band, S., Gissler, C., Ihmsen, M., Cornelis, J., Peer, A., and Teschner, M. (2018). Pressure boundaries for implicit incompressible sph. ACM Trans. Graph., 37(2):14:1– 14:11. [9] Batchelor, G. (1974). An Introduction to Fluid Dynamics. Cambridge Mathematical Library. [10] Becker, M. and Teschner, M. (2007). Weakly compressible sph for free surface flows. In Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation, SCA ’07, pages 209–217. [11] Becker, M., Tessendorf, H., and Teschner, M. (2009). Direct forcing for lagrangian rigidfluid coupling. IEEE Transactions on Visualization and Computer Graphics, 15(3):493– 503. [12] Belytschko, T., Guo, Y., Liu, K., W., and Xiao, P. (2000). A unified stability analysis of meshless particle methods. International Journal for Numerical Methods in Engineering, 48(9):1359–1400.
Referencias 113 [13] Bender, J., Müller, M., and Macklin, M. (2015). Position–based simulation methods in computer graphics. In Zwicker, M. and Soler, C., editors, EUROGRAPHICS 2015 Tutorials. Eurographics Association. [14] Cabezón, R. M., Garcia-Senz, D., and Relano, A. (2000). A one–parameter family of interpolating kernels for smoothed particle hydrodynamics studies. Journal of Computational Physics, 227(19):8523–8540. [15] Campbell, P. M. (1989). Some new algorithms for boundary value problems in smooth particle hydrodynamics. Mission Research Corp. [16] Capuzzo-Dolcetta, R. and Di Lisio, R. (2000). A criterion for the choice of the interpolation kernel in smoothed particle hydrodynamics. Applied Numerical Mathematics, 34(4):363–371. [17] Chow, A., Rogers, B., Lind, S., and Stansby, P. (2018). Incompressible sph (isph) with fast poisson solver on a gpu. Computer Physics Communications, 226:81–103. [18] Colagrossi, A. and Landrini, M. (2003). Numerical simulation of interfacial flows by smoothed particle hydrodynamics. Journal of Computational Physics, 191(2):448–475. [19] Cornelis, J., Ihmsen, M., Peer, A., and Teschner, M. (2015). Liquid boundaries for implicit incompressible sph. Computers & Graphics, 52:72–78. [20] DE, L., Le Touzé, D., and Alessandrini, B. (2009). Normal flux method at the boundary for sph. In Proceedings of the 4th SPHERIC Workshop. [21] Dehnen, W. and Aly, H. (2012). Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. Monthly Notices of the Royal Astronomical Society, 425(2):1068–1082. [22] Di Monaco, A., Manenti, S., Gallati, M., Sibilla, S., Agate, G., and Guandalini, R. (2011). Sph modeling of solid boundaries through a semi–analytic approach. Engineering Applications of Computational Fluid Mechanics, 5:1–15. [23] Dilts, G. A. (1999). Moving–least–squares–particle hydrodynamics-i.consistency and stability. International Journal for Numerical Methods in Engineering, 44(8):1115–1155. [24] Dyka, C. T., Randles, P. W., and Ingel, R. P. (1997). Stress points for tension instability in sph. International Journal for Numerical Methods in Engineering, 40(13):2325–2341. [25] Einfeldt, B., Munz, C., Roe, P., and Sjögreen, B. (1991). On godunov-type methods near low densities. Journal of Computational Physics, 92(2):279–295. [26] Erleben, K., Bender, J., and Teschner, M. (2010). Boundary handling and adaptive time– stepping for pcisph. In Workshop in Virtual Reality Interactions and Physical Simulation “VRIPHYS” (2010), pages 339–346. The Eurographics Association. [27] Feldman, J. and Bonet, J. (2007). Dynamic refinement and boundary contact forces in sph with applications in fluid flow problems. International Journal for Numerical Methods in Engineering, 72(3):295–324.
Referencias 114 [28] Ferrand, M., Laurence, D. R., Rogers, B. D., Violeau, D., and Kassiotis, C. (2012). Unified semi–analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless sph method. International Journal for Numerical Methods in Fluids, 71(4):446–472. [29] Ferrari, A., Dumbser, M., Toro, E. F., and Armanini, A. (2009). A new 3d parallel sph scheme for free surface flows. Computers and Fluids, 38(6):1203–1217. [30] Fourtakas, G., Stansby, P., Rogers, B., and Lind, S. (2018). An eulerian–lagrangian incompressible sph formulation (eli–sph) connected with a sharp interface. Computer Methods in Applied Mechanics and Engineering, 329(Supplement C):532–552. [31] Fulk, D. A. (1994). A numerical analysis of smoothed particle hydrodynamics. AIR FORCE INST OF TECH WRIGHT-PATTERSON AFB OH. [32] Gingold, R. and Monaghan, J. (1977). Smoothed particle hidrodynamics: Theory and application to non–spherical stars. Royal Astronomical Society, 181:375–389. [33] Harada, T., Koshizuka, S., and Kawaguchi, Y. (2007). Smoothed particle hydrodynamics on gpus. In Computer Graphics International, pages 63–70. [34] Hashemi, M., Fatehi, R., and Manzari, M. (2012). A modified sph method for simulating motion of rigid bodies in newtonian fluid flows. International Journal of Non–Linear Mechanics, 47(6):626–638. [35] He, X., Liu, N., Li, S., Wang, H., and Wang, G. (2012). Local poisson sph for viscous incompressible fluids. Computer Graphics Forum, 31(6):1948–1958. [36] Ihmsen, M., Cornelis, J., Solenthaler, B., Horvath, C., and Teschner, M. (2014a). Implicit incompressible sph. IEEE Transactions on Visualization and Computer Graphics, 20(3):426–435. [37] Ihmsen, M., Orthmann, J., Solenthaler, B., Kolb, A., and Teschner, M. (2014b). Sph fluids in computer graphics. In Eurographics 2014 - State of the Art Reports. The Eurographics Association. [38] J., C. A. (1997). A numerical method for solving incompressible viscous flow problems. Journal of Computational Physics, 135(2):118–125. [39] Kao, C., Osher, S., and Qian, J. (2004). Lax–friedrichs sweeping scheme for static hamilton–jacobi equations. Journal of Computational Physics, 196(1):367–391. [40] Kawata, D., Okamoto, T., Gibson, B., Barnes, D., and Cen, R. (2013). Calibrating an updated smoothed particle hydrodynamics scheme within gdc+. Monthly Notices of the Royal Astronomical Society, 428(3):1968–1979. [41] Kiara, A., Hendrickson, K., and Yue, D. (2013). Sph for incompressible free-surface flows. part i: Error analysis of the basic assumptions. Computers & Fluids, 86(5):611–624. [42] Koschier, D. and Bender, J. (2017). Density maps for improved sph boundary handling. In Spanish Computer Graphics Conference (CEIG), SCA ’17. ACM.
Referencias 115 [43] Leroy, A., Violeau, D., Ferrand, M., and Kassiotis, C. (2014). Unified semi–analytical wall boundary conditions applied to 2–d incompressible sph. Journal of Computational Physics, 261:106–129. [44] Lind, S., Xu, R., Stansby, P., and Rogers, B. (2012). Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves. Journal of Computational Physics, 231(4):1499–1523. [45] Liu, G. and Liu, M. (2003). Smoothed Particle Hydrodynamics:A Meshfree Particle Method. World Scientific. [46] Liu, M., Liu, G., and Lam, K. (2003). Constructing smoothing functions in smoothed particle hydrodynamics with applications. Journal of Computational and Applied Mathematics, 155(2):263–284. [47] Losasso, F., Talton, J., Kwatra, J., and Fedkiw, R. (2008). Two-way coupled sph and particle level set fluid simulation. IEEE Transactions on Visualization and Computer Graphics, 14(4):797–804. [48] Lucy, L. B. (1977). A numerical approach to the testing of the fission hypothesis. Astronomical Journal, 82(12):1013–1024. [49] Maciá, F., Antuono, M., González, L. M., and Colagrossi, A. (2011). Theoretical analysis of the no–slip boundary condition enforcement in sph methods. Progress of Theoretical Physics, 125(6):1091–1121. [50] Macklin, M. and Müller, M. (2013). Position based fluids. ACM Transactions on Graphics, 32(4):104:1–104:12. [51] Marrone, S., Antuono, M., Colagrossi, A., Colicchio, G., Le Touzé, D., and Graziani, G. (2011). δ-sph model for simulating violent impact flows. Computer Methods in Applied Mechanics and Engineering, 200(13–16):1526–1542. [52] Mayrhofer, A., Rogers, B., Violeau, D., and Ferrand, M. (2013). Investigation of wall bounded flows using sph and the unified semi-analytical wall boundary conditions. Computer Physics Communications, 184. [53] Modak, S. and Sotelino, E. D. (2002). The generalized method for structural dynamics applications. Advances in Engineering Software, 33(7):565–575. [54] Monaghan, J. (1992). Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics, 30:543–574. [55] Monaghan, J. (1994). Simulation free surface flows with sph. Journal of Computational Physics, 110(2):399–406. [56] Monaghan, J. (1997). Sph and riemann solvers. Journal of Computational Physics, 136(2):298–307. [57] Monaghan, J. (2005). Smoothed particle hydrodynamics. Rep. Prog. Phys., 68(8):1703– 1759.
Referencias 116 [58] Monaghan, J. (2011). A turbulence model for smoothed particle hydrodynamics. European Journal of Mechanics–B Fluids, 30(4):360–370. [59] Monaghan, J. and Kajtar, J. (2009). Sph particle boundary forces for arbitrary boundaries. Computer Physics Communications, 180(10):1811–1820. [60] Morris, J., , and Monaghan, J. (1997). A switch to reduce sph viscosity. Journal of Computational Physics, 136(1):41–50. [61] Morris, J. (1995). A study of the stability properties of smooth particle hydrodynamics. Publications Astronomical Society of Australia, 13(1):97–102. [62] Müller, M., Charypar, D., and Gross, M. (2003). Particle-based fluid simulation for interactive applications. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation, SCA ’03, pages 154–159. [63] Müller, M., Charypar, D., and Gross, M. (2003). Particle-based fluid simulation for interactive applications. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 154–159. Eurographics Association. [64] Müller, M., Schirm, S., Teschner, M., Heidelberger, B., and Gross, M. (2004). Interaction of fluids with deformable solids: research articles. Computational Animation Virtual Worlds, 15(3–4):159–171. [65] Newmark, N. M. (1959). A method of computation for structural dynamics. Journal of the engineering mechanics division, 85(3):67–94. [66] Oger, G., Doring, M., Alessandrini, B., and Ferrant, P. (2006). Two–dimensional sph simulations of wedge water entries. Journal of Computational Physics, 213(2):803–822. [67] Okamoto, T., Jenkins, A., Eke, V. R., Quilis, V., and Frenk, C. S. (2003). Momentum transfer across shear flows in smoothed particle hydrodynamic simulations of galaxy formation. Monthly Notices of the Royal Astronomical Society, 345(2):429–446. [68] Perea, J. J. and Cordero, J. M. (2016). A new potential-based sph method. In Spanish Computer Graphics Conference (CEIG). The Eurographics Association. [69] Perea, J. J. and Cordero, J. M. (2017). An improved parallel technique for neighbour search on cuda. In Spanish Computer Graphics Conference (CEIG). The Eurographics Association. [70] Price, D. J. (2008). Modelling discontinuities and kelvin–helmholtz instabilities in {SPH}. Journal of Computational Physics, 227(24):10040–10057. [71] Price, D. J. (2012). Smoothed particle hydrodynamics and magnetohydrodynamics. Journal of Computational Physics, 231(3):759–794. [72] Price, D. J. and Monaghan, J. J. (2004). Smoothed particle magnetohydrodynamics-i. algorithm and tests in one dimension. Monthly Notices of the Royal Astronomical Society, 348(1):123–138. [73] Puri, K. and Ramachandran, P. (2014a). Approximate riemann solvers for the godunov sph (gsph). Journal of Computational Physics, 270(Supplement C):432–458.
Referencias 117 [74] Puri, K. and Ramachandran, P. (2014b). A comparison of sph schemes for the compressible euler equations. Journal of Computational Physics, 256(0):308–333. [75] Raveendran, K., Wojtan, C., and Turk, G. (2011). Hybrid smoothed particle hydrodynamics. In Proceedings of the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pages 33–42. [76] S., M., Colagrossi, A., Antuono, M., Lugni, C., and Tulin, M. (2011). A 2d+t sph model to study the breaking wave pattern generated by fast ships. Journal of Hydrodynamics, Ser. B, 27(8):1199–1215. [77] Sigalotti, L. and López, H. (2008). Adaptive kernel estimation and sph tensile instability. Computers & Mathematics with Applications, 55(1):23–50. [78] Sigalotti, L., López, H., Donoso, A., Sira, E., and J., K. (2006). A shock-capturing sph scheme based on adaptive kernel estimation. Journal of Computational Physics, 212(1):124–149. [79] Sod, A. (1978). A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. Journal of Computational Physics, 27(1):1–31. [80] Solenthaler, B. and Pajarola, R. (2009). Predictive-corrective incompressible sph. ACM Transactions Graphics, 28(3):41–46. [81] Springel, V. (2010). Smoothed particle hydrodynamics in astrophysics. Annual Review of Astronomy and Astrophysics, 48(1):391–430. [82] Sun, P., Ming, F., and Zhang, A. (2015). Numerical simulation of interactions between free surface and rigid body using a robust sph method. Ocean Engineering, 98(4):32–49. [83] Swegle, J., Hicks, D., and S.W., A. (1995). Smoothed particle hydrodynamics stability analysis. Journal of Computational Physics, 116(1):123–134. [84] Tetsuya, T., Yoshinori, D., Tomoyuki, N., and C., L. M. (2017). An efficient hybrid incompressible sph solver with interface handling for boundary conditions. Computer Graphics Forum, 37(1):313–324. [85] Woodward, P. and Colella, P. (1984). The numerical simulation of two–dimensional fluid flow with strong shocks. Journal of Computational Physics, 54(1):115–173. [86] X., X. and D., X.-L. (2016). An improved weakly compressible sph method for simulating free surface flows of viscous and viscoelastic fluids". Computer Physics Communications, 201:43–62. [87] Xu, R., Stansby, P., and Laurence, D. (2009). Accuracy and stability in incompressible sph (isph) based on the projection method and a new approach. Journal of Computer Physics, 228(18):6703–6725. [88] Yang, X., Peng, S., and Liu, M. (2014). A new kernel function for {SPH} with applications to free surface flows. Applied Mathematical Modelling, 38(15-16):3822– 3833.