Aplicaciones del método smoothed particle hydrodinamics (SPH) en la dinámica de gases

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Aplicaciones del método smoothed particle hydrodinamics (SPH) en la dinámica de gases Autor: José Ignacio Pérez Puig Tutor: Miguel Pérez-Saborid Sánchez-Pastor Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 2 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Aplicaciones del método smoothed particle hydrodinamics (SPH) en la dinámica de gases Autor: José Ignacio Pérez Puig Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor titular Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 4 Trabajo Fin de Grado: Aplicaciones del método smoothed particle hydrodinamics (SPH) en la dinámica de gases Autor: José Ignacio Pérez Puig Tutor: Miguel Pérez-Saborid Sánchez-Pastor El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2017 El Secretario del Tribunal 6 A mi familia, amigos y profesores. 8 Agradecimientos Deseo expresar mi agradecimiento a todas las personas que me han acompañado durante estos años, en especial a mis padres por su apoyo incondicional y a mis amigos, los cuales han conseguido hacer el camino más fácil disfrutando de cada instante. También me gustaría agradecer a los profesores que me han formado durante estos años, especialmente a Don Miguel Pérez-Saborid Sánchez-Pastor por el gran apoyo recibido durante los últimos meses. José Ignacio Pérez Puig Sevilla, 2017 16 1 INTRODUCCIÓN as simulaciones mediante ordenadores se han convertido, durante estos años, en una herramienta muy importante en la resolución de problemas complicados dentro del ámbito de la ingeniería y la ciencia. La gran mayoría de los problemas fluido-mecánicos no tienen solución analítica y es necesario la utilización de métodos numéricos para su estudio. Estas simulaciones desempeñan un papel importante, sirviendo de prueba para numerosas teorías, que de otra manera serían muy difíciles de demostrar. Métodos numéricos basados en mallados tales como el método de las diferencias finitas, el método de los elementos finitos… han sido aplicados ampliamente en diferentes áreas de la dinámica de fluidos computacional a la vez que en la mecánica de sólidos. Una de las características más notables de los métodos basados en mallas, es la división de un medio continuo en pequeños subdominios discretos, mediante un proceso de discretización o mallado. De esta manera, las ecuaciones que gobiernan el problema pueden ser convertidas en un conjunto de ecuaciones algebraicas, donde las incógnitas corresponden a los valores de las variables en los nodos. A pesar de la gran utilidad, estos métodos numéricos presentan numerosas dificultades en otros aspectos, limitando así su aplicación en ciertos tipos de problemas. Las mayores dificultades se presentan en el uso del mallado, el cual puede llevar a lidiar con problemas en el uso de superficies libres, contorno deformable, interfaz móvil, grandes deformaciones… Por otro lado, en problemas con geometrías complicadas, la generación de una malla de calidad se ha convertido en una tarea complicada, costosa que requiere bastante tiempo. Figura (1-1) L “Lo que sabemos es una gota de agua; lo que ignoramos es el océano” -IsaacNewton - Estas dificultades que presentan los métodos son especialmente evidentes en la simulación de fenómenos hidrodinámicos tales como explosiones e impactos a grandes velocidades. Durante el proceso de la explosión existen grandes deformaciones, interfaces móviles, superficies libres… En los impactos a grandes velocidades, además, tiene lugar la propagación de ondas de choque. Por último, los métodos numéricos basados en mallas no son recomendables para situaciones donde el objeto de estudio es un conjunto discreto de partículas físicas en lugar de un medio continuo. Un ejemplo típico es el correspondiente a la interacción de estrellas en la astrofísica, movimiento de millones de átomos en estado de equilibrio o inequilibrio. Durante los últimos años, se ha investigado en una nueva generación de métodos computacionales, métodos libres de malla. La idea de estos métodos es proporcionar soluciones precisas y estables a las ecuaciones diferenciales que describen el movimiento, con cualquier tipo de condiciones de contorno. Uno de estos métodos es el Smoothed Particle Hydrodynamics (SPH), objeto de estudio en este proyecto. 1.1 Historia SPH El método SPH originalmente fue inventado para la resolución de problemas astrofísicos de tres dimensiones (Lucy 1977; Gingold and Monaghan,1977) en el espacio abierto. Dado que el movimiento colectivo de dichas partículas es similar al del movimiento de un líquido o un gas, se puede modelar mediante las ecuaciones clásicas de Newton para la hidrodinámica. Figura (1-2) Aunque los métodos tradicionales basados en mallas, como el método de las diferencias finitas o el método de los elementos finitos ya existían, y en general se encontraban más desarrollados que el SPH, presentaban dificultades a la hora de manejar complejos fenómenos como se ha visto antes. Esta circunstancia motivó a los investigadores a buscar alternativas para la resolución de este tipo de problemas, convirtiéndose el método SPH en una buena opción. A día de hoy, SPH es usado en una gran cantidad de áreas, desde simulaciones de estrellas binarias, supernovas, hasta el estudio de la evolución del universo. 18 1.2 Ventajas del SPH frente a otros métodos En el método SPH el sistema de estudio es representado por un conjunto de partículas, que tienen propiedades de un material e interactúan entre ellas dentro de un rango controlado por una función de peso o función kernel. La discretización de las ecuaciones que gobiernan el movimiento se basa en este sistema de partículas discretas. Se emplean numerosas formulaciones basadas en partículas para el cálculo de la densidad, velocidad y aceleración del fluido. La presión se obtiene mediante la ecuación de estado, de esta manera la aceleración se obtiene del gradiente de presiones y la densidad. Para el caso de fluidos viscosos, los efectos viscosos sobre la aceleración de la partícula también son incluidos. Como método Lagrangiano basado en partículas, SPH conserva exactamente la masa. El movimiento del fluido es representado por el movimiento de dichas partículas. El SPH presenta algunas ventajas especiales sobre los métodos numéricos tradicionales basados en mallas, anteriormente presentados. Estas ventajas se pueden resumir en los siguientes puntos: • Se trata de un método ideal para el modelado de superficies libres y problemas fluido dinámicos con interfaces. Las superficies libres, contornos móviles e interfaces de material puedes ser monitorizadas durante la simulación sin importar el complejo movimiento de las partículas. Figura (1-3) • SPH es un método como ya se ha comentado que no emplea un mallado. Este hecho permite un manejo directo y sencillo de las grandes deformaciones, ya que la conectividad entre las partículas es generada como parte de simulación y puede variar a lo largo del tiempo. • En el método, la partícula representa un volumen finito en la escala continua. Similar al método de dinámica molecular que emplea el uso de partículas para representar un átomo o molécula en nanoescalas. Por tanto, es natural extender el uso del SPH a escalas menores, especialmente en ámbitos como la biofísica o biomedicina. • SPH es apropiado para para problemas donde el objeto bajo consideración no es continuo. Por ejemplo, en bioingeniería, nano-ingeniería o astrofísica. Figura (1-4) • La implementación numérica es mucho más simple y el estudio de sistemas tridimensionales es mucho más natural que en el caso de los métodos tradicionales. 1.3 Objetivos y Estructura del trabajo El Objetivo de este Trabajo de fin de Grado es doble: presentar las bases teóricas del método SPH e implementar el método en MATLAB para resolución de diversas aplicaciones de dinámica de gases. Se trata del primer trabajo realizado de este método en la Escuela Superior de Ingeniería de Sevilla y pretende ser un punto de partida para futuros trabajos e aplicaciones del SPH. El método se aplicará en el estudio del problema unidimensional de ondas de choque en un tubo donde existe una discontinuidad en las propiedades físicas (Shock Tube). Además del estudio de la evolución con respecto al tiempo de un modelo simplificado de estrella y como alcanza el estado de equilibrio. Con estos objetivos por cumplir, el trabajo se estructurará en seis capítulos más, además del presente. En el siguiente capítulo se van a exponer los fundamentos teóricos sobre los que se basa el método SPH a implementar. Seguidamente, se explicará cómo se formulan las ecuaciones de Euler dentro del método. En los capítulos 4 y 5, se analizarán los problemas del tubo de choque y modelo simplificado de estrella, presentando los resultados obtenidos en las simulaciones junto a los códigos empleados. Finalmente, las conclusiones y futuros desarrollos serán expuestos en el capítulo 6. 20 2 FUNDAMENTOS TEÓRICOS “Las matemáticas son el lenguaje en el que Dios ha escrito el Universo” -Galileo Galilein el capítulo 1, se ha visto que el método Smoothed Particle Hydrodynamics (SPH) presenta características diferentes a los métodos tradicionales empleados en la resolución numérica de problemas dinámica de fluidos. En este capítulo se presentarán los conceptos y fundamentos teóricos esenciales para la formulación del método. Se analizará la estrategia del método SPH en la aproximación de variables mediante la representación integral de una función y la discretización en partículas. También se expondrán técnicas usadas en la derivación de las formulaciones SPH e interpolación en puntos diferentes a las partículas. 2.1 Aproximación SPH La gran mayoría de problemas hidrodinámicos se encuentran en forma de ecuación en derivadas parciales (EDP) de variables tales como densidad, velocidad, energía… La obtención de una solución analítica para este conjunto de ecuaciones raramente es posible, exceptuando casos simplificados. Para la resolución, primero se discretiza el dominio donde las EDPs están definidas. Una vez realizado, es necesario un método que facilite una aproximación de las funciones y sus derivadas en cualquier punto del dominio. Esta aproximación es aplicada al conjunto de EDPs llegando a ecuaciones diferenciales ordinarias (EDO) en forma discreta. Por último, se resuelven usando algún método de integración basado en diferencias finitas. En el caso del método SPH, el dominio se representa mediante un conjunto de partículas arbitrariamente distribuidas. Para la obtención de la formulación SPH existen básicamente dos pasos. El primer paso es representar una función y sus derivadas como una integral, este paso se suele denominar aproximación mediante la función kernel. El segundo paso consiste en la aproximación mediante partículas. En este paso, el conjunto de EDOs es resuelto mediante la discretización del dominio con una distribución inicial de partículas donde se aproximan las variables. 2.1.1 Aproximación mediante función kernel La representación integral de una función consiste en la aproximación mediante la suma de los valores correspondientes a las partículas cercanas, reduciendo la aproximación de la función a un punto o partícula. El concepto de la representación integral de una función 𝑓(𝒙) comienza con la siguiente identidad 𝑓(𝒙)=∫𝑓(𝒙′)𝛿(𝒙−𝒙′)𝑑𝒙′ Ω, (1) E donde 𝑓 es función del vector de posición tridimensional 𝒙. Por otro lado, 𝛿(𝒙−𝒙′) es la función delta de Dirac dada por la siguiente igualdad 𝛿(𝒙−𝒙′)={1 𝒙=𝒙′ 0 𝒙≠𝒙′. (2) La ecuación (1) implica que una función puede ser representada en forma integral. Debido al uso de la función delta de Dirac, la integral representada en la ecuación (1) es exacta, siempre y cuando 𝑓(𝒙) este definida y sea continua en el volumen Ω. Si remplazamos la función delta de Dirac por una función 𝑊(𝒙−𝒙′,ℎ), la representación integral de 𝑓(𝒙) viene dada por 𝑓(𝒙)=∫𝑓(𝒙′)𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′, Ω (3) donde W es la llamada función kernel. Dentro de la función, h es el parámetro llamado smoothing length, el cual define el área de influencia de la función kernel. Una vez se ha sustituido la función delta de Dirac, la representación integral de la ecuación (3) solo puede ser una aproximación. Por tanto, la ecuación se reescribe como <𝑓(𝒙)>=∫𝑓(𝒙′)𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′. Ω (4) La función 𝑊 normalmente es una función par, por razones que se darán posteriormente. Además, debe satisfacer un cierto número de condiciones, la primera de ellas es la condición de normalización ∫𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′ Ω=1. (5) Esta condición es también llamada la condición unidad, dado que la integración de la función kernel en el dominio ha de ser igual a la unidad. La segunda condición es la propiedad de delta de Dirac cuando el smoothing length se aproxima a cero lim ℎ→0𝑊(𝒙−𝒙′,ℎ)=𝛿(𝒙−𝒙′). (6) La tercera condición es la condición de compacidad 𝑊(𝒙−𝒙′,ℎ)=0 𝑠𝑖 |𝒙−𝒙′|>𝑘ℎ, (7) siendo 𝑘 es una constante relacionada con la función kernel para el punto 𝒙 y que define el área efectiva de la función. Esta área efectiva es llamada el dominio de soporte de la función kernel en el punto 𝒙. Usando esta condición de compacidad, la integración sobre el dominio completo se reduce a la integración sobre el dominio de soporte de la función kernel. 22 Esta aproximación tiene una precisión de orden ℎ2 o de segundo orden. Los errores en la representación integral se pueden estimar mediante el uso de una serie de Taylor de 𝑓(𝒙′) en torno a 𝒙. Si 𝑓(𝒙) es derivable y continua, tenemos <𝑓(𝒙)≥∫ [𝑓(𝒙)+𝑓′(𝒙)(𝒙′−𝒙)+𝑟((𝒙′−𝒙)2)]×𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′ Ω=𝑓(𝒙)∫𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′ Ω +𝑓′(𝒙)∫(𝒙′−𝒙)×𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′+𝑟(ℎ2). Ω (8) Si 𝑊 es una función par con respecto a 𝒙 y (𝒙′−𝒙), entonces 𝑊(𝒙−𝒙′,ℎ) debe ser una función impar. Por tanto, se obtiene que ∫(𝒙′−𝒙)𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′ Ω=0. (9) Usando la condición de unidad, y las anteriores ecuaciones llegamos a la siguiente expresión <𝑓(𝒙)>=𝑓(𝒙)+𝑟(ℎ2). (10) Quedando claro que la aproximación mediante el kernel en el método SPH de una función arbitraria tiene una precisión de segundo orden. La aproximación de una derivada espacial ∇·𝑓(𝒙) se obtiene simplemente sustituyendo 𝑓(𝒙) por ∇·𝑓(𝒙) en la ecuación (4), llegando a <∇·𝑓(𝒙)>=∫[∇·𝑓(𝒙′)]𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′. Ω (11) La divergencia en la integral es con respecto a la coordinada prima. Considerando la siguiente igualdad [∇·𝑓(𝒙)]𝑊(𝒙−𝒙′,ℎ)=∇·[𝑓(𝒙′)𝑊(𝒙−𝒙′,ℎ)]−𝑓(𝒙′)· ∇[𝑊(𝒙−𝒙′,ℎ)], (12) se obtiene la siguiente ecuación ∇·𝑓(𝒙)=∫∇·[𝑓(𝒙′)𝑊(𝒙−𝒙′,ℎ)]𝑑𝒙′ Ω−∫𝑓(𝒙′)·∇[𝑊(𝒙−𝒙′,ℎ)]𝑑𝒙′. Ω (13) La primera integral puede ser transformada mediante el teorema de la divergencia en una integral sobre la superficie 𝑆 del dominio de integración Ω <∇·𝑓(𝒙)>=∫𝑓(𝒙′)𝑊(𝒙−𝒙′,ℎ)·𝒏 󰇍 󰇍 𝑑𝒙′ S−∫𝑓(𝒙′)·∇𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′ Ω, (14) donde 𝒏 󰇍 󰇍 es el vector normal de la superficie 𝑆. Desde que la función kernel es normalmente definida para tener soporte compacto, el valor de 𝑊 en la superficie de la anterior integral es cero en SPH. Por consiguiente, la aproximación mediante el kernel de las derivadas se puede escribir como <∇·𝑓(𝒙)>=−∫𝑓(𝒙′)·∇𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′ Ω. (15) La aproximación de la derivada en el SPH permite que el gradiente sea determinado a partir de los valores de la función y de las derivadas de la función kernel, en lugar de los valores de la derivada de la propia función. La aproximación de derivadas de mayor orden se obtiene de una manera similar, sustituyendo 𝑓(𝒙) con la correspondiente derivada en la ecuación (4), usando integración por partes, el teorema de la divergencia y transformaciones triviales. Siguiendo análisis similares basados en series de Taylor, es fácil demostrar que la aproximación de las derivadas es también de segundo orden de precisión. Las aproximaciones mediante kernel en el método SPH son de segundo como ver en la siguiente figura (2-1). Figura (2-1) Sin embargo, para algunos escenarios no se cumple este hecho. Por ejemplo, problemas donde el dominio de soporte intersecta el dominio del problema. En estos casos la función 𝑊 es truncada por el contorno y la integración de la función kernel en el dominio ya no es igual a la unidad. Figura (2-2) 24 2.1.2 Aproximación mediante partículas En el método SPH, el sistema al completo es representado mediante un número finito de partículas, donde cada una de ellas tiene masa propia y ocupa una zona del espacio. La representación integral mediante la función kernel se puede convertir a una forma discreta siendo la suma del valor de las partículas en el dominio. Esta discretización es conocida como la aproximación mediante partículas en el método SPH. Considerando un dominio Ω formado por un conjunto de particulas, arbitrariamente distribuidas. El estado del sistema en cada instante es representado por estas particulas, cada una con sus propiedades. Estas particulas no son únicamente para la integración o interpolación, sino también para representar el material. Figura (2-3) La función 𝑊 es usada para la aproximación de las variables de la partícula 𝑖 usando la suma ponderada sobre las particulas 𝑗 ,dentro del dominio de soporte con una reducción de distancia 𝑘ℎ𝑖. El volumen infinitesimal 𝑑𝒙′ en las anteriores integrales de una partícula j, se remplaza por el volumen infinitesimal de la particula ∆𝑉𝑗, a través de la siguiente expresión: 𝑚𝑗=∆𝑉𝑗𝜌𝑗, (16) donde 𝜌𝑗 es la densidad de la partícula 𝑗(=1,2,…,𝑁) y 𝑁 es el número de partículas dentro del dominio de soporte de la partícula 𝑗. La representación integral de 𝑓(𝒙) puede escribirse ahora de forma discreta con la aproximación mediante partículas de la siguiente manera 𝑓(𝒙)=∫𝑓(𝒙′)𝑊(𝒙−𝒙′,ℎ)𝑑𝒙′ Ω≅∑𝑓(𝒙𝑗)𝑊(𝒙−𝒙𝒋,ℎ)∆𝑉𝑗 𝑁 𝑗=1 =∑𝑓(𝒙𝒋)𝑊(𝒙−𝒙𝒋,ℎ)𝑚𝑗 𝜌𝑗 𝑁 𝑗=1 . (17) Finalmente, la aproximación mediante partículas de una función en una partícula i puede ser escrita como: <𝑓(𝒙𝑖)>=∑𝑚𝑗 𝜌𝑗𝑓(𝒙𝒋)𝑊𝑖𝑗 𝑁 𝑗=1 , (18) donde 𝑊𝑖𝑗=𝑊(𝒙𝒊−𝒙𝒋,ℎ). (19) La ecuación (18) establece que el valor de una función en la posición de la partícula 𝑖 es aproximada usando la suma de los valores de la función en todas las particulas que se encuentran dentro del dominio de soporte de la partícula 𝑖, ponderadas por la función kernel. Siguiendo el mismo procedimiento, la aproximación mediante partículas de la derivada puede obtenerse como: <∇·𝑓(𝒙)>=−∑𝑚𝑗 𝜌𝑗𝑓(𝒙𝒋)·∇𝑊(𝒙𝒊−𝒙𝒋,ℎ) 𝑁 𝑗=1 , (20) donde el gradiente ∇𝑊 en la anterior ecuación es evaluado en la partícula 𝑗. La ecuación (20) establece que el valor del grandiente de una función en una partícula localizada en la posición 𝒙 puede ser aproximada usando la suma de los valores de la función en la posición de todas las partículas dentro del soporte de dominio, multiplicado al volumen de la partícula y ponderado por el grandiente de la función kernel. La aproximación mediante partículas convierte la forma continua de una función y sus derivadas en una suma discreta sobre un conjunto de partículas. El uso de estas sumas para la aproximación de la integral es un factor clave que hace el método SPH fácil de implementar. Otro aspecto importante en el método es que se introduce la masa y la densidad de las partículas dentro de las ecuaciones. Este hecho puede ser convenientemente aplicado a los problemas hidrodinámicos donde la densidad es una variable clave. Para la aplicación a la mecánica de sólidos, se necesitaría un tratamiento especial. 32 3.1.3 Viscosidad artificial A la hora de simular los problemas hidrodinámicos se requiere un tratamiento especial para permitir que los algoritmos sean capaces de modelar las ondas de choque. La aplicación de la conservación de la masa, momento y energía a través de la onda de choque requiere la simulación de la transformación de energía cinética en energía térmica. Físicamente, esta transformación de energía puede representarse como disipación viscosa. El modelo de viscosidad artificial desarrollado por Monaghan es el más usado en los algoritmos SPH. No solo proporciona la disipación necesaria para convertir la energía cinética en térmica, sino que también previene la penetración entre las partículas. A continuación, se describe el modelo de viscosidad utilizado. Se añade un termino Π𝑖𝑗 a la ecuación (39) 𝑑𝑣𝑖 󰇍 󰇍 󰇍 𝑑𝑡=−∑𝑚𝑗( 𝑗𝑃𝑗 𝜌𝑗2–𝑃𝑖 𝜌𝑖2+Π𝑖𝑗)·∇𝑖𝑊𝑖𝑗, (45) donde Π𝑖𝑗viene dado por: Π𝑖𝑗={−𝛼𝑐𝑖𝑗    𝜇𝑖𝑗+𝛽𝜇𝑖𝑗 2 𝜌𝑖𝑗     𝑣𝑖𝑗 󰇍 󰇍 󰇍 󰇍 ·𝑟𝑖𝑗 󰇍 󰇍 󰇍 <0 0 𝑣𝑖𝑗 󰇍 󰇍 󰇍 󰇍 ·𝑟𝑖𝑗 󰇍 󰇍 󰇍 >0. (46) Y la viscosidad 𝜇𝑖𝑗=ℎ𝑖𝑗𝑣𝑖𝑗 󰇍 󰇍 󰇍 󰇍 ·𝑟𝑖𝑗 󰇍 󰇍 󰇍 𝑟𝑖𝑗2+𝜂2. (47) Con 𝑣𝑖𝑗 󰇍 󰇍 󰇍 󰇍 =𝑣𝑖 󰇍 󰇍 󰇍 −𝑣𝑗 󰇍 󰇍 󰇍 y 𝑟𝑖𝑗 󰇍 󰇍 󰇍 =𝑟𝑖 󰇍 󰇍 −𝑟𝑗 󰇍 󰇍 . Las cantidades 𝜌𝑖𝑗     y 𝑐𝑖𝑗    hacen referencia a la media de la densidad y la velocidad del sonido y son calculadas de la misma manera que el smoothing length en la ecuación (41). Los parámetros 𝛼,𝛽,𝜂 son normalmente: 𝛼=1.0,𝛽=2.0,𝜂=0.1ℎ𝑖𝑗.. El termino asociado a 𝛼 representa la viscoisdad de volumen, mientras el segundo termino asociado con 𝛽 intenta impedir la penetración entre particulas con altos números de Mach. Por otro lado, el factor 𝜂 se introduce para prevenir divergencias numéricas cuando dos particulas se aproximan una a la otra. 3.1.4 Ecuación de la energía La ecuación de la energía (34) se convierte en formulación SPH aproximando la ecuación de la densidad con la ecuación (36) 𝑑𝑢𝑖 𝑑𝑡=𝑃 𝜌𝑖2∑𝑚𝑗(𝑣𝑖 󰇍 󰇍 󰇍 −𝑣𝑗 󰇍 󰇍 󰇍 )·∇𝑖𝑊𝑖𝑗 𝑗. (48) Si introducimos la viscosidad artificial (44), la ecuación de la energía se reescribe como: 𝑑𝑢𝑖 𝑑𝑡=𝑃 𝜌𝑖2∑𝑚𝑗(𝑣𝑖 󰇍 󰇍 󰇍 −𝑣𝑗 󰇍 󰇍 󰇍 )·∇𝑖𝑊𝑖𝑗 𝑗+12∑𝑚𝑗Π𝑖𝑗(𝑣𝑖 󰇍 󰇍 󰇍 −𝑣𝑗 󰇍 󰇍 󰇍 )·∇𝑖𝑊𝑖𝑗 𝑗. 3.2 Integración en el tiempo En las matemáticas, la condición de Courant es necesaria para la convergencia en la resolución numérica de ciertos tipos de problemas. El diferencial de tiempo entre cada iteración ha de ser lo suficientemente pequeño para que la simulación no produzca resultados incorrectos. Por ejemplo, en el caso de una onda moviéndose a través de una malla discreta si se quiere calcular la amplitud en periodos discretos de tiempo de igual duración, esta duración ha de ser menos que el tiempo que tarda la onda en viajar a la celda adyacente. Para la elección del diferencial de tiempo en cada iteración, Monaghan propone un método para asegurar la estabilidad del algoritmo. Es una combinación de dos condiciones, una condición de fuerza ∆𝑡𝑓=min 𝑖√ℎ |𝑎 𝑖|, (1.49) junto con la condición de Courant, la cual como se ha visto incluye condiciones de estabilidad ∆𝑡𝑐=min 𝑖ℎ 𝑐𝑖+0.6(𝛼𝑐𝑖+𝛽max 𝑗𝜇𝑖𝑗). (50) La viscosidad 𝜇𝑖𝑗 es dada por la ecuación (45) y los parámetros 𝛼,𝛽 son los mismos que en la ecuación (44). El diferencial de tiempo es tomado por la siguiente expresión ∆𝑡=0.25min(∆𝑡𝑓,∆𝑡𝑐). 34 4 TUBO DE CHOQUE “Si buscas resultados distintos, no hagas siempre lo mismo” -Albert Einstein – l problema de tubo de choque es una referencia muy utilizada para la evaluación de métodos numéricos. Muchos investigadores del SPH han estudiado este problema (Monaghan 1983, Hernquist 1989…), en el cual emplean el spline cúbico como función kernel. 4.1 Descripción del problema El tubo de choque consiste en un tubo largo y recto relleno de gas, que es separado por una membrana en dos zonas con diferente densidad y presión. El gas en cada zona se encuentra en equilibrio con presión, densidad y temperatura constante. Cuando la membrana es retirada, se produce una onda de choque, una onda de expansión o rarefacción en la zona de discontinuidad. La onda de choque se mueve hacia la región con menor densidad, mientras que la onda de expansión viaja a través de la zona con mayor densidad. La superficie de contacto se forma cerca del centro y se mueve dentro de la región de baja densidad detrás de la onda de choque. En las siguientes figuras se puedes apreciar esquemas de la evolución del fluido Figura (4-1) E Figura (4-2) Existen 4 zonas claramente diferenciadas como se puede ver en la figura (4-2). La zona 1 corresponde al gas de baja presión que no ha sido modificado por la onda de choque. Dentro de la zona 2, distinguimos entre 2 y 2′ por la superficie de contacto. Esta zona es la situada inmediatamente después de la onda de choque viajando a velocidad constante. A través de la superficie de contacto la densidad y la temperatura son discontinuas. La zona ubicada entre el inicio y el final de la onda de expansión corresponde a la zona 3, dentro de la cual las propiedades del fluido cambian gradualmente debido a que se trata de un proceso isentrópico. Por último, la zona 4 es la zona sin perturbar del gas a alta densidad. 4.2 Resultados En esta simulación las condiciones iniciales vienen dado por el siguiente conjunto de valores 𝑥≤0 𝜌=1.0 ,𝑃=1.0 ,𝑣=0 ,𝑢=2.5, (51) 𝑥>0 𝜌=0.125 ,𝑃=0.1 ,𝑣=0 ,𝑢=2.0. (52) Para la simulación del problema se utilizan 𝑁=400 particulas, donde 0.8∗𝑁 particulas se distribuyen uniformemente en la zona de alta densidad 𝑥∈[−0.6,0] mientras que el resto de particulas se ubican en la zona de baja densidad. Esta distribución de partículas se debe a que la densidad es ocho veces mayor en la zona izquierda que en la derecha de la membrana, por consiguiente, el espaciado debe ser 8 veces menor en la zona de alta densidad. La simulación del tubo de choque se resuelve mediante las ecuaciones de la densidad (37), aceleración (43) y la energía interna (47). La presión se obtiene mediante la ecuación de estado en cada intervalo de tiempo. 36 En la figura (4-3) se representa la evolución de la densidad vs la posición para el caso sin viscosidad artificial. Como era de esperar, se observan oscilaciones en la solución que tiende a ser inestable en la simulación de ondas de choque. Figura (4-3) Incluyendo los términos de viscosidad artificial se obtienen los siguientes perfiles: Figura (4-4) Figura (4-5) Figura (4-6) 38 Figura (4-7) 4.3 Códigos MATLAB Para acabar, se presentan los códigos empleados para la simulación del tubo de choque. En los códigos existe una función principal donde se realiza la simulación y que contiene las ecuaciones de la aceleración, energía, densidad y presión incluyendo los términos viscosos. Dentro de esta función se llama a las funciones kernel empleadas, en este caso Spline kernel. function shock_tube %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %------------------------FUNCIÖN PRINCIPAL TUBO DE CHOQUE------------------ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; clc; close all; %% Datos y configuración inicial Gam = 1.4; N = 400; %Número de particulas N1 = ceil(4*N/5); N2 = N-N1; %Zona de alta presión x(1:N1,1) = linspace(-0.6,0,N1); rho(1:N1,1) = 1; m(1:N1,1) = 0.6*rho(1:N1,1)/N1; v(1:N1,1) = 0; u(1:N1,1) = 1/1/(gam-1); %Zona de baja presión x((N1+1):N,1) = linspace(0,0.6,N2); rho((N1+1):N,1) = 0.125; m((N1+1):N,1) = 0.6*rho((N1+1):N,1)/N2; v((N1+1):N,1) = 0; u((N1+1):N,1) = 0.1/0.125/(gam-1); %Otros parámetros H = 1*(x(N,1)-x(N-1,1)); Alpha = 1.0; Beta = 2.0; Eta = 0.1*h; %% Itegración en el tiempo Ntime = 200; %Número máximo de iteraciones T = 0; for itime=1:Ntime %Primero calculamos las distancias relativas y funciones kernel for i=1:N for j=1:N, r(i,j) = x(i,1)-x(j,1); w(i,j) = spline_kernel(r(i,j),h,1); dw_dx(i,j) = spline_kernelp(r(i,j),h,1); end end % Calculo de la densidad rho(:,1)=w*m; % Calculo de la presion y velocidad del sonido p(:,1)=(gam-1)*rho.*u; c(:,1)=sqrt(gam*p./rho); %Calculo del gradiente de presiones y viscosidad for i=1:N, for j=1:N, f_p(i,j) = (p(i,1)/rho(i,1)^2+p(j,1)/rho(j,1)^2); cb(i,j) = 0.5*(c(i,1)+c(j,1)); rhob(i,j)= 0.5*(rho(i,1)+rho(j,1)); vdif(i,j)= v(i,1)-v(j,1); xdif(i,j)= x(i,1)-x(j,1); %viscosidad mu(i,j) = h*vdif(i,j)*xdif(i,j)/(xdif(i,j)^2+eta^2); Pi(i,j) = 0; if vdif(i,j)*xdif(i,j) < 0, Pi(i,j)=(-alpha*cb(i,j)*mu(i,j)+beta*mu(i,j)^2)/rhob(i,j); end end end 40 aux1(:,1) = (vdif.*dw_dx)*m; aux2(:,1) = 0.5*(Pi.*vdif.*dw_dx)*m; %Condición de courant Dtf = min(sqrt(h./c)); mu_max(:,1) = max(mu')'; dtc = min(h./(c+0.6*(alpha*c+beta*mu_max))); dt = 0.25*min(dtf,dtc); %Actualización de las variables T = T+dt; vt(:,1) = v(:,1)-dt*((f_p+Pi).*dw_dx)*m; ut(:,1) = u(:,1)+dt*((p./rho.^2).*aux1+aux2); xt = x+0.5*dt*(vt+v); rhot(:,1) = w*m; % Promediado Nr = 1000; r = zeros(Nr,1); r(:,1) = linspace(-0.6,0.6,Nr); rho_r = zeros(Nr,1); for i=1:Nr sum = 0; sumv = 0; sump = 0; sumu = 0; for j=1:N, uij = r(i,1)-x(j,1); sum = sum+m(j,1)*spline_kernel(uij,h,1); sumv = sumv+vt(j,1)*m(j,1)/rho(j,1)*spline_kernel(uij,h,1); sump = sump+p(j,1)*m(j,1)/rho(j,1)*spline_kernel(uij,h,1); sumu = sumu+ut(j,1)*m(j,1)/rho(j,1)*spline_kernel(uij,h,1); end rho_r(i,1)=sum; v_r(i,1) =sumv; p_r(i,1) =sump; u_r(i,1) = sumu; end % Representacion figure(1) plot (xt,rhot,'o') hold on plot(r(:,1),rho_r,'r') axis([-0.4 0.4 0 1.2]) pause(0.001) title('Densidad vs Posición') xlabel('x') ylabel('\rho') hold off %Actualizando valores t-1 v=vt; u=ut; rho=rhot; x=xt; end figure(2) plot(xt,vt,'o') hold on plot(r(:,1),v_r,'r') axis([-0.4 0.4 0 1.2]) title('Velocidad vs Posición') xlabel('x') ylabel('v') figure(3) plot(xt,p,'o') hold on plot(r(:,1),p_r,'r') axis([-0.4 0.4 0 1.2]) title('Presión vs Posición') xlabel('x') ylabel('P') figure(4) plot(xt,ut,'o') hold on plot(r(:,1),u_r,'r') axis([-0.4 0.4 1.2 3.2]) title('Energia interna vs Posición') xlabel('x') ylabel('u') % % function w=spline_kernel(r,h,d) %SPLINE_KERNEL- Smoothing kernel function, reduces computation time. % Format:w=gaussian_kernel(r,h,d) % d-dimension of the problem % h-kernel length q = norm(r)/h; ch = [1/(h),15/(7*pi*h^2),3/(2*h^3*pi)]; if 0<=q && q<1 w = 2/3-q^2+q^3/2; elseif 1<=q && q<2 w = 1/6*(2-q)^3; else w = 0; end w = ch(d)*w; 48 Figura (5-7) Figura (5-8) Figura (5-9) Figura (5-10) 50 Figura (5-11) De las primeras simulaciones se observa una serie de aspectos claves del método SPH en su comparación con los resultados analíticos. Ante todo, la reproducción de la solución analítica es bastante buena. La densidad central 𝜌(0) se aproxima mejor cuanto mayor sea el número de partículas 𝑁 empleado. Usando únicamente 𝑁=50 partículas la densidad obtenida es ligeramente inferior a la obtenida por la expresión analítica. Para corregir este hecho, se emplean más partículas y se puede observar que el valor obtenido con 𝑁= 100 y 𝑁=300 es muy cercano al valor buscado. Otro detalle a comentar es cómo en la representación del perfil de densidades, la aproximación se desvía de la solución analítica cuando 𝑟→𝑅 . En el modelo analítico la densidad converge a cero. En este caso la forma de la función kernel altera el perfil de densidad siendo algo difícilmente corregible en el método SPH. Variando el smoothing length podemos llegar a mejores resultados. Generalmente, una reducción de ℎ proporciona una mejora en la respuesta siempre y cuando existan suficientes partículas en los alrededores de cada partícula. Reducir h para obtener mejores respuestas implica aumentar el número de partículas. A continuación, se varía las funciones kernel presentadas en el trabajo, analizando su efecto en la evolución de la densidad central y el perfil de densidad. Cambiando entre Gaussian kernel y Spline kernel se obtienen las siguientes gráficas: Figura (5-12) Figura (5-13) 52 Aparentemente no se aprecia variación significativa entre la utilización del Gaussian kernel o el spline kernel en el perfil de densidad, como podemos ver en la Figura (5-13). Gingold & Monaghan recomiendan la utilización del Gaussian kernel para obtener mejores respuestas físicas, sin embargo, el spline kernel también es efectivo y realiza los cálculos más eficientemente como ya se ha explicado. Por último, se varía el valor empleado de la amortiguación: Figura (5-14) El incremento de 𝜈 reduce el tiempo que emplea la estrella en llegar al equilibrio, Figura (5-14). Sin amortiguación la estrella continuaría oscilando indefinidamente. Además, también se ha simulado la colisión de dos estrellas con una velocidad relativa de 𝑣=0.4, representando la posición de las partículas y el valor de la presión y densidad máximas como función del tiempo en las figuras (5-18) y (5-19). Ejemplo claro de como el método SPH puede ser usado para monitorizar la evolución de parámetros físicos en el tiempo que de otra manera serían muy difícil o imposible de calcular analíticamente. Figura (5-15) Figura (5-16) 54 Figura (5-17) Figura (5-18) Figura (5-19) 5.3 Códigos Matlab Para acabar, se presentan los códigos empleados para la simulación. En primer lugar, se presenta el programa utilizado para la simular como alcanza el equilibrio una estrella a partir de una configuración aleatoria de radio 𝑅. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %----------------------------Programa principal---------------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Modelo simplificado de Estrella 2D clear all close all rng('default') %% Datos N = 300; %Numero de particulas d = 2; %Dimension del problema h = 0.06/sqrt(N/1000); %Smoothing length dt = 0.04; %Time step nu = 3; %Coeficiente de viscosidad k = 0.1; %Constante de presion n = 1; %Indice politropico M = 2; %Masa de la estrella m = M/N; %Masa de cada partícula R = 0.75; %Radio de la estrella max_time_step = 400; 56 %Termino gravitacional lambda = 2*k*(pi^(-1/n))*((M*(1+n)/(R^2))^(1+1/n))/M; %% Creamos puntos random dentro de un circulo de radio R x = zeros(1,d); % Primer punto en el centro ii = 2; % Iniciamos el indice while length(x(:,1))<N y=-R+(R+R).*rand(1,d); if norm(y)^2<=R^2 x(ii,1) = y(1); x(ii,2) = y(2); ii = ii+1; end end %% Bucle %Inicializamos las variables v_mhalf = zeros(N,d); a = zeros(N,d); rho0 = zeros(N,1); %vector para almacenar el maximo valor de la densidad time = 0; for i=1:max_time_step disp(i) %Actualiacion posicion mediante el algorítmo Leap frog v_phalf = v_mhalf+a*dt; x = x+v_phalf*dt; v = 0.5*(v_mhalf+v_phalf); v_mhalf = v_phalf; %Actualizamos las variables rho = calculate_density(x,m,h,d,N); P = k*rho.^(1+1/n); %relacion de fluido politropico a = calculate_acceleration(x,v,m,rho,P,nu,lambda,h,d,N); rho0(i) = max(rho); %Animación %2D figure(1) scatter(x(:,1),x(:,2),[],rho,'filled') axis([-1 1 -1 1 ]) hold on colorbar plot(xc,yc,'r') pause(0.001) hold off end %% Evolución de rho_max figure(2) hold on plot(linspace(0,16,400),rho0,'r') xlabel('t') ylabel('\rho(0)') title('Evolución de la densidaad central como función de t') %legend('\nu=1','\nu=3') %% Representar rho(r) en el equilibrio mediante un promediado Nr = 1000; r = zeros(Nr,d); r(:,1) = linspace(0,R,Nr)'; rho_r = zeros(Nr,1); for i=1:Nr sum=0; for j=1:N, uij = r(i,:)-x(j,:); sum = sum+m*spline_kernel(uij,h,d); end rho_r(i,1)=sum; end %Curva analitica d_ana = (lambda/(2*k*(1+n))*(R^2-r.^2)).^n; %Representación de ambas curvas figure(3) hold on plot(r(2:end,1),rho_r(2:end),'m') hold on plot(r(2:end,1),d_ana(2:end,1),'--r') xlabel('r') ylabel('\rho(r)') title('Perfil de densidad en el equilibrio') legend('h=0.04/sqrt(N/1000)','Sol analítica','h=0.03/sqrt(N/1000)') %% Comparar con valor analitico en rho(0) rho_maxana = lambda/(4*k)*R^2; %Densidad en el centro de la extrella, valido para n politropico igual a 1 disp([rho_maxana max(rho_r)])