More Stable SPH-based Fluid Simulation

CEIG - Spanish Computer Graphics Conference (2012) Isabel Navazo and Gustavo Patow (Editors) More Stable SPH-based Fluid Simulation J.J. Perea and J.M. Cordero University of Seville, Spain Abstract The simulation of fluid is a field largely discussed, from the viewpoint of engineering, in the so-called computational fluid dynamics (CFD). A problem that occurs repeatedly in the literature of CFD is the stability of the simulation. One of the elements that introduce instability in any CFD model is the temporal integration method. Among the most commonly used CFD techniques, both in the field of engineering and in the computer graphics, is SPH, which is particularly relevant to its computational simplicity. However, this technique adds stability problems due to the treatment of viscous and pressure terms arising from the formulation of Navier–Stokes. SPH also introduces greater instability by the need to use interpolation functions appropriate to the phenomena that occur in fluid dynamics. We present solutions to these four problems of inestability mentioned: treatment of the terms of pressure and viscosity that appear in the Navier–Stokes, the formulation of the characteristic interpolation function of SPH and treatment of temporal integration. Categories and Subject Descriptors (according to ACM CCS): I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling—Physically based modeling; I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—Animation 1. Introducción La animación del flujo de un fluido sigue siendo actualmente un campo de investigación abierto. Desde un punto de vista mecánico los fluidos son objeto de análisis general por la Mecánica de Fluidos y, desde un punto de vista numérico, son objeto de estudio particular por la Dinámica de Fluidos Computacional (Computational Fluid Dinamics, CFD), que describe métodos numéricos para resolver y analizar los problemas implicados en el flujo de un fluido. El comportamiento dinámico del flujo de un fluido se describe comunmente a través de la ecuación de Navier– Stokes, que caracterizada para un fluido incompresible e isotérmico podemos expresarla en términos de la formulación Larangiana como la ecuación: ρDv Dt =−∇P+ρg+µ∇2v,(1) En esta ecuación ρes la densidad de masa, el término D/Dt es la derivada sustancial, ves la velocidad, Pes la presión, ges la aceleración de la gravedad y µes la constante de viscosidad del fluido. El término −∇P, denominado de presión, es la fuerza que surge en el fluido como oposición a la variación de presión. En consecuencia, en zonas del fluido donde la presión es constante no aparecerán fuerzas asociadas al término de presión. El término ρges la fuerza que surge por estar sometido el fluido a la gravedad y µ∇2ves el denominado término viscoso, de naturaleza disipativa y que depende, a través µ, de las propiedades constitutivas del fluido. No es posible obtener, de manera general, una solución analítica a la ecuación diferencial (1) de Navier–Stokes, por lo que en CFD se utilizan métodos numéricos. Una de las técnicas numéricas más utilizadas en CFD y en animación por ordenador (Computer Graphics, CG) es Smoothed Particle Hydrodynamics (SPH) [Mon05]. SPH transforma la ecuación de Navier–Stokes, definida en un dominio continuo, en un sistema de ecuaciones diferenciales en el tiempo, resoluble a través de métodos basados en el desarrollo en serie de Taylor. La simulación de fluidos con SPH presenta dos grupos de problemas [TY09]: los relativos a la incompresibilidad del fluido y los relativos a la estabilidad. En este trabajo sólo tratamos de resolver algunos de los problemas relacionados con la estabilidad en la simulación de fluidos basada en el método SPH. Algunas de las causas de inestabilidad en SPH son: Fluctuaciones en el flujo del fluido debido al término de presión [Mon92,DC96] y al término de viscosic The Eurographics Association 2012. DOI: 10.2312/LocalChapterEvents/CEIG/CEIG12/011-020 J.J. Perea & J.M. Cordero / More Stable SPH-based Fluid Simulation dad [Mon92], y determinación incorrecta de la función de interpolación [Mor95,SHS95]. En este trabajo presentamos propuestas para aumentar la estabilidad de la simulación mediante SPH para los problemas mencionados: 1. En la sección 4definimos un nuevo término de presión para mejorar la estabilidad en las regiones isobáricas. 2. En la sección 5proponemos un término de amortiguamiento para mejorar la estabilidad generada por el término de viscosidad. 3. En la sección 6presentamos una nueva función de interpolación. 4. En la sección 7describimos un método de integración temporal, con un orden de error adecuado a las necesidades de estabilidad de SPH. 2. Trabajo Relacionado El uso de SPH en su formulación clásica [GM77] plantea algunos problemas desde el punto de vista de la estabilidad [AS05], que aparecen reflejados como oscilaciones irreales de las partículas durante la simulación. Podemos agrupar estos problemas en tres fuentes de inestabilidad: el término de presión, el término de viscosidad y la función de interpolación. Para solventar los problemas de estabilidad asociados al término de presión, algunos autores proponen introducir un término de amortiguamiento viscoso ficticio. Uno de los primeros en introducir dicho término fue Monaghan [Mon92], haciéndolo depender de la velocidad, el tipo y de las dimensiones características del fluido, entre otras magnitudes. Otros autores, como Fangming et al. [FOC07], calculan el término de amortiguamiento a partir de un potencial escalar de tipo Lenard-Jones. El término obtenido es relativamente sencillo y favorece la estabilidad, si bien requiere de una distribución simétrica de partículas. En caso contrario se produce un agrupamiento irreal de partículas y un aumento en su comportamiento errático, por lo que el sistema se vuelve inestable. Toro [Tor01] afronta este problema de estabilidad estableciendo una formulación variacional del término de presión, de este modo no es necesario considerar el término de amortiguamiento de Monaghan. Con este método se consigue reducir las oscilaciones, aunque las expresiones obtenidas son complejas de resolver. Riadh et al. [AS05], en una linea intermedia, utiliza el método de resolución de Riemann, para calcular un término de amortiguamiento más sencillo que el de Monaghan con una menor dependencia de los parámetros, con lo que se gana versatilidad y estabilidad. En una línea de investigación similar Vacondio et al. [VRSM12] vuelve a utilizar el método de resolución de Riemann, pero imponiendo el flujo de Lax–Friedrichs. La expresión que obtiene garantiza un uso aún más generalizado del término de amortiguamiento y una mayor estabilidad. Sigalotti et al. [SLT09] propone el uso de un método SPH adaptativo para calcular la presión en las zonas de contorno. El problema que plantea esta técnica es que depende de su forma, con lo que el proceso se complica en simulaciones con contornos complejos. Jiang et al. [JOR∗12] corrige la formulación fundamental de SPH, para ello efectúa un desarrollo en serie Taylor de dicha formulación, a partir de la cual, rectifica la expresión de las magnitudes dinámicas que intervienen en la ecuación de Navies–Stokes. De este modo evita el uso del término de amortiguamiento propuesto por Monaghan [Mon92]. Hashemi et al. [HFM12] identifica que el problema de la inestabilidad se encuentra en el agrupamiento erróneo de las partículas producido en las zonas del contorno y lo rectifica mediante un proceso predictor– corrector. Para corregir las posiciones de las partículas agrupadas, emplea un desarrollo en serie Taylor de las magnitudes que intervienen en la ecuación de Navier–Stokes. Para evitar la inestabilidad en las zonas del contorno otros autores como Safdari et al. [SZSY11] y Schechter et al. [SB12] proponen el uso de partículas virtuales o fantasmas, situadas en las proximidades de las zonas de interacción sólido-líquido. De este modo se inducen un efecto de campo local en la presión, favoreciendo la estabilidad. Una línea diferente es la seguida por Bodin et al. [BLS11] que utiliza formulación de mecánica analítica para evitar los problemas de inestabilidad en el contorno. Para ello, considera que las interacciones que presenta el fluido en el contorno son formuladas, en la ecuación dinámica, como ligaduras holónomas y las fuerzas de ligadura, mediante multiplicadores de Lagrange. El método obtenido es de gran precisión y reduce la inestabilidad, si bien el uso de multiplicadores de Lagrange aumenta la complejidad y exige del uso de metodología implícita, con la consecuente pérdida de precisión debida a la viscosidad numérica [TY09]. Según Müller [MCG03] el término de viscosidad produce inestabilidad. Para solucionarlo formula funciones de interpolación ad hoc [MCG03,MSKG05]. Otra formulación fue presentada por Amini et al. [AEF11], que consiste en una reformulación del término viscoso, en el que se fijan parámetros dependientes de las propiedades constitutivas del fluido. Con esta formulación no se hace necesario el uso de una función específica, como ocurría con la propuesta de Müller. Sin embargo, el estudio de Amini presenta dos incovenientes: el primero la elección precisa de los valores de los parámetros, de modo que diferentes condiciones en el fluido pueden afectar a la estabilidad de la simulación; el segundo es que la expresión obtenida no siempre satisface la conservación del momento lineal. Una elección arbitraria de la función interpolación produce inestabilidades [Mon92,LXSR12]. Diversos estudios se centran en las condiciones que deben cumplir las funciones de interpolación para que se reduzca la inestabilidad. Swegle [SHS95] y Fulk [Ful94] se centran en determinar las condiciones de continuidad que tienen que satisfacer dicha función para que garantice estabilidad. Otros autores, como Liu et al. [LLL03] y Müller et al. [MSKG05], proponen el c The Eurographics Association 2012. 12 J.J. Perea & J.M. Cordero / More Stable SPH-based Fluid Simulation uso de funciones polinímicas, y fijan los criterios que deben cumplir sus derivadas para garantizar la estabilidad. Concluyen que es necesario el uso de distintas funciones de interpolación para cada uno de los términos de la ecuación de Navier–Stokes. En esta misma línea se sitúan los estudios de Solenthaler et al. [SBC∗11] y Lind et al. [LXSR12], si bien la elección ad hoc que ellos proponen no siempre opera de manera adecuada en problemas distintos para los que han sido diseñado. 3. El Método SPH En esta sección vamos a definir los términos relacionados con SPH que utilizaremos en secciones posteriores. El método SPH trata el medio continuo como un sistema de partículas que interaccionan. Estas partículas representan fracciones discretas del volumen del fluido. Es por tanto, una descripción Lagrangiana del medio. Analíticamente, el método SPH se basa en la expresión integral de cualquier función real obtenida a partir la función δde Dirac [Mon92]. Una generalización de dicha formulación sustituye a la función δde Dirac por una función de interpolación. Cualquier función f(r), definida sobre un dominio Ω, se puede expresar en términos de una función de interpolación como: hf(r)i=ZΩ f(r0)W(r−r0,h)dr0,(2) donde hi indica que es una aproximación, Ωes el domino de definición de f(r), y W(r−r0,h)es la función de interpolación denominada también función kernel ofunción de suavizado, que es no nula en una región del espacio llamada dominio soportado, con centro en el punto ry radio h. La función kernel desempeña un papel fundamental en el método SPH. Se trata de una función monótona decreciente con la distancia, que debe satisfacer la condición de normalización: ZΩ W(r−r0,h)dr0=1 (3) y al menos dos condiciones más: debe ser compacta en el dominio soportado y debe ser de paridad definida. Cualquier función que cumpla estos tres requisitos podría ser utilizada como función kernel. Con la ecuación (2) se tiene una expresión aproximada de cualquier función desconocida f(r), a partir de valores f(r0), en puntos de posición r0del dominio soportado, interpolados por la función kernel W(r−r0,h). Al aplicar la ecuación (2) sobre un dominio discretizado por un sistema de partículas, las integrales deben transformarse en sumatorios, en la forma expresada por la ecuación: fi=f(ri) = N(i) ∑ j mj fj ρj W(ri−rj,h),(4) Magnitudes como la aceleración, la velocidad, la posición, la densidad de masa o las fuerzas asociadas a una partícula i, pueden ser calculadas mediante la ecuación (4), donde f(ri) hace referencia a la magnitud fdefinida en la posición rique ocupa la partícula i,N(i)es el conjunto de las jpartículas vecinas a ique se encuentran dentro del dominio soportado en la partícula i, y donde mj,rjyρjson la masa, la posición y la densidad de masa respectivamente de la partícula j. Utilizar SPH simplifica las expresiones de las derivadas espaciales presentes en la ecuación (1) de Navier–Stokes, ya que se reducen a las derivadas espaciales de la función kernel, como vemos en las siguientes ecuaciones: ∇f(ri) = N(i) ∑ j mj fj ρj ∇W(ri−rj,h)(5) ∇2f(ri) = N(i) ∑ j mj fj ρj ∇2W(ri−rj,h).(6) Aunque con el método SPH se garantiza la conservación de la masa, no ocurre así con la conservación del momento lineal, por lo que se hace necesario imponer ciertas condiciones de simetría para los términos de presión y viscosidad de la ecuación (1) de Navier–Stokes, como veremos en las secciones 4y5. 4. Tratamiento del Término de Presión En esta sección vamos a describir el proceso seguido para obtener un término de presión para SPH que mejore la estabilidad de la simulación. En primer lugar vemos la necesidad de una expresión simétrica en el término de presión, y a continuación describimos un nuevo término de presión simétrico que contempla las regiones isobáricas. Del análisis realizado por Monaghan [Mon92] se puede deducir que dos son las condiciones que debe cumplir el término de presión, −∇P, para que la simulación sea estable: la primera es que cumpla la conservación del momento lineal, y la segunda es que no genere comportamientos oscilatorios. Para obtener una expresión del término de presión que cumpla la conservación del momento, Monaghan [Mon92] presenta un método basado en relaciones de identidad asociadas al cálculo diferencial, expresadas en las ecuaciones: ∇f(r) = 1 ρ(∇(ρf(r))−f(r)∇ρ)(7) ∇f(r) = ρ∇f(r) ρ+f(r) ρ2∇ρ(8) Al aplicar estas ecuaciones en el cálculo del término de presión, obtenemos expresiones en las que intervienen las presiones por pares de partículas, expresiones simétricas, que aseguran la conservación del momento lineal. En este sentido Müller et al. [MCG03] utilizan la presión promedio entre cada par de partículas, que a diferencia de la propuesta de Monaghan, no se fundamenta en relaciones de identidad. Aunque la expresión obtenida por Müller es más simple que las obtenidas por Monaghan, requiere de una elección meticulosa de la función kernel para evitar los comportamientos c The Eurographics Association 2012. 13 J.J. Perea & J.M. Cordero / More Stable SPH-based Fluid Simulation oscilatorios y los agrupamientos inapropiados de las partículas [MCG03]. Las expresiones simétricas que se pueden obtener, bien mediante el método de Monaghan, bien con la propuesta de Müller, no evitan la aparición de inestabilides en el campo de velocidades. Para evitar este problema Monaghan [Mon92] propone el uso de un término artificial de amortiguamiento viscoso (viscosidad artificial). Se trata de un coeficiente cuadrático de la velocidad que se ha mostrado muy eficaz en simulación de fluidos a alta velocidad. Una de las ventajas que presenta este término de viscosidad es que no altera la conservación del momento en aquellas expresiones simétricas que previamente lo cumplían. Sin embargo, por depender de la velocidad al cuadrado, en simulaciones de fluidos a velocidades muy inferiores a la velocidad del sonido, el término no se muestra tan eficaz. Además, para que este término no induzca inestabilidades, se requiere de una correcta selección de dos de los parámetros que en él intervienen [FDTA09]. Un efecto importante, que no se ha considerado en estos estudios, es el comportamiento del fluido en zonas con igual presión, las regiones isobáricas. En estas formulaciones, al utilizar la suma de presión entre partículas, el gradiente no se anula en las regiones isobáricas, como cabría esperar, lo que genera la aparición de fuerzas no nulas donde deberían serlo. Estas fuerzas inducirían un aumento de la velocidad no contemplado por la ecuación de Navier–Stokes, lo que resulta en una simulación menos estable. Este hecho se hace aún más evidente en situaciones en las que el fluido tiende a estancarse, ya que la viscosidad artificial propuesta por Monaghan casi no interviene. Una forma simple de corregir este problema sería considerar la diferencia de presiones obtenida a partir de la ecuación (7). Sin embargo, utilizar esta expresión en todo el fluido genera inestabilidad fuera de las regiones isobáricas. Veamos como obtener el término de presión propuesto. A partir de las relaciones de identidad (7)y(8) podemos obtener dos expresiones para el gradiente de presión: ∇Pi=1 ρi"N(i) ∑ j mjPj−Pi∇W(rj−ri,h)#(9) ∇Pi=ρi"N(i) ∑ b mj Pj ρ2 j +Pi ρ2 i!∇W(rj−ri,h)#(10) Sumando estas ecuaciones obtenemos: ∇Pi=1 2ρi"N(i) ∑ j mj∇W(rj−ri,h)Pm#(11) Pm="Pj−Pi+ρ2 i Pj ρ2 j +Pi ρ2 i!# (12) que resulta en una expresión con sumandos simétricos, por lo que satisface la conservación del momento. Sin embargo, sigue sin anularse en las zonas en las que el gradiente de presión es nulo. Para solventar este problema reformulamos la ecuación (11) como una combinación lineal de ambos sumandos, pero ponderada mediante dos coeficientes de control dependientes de la presión relativa, obteniendo finalmente las ecuaciones (13)y(14), que representan el término de presión propuesto. De este modo, aplicamos la diferencia de presión en las regiones isobáricas y la suma fuera de ella. ∇Pi=1 2ρi"N(i) ∑ j mj∇iW(rj−ri,h)Pp#(13) Pp="αPj−Pi+βρ2 i Pj ρ2 j +Pi ρ2 i!# (14) En (14), αyβson los coeficientes de control. Estos coeficientes deben cumplir dos criterios, el primero es que 0 ≤ α≤2y0≤β≤2, y el segundo es que α+β=2. Con estas restricciones garantizamos que cuando Pi∼Pj, entonces se cumpla que α∼2 y β∼0. Contrariamente, si Pi>> Pj oPi<< Pjentonces α∼0yβ∼2. Para valores de presión intermedios utilizamos la siguiente relación cuadrática, dependiente de la presión relativa: α=−1+1 4δP2+1 δP+2 (15) donde δes un escalar que denominaremos factor isobárico, yPes un coeficiente de presión relativa dado por la ecuación: P=|Pi−Pj| (max{Pi,Pj}).(16) A través del factor isobárico controlamos el peso entre los sumandos de la ecuación (13). Este factor ofrece buenos resultados para valores en el rango 0.05 ≤δ≤0.25. Presenta dos ventajas. La primera es que establecemos un umbral para un mejor control de las regiones isobáricas, donde a mayor δ menor precisión en las zonas con igual presión. La segunda es que homogeneizamos el campo de presión entre partículas con valores muy próximos. La consecuencia es una mejora en la estabilidad, ya que inducimos un efecto de campo local entre las partículas que se encuentran muy próximas con valores de presión similares, sin que ello genere agrupamientos irreales. Cuando las presiones superan el valor acotado por δ, comienza a actuar el término de la suma de presiones de la ecuación (13), aumentando el gradiente de presión. Este control de estabilidad se hace más relevante en simulaciones donde el fluido tienda a estancarse, ya que las presiones de las partículas tienen valores muy similares. Además, con el uso de la ecuación (13), conseguimos simulaciones estables sin necesidad de utilizar el término de viscosidad artificial propuesto por Monaghan con lo que evitamos los problemas, ya mencionados, asociados a una mala selección de los parámetros que intervienen. 5. Tratamiento del Término de Viscosidad Una incorrecta formulación del término de viscosidad puede introducir inestabilidad en la simulación debido a las flucc The Eurographics Association 2012. 14 J.J. Perea & J.M. Cordero / More Stable SPH-based Fluid Simulation tuaciones u oscilaciones irreales de las partículas. Veamos inicialmente las principales propuestas en este campo, para a continuación describir el término de viscosidad que proponemos en comparación con estas. El término de viscosidad debe cumplir, al igual que el término de presión de la sección 4, la condición de simetría, para satisfacer la conservación del momento lineal. A este respecto se han propuesto distintas formulaciones en las que intervienen la velocidad de las partículas por pares. Uno de los más empleados en CG es el propuesto por Müller et al. [MCG03]. Proponen un término basado en el promedio de la velocidad relativa de cada par de partículas. Para evitar oscilaciones irreales de las partículas, requiere del uso de funciones kernel cuyo laplaciano sea positivo, que a su vez conlleva a que la viscosidad adquiera también valores positivos, lo cual afecta al caracter disipativo del término µ∇2v. Para solventar el problema de las oscilaciones irreales, otras investigaciones proponen un factor de amortiguamiento, que no debe alterar la simetría del término viscoso para cumplir con la conservación del momento lineal, ni tampoco debe alterar los resultados dinámicos especialmente cuando las distancias entre las partículas es muy pequeña. En este sentido Zhou et al. [ZCGL10] añade un factor de viscosidad dependiente de la densidad, la velocidad del fluido y de coeficientes obtenidos a partir de las propiedades físicas del fluido. De este modo obtienen resultados estables. Amini et al. [AEF11] presenta un término de viscosidad que depende de la velocidad relativa entre las partículas, de manera que a mayor velocidad mayor amortiguamiento. Evitan fluctuaciones irreales de las partículas, pero afecta a la simetría, y por tanto a la conservación del momento lineal, del término de viscosidad. Para evitar los problemas de fluctuaciones y conservar el momento lineal, y siguiendo la propuesta de Zhou et al. [ZCGL10] y de Amini et al. [AEF11], presentamos nuestra formulación del término de viscosidad: µ∇2vi=ρi"N(i) ∑ b mj µjvj ρ2 j +µivi ρ2 i!∇2W(rj−ri,h)#Γ (17) Γ=exp −2|vj−vi| hiρiρj|rj−ri|!(18) Esta ecuación está calculada a partir de la ecuación (8), con lo que se garantiza la conservación del momento lineal. Le añadimos un factor de amortiguamiento exponencial, ecuación (18), dependiente de la velocidad relativa entre partículas, su distancia y el radio del dominio soportado. De este modo, obtenemos un factor simple pero eficiente que aumenta el amortiguamiento conforme crece la distancia y la velocidad relativa. A diferencia de la propuesta de Amini et al. [AEF11], este factor mantiene la conservación del momento lineal. Además no requiere de funciones kernel con un laplaciano positivo, y por tanto evita los problemas de la solución de Müller et al. [MCG03]. En la ecuación (17), µiyµjson las constantes de amortiguamiento viscoso asociadas a las posiciones de las partículas iyjrespectivamente, y hies el dominio soportado asociado a la partícula i. Como podemos observar al analizar la ecuación (17), al intervenir hien el denominador de la exponencial favorecemos la interacción entre partículas dentro del dominio soportado, y aumentamos el amortiguamiento cuando se aproximan al límite, es decir cuando |ri−rj| → hi. Además, formulamos un amortiguamiento selectivo con el que permitimos la interacción a corta distancia, y controlamos los desplazamientos irregulares a larga distancias, amortiguándolos. De este modo favorecemos la cohesión propia del fluido y evitamos fluctuaciones, favoreciendo todo ello la estabilidad. 6. Nueva Función Kernel Veamos en primer lugar, según los trabajos relacionados, qué características debe cumplir una función kernel. Destacaremos aquellas características que garantizan o mejoran la estabilidad del sistema simulado. Finalmente proponemos una nueva función kernel. En la sección 3, mencionamos los requisitos básicos que debe cumplir una función kernel: normalización, compacidad y paridad. Fulk [Ful94] afirma que las funciones kernel cuyas derivadas sean continuas y suficientemente suaves mejoran la estabilidad. De lo que se deduce que una función cuya derivada tenga discontinuidad de salto puede inducir inestabilidad. Swegle et al. [SHS95] demostraron que la inestabilidad estaba íntimamente relacionada con las características de la función kernel. Morris [Mor95] estudió la precisión que se obtenía con distintas funciones kernel, y demostró que una correcta elección podía mejorar la estabilidad del sistema. Capuzzo–Di Lisio [CDDL00], en una línea analítica similar a la de Morris, estableció las condiciones que debe cumplir la función kernel para obtener estabilidad. Belytschko et al. [BGL∗00], efectúan un análisis sobre el grado polinómico de las funciones kernel. Liu et al. [LLL03] relacionan la estabilidad de la función kernel con su continuidad y la de sus derivadas. Demuestran que cuanto mayor sea la clase de la función utilizada mejor será la estabilidad obtenida, y exigen que una función kernel, para que garantice estabilidad, debe ser al menos de clase C2. Hongbin–Xin [HX05] establece tres criterios para demostrar si una función kernel es estable: primero, las magnitudes dinámicas de sus partículas son iguales; segundo, la función kernel da como resultado que cualquier magnitud dinámica adimensional vale 1 y la de su derivada vale 0; y tercero, que la derivada de la función kernel es nula cuando la distancia tiende a cero. Establecemos por tanto que según los trabajos mencionados una función kernel es estable si: •Es simétrica respecto del origen. •Tanto ella como su primera derivada se anulan en el límite del dominio soportado, con lo que se evita la aparición de inestabilidades en el término de presión. •Es definida positiva, para impedir que la densidad de masa sea negativa. c The Eurographics Association 2012. 15 J.J. Perea & J.M. Cordero / More Stable SPH-based Fluid Simulation •Cuanto mayor sea el grado polinómico de la función, mejores serán los resultados obtenidos. Este hecho está relacionado con la clase Cnde la función kernel, por lo que una función de clase C∞, sería muy adecuada desde la pespectiva de la estabilidad. •Cumple con los mencionados criterios establecidos por Hongbin–Xin. Veamos algunos trabajos que proponen funciones kernel que cumplen algunas de estos criterios de estabilidad. Monaghan [Mon92] propuso una función gaussiana que no satisface la condición de compacidad en el dominio soportado ya que presenta una extensión espacial infinita. Morris [Mor96] solventó ese problema de compacidad mediante el desarrollo de funciones Spline de quinto grado que presentan un comportamiento óptimo en simulaciones a alta velocidad, lo que las hace ideales en simulaciones de gases, no obstante no ofrecen un comportamiento tan óptimo en simulaciones de líquidos a bajas velocidades. Dresbum y Cani [DC96], proponen un polinomio de tercer grado, con el que se obtienen buenos resultados en situaciones de baja densidad nodal, sin embargo, presenta dos inconvenientes: el primero es que no es adecuado para aproximar el término de viscosidad de la ecuación de Navier–Stokes, [MCG03], el segundo es que la primera derivada presenta una discontinuidad de salto, lo cual, según Fulk [Ful94], afecta a la estabilidad. Müller et al. [MCG03] resuelven este último problema planteado definiendo una función cúbica específica para el término viscoso. Pero Johnson et al. [JSB96] demuestran que las funciones cúbicas producen inestabilidades. Liu et al. [LLL03] proponen una función polinómica de cuarto grado que ofrece buenos resultados si se operaba con una alta densidad de partículas. De todo lo afirmado hasta ahora, se deduce que la mejor función kernel será aquella que sea de clase C∞, con una forma similar a la función gaussiana, pero compacta, simétrica, definida positiva y que se anule en el límite del dominio soportado. Lo mismo exigiremos para su derivada. Con estas premisas calculamos nuestra función kernel. Se trata de una función sinusoidal elevada al cuadrado, con lo que garantizamos que es definida positiva en el dominio soportado, con una forma muy similar a la función de Gauss, compacta, simétrica, y que se anula en el límite del dominio soportado, tanto ella como su derivada. La expresión de nuestra función kernel es: W(R,h) =              A0   sen2π hR π hR2  −h≤R≤h 0R<−h,R>h (19) donde Res la distancia relativa entre partículas, hes el radio del dominio soportado y A0es la constante de normalización. Para calcular A0interpolamos nuestra función con una función Q–Spline, obteniendo un valor de A0=5.471870115/(πh3) Figure 1: Función kernel propuesta, su gradiente y su laplaciano. Hemos considerado R=πen la representación. La forma que presenta nuestra función kernel, así como sus derivadas, se muestra en la figura 1, donde se evidencian algunas de sus características relevantes, como son la compacidad, la simetría, el contorno suave y la nulidad en el límite del domino soportado tanto de ella como de su derivada. Además de las características analíticas de nuestra función kernel, que ayudan a la estabilidad del sistema, presenta la ventaja adicional de poder utilizarla como única función kernel para calcular la densidad de masa, el término presión y el término viscoso, a diferencia de otras formulaciones [CBP05,MCG03] que requieren de distintas funciones kernel para cada uno de los términos de la ecuación de Navier–Stokes. Como vemos por la forma del gradiente en la figura 1, en los límites del dominio soportado, distancias próximas a h, la forma sinusiodal induce la aparición de interacciones repulsivas y atractivas. Esto evita los problemas de uso de fuerzas ficticias sinusoidales, como las propuestas por Tartakousky et al. [TMW09], para conseguir estabilidad en simulaciones en las que se analiza la formación e interacción de gotas de líquidos. A pesar de que en [MCG03] se afirma que una función kernel cuyo laplaciano presente valores negativos puede generar inestabilidades, en nuestro caso tal fenómeno no aparece, como demostraremos más adelante. Para evidenciar la estabilidad de nuestra función kernel utilizaremos los criterios de estabilidad fijados en el estudio de Hongbin–Xin [HX05], según el cual a partir de un sistema en situación estable con magnitudes dinámicas adimensionales se cumplirán las ecuaciones: MR d,α=WR αd+M0R d,α(20) M0R d,α= n ∑ s=1Ws α−R αd+Ws α+R αd (21) ∂W ∂rr=0 =0 (22) c The Eurographics Association 2012. 16 J.J. Perea & J.M. Cordero / More Stable SPH-based Fluid Simulation Compacity Simetry C1 C∞ Null in R=h Spiky X X 0,1,2 Exponencial X X X Quartic X X X 0,1 Quintic X X X 0,1 Q-Spline X X X 0,1,2 Proposed X X X X 0,1 Table 1: Comparativa entre las características de la función propuesta y otras funciones kernel. NR d,α=∂W ∂rR αd+N0R d,α(23) N0R d,α= n ∑ s=1∂W ∂rs α−R αd+∂W ∂rs α+R αd (24) donde Mes una relación obtenida a partir de la ecuación (4), Nes la derivada de la expresión M,des la distancia entre las partículas en la configuración inicial, αes el radio soportado adimensional y nhace referencia al número de términos que se desea considerar en el sumatorio. Hongbin–Xin afirma que las funciones kernel que garantizan la estabilidad son aquellas que cumplen que el valor de la expresión (20) es 1 y el de la expresión (23) es 0, independientemente del valor R/d. En la figura 2mostramos los valores que adquieren las magnitudes MyNpara diferentes valores de R/d. Podemos observar en la figura que la función kernel que proponemos cumple con los criterios de estabilidad fijados por Hongbin– Xin, ya que: M=1; N=0; y∂W ∂rr=0 =0 (25) independientemente del valor R/d. Además de la probada estabilidad de nuestra función kernel, vamos a evidenciar sus buenas características realizando una comparativa con otras funciones cuyo uso en SPH es muy extendido. Dicha comparativa se muestra en la tabla 1. En la tabla 1el campo Null in R =hhace referencia al orden diferencial que se anula en los límites del dominio soportado, 0 indica que se anula la función kernel, 1 se anula la primera derivada, 2 se anula la segunda derivada. De los resultados obtenidos podemos afirmar que nuestra función kernel cumple todas las características de estabilidad. 7. Nuevo Método de Integración Temporal Resolver la ecuación diferencial (1) de Navier–Stokes permite conocer el comportamiento dinámico del fluido. No es posible obtener, de manera general, una solución analítica a esta ecuación, por lo que en CFD se utilizan métodos numéricos. Uno de los problemas que presentan estos métodos, es el de la estabilidad en la simulación, especialmente en la integración temporal. Para mejorar la estabilidad es posible utilizar métodos implícitos de integración temporal, que transforman la ecuación de Navier–Stokes en un sistema de ecuaciones cuya dimensión depende del número de partículas. Para obtener simulaciones precisas requerimos de miles de partículas, lo que se traduce en un elevado coste computacional. Esta desventaja es aún mayor cuando aplicamos métodos implícitos a sistemas dinámicos operados mediante una formulación Lagrangiana de la ecuación del fluido, ya que el sistema obtenido aparece fuertemente acoplado, lo cual ralentiza aún más el proceso de resolución. Estas desventajas han llevado a que los métodos implícitos de integración temporal suelan emplearse con formulaciones Eulerianas [Lan11]. Un modo de simplificar el sistema obtenido por la aplicación de una metodología implícita, es mediante el uso de métodos aproximativos, si bien inducen amortiguamientos ficticios que afectan a la precisión y estabilidad de los resultados [SCJ12]. La utilización de métodos explícitos mejora los costes computacionales y reduce la complejidad de la resolución del problema, si bien para obtener resultados estables debemos utilizar pasos de tiempo inferiores a los utilizados en los métodos implícitos, o bien debemos emplear métodos multipaso, que redundarían en el aumento del coste computacional. Por lo general, en sistemas operados mediante SPH se utiliza el método de Verlet [CM09]. Si bien se obtienen buenos resultados visuales, el error cometido es del orden O(n4)como máximo, lo cual exige que el paso de integración sea muy pequeño para que los resultados obtenidos sean precisos. Proponemos la utilización de un método explícito con un orden de error muy bajo, que nos permita aumentar el paso de tiempo sin inducir inestabilidad. El problema de la imprecisión es más evidente en situaciones en las que el sistema comienza a ser inestable, ya que ello implica que las fuerzas presentes en la integración cambian apreciablemente dentro de cada paso de integración, lo que favorece aún más la inestabilidad. Blanc y Pastor [BP12] reducen este problema utilizando el método de Runge–Kutta, aumentando el orden de error hasta como máximo O(n5). Así se mejora la precisión de los resultados, pero el proceso se torna más complejo, lo que implica un aumento del coste computacional. Para obtener unos resultados precisos que mejoren la estabilidad, sin que ello implique un aumento en el tiempo de cálculo, presentamos nuestro método de integración temporal. Básicamente, para desarrollarlo seguimos un proceso de cálculo similar al descrito por Modak y Sotelino [MS02], pero adaptado para obtener un método explícito. Partiendo del desarrollo en serie de Taylor para las magnitudes posición, velocidad y aceleración, y considerando el teorema del valor medio, obtenemos: un+1=un+∆tvn+∆t2 2an+β2∆t2 3! (2an−2an−1+an−2), (26) c The Eurographics Association 2012. 17 J.J. Perea & J.M. Cordero / More Stable SPH-based Fluid Simulation Figure 2: Comparativa de funciones kernel según criterios de estabilidad. vn+1=vn+∆tan+∆t 2β1(2an−2an−1+an−2)(27) an+1=β0(3an−2an−1+an−2)∆t(28) donde an+1,vn+1,un+1hacen referencia a la aceleración, velocidad y posición en el instante tn+1, y donde β0,β1yβ2 son parámetros de estabilidad que cumplen las relaciones β0=1 β2+β1−1.082β2+1.19β1=0 El orden de error cometido para la posición es de O(n6), muy superior al ofrecido por el método de Verlet o el de Runge–Kutta. Sin embargo, esta mejora sólo exige, si se compara con el método de Verlet, almacenar un valor más, el asociado a an−2, con lo que ganamos precisión sin un apreciable coste de almacenamiento. El paso de integración tiene que satisfacer la condición de Courant–Friedrich–Levy [HK89]. 8. Resultados Para validar todas y cada una de las propuestas presentadas en este trabajo, hemos desarrollado un simulador de fluidos basándonos en la propuesta WCSPH [BTT09]. Aunque la ecuación (1) se obtiene bajo la hipótesis de fluido incompresible, es más adecuado, sobre todo cuando se trata con líquidos, considerar que estos poseen cierto grado de compresibilidad. Analíticamente, que un fluido sea ligeramente compresible implica que la densidad de masa y la presión no son independientes, sino que están relacionadas a través de una expresión denominada ecuación de estado. Existen diversas formulaciones para establecer la ecuación de estado de un líquido, una de la más utilizadas es la ecuación de Tait [Bat74]. Para esta ecuación hemos utilizado un valor del coeficiente empírico γinferior a 7 como en [AEF11]. Para la detección de partículas vecinas utilizamos una técnica basada en un algoritmo NNS, según se describe en Teschner et al. [THM∗03]. Para la respuesta consideramos el modelo con pérdida de energía, en la interacción sólido– fluido, descrito por [SBC∗11]. Para considerar la tensión superficial hemos seguido el método descrito por Müller et al. [MCG03]. Los resultados gráficos obtenidos se muestran en la figura 3. En la parte izquierda de la figura mostramos una captura de la simulación con partículas, donde las partículas han sido coloreadas en base a la densidad asociada a las mismas. En la parte derecha mostramos el mismo estado de simulación que aparece en la figura de la izquierda, pero renderizado. Para obtener estos resultados hemos utilizado los siguientes valores para los parámetros, en unidades internacionales: 0.001 para ∆t, 7018 como número de partículas, 0.01 para h, 985.63 para ρ, 0.26 para µ, 25·10−3como masa de la partícula, 5.65 para γ, 4.5·10−3para el coeficiente isobárico δ, 1 para β0, 0.7586 para β1y 1.101 para β2. En las simulaciones generadas observamos una mayor cohesión entre capas, como se evidencia por la forma de la cresta de la ola en la figura 3. Esto se debe a que en zonas de presión muy próximas, regiones isobáricas, la componente de las fuerzas asociadas a la presión desaparece, apareciendo un desplazamiento laminar más suave, lo que también favorece una forma más real del frente de la ola. 9. Conclusiones Hemos revisado en este trabajo algunos elementos utilizados en la simulación de fluidos mediante SPH, con la intención de aumentar la estabilidad de la misma. Hemos operado sobre el término de presión de la ecuación de Navier–Stokes, en el que hemos introducido un coeficiente de presión para simular adecuadamente y con precisión las zonas del fluido con igual presión. La expresión obtenida nos da estabilidad sin necesidad de introducir un término de amortiguamiento c The Eurographics Association 2012. 18 J.J. Perea & J.M. Cordero / More Stable SPH-based Fluid Simulation Figure 3: Una ola simulada haciendo uso de los mecanismos de mejora de la estabilidad descritos. A la izquierda la simulación generada por el simulador basado en SPH. A la derecha el mismo estado de simulación pero renderizado. ficticio. Hemos compensado las posibles oscilaciones aleatorias introduciendo un factor de amortiguamiento viscoso, simple pero efectivo, dependiente de la velocidad relativa y de la distancia entre las partículas. De este modo favorecemos los movimientos a corta distacia, pero amortiguamos conforme ésta aumenta, favoreciendo la cohesión propia de los fluidos reales. Hemos propuesto una función kernel cuyas características mejoran la estabilidad ya conseguida con los términos de presión y viscosidad. Una función kernel que puede ser utilizada para calcular todos lo términos de la ecuación de Navier–Stokes. Finalmente describimos un método de integración temporal explícito, desarrollado a partir de un método implícito y con una precisión superior a la obtenida por otros métodos de integración explícitos. La conjunción de estas cuatro mejoras permite obtener resultados satisfactorios incluso con un número de partículas relativamente bajo. Para validar todas estas conclusiones se ha desarrollado un simulador de SPH en el que se han implementado todas las propuestas. No hemos entrado en este trabajo en las soluciones propuestas para paralelizar el método, que sin duda mejorarían el coste computacional. Las propuestas presentadas pueden ser aplicadas tanto a soluciones secuenciales como la utilizada en el apartado de resultados y a soluciones paralelizables. 10. Agradecimientos Esta investigación ha sido financiado a través de los proyectos QSimTest (TIN2007-67843-C06_03) y Tempros (TIN2010-20057-C03-02) del Ministerio de Educación y Ciencia de España. References [AEF11] AMINI Y., EMDAD H., FARID M.: 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 (Oct 2011), 184–194. doi:10.1016/ j.euromechflu.2010.09.010.2,5,8 [AS05] ATA R., SOULAÏMANI A.: A stabilized sph method for inviscid shallow water flows. International Journal for Numerical Methods in Fluids 47 (Jan 2005), 139–159. doi:10.1002/ fld.801.2 [Bat74] BATCHELOR G.: An Introduction to Fluid Dynamics. Cambridge Mathematical Library, 1974. 8 [BGL∗00] BELYTSCHKO T., GUO Y., LIU K., W., XIAO P.: A unified stability analysis of meshless particle methods. International Journal for Numerical Methods in Engineering 48, 9 (May 2000), 1359–1400. doi:10.1002/ 1097-0207(20000730)48.5 [BLS11] BODIN K., LACOURSIÈRE C., SERVIN M.: Constraint fluids. IEEE Transactions on Visualization and Computer Graphics 18, 3 (Feb 2011), 516–526. doi:10.1109/TVCG. 2011.29. 2 [BP12] BLANC T., PASTOR M.: A stabilized fractional step, runge-kutta taylor sph algorithm for coupled problems in geomechanics. Computer Methods in Applied Mechanics and Engineering 221-222 (May 2012), 41–53. doi:10.1016/j.cma. 2012.02.006. 7 [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.8 [CBP05] CLAVET S., BEAUDOIN P., POULIN P.: Particlebased viscoelastic fluid simulation. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation (2005), SCA ’05, ACM, pp. 219–228. doi:10.1145/ 1073368.1073400.6 [CDDL00] CAPUZZO-DOLCETTA R., DILISIO R.: A criterion for the choice of the interpolation kernel in smoothed particle hydrodynamics. Applied Numerical Mathematics 34, 4 (Aug 2000), 363–371. doi:10.1016/S0168-9274(99)00036-7.5 [CM09] CSABA P., MATEUSZ G.: Validation of an sph sloshing simulation by experiments. ASME Conference Proceedings 2009, 5 (Jun 2009), 729–739. doi:10.1115/ OMAE2009-79792.7 [DC96] DESBRUN M., CANI M.: Smoothed particles: A new paradigm for animating highly deformable bodies. In Eurographics Workshop on Computer Animation and Simulation, EGCAS ’96, August, 1996 (Aug 1996), Boulic R., Hégron G., (Eds.), Springer, pp. 61–76. 1,6 [FDTA09] FERRARI A., DUMBSER M., TORO E. F., ARMANINI A.: A new 3d parallel sph scheme for free surface flows. Computers and Fluids 38, 6 (Jun 2009), 1203–1217. doi:10.1016/ j.compfluid.2008.11.012. 4 c The Eurographics Association 2012. 19