Modelling the Fluid–Boundary Interaction in SPH

CEIG - Spanish Computer Graphics Conference (2018) I. García-Fernández and C. Ureña (Editors) Modelling the Fluid–Boundary Interaction in SPH Juan J. Perea and Juan M. Cordero University of Seville, Spain Abstract Smoothed Particle Hydrodynamics (SPH) is a numerical method based on mutually interacting meshfree particles, and has been widely applied to fluid simulation in Computer Graphics. Originally SPH does not define the behaviour of the particle system in the contour, so the different variants of SPH have been solving this deficiency with different techniques. Some of these techniques are based on fictitious forces, specular particles or semi–analytic fields. However, all these proposals present a drawback, that are may introduce additional inaccuracy as a divergent behaviour of the particle dynamics or an artificial separation between the fluid limits and the contour. To solve these limitations at this paper presents a new technique based on contour particles that are used during simulation to model the interaction with the fluid. The use of contour particles had already been used in other works to construct the contour like a particle layer. That solution presents problems especially when increasing the complexity of the contour shape. In addition, unlike other techniques, this paper presents an additional advantage, the possibility of obtaining all the dynamic magnitudes for improving efficiency and versatility. CCS Concepts •Computing methodologies →Collision detection; Physical simulation; 1. Introducción La complejidad de los fenómenos que aparecen en el flujo de un fluido, se incrementa cuando éste interactúa con el contorno. Esta interacción suele inducir la aparición de fenómenos como la ruptura del flujo o la formación de turbulencia, que sólo se pueden cuantificar mediante el uso de métodos numéricos, siendo el método SPH uno de los más utilizados [GM77,Luc77]. Descriptivamente, se trata de un método Lagrangiano en el que el fluido es modelado mediante un sistema de partículas que interactúan mutuamente, interacción que es ponderada mediante una función, monótona decreciente con la distancia, denominada función kernel. La formulación compacta que se obtiene al utilizar SPH, en la que se garantiza la conservación de la masa, junto con su adaptabilidad y versatilidad, al operar con modelados mediante partículas, ha convertido a SPH en uno de los métodos numéricos más efectivos para la simulación de fluidos. Sin embargo, a pesar de las virtudes que posee, este método presenta importantes limitaciones cuando se modela la interacción fluido–contorno [Cam89,EBT10]. La consecuencia de estas limitaciones, es la aparición de problemas de consistencia y estabilidad que se traducen en errores de simulación [EBT10,MAGC11]. Además, por lo general, sólo es posible la implementación de condiciones de contorno relativamente sencillas como la de Dirichlet [HFM12]. Analíticamente, estas limitaciones se suelen relacionar con el truncamiento de la función kernel en las proximidades del contorno [MK09,KHY13,SMZ15]. La consecuencia, es que se generan aceleraciones espurias que inducen dispersiones erróneas de las partículas [AIA∗12,KHY13,SMZ15]. Para solventar este problema, se han propuesto diferentes técnicas, cuya capacidad está condicionada por dos restricciones: la primera es la complejidad de la forma del contorno y la segunda es la imposición de las condiciones de contorno. Si la forma del contorno es angulosa se pueden producir penetraciones erróneas de las partículas dentro del contorno, lo que afecta a la simulación. Para evitar este problema se han desarrollado técnicas basadas en la ubicación de partículas en el contorno, dispuestas en un único estrato [MAC∗11], o en varios estratos dentro de él [CL03]. De este modo, se optimiza la forma del perfil del contorno y se reduce la posibilidad de que las partículas fluidas atraviesen el contorno. Sin embargo, aunque su uso permite obtener resultados realistas, el espectro de condiciones de contorno simulables es limitado. Además, en el caso de una distribución de partículas en varios estratos, se pueden obtener resultados sobrevaluados cuando el contorno es complejo, lo que afecta al realismo de la simulación [KHY13]. Para evitar estas limitaciones se han desarrollado técnicas basadas en la renormalización del contorno [FLR∗12] o en su modelado semi–analítico [MRVF13]. Con estos métodos, se mejora la simulación en contornos complejos e incluso permiten imponer condiciones de contorno complejas. No obstante, la complejidad de su formulación los hacen poco adecuado para su uso en el ámbito de Computer Graphics (CG). Considerando las ventajas de estas técnicas y con el objetivo de evitar sus inconvenientes, hemos desarrollado nuestra técnica c 2018 The Author(s) Eurographics Proceedings c 2018 The Eurographics Association. https://diglib.eg.org https://www.eg.org DOI: 10.2312/ceig.20181147 Juan J. Perea & Juan M. Cordero / Modelling the Fluid–Boundary Interaction in SPH para modelar la interacción fluido–contorno. Se trata de una técnica que aúna las ventajas del método de partículas de contorno y del método del contorno semi–analítico. Hemos seleccionado el método de partículas de contorno por las ventajas en lo referente a versatilidad y adaptabilidad a contornos complejos, como se demuestra en los estudios de [BTT09,AIA∗12]. Pero en nuestra técnica vamos a obtener una formulación en la que no es necesario introducir magnitudes dinámicas artificiosas en las partículas del contorno, ni adaptar el paso de tiempo. En lo referente a considerar el método del contorno semi–analítico, nuestro objetivo es proporcionar un método con el que imponer condiciones de contorno complejas. Condiciones que no son modeladas en los métodos partículas de contorno existentes. Además, se conseguirá un entorno que favorezca simulaciones estables y precisas, como garantizan los métodos semi–analíticos, pero con una formulación relativamente menos compleja que facilite su uso en CG. Tras esta introducción el resto del artículo se organiza del siguiente modo. En la sección 2describimos las técnicas más destacables en el modelado de la interacción fluido–contorno, nos centraremos en sus características más notables, así como en las limitaciones que puedan presentar. En la sección 3describiremos la técnica que proponemos, desglosaremos el proceso en dos etapas y evidenciaremos las ventajas que ofrece. Por otro lado, en la sección 4mostraremos los resultados que se obtienen al implementar nuestro modelo propuesto. Finalmente, en la sección 5expondremos las conclusiones que se pueden deducir de los resultados obtenidos. 2. Trabajo Relacionado Para modelar la interacción fluido–contorno, se han desarrollado diversas técnicas orientadas a mejorar el realismo de la simulación, evitar los problemas relacionados con la forma del contorno e imponer condiciones de contorno descritas mediante ecuaciones diferenciales. Entre todas ellas destacan: el método del flujo normal [ZSMC17], la renormalización del contorno [FLR∗12], el metodo de contorno semi–analítico [DMMG∗11,MRVF13], la implementación de fuerzas límites [MK09,AIA∗12], el uso de partículas fijas en el contorno [MAC∗11] y la distribución de partículas especulares o fantasmas [CL03]. El método de flujo normal, como se describe en [ZSMC17], se basa en la formulación integral de SPH para resolver el problema de truncamiento de la función kernel en las proximidades del contorno. En síntesis, esta técnica opera sobre la formulación de una integral de superficie definida sobre el contorno. De este modo, evita la ruptura de simetría de la función kernel, aunque plantea otros problemas relacionados con la resolución de la integral de superficie. En especial, cuando la geometría del contorno presenta perfiles angulosos [LVFK14]. Además, la complejidad de la formulación definida en este método, lo hace poco adecuado para su uso en CG. En este mismo contexto, es decir, evitar el problema de truncamiento de la función kernel, se ha desarrollado la técnica de renormalización del contorno [FLR∗12,HFM12,LVFK14]. En síntesis, consiste en renormalizar la función kernel considerando la posición relativa de las partículas y el contorno. Esta renormalización se puede realizar mediante métodos aproximativos, cómo propusieron DeLeffe et al. [DLTA09], o analíticos, según la descripción de Feldman y Bonet [FB07]. De este modo, la función renormalizada corrige la ausencia de partículas más allá los límites del contorno. No obstante, a pesar de que resuelve, en parte, el problema de truncamiento de la función kernel, y que la formulación obtenida es menos compleja que la del método del flujo normal, sólo es efectiva para contorno suaves. En el caso de contornos complejos, el método requiere del uso de diferentes funciones kernel normalizadas, lo que suele inducir incongruencias numéricas en los resultados, como se destaca en [MK09,LVFK14]. Desde otro punto de vista y con el objetivo de poder modelar cualquier condición de contorno se ha desarrollado el método semi–analítico [DMMG∗11,MRVF13]. En esencia, el método consiste modelar el contorno como un continuo sobre el que se define la formulación integral, no sólo de las condiciones de contorno, sino también la conservación de las magnitudes tanto dinámicas como termodinámicas. De este modo, se desarrolla un entorno matemático que permite manejar contorno complejos y cualquier condición de contorno, independientemente de su orden diferencial. En la investigación desarrollada por [MRVF13] se lleva a cabo una descripción exhaustiva de este método. La clave es la formulación integral del denominado factor de corrección definido sobre el contorno. Este factor de corrección ponderará las magnitudes en las proximidades del contorno y permitirá modelar las ecuaciones diferenciales del contorno. De este modo, se consigue mejorar la simulación del fluido y ampliar el espectro de condiciones de contorno simulables. Sin embargo, presenta el problema de la complejidad de las ecuaciones que se obtienen. Complejidad que se incrementa para contornos con formas relativamente poco complejas. La consecuencia es que, con su formulación completa, el método no es apropiado para su uso en CG. En otro orden, se han desarrollado técnicas que obvian el problema del truncamiento de la función kernel y se han centrado en modelar la dinámica específicamente en las regiones de contorno. Las más extendidas se basan en la implementación de fuerzas ficticias, cuyo módulo depende de la distancia entre las partículas y el contorno. Por depender de la distancia relativa, suelen ser formuladas a partir de potenciales conservativos, como la propuesta de Müller et al. [MST∗04] cuyas fuerzas ficticias se derivan del potencial de Lennard–Jones. La principal ventaja que ofrecen estas técnicas es su fácil implementación, con un coste computacional relativamente bajo, como se destaca en [MST∗04,Mon05,MK09]. Sin embargo, para que los resultados sean estables, se requiere de una elección específica de las constantes que controlan las fuerzas ficticias [EBT10] y de pasos de tiempo muy cortos [AIA∗12]. Además, suelen introducir un ruido numérico, no físico, en la velocidad y en la presión de las partículas, como analizaron [SCA∗11,FLR∗12]. El resultado, son dispersiones erróneas de las partículas o su separación irreal del contorno [LVFK14]. Otra limitación que presentan, es que no permiten un control específico del campo de presiones ni de velocidades, según es destacado por [HFM12,ZSMC17]. Otras técnicas alternativas, desarrolladas para evitar los problemas de las fuerzas límites, se basan en el modelado del contorno mediante partículas ficticias, bien fijando un único estrato fijo en el contorno [MAC∗11,AIA∗12], comunmente denominada partículas de contorno, bien generando diversos estratos en el contorno [CL03,HKK07], conocida como técnica de partícuc 2018 The Author(s) Eurographics Proceedings c 2018 The Eurographics Association. 12 Juan J. Perea & Juan M. Cordero / Modelling the Fluid–Boundary Interaction in SPH las especulares o fantasmas. En referencia a la técnica de partículas de contorno, su uso simplifica la detección del contorno, así como la respuesta a la interacción fluido–contorno. Sin embargo, su aplicación está restringida a condiciones de contorno relativamente sencillas [ZSMC17] y puede presentar problemas a la hora de evaluar el gradiente de presión [AIA∗12]. Aunque investigaciones recientes han intentado mejorar la versatilidad y precisión de los métodos basados en partículas de contorno. Especial mención merecen las técnicas propuestas por Cornelis et al. [CIPT15] y por Tetsuya et al. [TYTC17]. La técnica propuesta por [CIPT15] presenta dos importantes mejoras: la primera es que no requiere de un procesamiento por separado de las partículas fluidas y las de contorno, lo que mejora la eficiencia en el procesamiento, la segunda es que permite la simulación de transiciones de fase, esto es, permite la simulación de fenómenos como es el de la licuefacción. No obstante, a pesar de estas ventajas, puede presentar dos limitaciones. La primera es que se requiere del uso de mallas como paso previo a la distribución de las partículas de contorno, lo que afecta a los resultado en simulaciones con contornos complejos. El segundo es su formulación, que no permite simular condiciones de contorno complejas, ni modelar las ecuaciones de termodinámicas que rigen los fenómenos de transición. En referencia a la técnica desarrollada por [TYTC17], su propuesta distingue entre entre las partículas del fluido y las del contorno sobre, las que define funciones kernel distintas. La densidad total se calcula como una composición de las densidades asociadas a las partículas fluidas y a las de contorno. Con esta formulación es posible imponer condiciones de Dirichlet, para la superficie libre, y de Neumann para la interacción fluido–contorno. Sin embargo, a pesar de las mejoras que se obtienen, el método presenta dos limitaciones: la primera es que requiere una reformulación completa para imponer condiciones de contorno distintas de las de Neumann, la segunda es que requiere dotar de magnitudes dinámicas a las partículas de contorno, lo cual reduce la eficiencia del método. Además, sólo permite el control de la densidad y de la presión, lo cual impide la implementación de condiciones de contorno en las que interviene la velocidad. Con respecto a las partículas especulares, o fantasmas [CL03], el método consiste en generar tantos estratos, en el contorno, como sean necesarios para evitar el problema del truncamiento de la función kernel. A través de las partículas fantasmas se obtienen resultados estables, el problema es que requiere el uso de un gran número de partículas fantasmas, lo que reduce la eficiencia. Además, en contornos complejos, con bordes angulosos, se produce un solapamiento de partículas fantasma lo que induce compartamientos erráticos de las partículas [KHY13]. Para evitar este problema de solapamiento, Koschier et al. [KB17], desarrollan una técnica basada en la formulación de un mapa de densidad en las regiones del contorno. En esencia, la densidad en las proximidades del contorno, es modelada como la suma de la densidad de las partículas del fluido, establecida mediante el método SPH, y de la densidad de las del contorno, formulada mediante una función densidad continua. La integral resultante de la composición, es aproximada mediante funciones de interpolación en un proceso similar al Método de Elementos Finitos (FEM) [ZTN14]. Aunque esta propuesta representa una importante mejora respecto de los métodos previos, presenta ciertas limitaciones relacionadas con la elección de las funciones de interpolación y con la distorsión, problemas descritos en la literatura del FEM [ZTN14]. 3. Modelo Propuesto Para la simulación de la interacción fluido-contorno hemos desarrollado un modelo con el que se evitan las limitaciones de los modelos basados en partículas de contorno y la complejidad de los métodos semi–analíticos, que restringe su uso en CG. Comparte similitudes con los métodos basados en partículas de contorno como los desarrollados por [EBT10,AIA∗12]. Sin embargo, nuestra propuesta no requiere modificar el paso de tiempo, lo que representa una ventaja desde el punto de vista de la eficiencia. Además, tampoco requiere del uso de magnitudes artificiosas definidas sobre las partículas de contorno. Por otro lado, vamos a desarrollar una formulación, que describiremos en la sección 3.2, con la que imponer cualquier condición de contorno, pero evitando la complejidad de los métodos semi–analíticos. Como punto de partida, consideramos un conjunto de partículas posicionadas en el contorno, en lo sucesivo referidas como partículas de contorno, con las que interactúan las partículas fluidas que representan al fluido. A través de la cinemática de las partículas fluidas y de su posición relativa con las de contorno, conoceremos no sólo si existe interacción fluido–contorno, sino cómo se modela dicha interacción a partir de las condiciones de contorno. En el algoritmo 1describimos los pasos seguidos para detectar la interacción e imponer las condiciones de contorno. Algorithm 1 Proceso de interacción propuesto. Require: El contorno esté modelado mediante partículas de contorno. Ensure: Las partículas fluidas cumplen las condiciones de contorno. Input: Lista partículas fluidas, Lista de partículas de contorno. 1: for Cada partícula fluidas do 2: Determinar las partículas fluidas vecinas 3: Determinar las partículas de contorno vecinas 4: Calcular dinámica mediante SPH según algoritmo 2 5: if 6 ∃ interacción con partícula de contorno then 6: Integrar nueva posición. 7: Detectar colisiones. 8: if ∃colisión then 9: Determinar las partículas de contorno con la que interactúa. 10: Ampliar la lista de partículas de contorno vecinas. 11: Modificar dinámica de partículas fluidas según algoritmo 3 12: end if 13: end if 14: Integrar nueva posición, (ec. 16) 15: end for Para implementar el proceso descrito en el algoritmo 1, vamos a distinguir dos etapas: la etapa de interacción y la etapa de restricción. En la etapa de interacción, evaluamos si existe interacción fluido–contorno. En la etapa de restricción, imponemos las condiciones de contorno sobre las partículas fluidas que interactúan con c 2018 The Author(s) Eurographics Proceedings c 2018 The Eurographics Association. 13 Juan J. Perea & Juan M. Cordero / Modelling the Fluid–Boundary Interaction in SPH i (a) i (b) i (c) Figura 1: Casuística en la etapa de interacción. En todas ellas, las partículas de contorno vecinas están coloreadas en verde y la partícula de interacción se muestra con núcleo blanco. La partícula roja, Figuras 1b y1c, representa la colisión de la partícula fluida con el contorno, lo que requiere el uso del proceso alternativo. Algorithm 2 Ensure: La dinámica de la partícula fluida queda actualizada. Ensure: Si existe interacción Input: Partícula fluida, lista de partículas fluidas vecinas, lista de partículas de contorno vecinas. 1: Calcular dinámica mediante SPH: ρi,vi,Pi 2: for Cada partícula de contorno vecina do 3: Calcular orientación relativa, parámetro θk(ec. 1). 4: if ∃θk≃1then 5: Se detecta interacción. 6: Calcular velocidad restringida (ec. 9). 7: Calcular la tasa de avance (ec. 17). 8: Calcular presión restringida (ec. 20). 9: Calcular densidad restringida (ec. 22). 10: end if 11: end for Algorithm 3 Ensure: La dinámica de la partícula fluida queda actualizada. Input: Partícula fluida, lista de partículas fluidas vecinas, lista de partículas de contorno vecinas. 1: for Cada partícula de contorno vecina do 2: Calcular término de deslizamiento (ec. 8). 3: Calcular la tasa de avance (ec. 17). 4: Calcular presión restringida (ec. 20). 5: Calcular densidad restringida (ec. 22). 6: end for el contorno. Tanto el proceso seguido para detectar la interacción, como el modelado de las ecuaciones de contorno para imponer las restricciones, se describirán a continuación. 3.1. Etapa de Interacción Para detectar si existe interacción fluido–contorno, consideramos la posición relativa entre las partículas fluidas y las de contorno, así como la velocidad de las partículas fluidas. Con la primera detectamos la proximidad entre ambas, con la segunda si esa proximidad es compatible con la condición de interacción. De modo análogo a otros estudios [BTT09,EBT10], vamos a considerar que el contorno se encuentra inmóvil. Ello no representa una pérdida de generalidad, ya que las ecuaciones que se obtendrán serán válidas para contornos que se desplazan, simplemente considerando la velocidad relativa, en lugar de la velocidad de las partículas fluidas exclusivamente. Dependiendo de la velocidad y de la posición relativa entre las partículas fluidas y las de contorno, pueden darse tres situaciones: la primera, es cuando la partícula fluida tiene partículas de contorno vecinas e interactúa con algunas de ellas, Figura 1a, la segunda es cuando, aún teniendo partículas de contorno vecinas, no interactúa con ninguna de ellas, Figura 1b, y la tercera es cuando no tiene partículas de contorno vecinas, pero la partícula fluida acaba colisionando con el contorno, Figura 1c. La condición necesaria para que una partícula fluida interactúe con el contorno, es que su velocidad y su posición relativa, respecto a las de contorno, sean paralelas. Cuantitativamente, esta condición la definimos mediante la ecuación 1: θk=~vi·(~rck−~ri) |~rck−~ri||~vi|≃1 (1) donde θklo denominamos parámetro de orientación,~vies la velocidad de la partícula fluida i,~ries su posición y~rckes la posición de la partícula de contorno k. En la interacción fluido–contorno, el caso más habitual es el primero (Figura 1a) esto es, la partícula fluida tiene partículas de contorno dentro de su dominio soportado con las que interactúa. En esta situación, todas las partículas de contorno vecinas participarán a la hora de imponer la condición de contorno. En este caso, el proceso seguido para la calcular la dinámica de las partículas fluidas es el descrito entre los puntos 4y10 del algoritmo 2. Para las otras dos situaciones en las que existe interacción, pero ésta no se detecta; bien porque la partícula fluida tiene partículas de contorno vecinas, pero con ninguna se cumple la condición 1 (Figura 1b), o bien porque la partícula fluida no tiene partículas de contorno vecinas (Figura 1c); tenemos que utilizar un proceso alternativo ya que la partícula fluida no “conocerá” la ubicación del contorno y acabará colisionando. Este proceso alternativo consiste en localizar entre todas las partículas de contorno, cuales satisfacen la condición 1. Una vez encontradas, se las consideran partículas de contorno vecinas aunque no se encuentre dentro del dominio soportado de la partícula fluida. El proceso seguido para calcular la dinámica de estas partículas fluidas es el descrito entre los puntos 5y13 del algoritmo 1. c 2018 The Author(s) Eurographics Proceedings c 2018 The Eurographics Association. 14 Juan J. Perea & Juan M. Cordero / Modelling the Fluid–Boundary Interaction in SPH Tras la etapa de interacción, conocemos cuales son las partículas de contorno vecinas de las partículas fluidas. A partir de esta relación, desarrollaremos la etapa de restricción que describiremos en la sección 3.2. 3.2. Etapa de Restricción Una vez determinado que el fluido interactúa con el contorno, el siguiente paso es imponer las condiciones de contorno que restringen su dinámica. Nuestra propuesta tiene tres objetivos principales: 1. Desarrollar una formulación general, con la que se pueda imponer cualquier condición de contorno, independientemente de su complejidad. 2. Satisfacer la conservación del momento. 3. Optimizar las ecuaciones obtenidas para su uso en modelados de contornos mediante partículas. A lo largo de la descripción, vamos a tener en cuenta que el contorno es rígido. Como punto de partida, consideramos la ecuación que relaciona la velocidad antes y después de la interacción, ecuación 2, y, a partir de ella, vamos a obtener la ecuación propuesta para calcular la velocidad restringida. ~vi 0=−α~vni+β~vτi(2) donde el superíndice prima hace referencia a la velocidad restringida por la condición de contorno, αyβson parámetros de interacción y los subíndices nyτhacen referencia a las direcciones normal y tangencial de la superficie del contorno. Cabe destacar que el signo negativo con el que aparece α, es el modo habitual de compactar la condición cinemática normal y la impenetrabilidad del contorno. El primer problema que plantea la ecuación 2es que se formula sobre la premisa de una dirección normal y una tangencial, y en nuestro modelo basado en partículas de contorno, es difícil establecer una dirección normal y una tangencial, ya que de ellas sólo conocemos su posición. Diversas propuestas, como las desarrolladas por [BTT09,AIA∗12], asignan las direcciones normal y tangencial partir de la discretización triangular del contorno. El problema que ello plantea, es que se requiere de esa malla de triángulos. Además, el método pierde capacidad ya que las partículas se adaptan mejor al perfil del contorno. Para evitar el uso de direcciones normal y tangencial, vamos a reformular la ecuación 2de modo que sólo dependa de la velocidad de las partículas fluidas y de la posición relativa entre éstas y las de contorno. Para ello, en primer lugar, eliminamos la dependencia, en la ecuación 2, de la componente normal, con lo que obtenemos: ~vi 0=−α~vi+ (α+β)~vτi(3) Aunque en la ecuación 3se ha eliminado la componente normal, todavía aparece la tangencial. Así pues, el siguiente paso es expresar la componente tangencial, también llamada de deslizamiento, en función de la velocidad ~viy de las posiciones relativas. Para llevar a cabo esta tarea, vamos a considerar la formulación establecida por SPH. No obstante, la formulación estándar de SPH para la velocidad, viola la conservación del momento [Mon92]. Para evitar este problema, Monaghan [Mon92] desarrolla un método basado en relaciones de identidad asociadas al cálculo diferencial. Con esta formulación, se pueden obtener ecuaciones simétricas en las que intervienen los valores relativos de las magnitudes dinámicas. Para ello, la ecuación general de partida es: ∇f(r) = 1 ρ(∇(ρf(r))−f(r)) Particularizando para el caso de la velocidad se obtiene: ∇~vi=1 ρi∑ j~vj−~vi∇W(~rj−~ri,h) A partir de esta ecuación, se puede obtener una expresión simétrica en la que también intervengan las densidades relativas. La expresión que se obtiene es: ∇~vi=ρi∑ j ~vj ρ2 j −~vi ρ2 i!∇W(~rj−~ri,h)(4) donde ρi,ρj,~viy~vjson las densidades y velocidades de las partículas fluidas iyjrespectivamente. La ecuación 4, permitiría calcular el gradiente de la velocidad. Sin embargo, nuestro objetivo no es el gradiente sino la velocidad. Para ello, utilizaremos la formulación general de SPH [Mon92]. Según esta formulación, la derivada de cualquier magnitud se puede aproximar cómo: ∇lf(ri) = ∑ j mj fj ρj ∇lW(rj−ri,h),(5) donde lhace referencia al orden diferencial. A partir de la ecuación 5, para el caso de l=0, la ecuación 4 puede expresarse en términos de términos: ~vi=ρi∑ j mj ~vj ρ2 j −~vi ρ2 i!W(~rj−~ri,h)(6) Con lo que hemos obtenido una expresión simétrica que satisfará la conservación del momento. Seguidamente, a partir de 6, formulamos nuestra propuesta para la velocidad de deslizamiento,~vτi. Dicha ecuación es: ~vτi=mi 2ρi∑ s (~vi−~vs)W(~rs−~ri,h)(7) donde ~vsjugaría el papel de la velocidad de deslizamiento de las partículas de contorno. Sin embargo, esta consideración es artificiosa, ya que las partículas de contorno no pueden tener una velocidad de deslizamiento sobre sí mismas. Para evitar esta incongruencia vamos a a modificar la ecuación 7, de modo que no intervenga~vs. Esto lo logramos utilizando la relación vectorial entre c 2018 The Author(s) Eurographics Proceedings c 2018 The Eurographics Association. 15 Juan J. Perea & Juan M. Cordero / Modelling the Fluid–Boundary Interaction in SPH velocidades, ~vs=~v−~vn, y considerando que ~vn≃proy~vi(~rck−~ri). La expresión obtenida es: ~vτi=mi 2ρi∑ kvi−~vi·(~rck−~ri) |~rck−~ri|(~rck−~ri)W(~rck−~ri,h)(8) donde hemos conseguido que la velocidad de deslizamiento sólo dependa de la velocidad de las partículas fluidas y de sus posiciones relativas respecto del contorno. Según la ecuación 8, esta componente de la velocidad será máxima cuando posición relativa partículas fluida–contorno y velocidad sean paralelas, lo cual es coherente con la condición impuesta por 1. Finalmente, sustituyendo las ecuaciones 6y8en la expresión 2, obtenemos la formulación para la velocidad restringida por el contorno, en la que no interviene, explícitamente, ni la dirección normal ni la tangencial. ~vi 0=−αρi∑ j mj ~vj ρ2 j −~vi ρ2 i!W(~rj−~ri,h) +(α+β)mi 2ρi∑ k~vi−~vi·(~rck−~ri) |~rck−~ri|(~rck−~ri)W(~rck−~ri,h)(9) Con la ecuación 9, hemos logrado dos de los tres objetivos iniciales, el primero es la conservación del momento y el segundo es que que sólo dependa de la velocidad de la partícula fluida y de la posición relativa entre las partículas fluidas y las de contorno. Además, evitamos asignar magnitudes dinámicas artificiales sobre las partículas de contorno, como se opera en metodos de similares característica [AIA∗12,ZSMC17]. También, aseguramos que cuanto más paralela sea la velocidad de la partícula fluida a la posición relativa,~rck−~ri, mayor será la participación del contorno sobre el fluido, lo cual es coherente con las restricciones dinámicas asociadas a la interacción fluido–contorno, según se destaca en [ZTN14]. 3.2.1. Imposición de la condición de contorno Conseguido dos de los objetivos planteados, queda generalizar la ecuación 9para poder imponer condiciones de contorno más complejas. Nuestro objetivo es establecer una relación entre los parámetros αyβde la ecuación 9, que habitualmente tienen valores constantes, y las ecuaciones en derivadas parciales con las que se describen las condiciones de contorno. En nuestra descripción, nos vamos a centrar en el parámetro α, que en dinámica de fluidos suele relacionarse con la componente normal y ésta es relevante en la interacción con el contorno [ZTN14]. No obstante, el proceso que describiremos es extrapolable para el parámetro β, simplemente considerando la dirección tangencial. La condición de Robin es la formulación más general de la condición de contorno [ZTN14]. Se trata de la composición de las condiciones de contorno de Dirichlet y de Neumann. En términos generales, la condición de Robin se formula como: γ∂l ∂χlζ+σζ =ζδΩ (10) donde γ,σson contantes definidas en la formulación, ζes la magnitud dinámica que debe satisfacer la condición de contorno, ζδΩ es el valor de dicha magnitud en el contorno δΩ,les el orden diferencial, que siempre debe ser inferior al de la ecuación dinámica y χes la dirección del espacio donde se aplica la restricción. Analíticamente, para imponer esta restricción, se sustituye la solución general, asociada al flujo libre, en las ecuaciones de contorno y se resuelven. La solución obtenida, satisfará las condiciones de contorno [ZTN14]. Tomando como punto de partida este proceso, desarrollamos nuestra propuesta. Para ello, fijamos la hipótesis de que las magnitudes dinámicas son multiplicadas por un parámetro, que posteriormente denominaremos cómo parámetro de interacción, a través del cual, la condición de contorno es modelada. Esta hipótesis no es restrictiva, ya que el parámetro tendrá un valor de 1 para las partículas fluidas que no interactúan con el contorno y distinto de 1 en caso contrario. Además, esta hipótesis está acorde con los fundamentos de la dinámica de fluidos [ZTN14]. En base a esta hipótesis, la velocidad de la partícula fluida quedaría definida cómo: ~v0 i=AδΩ~vi(11) donde~v0 isería la velocidad restringida por el contorno,~vies la velocidad de flujo libre y AδΩ es el parámetro de interacción. Este parámetro será el responsable de “modificar” la velocidad del fluido libre para que se satisfaga la condición de contorno. Cómo ~vies la velocidad de flujo libre, esto es, la velocidad obtenida directamente de la ecuación dinámica, podemos de nuevo utilizar la formulación de la velocidad obtenida a partir de la propuesta por Monaghan [Mon92]. Según esta afirmación, la velocidad restringida se definiría a partir de la 6. Así, la ecuación 11 se expresa en términos de: ~v0 i=AδΩρi∑ j mj ~vj ρ2 j −~vi ρ2 i!W(~rj−~ri,h)(12) El siguiente paso es relacionar la ecuación 12 con la formulación de la condición de contorno, ecuación 10. Dicha ecuación, particularizada para el campo de velocidades, se expresa en términos de la ecuación γ∇n~v0 i+σ~v0 i=~viδΩ (13) donde ∇nse refiere a la componente normal del gradiente y~viδΩ es el valor de la velocidad en el contorno. c 2018 The Author(s) Eurographics Proceedings c 2018 The Eurographics Association. 16 Juan J. Perea & Juan M. Cordero / Modelling the Fluid–Boundary Interaction in SPH Así pues, al sustituir la ecuación 12, en la ecuación 13 se tiene: AδΩρi∑ j mj~vrel [γ∇rck riW(~rck−~ri,h)+σW(~rck−~ri,h)] =~viδΩ (14) donde ~vrel hace referencia al término de velocidad relativa de la ecuación 12, y donde hemos consideramos que ∇n≃ ∇rck ridonde ∇rck ries la proyección del operador ∇en la dirección~rck−~ri. Despejando el parámetro de interacción, AδΩ, de la ecuación 14, obtenemos: Avi δΩ =~viδΩ  ρi∑jmj|~vrel|(γ|∇rck riW(~rck−~ri,h)|+σW(~rck−~ri,h)) (15) donde se ha considerado el módulo de las magnitudes vectoriales, ya que Avi δΩ es una magnitud escalar. Sustituyendo la ecuación 15, en la ecuación 12, y comparando con la ecuación 9, se puede concluir que AδΩ ≡α. Así somos capaces de imponer cualquier condición de contorno, sin que se viole la conservación del momento y evitando las complejidades de los métodos semi–analíticos. 3.2.2. Posición restringida Conocida la velocidad que cumple la condición de contorno, el siguiente paso es conocer la posición restringida de la partícula. La hipótesis fundamental que consideramos es que, independiente del momento en el que se produzca la colisión, dentro del paso de tiempo, la colisión se producirá al final del paso de tiempo. No obstante, vamos a evitar que el tiempo aparezca explícitamente en nuestra ecuación. Por otro lado, en la ecuación de la posición que proponemos evitamos que dependa de las direcciones normal o tangencial, ya que, como mencionamos con anterioridad, el objetivo en toda la formulación de nuestro modelo, es que sólo intervengan las posiciones relativas entre partículas de contorno y las fluidas y la velocidad de éstas. Con estas premisas, la ecuación que formulamos es: ~ri 0=~ri+ˆ ~rc0 ki(1−λ)(16) donde ˆ ~rc0 ki, que designamos como tasa de avance, satisface: ˆ ~rc0 ki=1 j0 ∑ ck0~rck0−~ri!(17) donde el subíndice ck0hace referencia a las partículas de contorno que satisfacen la condición 1,j0es el número de partículas de contornos vecinas de la partícula fluida ique cumplen dicha condición yλ, que denominamos tasa de interacción, es λ=Rf+Rc+δcon RfyRcel radio de la partícula fluida y de la partícula de contorno respectivamente, y δel espesor de contacto que permitirá controlar la aproximación al contorno de la partícula fluida. Por lo general, el valor de esta magnitud es inferior al 10% del radio de la partícula fluida. 3.2.3. Presión y densidad restringida Obtenida la velocidad restringida por la condición de contorno, según ecuaciones 9y15, es necesario calcular la nueva presión y densidad. Ambas magnitudes se calcularán a partir de la velocidad. Para el caso de la presión, nos basaremos en la hipótesis fundamental sobre la que se desarrolla el método PCISPH, ideada por Solenthaler et al. [SP09]. Por otro lado, para calcular la densidad nos basaremos en la ecuación de Venturi, en la que se relacionan las magnitudes presión, velocidad y densidad. En primer lugar, vamos a formular la ecuación para calcular la presión asociada a la velocidad restringida. Para ello consideramos la hipótesis establecida por Solenthaler et al. [SP09] en su modelo PCISPH. Según esta hipótesis, en el cambio de la posición de la partícula fluida, la fuerza más relevante es la fuerza de presión, ~ Fp i, por lo que el resto de fuerzas pueden despreciarse. Cuantitativamente, esta afirmación indica que: ∆~ri=∆t2~ Fp i mi (18) Expresando la ecuación 18, en términos de la velocidad, obtenemos: ∆~vi=∆t ~ Fp i mi (19) donde ∆~vies el cambio de la velocidad debido a la interacción con el contorno, ∆tes el paso de tiempo de integración y ~ Fp ies la fuerza de presión. Para calcular la presión, es necesario relacionarla con la fuerza de presión. Para ello, volvemos a considerar la hipótesis de Solenthaler et al. [SP09], según la cual, la fuerza de presión de la partícula fluida y la de sus vecinas es muy similar, con lo que, junto con la formulación conservativa de la fuerza de presión, desarrollada por Monaghan [Mon92], permite formular: |∆~vi| ∆t=mi 2P0 i ρ2 0 ∑ j ∇W(~rj−~ri,h) ⇒ P0 i=ρ2 0|∆~vi| 2mi∑j∇W(~rj−~ri,h) (20) Finalmente, conocida la velocidad,~v0 i, y la presión, P0 i, restringidas, queda conocer la nueva densidad ρ0 i. El punto de partida que vamos a considerar es la ecuación de Venturi. En síntesis, esta ecuación permite afirmar que la velocidad, la presión y densidad en un fluido están relacionadas, de modo que un cambio en alguna de ellas, en dos regiones distintas del fluido, afecta a las otras. En forma compacta, la ecuación de Venturi es: Pi ρig+~ri+v2 i g=P0 i ρ0 ig+~r0 i+v0 i2 g(21) c 2018 The Author(s) Eurographics Proceedings c 2018 The Eurographics Association. 17 Juan J. Perea & Juan M. Cordero / Modelling the Fluid–Boundary Interaction in SPH Operando sobre la ecuación 21, obtenemos la densidad restringida. La ecuación obtenida es: ρ0 i=ρiP0 i Pi+ρig+∆~ri+∆v2 i(22) Así pues, con las ecuaciones 9,16,20 y22 podemos calcular las magnitudes del fluido al interactuar con el contorno. 4. Resultados y Discusión En esta sección, vamos llevar a cabo la implementación de nuestra propuesta junto con otras tres técnicas basadas en la discretización del contorno mediante partículas. De este modo, será posible desarrollar un estudio comparativo entre cada una de ellas. Así, se podrá visualizar las capacidades que nuestra propuesta ofrece. Las tres técnicas seleccionadas son: la desarrollada por Oger et al. [ODAF06], la propuesta por Becker et al. [BTT09] y la técnica de Akinci et al. [AIA∗12]. Hemos seleccionado estas técnicas porque presentan ciertas semejanzas con nuestro modelo, como son el tratamiento del contorno mediante partículas y la detección de la interacción fluido-contorno a través de técnicas geométricas. 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 [BTT09] o la función kernel utilizada por [AIA∗12] para calcular la densidad de las partículas de contorno. Para llevar a cabo la implementación de estas diferencias, hemos utilizado los valores y recomendaciones dada por los autores de cada técnica. Para las simulaciones hemos utilizado un procesador Intel Core i3-4130, con 16 GB de memoria RAM y una tarjeta gráfica NVIDIA GTX780. En cada implementación, vamos a considerar un fluido incompresible con una tasa de incompresibilidad de η= 1% y con una densidad estática ρ0=1000. Para ello vamos a utilizar la propuesta desarrollada por [PC16], que evita los problemas de inestabilidad asociados al gradiente de presión y permite imponer la condición de incompresibilidad. Como método de integración temporal hemos seleccionado un método semi–implícito desarrollado por Modak y Sotelino [MS02] que ofrece un orden de precisión O(n5), sin un aumento significativo en los costes de computación. Los valores de los parámetros de integración que intervienen en el método son β0=1, β1=0.757, β2=1.101. El paso de tiempo que utilizaremos es ∆t=0.025. En cuanto a la discretización de nuestro contorno, hemos considerado una distribución partículas, cuyo radio es el doble del que poseen las partículas fluidas. Además, hemos impuesto la condición de que estén solapadas un 5%, es decir, que la distancia entre el centro de dos partículas de contorno contiguas es 1.9Rf, con Rf el radio de la partícula fluida. Así, conseguimos un recubrimiento completo del contorno. En referencia a la busqueda de partículas vecinas, tanto fluidas como de contorno, hemos utilizado la técnica basada en el etiquetado mediante función hash y de segmentación del espacio desarrollada por [PC17]. Para la simulación hemos utilizado 50k partículas fluidas todas ellas de igual masa m=0.03, en cada una de las cuales hemos utilizado una spline cúbico [Mon05] como función kernel, ya que diversos estudios [BTT09,AIA∗12] demuestran su idoneidad en la simulación de la interacción fluido–contorno. El radio de dominio soportado para esa función kernel es h=0.38. Con estas condiciones hemos obtenido los resultados, para cada una de las técnicas implementadas, mostrados en la Figura 2, 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 2d, ofrece unos buenos resultados desde el punto de vista de estabilidada y adaptabilidad al contorno. Analizando cada imagen, se puede afirmar que el resultado obtenido a partir de la técnica propuesta por [ODAF06], Figura 2a, 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 sobre–amortiguamiento artificial, la segunda es que se produce una ligera penetración en el contorno. Ello se debe a que el modelo de [ODAF06] 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. [BTT09], Figura 2b, se solventa este problema, es decir, no se muestra sobre–amortiguamiento, 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 a la dispersion. Con respecto al resultado obtenido a partir de la técnica propuesta por Akinci et al. [AIA∗12], Figura 2c, 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 2d. 5. Conclusiones En este artículo hemos llevado a cabo una revisión de las técnicas más importantes en la simulación de la interacción fluido–contorno. Nos hemos centrado en las características que influyen en la versatilidad y eficiencia. A partir de estas características, hemos propues- to una técnica que aúna las virtudes de dos de los mejores métodos, estos son, el método de partículas de contorno y el método de contorno semi–analítico, pero evitando sus limitaciones. Particularmente hemos logrado: •Un método optimizado para contornos modelados mediante partículas. Esta afirmación queda confirmada por las ecuaciones 9,20 y22 que sólo dependen de la velocidad de la partícula fluida y de la posición relativa entre los dos tipos de partículas. •Obtener un conjunto de ecuaciones que permiten calcular no sólo la velocidad restringida por el contorno, sino también las otras magnitudes que intervienen en la dinámica del fluido, estas son: la presión, ecuación 20, y la densidad, ecuación 22. Especialmente destacable es la ecuación 16, con la que se obtiene la posición, ya que permite controlar, eficientemente, la posición gracias a la tasa de interacción. c 2018 The Author(s) Eurographics Proceedings c 2018 The Eurographics Association. 18 Juan J. Perea & Juan M. Cordero / Modelling the Fluid–Boundary Interaction in SPH (a) (b) (c) (d) Figura 2: Comparativa de los resultados obtenidos con diferentes técnicas basadas en partículas de contorno. En la Figura 2a, se muestra el resultado de la técnica propuesta por [BTT09]. En la Figura 2b, el obtenido a partir del modelo de [ODAF06]. En la Figura 2c, el de la propuesta de [AIA∗12]. Finalmente, en la Figura 2d, 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. •Desarrollar un método que mejora la respuesta a la interacción fluido–contorno. Esta afirmación esta fundamentada en las imágenes de la Figura 2, donde los mejores resultados se obtienen a partir de nuestra técnica. •Un método que no requiere adaptar el paso de tiempo, lo que favorece su eficiencia ni del uso de procesos predicción– corrección. •Un método estable que permite obtener simulaciones realistas con un alto control de la interacción fluido–contorno, como se puede deducir de la Figura 2. 6. Agradecimientos Esta investigación ha sido apoyada por el proyecto Pololas (TIN2016-76953-C3-2-R) del Ministerio de Economía y Competitividad de España. References [AIA∗12] AKINCI N., IHMSEN M., AKINCI G., SOLENTHALER B., TESCHNER M.: Versatile rigid–fluid coupling for incompressible sph. ACM Trans. Graph. 31, 4 (Jul 2012), 62:1–62:8. doi:10.1145/ 2185520.2185558.1,2,3,5,6,8,9 [BTT09] BECKER M., TESSENDORF H., TESCHNER M.: Direct forcing for lagrangian rigid–fluid coupling. IEEE Transactions on Visualization and Computer Graphics 15, 3 (May 2009), 493–503. doi:10.1109/ TVCG.2008.107.2,4,5,8,9 c 2018 The Author(s) Eurographics Proceedings c 2018 The Eurographics Association. 19