A new Potential-based SPH Method

CEIG - Spanish Computer Graphics Conference (2016) Alejandro García-Alonso and Belen Masia (Editors) A new Potential-based SPH Method J.J. Perea and J.M. Cordero University of Seville, Spain Abstract In Computer Graphics the fluid simulation plays a relevant role. Amongst the existing methods, the SPH method is highlighted by its simplicity, versatility and adaptability. Nevertheless, several research show some limitations that affect both the stability and the realism of simulation. These limitations are related with the calculation of the mass density and the pressure gradient, mainly when the number of neighbouring particles is low as boundary and discontinuities. Also, the simulation realism is compromised when the incompressibility constraint is enforced and the fluid heating is not dissipated. To overcome these limitations, we present a new technique to simulate the fluid flow. It involves of a hybrid technique that allow us to calculate a stable pressure gradient and the mass density even when the number of neighbouring particle is low. Also, we directly control the incompressibility and dissipate the fluid heating with accurate. Thus, we obtain stable and realistic fluid simulation, even with higher time steps. Categories and Subject Descriptors (according to ACM CCS): I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—Animation 1. Introducción La simulación dinámica de fluidos presenta un especial interés en el ámbito de Computer Graphics (CG). Desde un punto de vista analítico, es descrita por la Mecánica de Fluidos (FM). Desde una perspectiva numérica, es objeto de estudio a través la Dinámica de Fluidos Computacional (CFD). Mediante la CFD, las ecuaciones en derivadas parciales, que describen la dinámica del fluido, son transformadas en un conjunto de ecuaciones algebraicas, cuya resolución permite simular el flujo. En la FM se formulan hasta siete ecuaciones en derivadas parciales para describir el flujo [ZTN14]. No obstante, en la mayoría de las simulaciones llevadas a cabo en CG, sólo se considera la ecuación de Navier–Stokes, sobre la que se suele imponer algunas restricciones, como es el caso de la incompresibilidad [SP09]. Según la formulación Lagrangiana, la ecuación de Navier–Stokes es: ρ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 fuerza de presión, es la fuerza que surge en el fluido como oposición a la variación de presión entre diferentes áreas del fluido. El término ρges la fuerza que surge por estar sometido el fluido a la gravedad y µ∇2ves la denominada fuerza de viscosidad, de naturaleza disipativa y que depende, a través µ, de las propiedades constitutivas del fluido. Por lo general, la ecuación Navier–Stokes sólo puede resolverse numéricamente. Uno de los métodos numéricos más utilizados en CG es SPH [Luc77,GM77]. Descriptivamente, mediante SPH el fluido es modelado por un sistema de partículas, que representan fracciones discretas del volumen total del fluido. La interacción mutua entre partículas, es ponderada mediante una función monótona y decreciente con la distancia, la función kernel, cuyo radio de interacción, el radio del dominio soportado, debe ser limitado [Pri12]. Este radio, acota el número que partículas que interactúan mutuamente, es decir, restringe el número de partículas vecinas. El método SPH es especialmente atractivo por su versatilidad y adaptabilidad en complejos entornos de simulación [MCG03], por satisfacer la conservación de la masa [MCG03,Mon92], así como por reducir la complejidad de las ecuaciones a resolver, al basarse en la formulación Lagrangiana [Mon92,LL03]. Sin embargo, a pesar de estas ventajas, diversos estudios evidencian problemas de estabilidad que afectan al realismo de las simulaciones. Estos problemas están relacionados con: •La imprecisión en la fuerza de presión y la densidad, especialmente cuando el número de partículas vecinas es bajo, como ocurre en las regiones donde el fluido interactúa con el contorno o en las discontinuidades [SP09,Pri12,PR14]. c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. DOI: 10.2312/ceig.20161315 J.J. Perea & J.M. Cordero / A new Potential-based SPH Method •La imposición de la incompresibilidad, especialmente en los líquidos [SP09,MM13]. •La disipación del calentamiento que sufre el fluido, principalmente en las zonas de contorno y en las discontinuidades [Pri12,DA12]. La consecuencia, es que se generan simulaciones con dispersiones espurias y agrupaciones irreales de partículas en diversas regiones del fluido [SP09,Pri12,DA12]. En este artículo proponemos una nueva técnica para resolver estos problemas de estabilidad. Se trata de una técnica híbrida SPH- potencial de presión con la que obtenemos estabilidad, eficiencia y versatilidad. En síntesis, formulamos la fuerza de presión a partir de un potencial conservativo, mientras que el resto de magnitudes de la ecuación (1) son calculadas mediante SPH. Así eliminamos la dependencia de la fuerza de presión respecto de la densidad que suele ser la causa de inestabilidad, especialmente en regiones donde el número de partículas es bajo [SP09,Pri12]. Además, por la propia formulación de esta fuerza, controlamos la condición de incompresibilidad y calculamos el calentamiento del fluido que disipamos. Para calcular la densidad, utilizamos la formulación de SPH pero relacionándola con la fuerza de presión. De este modo, la densidad “hereda” la estabilidad y la incompresibilidad que se consigue en la fuerza de presión. Finalmente, proponemos una fuerza de viscosidad que favorece la cohesión del fluido sin que ello implique sobreamortiguamiento. Tras esta introducción el resto del artículo se organiza del siguiente modo. En la sección 2revisamos las investigaciones más relevantes desde la perspectiva de la estabilidad y realismo. En la sección 3describimos los fundamentos del método SPH necesario para entender la formulación de la densidad y de la fuerza viscosidad desarrolladas. Después de esto, en la sección 4, describimos el modelo propuesto con el que se mejora la estabilidad y realismo. Dentro de él, en la sección 4.1, fijamos las condiciones que debe satisfacer la fuerza de presión para conseguir estabilidad y formulamos dicha fuerza. Por otro lado, en la sección 4.2 describimos el proceso propuesto para imponer la incompresibilidad. Seguidamente, en la sección 4.3, mostramos la técnica seguida para calcular la densidad. En la sección 4.4 proponemos la formulación de la fuerza de viscosidad y analizamos sus características analíticas. En otro orden, en la sección 5, mostramos los resultados obtenidos al implementar nuestra propuesta, evidenciamos la mejora en la estabilidad y el realismo, así como los rangos de valores en los que el modelo es estable. Finalmente, en la sección 6describimos los objetivos logrados con nuestra técnica. 2. Trabajo Relacionado Diversas investigaciones describen los problemas de estabilidad que aparecen en simulaciones mediante SPH. Entre ellas, las tres causas más destacables son la imprecisión en la fuerza de presión, especialmente cuando el número de partículas vecinas es bajo; la imposición de la incompresibilidad en los líquidos y la disipación del calentamiento que se produce en el fluido. En referencia a la imprecisión en la fuerza de presión, existe un consenso en relacionarlo con el número de partículas vecinas [LLL03,SP09,Pri12,DA12,PR14]. Conforme el número de partículas vecinas disminuye, la imprecisión de la fuerza de presión aumenta. La consecuencia, es que las partículas se dispersan y siguen trayectorias erráticas. Usualmente, se ha considerado que incrementando el número de partículas vecinas, aumentando el radio del dominio soportado, se reduce la imprecisión. El problema es que existen regiones en el fluido, como en el contorno o en las discontinuidades, donde el incremento del número de partículas no resuelve el problema [SP09,DA12,PR14]. Además, investigaciones cómo las llevada a cabo por [YPL14,Pri12,DA12] cuestionan que el incremento del número de partículas vecinas resuelva este problema y advierten, que aumentar el radio del dominio soportado, disminuye la resolución de los resultados, favoreciendo el comportamiento espurio de las partículas. Otros estudios, relacionan la imprecisión en la fuerza de presión, no sólo con el número de partículas vecinas, sino también con las características analíticas del gradiente de la función kernel, el cual interviene en la fuerza de presión. En este ámbito, [AS05,CGSR00,DA12] demuestran que la presencia de un punto de inflexión en el gradiente de la función kernel, y por lo tanto en la fuerza de presión, favorece la estabilidad. A partir de estas causas de inestabilidad en la fuerza de presión, se han desarrollado diversas técnicas para evitarlas. [DRI97] propone la introducción de partículas de compensación, su función sería amortiguar la dispersión errática en los contornos del fluido o evitar el sobreamortiguamiento fuera ellos. En una linea similar, [YMSH04] propone añadir fuerzas de control que inducen un efecto similar al de las partículas de compensación [DA12]. No obstante, aunque se consigue reducir, en cierta medida, el sobreamortiguamiento no ocurre lo mismo con la dispersión de partículas en los contornos. Otros autores, [Mor96,Abe11] proponen la reformulación de la discretización del gradiente de la función kernel, de manera que varíe dependiendo de las regiones en las que opere. Aunque con esta propuesta se logran resultados satisfactorios, para condiciones de simulación muy específicas, la variación del gradiente, entre diferentes regiones del fluido, puede violar la consevación del momento [Mon11,DA12]. Para evitar los problemas asociados al bajo número de partículas vecinas y prescindir del uso del gradiente de presión, se ha desarrollado un conjunto de técnicas alternativas [CBP05,ATO15]. Esencialmente, en estas técnicas se “acomodan” las posiciones de las partículas para satisfacer las restricciones impuestas sobre el fluido, de ahí su designación cómo basada en la posición de partículas (DBP). Esta colocación de las partículas se fija considerando potenciales armónicos [CBP05], o imponiendo formulación de extremos condicionados [MM13]. Aunque los resultados que se obtienen son estables y realistas [ATO15], presentan algunas limitaciones ya que se pierde información dinámica. Además, diversos estudios alertan de que inducen un sobreamortiguamiento que puede afectar al realismo de la simulación [FSJ01,MM13]. Respecto al segundo factor que condiciona la estabilidad de la simulación, la imposición de la incompresibilidad, se han desarrollado diversas propuestas [CM99,SP09,MM13]. En [IOS∗14] se puede encontrar una descripción precisa sobre estas propuestas. En base a su formulación se pueden clasificar en dos categorías: las basadas en la ecuación de estado (EOS) [Mon92,BT07,SP09] y las basadas en métodos de proyección [MM13]. Dentro de los algoritmos basados en la EOS, la manera más básica de imponer la incompresibilidad del fluido es añadir fuerzas ficticias que controlen c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. 60 J.J. Perea & J.M. Cordero / A new Potential-based SPH Method el gradiente de presión, de manera que las partículas no puedan aproximarse a una distancia incompatible con la incompresibilidad [Mon92,Mon94]. Con estas técnicas se obtienen buenos resultados, tanto en fluidos compresibles [Mon92], cómo débilmente compresibles [BT07]. Sin embargo, para tasas de compresibilidad bajas, inferiores al 5% se requiere de incrementos de tiempo relativamente bajos [SP09], del orden de 10−6s–10−5s. Esta limitación es consecuencia de no tener un control implícito de la incompresibilidad sobre la fuerza de presión. Para solventar esta limitación y mejorar el rendimiento de las simulaciones del fluido, [SP09] desarrolla una técnica basada en un proceso predictor–corrector en el que, en cada paso de tiempo, la fuerza de presión es rectificada, desde un valor inicial, hasta un valor acorde con la condición de incompresibilidad. En una línea de investigación similar, [HLL∗12] desarrolla un modelado local de la ecuación de presión de Poisson (PPE), con el objetivo de mejorar el rendimiento del proceso predictor–corrector. Cabe destacar que, si bien estas técnicas mejoran los resultados obtenidos para fluidos débilmente compresibles, inducen un incremento en el coste computacional, debido al proceso predictor–corrector. Además, son susceptibles de mostrar sobreamortiguamiento [MM13] y pueden volverse imprecisas cuando el número de partículas vecinas es bajo, como ocurre en los contornos del fluido [SP09]. Como alternativa a las metodologías basadas en la ecuación de estado, se han desarrollado los denominados métodos basados en proyección. Esencialmente, en estos métodos, se calculan las velocidades intermedias de las partículas que son predichas sin considerar las fuerzas de presión. Determinadas las velocidades intermedias, se calcula la fuerza de presión a partir de la PPE. Con esta fuerza de presión se corrigen las velocidades intermedias de manera que cumplan la ecuación dinámica. El problema que presentan estas técnicas se debe a que fueron formuladas para métodos basados en la discretización del espacio mediante mallas y, su adaptación a SPH, suele implicar una pérdida de precisión en la fuerza de presión, lo que induce inestabilidad. Para evitar este problema se suele utilizar un modelo dual malla–partículas [LTKF08] que, por lo general, induce una pérdida de “información” entre los puntos de muestreo de la malla y las partículas [RWT11]. Para evitar esta pérdida, [RWT11] propone el uso de la EOS para rectificar los resultados asociados a las partículas. No obstante, aunque se obtiene una importante mejora en términos de precisión, su complejidad y pérdida de eficiencia lo hace poco adecuado para los resultados requeridos en CG [LXSR12]. Esto a hecho que las técnicas que imponen la incompresibilidad explícitamente, sin necesidad de utilizar la ecuación PPE sean muy empleadas [SP09,MM13]. Finalmente, en referencia a la tercera causa de inestabilidad, es sabido que durante la evolución del flujo aparece un calentamiento que debe ser disipado para conseguir simulaciones estables [Mon97,PR14]. Este fenómeno es especialmente relevante en las discontinuidades, en las interacciones del contorno y en la simulación de la incompresibilidad [Pri08,KOG∗13]. Investigaciones llevadas a cabo por [Pri12,Mon97,Pri08,OJE∗03,KOG∗13] describen varias propuestas para cuantificar y disipar este calentamiento. Generalmente, esto se logra añadiendo un término disipativo a la ecuación de Navier–Stokes. El problema que se plantea es cuantificar el calentamiento en las regiones donde se produce, para una correcta disipación, ya que, en caso contrario, la simulación se vuelve imprecisa, sobreamortiguada e inestable [FDTA09,PR14]. 3. El Método SPH En esta sección, vamos a describir los fundamentos del método SPH. Esta descripción constituirá el marco numérico para las ecuaciones que empleamos en nuestro modelo. El método SPH, inicialmente desarrollado por Gingold y Monaghan [GM77] y Lucy [Luc77], suele formularse a partir de la formulación integral de la función delta de Dirac: f(r) = ∫f(r′)δ(r−r′)dr′,(2) donde frepresenta a cualquier función escalar y δes la función delta de Dirac. Considerando la teoría de interpolación de funciones, la función δpuede ser interpolada, o sustituida, por cualquier función con dominio espacial. Según este criterio, la ecuación (2) puede generalizarse en términos de la ecuación: f(r) = ∫Ω f(r′)W(r−r′,h)dr′+O(h2),(3) dondeWes la función kernel, Ωrepresenta el volumen del dominio de definición Wyhes el radio del dominio soportado. Analíticamente, la función kernel tiene que satisfacer, al menos, las siguientes características: 1. Debe estar normalizada y tender a δcuando htiende a cero, es decir: lim h→0W(r−r′,h) = δ(r−r′)y∫V WdV′=1.(4) 2. Debe ser definida positiva y continuamente decreciente con (r− r′). 3. Debe ser simétrica respecto del a (r−r′). Cuantitativamente, modelar un medio continuo mediante un sistema de partículas, requiere transformar la expresión integral, ecuación (3), en un sumatorio donde cada partícula representa a un elemento de masa ρdV, es decir: ⟨f(r)⟩=∫Ω f(r′) ρ(r′)W(r−r′,h)ρ(r′)dr′≈∑ j∈N(i) mj fj ρj W(rj−ri,h), (5) donde ⟨⟩hace referencia a la función aproximada de f, ya que solo se consideran hasta el término de segundo orden en el desarrollo en serie, N(i)representan al conjunto de partículas jvecinas de iyρj es la densidad de masa asociada a cada partícula j. A partir de la propiedad de simetría, (punto 3.), el gradiente de cualquier magnitud escalar f(r)se puede simplificar cómo: ⟨∇f(r)⟩=∂ ∂r∫f(r′) ρ(r′)W(r−r′,h)ρ(r′)dr′≈ c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. 61 J.J. Perea & J.M. Cordero / A new Potential-based SPH Method ≈∑ j∈N(i) mj fj ρj∇W(rj−ri,h).(6) Generalizando la ecuación (6), se puede obtener una ecuación válida para cualquier orden diferencial. Esta ecuación formal es: ⟨∇lf(ri)⟩=∑ j∈N(i) mj fj ρj∇lW(rj−ri,h),(7) donde lhace referencia al orden diferencial. Analíticamente, aunque la ecuación (7) cumple las premisas de SPH, ésta presenta dos limitaciones relevantes. La primera es que no satisface la conservación del momento, por no tratarse de una expresión simétrica. La segunda es que no considera términos de segundo orden, esta segunda limitación afecta, especialmente cuando se aplica a la (1), a la fuerza de presión, en los límites del contorno y en las discontinuidades [DA12,Pri12,PR14]. Para evitar la primera limitación asociada a la ecuación (7), Monaghan [Mon92] propone una relación de simetrización que con la que se pueden obtener términos de la ecuación de Navier–Stokes que cumplen la conservación del momento. Esta expresión es: ∇lf(r) = ρ(∇l(f(r) ρ)+f(r) ρl+1∇ρ).(8) 4. Modelo Propuesto En esta sección vamos a describir el modelo propuesto. En primer lugar vamos a establecer, y justificar, la hipótesis fundamental en la que se basa. Seguidamente, estableceremos las restricciones que debe satisfacer la fuerza de presión, para garantizar estabilidad en cualquier condición. Considerando estas restricciones calcularemos la fuerza de presión. A continuación, con la fuerza de presión, y gracias a sus características, impondremos la incompresibilidad y evaluaremos el calentamiento. Una vez que la incompresibilidad es impuesta, vamos a describir el proceso por el que calculamos la densidad. En este proceso aunamos la formulación fundamental de SPH y la fuerza de presión calculada. Finalmente, vamos describir nuestra propuesta de fuerza de viscosidad. Justificaremos porqué esta fuerza favorece la estabilidad y la cohesión sin introducir sobreamortiguamiento. Analíticamente, la fuerza de presión permite cuantificar dinámicamente la fluidez del fluido. Asociada a un medio continuo, sólo aparece cuando la distancia relativa, entre secciones contiguas, varía respecto de una situación de equilibrio [ZTN14]. A partir de este hecho, establecemos la hipótesis fundamental de nuestro modelo: la fuerza de presión se debe al cambio de la distancia relativa entre partículas. Con esta premisa, eludimos la intervención de la densidad en la fuerza de presión, con lo que evitamos la inestabilidad que esta magnitud induce cuando el número de partículas es bajo [PM04,DA12]. Por otro lado, al depender la fuerza de presión de la distancia relativa, podemos formular un potencial escalar conservativo –en lo sucesivo referido como potencial de presión– del que se deriva la fuerza, es decir:  Fpi=−∇Vpi(R),(9) donde  Fpies la fuerza de presión asociada a la partícula i,Vpi(R)es el potencial de presión y R=|rj−ri|es la distancia relativa entre las partículas. Con esta formulación evitamos el uso de la teoría de interpolación, que es la base analítica del método SPH; por lo que la función potencial, y en consecuencia la fuerza de presión que se deriva de ella, no requiere estar normalizada, lo que simplifican los cálculos. Establecida la hipótesis fundamental, seguidamente vamos a fijar las condiciones que debe cumplir el potencial de presión, del que posteriormente obtendremos la fuerza de presión. 4.1. Potencial y Fuerza de Presión Diversos estudios relacionan la fuerza de presión con la estabilidad y el realismo en las simulaciones. A partir de estas investigaciones, vamos a fijar las características que debe cumplir nuestro potencial de presión del que derivamos la fuerza de presión: 1. Debe ser independiente del número de partículas vecinas, de este modo, se garantizaría la estabilidad incluso en lugares donde éste sea muy bajo. Deducimos esta restricción de los estudios desarrollados por [SP09,Pri12,DA12,MM13], donde se indica que, un bajo número de partículas vecinas, induce inestabilidad no sólo en la fuerza de presión, sino también en la densidad. 2. Debe presentar concavidad con respecto del eje de abscisas, con un valor mínimo de potencial. Esta condición es una extrapolación de las características que debe satisfacer el gradiente de la función kernel en la formulación SPH [CGSR00,DA12,PR14]. Según estas investigaciones, un gradiente cóncavo y con un punto de inflexión favorece la estabilidad. 3. Su dependencia con la distancia debe ser asimétrica, respecto del valor mínimo de potencial. De este modo, se satisfacen dos requisitos: el primero es que se incremente el valor del potencial conforme la distancia relativa disminuye, como se deduce del estudio de [LL03]; el segundo es ponderar este potencial con la propia distancia, como puede concluirse de los estudios realizados por [Pri12,DA12]. Un potencial armónico satisfaría la primera condición pero no la segunda, lo que no garantiza una distribución ponderada de las magnitudes dinámicas, afectando a la estabilidad [Pri12]. 4. Debe considerar la restricción de incompresibilidad [SP09, MM13], por lo tanto, debe limitar la distancia mínima a la que las partículas se pueden aproximar. Por ello, el potencial debe presentar un comportamiento asintótico con la distancia, para impedir que las partículas se aproximen a una distancia incompatible con la incompresibilidad. c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. 62 J.J. Perea & J.M. Cordero / A new Potential-based SPH Method A partir de estas características, formulamos nuestra función potencial: V(ri) = 5.25 6πh∑ j∈N(i)(2 (R−0.9ε)2−3.5π2 ζ(R−(0.85ε)) +(3π2α)), (10) donde εlo denominamos factor asintótico, con el que podemos tasar la incompresibilidad, ζes el factor de profundidad, con el que controlamos el agrupamiento de las partículas y αes el coeficiente de dispersión, a través del cual podemos cuantificar el valor del potencial en el limite de interacción. El rango en el que estos valores favorecen la estabilidad son: 0.3≤ε≤h, 0.5≤ζ≤2.0 y 0.8≤α≤2.25. En la Figura 1se muestra el perfil del potencial de presión formulado. Analizando la fórmula del potencial, ecuación (10), podemos afirmar que, la interacción entre cada dupla de partículas sólo depende de su distancia relativa. Así, evitamos los problemas inducidos por la densidad cuando el número de partículas de vecinas es bajo [DA12,PR14] y cumplimos con la exigencia impuesta en el punto 1. Figure 1: Grafica del potencial. Los valores que hemos considerado son h=3, ε=0.6, ζ=1.25 y α=1.0. Analizando la gráfica del potencial, Figura 1, podemos comprobar que es cóncavo respecto de la abscisa, como se exigía en el punto 2. Támbién, que presenta un perfil asimétrico respecto del mínimo del potencial, acorde con el punto 3. Además, por poseer un límite asintótico vamos a poder controlar la incompresibilidad, como se exigía en el punto 4. Una vez que hemos formulado el potencial de presión vamos a obtener la fuerza de presión, según la ecuación (9). La fuerza de presión obtenida es:  Fpi=−5.25 6πh∑ j∈N(i)(3.5π2 ζ(R−(0.85ε))2−4 (R−0.9ε)3)Γ R R, (11) donde Γ=h−R. Hemos introducido este término para considerar la restricción de que en el fluido la fuerza de presión se satura con la distancia [ZTN14]. Esta restricción es destacada en investigaciones desarrolladas por [Mon92,Mor95,DA12,YPL14]. Cabe destacar que la fuerza de presión, además de las características “heredas” del potencial, cumple la conservación del momento, condición esencial en la dinámica. Seguidamente, establecida la fuerza de presión, vamos a describir cómo controlamos de la incompresibilidad aprovechando su límite asintótico. 4.2. Imposición de la Incompresibilidad y Disipación del Calentamiento Desde un punto de vista dinámico, la incompresibilidad se relaciona con la fuerza maxima que es capaz de “soportar” el fluido sin modificar su volumen local [ZTN14]. En base a este hecho, cuando el fluido se modela mediante partículas, la incompresibilidad se puede explicar como la fuerza de presión máxima, que es capaz de soportar cada partícula, sin modificar el volumen local asociado. En este contexto, la fuerza que hemos calculado, ecuación (11), permite modelar directamente la incompresibilidad, al presentar un valor asintótico que limitará la fuerza de presión máxima. Sin embargo, para tasas de incompresibilidad inferiores al 2%, tenemos que introducir un proceso para controlar los posibles errores numéricos acumulados en cada paso de tiempo. Estos errores no son relevantes para tasas de incompresibilidad superiores. Para tasas de incompresibilidad inferiores al 2%, desarrollamos un proceso de control dividido en dos fases. En la primera calculamos la fuerza de presión, según la ecuación (11). En la segunda, para la vecindad de cada partícula, comprobamos cuales son las partículas cuya fuerza de presión viola la condición de incompresibilidad, se cuantifica el exceso de fuerza y rectificamos su dinámica y su posición. Además, también en esta segunda etapa, disipamos el calentamiento inducido por el exceso de fuerza. El procedimiento seguido se describe en el algoritmo 1. Algorithm 1 Proceso de control incompresibilidad. Require: Conocer todas las vecindades N(i). Ensure: La condición de incompresibilidad queda impuesta y el calentamiento disipado. Input: Conocer  FC j|i, ecuación (12) y ξ, ecuación (13). 1: for Cada partícula ido 2: Calcular cada  Fpi. 3: end for 4: for Cada partícula ido 5: for ∀j∈N(i)do 6: if (| Fpj|>| FC j|i|)then 7: Se asigna  Fpj← FC j|i 8: Se rectificarj←ξ( R/R) 9: Se calcula el calentamiento ∆Qj|i 10: end if 11: end for 12: end for La primera magnitud que calculamos es la fuerza máxima de compresión, esto es, la fuerza a partir de la cual se viola la condición de incompresibilidad. Para ello, particularizamos la ecuación c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. 63 J.J. Perea & J.M. Cordero / A new Potential-based SPH Method (11) con una distancia acorde con la densidad estática del fluido. La ecuación obtenida es:  FC j|i=−5.25 6πh(3.5π2 ζ(ξ−(0.85ε))2−4 (ξ−0.9ε)3) R |R|,(12) donde  FC j|ies la fuerza de compresión máxima que puede soportar la partícula jrespecto de iyξes la distancia máxima a la que se puede aproximar la partícula ja la partícula i, sin que se viole la incompresibilidad. La ecuación que satisface este criterio es: ξ=3 √3m 4π(1−η)ρ0,(13) donde ηes la tasa de compresibilidad y ρ0es la densidad estática. En el caso de que | FC j|i|>| Fpj|i|, tenemos que controlar la incompresibilidad. Para ello, modificamos la fuerza de la partícula j y su posición, hasta los valores máximos tolerados, es decir:  Fpj|i= FC pj|iyr′ j=ξ R R,(14) donde  Fpj|ies la fuerza de presión que tiene la partícula jrespecto de la partícula iyr′ jes su nueva posición. Esta rectificación es válida siempre que la diferencia entre |FC pj|i| y|FC pj|i|no supere el 10%. Aunque esta limitación parezca restrictiva, el hecho de que la propia fuerza de presión controle la incompresibilidad, permite que la diferencia entre las dos fuerzas no alcance a esta cifra. Por otro lado, el control establecido en las expresiones (14), genera un exceso de carga sobre el fluido, lo que induce un calentamiento que debe ser disipado. En nuestro modelo, podemos cuantificar este calentamiento y disiparlo en forma de calor. De este modo, ganamos estabilidad principalmente en las discontinuidades y satisfacemos la conservación de la energía y el momento. La ecuación propuesta para disipar el calentamiento inducido por el exceso de carga es: ∆Qj|i=(Fpj|i−FC pj|i)R.(15) Aunque el proceso que acabamos de describir muestra ciertas similitudes con otras técnicas que imponen explícitamente la incompresibilidad, como es el caso de PCISPH [SP09], presenta tres importantes diferencias. La primera es que es un proceso que se realiza en sólo dos fases y en la segunda fase sólo se opera sobre un 20% del total de las partículas, como máximo. Esto representa una mejora en la eficiencia respecto de los modelos basados en predicción–corrección ya que, además de requerir más de dos fases, el procesamiento se efectúa sobre el 100% de las partículas. La segunda diferencia es que no necesitamos utilizar la EOS, con lo que evitamos los problemas de eficiencia que su uso implica [ICS∗13]. La tercera diferencia es que podemos cuantificar el calentamiento que sufre el fluido durante la compresión y disiparlo, mejorando la estabilidad como se indica en [Sod78,SL08,PR14]. 4.3. Cálculo de la Densidad Una vez que todas las partículas satisfacen la condición de incompresibilidad, calculamos la densidad a partir de la fuerza de presión. Para ello, nos basamos en la formulación fundamental de SPH, ecuación (5), particularizada para la densidad. En esencia, partimos de la misma premisa establecida en el método PCISPH por [SP09], según la cual, los incrementos en la densidad sólo están relacionados con los cambios en la posición de las partículas y estos cambios son debidos a la fuerza de presión. Sin embargo, nuestro planteamiento presenta dos diferencias respecto del PCISPH: la primera es que relacionamos la fuerza y los cambios de posición mediante un truncado basado en el método de Newmark [New59], con el conseguimos una precisión del orden de O(n5); la segunda es que no necesitamos un proceso predictor–corrector para calcular la densidad, ya que la fuerza de presión de ninguna partícula violará la incompresibilidad. En base a este planteamiento, la nueva densidad en cada paso de tiempo, ρn+1 i, cumplirá: ρn+1 i=m∑ j∈N(i) W((rn i+∆ri)−(rn j+∆rj)),h),(16) donde mes la masa de las partículas y ∆ri,∆rjson los desplazamiento de la partícula iyj, respectivamente. Según la hipótesis inicial, los desplazamientos de la ecuación (16) están relacionados con la fuerza de presión. Para fijar esta relación, utilizamos el desarrollo en serie de Taylor para la aceleración y la posición, que truncaremos basándonos en el método de Newmark, es decir: ¨rn+1=¨rn+... xn+˙ ... rn 2! ∆t2+¨ ... rn 3! ∆t3+... ⇒ ¨rn+1−¨rn ∆t=... rn+˙ ... rn 2! ∆t+3! 5! ¨ ... rn∆t2+... (17) rn+1=rn+˙rn∆t+¨rn 2!∆t2+ ... rn 3! ∆t3+˙ ... rn 4! ∆t4+¨ ... rn 5! ∆t5+... ⇒ rn+1−rn=˙rn∆t+¨rn 2!∆t2+∆t3 3! (... rn+3! 4! ˙ ... rn∆t+3! 5! ¨ ... rn∆t2+...) (18) Comparando las ecuaciones 17 y18, podemos aproximar el incremento cómo: ∆r=rn+1−rn=˙rn∆t+¨rn 2!∆t2+∆t3 3! β1(¨rn+1−¨rn ∆t)⇒ ∆r=˙rn∆t+(3−β1 6)¨rn∆t2+β1∆t3Fp,(19) donde β1es el parámetro de interpolación y Fpes el módulo de la fuerza de presión modificada para imponer la compresión. c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. 64 J.J. Perea & J.M. Cordero / A new Potential-based SPH Method Particularizando la ecuación (19) para cada partícula obtenemos ∆ri. Sustituyendo este valor en la ecuación (16) es posible calcular la densidad que, al obtenerla a partir de la fuerza de presión, cumple la condición de incompresibilidad. 4.4. Fuerza de Viscosidad El principal problema asociado a la fuerza de viscosidad, desde el punto de vista del realismo de la simulación, es el sobreamortiguamiento inducido para de evitar oscilaciones aleatorias de las partículas. A ello hay que añadir que, por tratarse de una fuerza, debe satisfacer la conservación del momento. Existen diversos estudios que proponen soluciones a este problema. Cabe destacar la propuesta realizada por [ZCGL10], según la cual, el laplaciano de la velocidad se calcula partiendo de la formulación simétrica establecida en la ecuación (8). A esta formulación [ZCGL10] añade un factor de viscosidad artificial con el que se reducen las oscilaciones aleatorias de la partículas. Esta formulación consigue simulaciones estables, sin embargo, la complejidad del término de viscosidad la hace poco adecuada para simulaciones en CG. Además, si la velocidad de las partículas es baja, no se consigue mitigar el sobreamortiguamiento. Una propuesta similar es realizada por [AEF11]. Presenta un término de viscosidad dependiente de las velocidades relativas de las partículas, con lo que se optimiza la actuación del término de amortiguamiento sólo cuando la velocidad relativa es alta. Aunque los resultados obtenidos son realistas, la ecuación obtenida no es simétrica, con lo que no se cumple la conservación del momento. Además, se muestra ineficiente cuando la velocidad relativa de las partículas es baja. Con el objetivo de obtener una fuerza de viscosidad que impida el sobreamortiguamiento y satisfaga la conservación del momento presentamos nuestra formulación. Se trata de una ecuación simétrica obtenida a partir de la ecuación (8), sobre la que imponemos una relación de amortiguamiento. Conceptualmente, presenta ciertas similitudes con la propuesta de [ZCGL10] pero evitando las carencias que presenta cuando las partículas se desplazan a baja velocidad. Además, vamos a eliminar las posibles oscilaciones irreales de las partículas, sin que se vea afectada la conservación del momento. Nuestra expresión de la fuerza de viscosidad es: µ∇2vi=ρi[∑ j∈N(i) mjµ(vj ρ2 j +vi ρ2 i)∇2W(R,h)]χ(20) χ=exp −2|vj−vi| 3 √h√ρiρjR .(21) Analizando la ecuación (20), podemos comprobar que se cumple la conservación del momento, ya que la ecuación es simétrica y el término de amortiguamiento interviene multiplicando, lo que no afecta a dicha simetría. Además, centrándonos en la formulación del exponente, ecuación (21) podemos destacar que actúa aumentando el amortiguamiento cuando la velocidad relativa y la distancia crece. Por otro lado, como en el exponente se compensan mutuamente la distancia y la velocidad relativa, su “acción” se reduce al disminuir la distancia, evitando el sobreamortiguamiento. Otra característica relevante es el hecho de que haparezca en el denominador, lo que favorece la interacción entre partículas dentro del dominio soportado y aumenta el amortiguamiento cuando se aproximan al limite, es decir cuando R→h, pero este incremento se realiza en orden de 3 √h, con lo que reducimos la posibilidad del sobreamortiguamiento. Así pues, en base a este análisis, podemos afirmar que hemos formulado un amortiguamiento selectivo que permite la interacción a corta distancia y controla las dispersiones en distancias del rango de h, amortiguándolas. De este modo, favorecemos la cohesión propia del fluido y evitamos inestabilidades. 5. Resultados En esta sección mostraremos los resultados obtenidos utilizando el modelo propuesto. Desarrollaremos una comparativa entre estos resultados y los obtenidos sin implementar nuestras mejoras. Con ello, buscamos establecer un marco comparativo en el que analizar las mejoras que se obtienen al implementar nuestro modelo. Para cuantificar la estabilidad, hemos implementado el experimento de Sod [Sod78], ya que existe un consenso en considerarla como una prueba esencial para medir la estabilidad que se obtiene con un método de simulación de fluidos [Sod78,SL08,PR14]. Por otro lado, para mostrar el realismo conseguido, hemos implementado un líquido cuya tasa de compresibilidad no excede el 1%, por lo que puede considerarse un fluido incompresible. Finalmente, hemos desarrollado un conjunto de pruebas con diferentes valores en los parámetros de simulación. Con ello queremos mostrar la versatilidad que se consigue con nuestra propuesta. En cada una de las pruebas hemos utilizado el método semi– implícito de Euler para la integración temporal. Para evaluar la vecindad de cada partícula, hemos utilizado una técnica basada en el algoritmo NNS, según se describe en [THM∗03]. Para la interacción fluido–sólido hemos utilizado el modelo descrito por [SBC∗11]. La función kernel que hemos utilizado es el Q–spline, aunque puede utilizarse cualquier otra función kernel. Para llevar a cabo cada una de las simulaciones, hemos utilizado un equipo con un procesador Intel Core i3-4130, con 4GB de memoria RAM y equipado con una tarjeta gráfica NVIDIA GTX780. La primera prueba que hemos implementado es el experimento de Sod. Consiste en la expansión de un fluido sin viscosidad desde una región de alta a otra de baja presión que inicialmente están separadas por una membrana no porosa. Cuando la membrana es retirada se genera una onda de choque que evoluciona desde la región de alta hacia la de baja presión. Durante el proceso de expansión, se generan cuatro regiones con diferentes valores de presión, densidad y energía térmica, aunque las zonas de mayor relevancia, desde el punto de vista de cuantificar la estabilidad, son la localización de la onda de choque y de la cola de rarefacción. Una descripción más precisa de esta fenomenología se desarrolla en [SL08,PR14]. Para llevar a cabo la prueba, hemos utilizado un dominio [-0.5,0.5]. Inicialmente la región de alta presión está comprendida entre [-0.5,0.0)y la de baja presión (0.0,0.5]. Por otro lado, los valores iniciales de densidad, presión y velocidad son (ρh,ph,vh) = (1.0,1.0,0.0) y (ρl,pl,vl) = (0.125,0.1,0.0), donde los subíndices hylhacen referencia a las regiones de alta y baja presión, respectivamente. El número de partículas utic 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. 65 J.J. Perea & J.M. Cordero / A new Potential-based SPH Method (a) (b) (c) Figure 2: Graficas de la solución exacta y de los resultados numéricos de la densidad, presión y energía térmica obtenidos en el experimento de Sod. Los resultados numéricos se han obtenido implementando el método SPH clásico y el método híbrido potencial–SPH propuesto. lizado ha sido nh=1200 y nl=300, la masa de las partículas m=4.3910−3, el coeficiente de viscosidad µ=0.0 y el paso de tiempo ∆t=4.010−4, según se sugiere en [PR14]. A partir de estos valores iniciales hemos obtenido los resultados mostrados en la Figura 2en el instante t=0.15s. Por otro lado, para mostrar el realismo hemos implementado un fluido mediante 38500 partículas, en cada una de las cuales m= 0.03, h=0.18. Los otros valores utilizados son: ρ0=1000, ε=0.6 yη=1%. Los resultados obtenidos los mostramos en la Figura 3. Finalmente, para cuantificar los rangos de estabilidad que se obtinen con nuestra propuesta hemos llevado a cabo diversas pruebas. Los valores que hemos obtenido se muestran en la Tabla 1. 6. Conclusiones En este artículo hemos llevado a cabo una revision de los factores que afectan a la estabilidad y al realismo en la simulación de fluidos mediante sistema de partículas. A partir de estos factores, hemos propuesto un modelo híbrido potencial–SPH con el que evitamos los problemas de estabilidad que presentan las simulaciones mediante SPH. Particularmente hemos conseguido: 1. Mejorar la estabilidad en la fuerza de presión, formulando una fuerza independiente de la densidad. Así, hemos obtenido resultados estables en todas las regiones del fluido, incluyendo a los contornos. Este hecho se puede deducir a partir de los resultados mostrados en las gráficas de la Figura 2. Concretamente, en la discontinuidad de contacto, situada entorno a x=0.12 y en la cola de rarefacción, localizada alrededor de x=−0.25. En esta posición, la menor amplitud y longitud de la oscilación la obtenemos al implementar nuestro modelo. Esta mejora en la estabilidad, también puede visualizarse en los resultados de la simulación del fluido, Figura 3, principalmente en los contornos. 2. Imponer eficientemente la incompresibilidad y determinar el calentamiento en el fluido. Los resultados de esta imposición quedan reflejados en las imágenes de simulación del fluido, Figura 3, especialmente en los contornos. En referencia al calentamiento disipado, la mejora obtenida puede cuantificarse en las gráficas de la Figura 2y más concretamente en la cola de rarefacción, entorno a x=−0.25. En esta zona, la amplitud de la oscilación se asocia con el calentamiento, que es disipado en nuestro modelo. 3. Una fuerza de viscosidad que favorece la cohesión sin que ello induzca un sobreamortiguamiento. Esto se puede comprobar en las imágenes de simulación de la Figura 3. 4. Desarrollar un modelo estable y eficiente para un amplio rango de valores. Esta afirmación queda confirmada por los resultados mostrados en la Tabla 1. Especial interés tiene el rango del número de partículas vecinas en la que podemos operar y los pasos de tiempo en los que el sistema es estable. Agradecimientos Esta investigación ha sido apoyada por loa proyectos Megus (TIN2013-46928-C3-3-R) y el SoftPLM Network (TIN2015- 71938-REDT) del Ministerio de Econocmía y Competitividad de España. References [Abe11] ABEL T.: rpsph: A novel smoothed particle hydrodynamics algorithm. Monthly Notices of the Royal Astronomical Society 413, 1 (May 2011), 271–285. doi:10.1111/j.1365-2966.2010.18133.x. 2 [AEF11] AMINI Y., EMDAD H., FARID M.: A new model to solve fluidhypo-elastic solid interaction using the smoothed particle hydrodynamics (sph) method. European Journal of Mechanics - B/Fluids 30, 2 (Oct 2011), 184–194. 7 [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 [ATO15] ALDUÁN I., TENA A., OTADUY M. A.: Efficient and robust position-based fluids for vfx. In Proc. of Congreso Español de Informática Gráfica (2015). 2 [BT07] BECKER M., TESCHNER M.: Weakly compressible sph for free surface flows. In Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation (2007), SCA ’07, pp. 209–217. 2,3 [CBP05] CLAVET S., BEAUDOIN P., POULIN P.: Particle-based 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.2 c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. 66 J.J. Perea & J.M. Cordero / A new Potential-based SPH Method (a) (b) (c) (d) Figure 3: Secuencias de simulación en el instante t=1.25spara distintos pasos de tiempo: (a) ∆t=0.001s, (b) ∆t=0.008s, (c) ∆t=0.02s y (d) ∆t=0.06s. Potencial–SPH Hibrido Sin Compresibilidad SPH Radio Dominio Soportado (h) 0.06 – 0.8 0.2 – 0.45 0.2 – 0.35 Densidad (ρ) 125 – 2500 800 – 1800 800 – 1200 Numero de Particulas Vecinas (N(i)) 6 – 85 20 – 85 50 – 85 Viscosidad (µ) 2·10−3– 0.8·1032·10−2– 0.8·1032·10−2– 0.6·103 Masa de la Particula (m) 1·10−4– 1·10−12·10−3– 1·10−20.02 – 0.05 Paso de Tiempo (∆t) 6·10−2– 1·10−51·10−6– 1·10−51·10−6– 1·10−5 Table 1: Rangos de valores de estabilidad. Las magnitudes se expresan en unidades del sistema internacional. [CGSR00] CABEZÓN R. M., GARCIA-SENZ D., RELANO A.: A one– parameter family of interpolating kernels for smoothed particle hydrodynamics studies. Journal of Computational Physics 227, 19 (Oct 2000), 8523–8540. doi:10.1016/j.jcp.2008.06.014.2,4 [CM99] CUMMINS S. J., M. R.: An sph projection method. Journal of Computational Physics 152 (Jul 1999), 584–607. doi:10.1006/jcph.1999.6246.2 [DA12] DEHNEN W., ALY H.: Improving convergence in smoothed particle hydrodynamics simulations without pairing instability. Monthly Notices of the Royal Astronomical Society 425, 2 (Apr 2012), 1068–1082. doi:10.1111/j.1365-2966.2012.21439.x.2,4,5 [DRI97] DYKA C. T., RANDLES P. W., INGEL R. P.: Stress points for tension instability in sph. International Journal for Numerical Methods in Engineering 40, 13 (Jul 1997), 2325–2341. doi:10.1002/(SICI)1097-0207(19970715)40.2 [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. 3 [FSJ01] FEDKIW R., STAM J., JENSEN H.: Visual simulation of smoke. c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. 67