You are currently browsing the tag archive for the ‘SPH’ tag.

Definimos a los kernels como funciones del tipo:

W_{ab}=W(\boldsymbol{r}_a - \boldsymbol{r}_b,h),

donde a es la partícula en la que está centrada la función y b es una partícula dentro del soporte compacto de la función kernel, controlado éste último por h, la smoothing length (longitud de suavizado).

En este post básicamente pretendo aclarar lo que significa \nabla_a W_{ab} cuando, por ejemplo, tenemos definido el kernel como:

W(q) = \alpha_D \exp (-q^2) con 0 \leq q \leq 2.

Para empezar, \alpha_D es una constante de dimensionalidad, por lo que la fórmula está escrita de manera compacta y sirve para cualquier dimensión. Además, tenemos que q = \frac{r}{h}, siendo r la distancia ente las partículas, por lo que:

r = |\boldsymbol{r}_a - \boldsymbol{r}_b| =^{(3D)} \sqrt{(x_a-x_b)^2 + (y_a-y_b)^2 + (z_a-z_b)^2}.

Si fijamos la posición de la partícula a, la función que nos da la distancia de esta a cualquier punto dentro del soporte compacto es:

r_a (\boldsymbol{r}) = |\boldsymbol{r}_a-\boldsymbol{r}| =^{(3D)} \sqrt{(x_a-x)^2 + (y_a-y)^2 + (z_a-z)^2},

siendo q_a lo mismo añadiendo el factor h.

Por lo tanto, en este caso tenemos, en tres dimensiones y donde b es una partícula en una posición arbitraria (x,y,z):

\nabla_a W_{ab}(q) =^{(3D)} (\partial_x W_{ab}(q_a), \partial_y W_{ab}(q_a), \partial_z W_{ab}(q_a)) =

= \alpha_D \exp(-q^2) (-2q) (\partial_x (q_a), \partial_y (q_a), \partial_z (q_a)) = \alpha_D \exp(-q^2) (-2q) \nabla_a q_a

donde:

\nabla_a q_a = \frac{-1}{h r_a} (x_a-x,y_a-y,z_a-z).

De la misma manera, si tenemos:

W(q) = \alpha_D \begin{cases} 1-\frac{3}{2} q^2 + \frac{3}{4} q^3, 0 \leq q < 1\\ \frac{1}{4} (2-q)^3, 1 \leq q < 2 \\ 0, q \geq 2 \end{cases}

entonces:

\nabla_a W_{ab}(q) = \alpha_D \begin{cases} (-3q + \frac{9}{4}q^2) \nabla_a q_a, 0 \leq q < 1 \\ -\frac{3}{4} (2-q)^2 \nabla_a q_a, 1 \leq q < 2 \\ 0, q \geq 2 \end{cases}

Así pues:

\nabla W(q) = \frac{d}{dq} W(q) \nabla q.

 

Anuncios

Construimos una malla de 100\times 100 puntos de nuestros kernels en 2D. La primera animación contiene el kernel Gaussiano en el primer frame, el cúbico en el segundo, el cuártico en el tercero y el quíntico en el último:

En la siguiente animación mostramos los puntos de la primera derivada de los kernels:

Finalmente, la animación con puntos de las segundas derivadas de los kernels de nuestra aplicación:

En 2D nos encontramos las siguientes gráficas.

En primer lugar, dos imagenes, cualitativas y parciales, del Kernel Gaussiano y su derivada (en blanco opaco y hueco respectivamente) junto al Kernel cúbico y su derivada (en azul), y a continuación dos del Kernel Gaussiano y su derivada junto al Kernel quíntico y su derivada (en rojo):

A continuación tenemos diferentes funciones kernel, el linea contínua, junto a su primera y segunda derivada, en lineas discontinuas de mayor y menor longitud respectivamente.

En color negro el kernel gaussiano, su primera derivada y su segunda derivada:

en azul el cúbico:

en verde el cuártico:

y en rojo el quíntico:

y representa f(q), x representa q con r/h y estamos en 1D.

Actualizamos la minianimación que teniamos incluyendo los algoritmos y las estructuras de datos necesarias para la NNPS.

Con el fin de poder comprobar la determinación de vecinos, el color de cada péndulo depende, en un instante dado, del número de vecinos que tiene a una distancia menor o igual a h.

En la siguiente animación se utiliza una LinkedList:

En el método SPH necesitamos conocer cual es su conjunto de vecinos de una partícula a. Este viene determinado por la smoothing length h, uno de los parámetros de las funciones kernel.

Supongamos que tenemos N partículas. Si N es lo suficientemente pequeño, lo único que tenemos que hacer es, para cada partícula, recorrer el resto de partículas calculando su distancia a nuestra partícula y tener en consideración solo aquellas que esten a distancia menor que h (all-pair search). La complejidad del algorítmo es O(N^2).

Sin embargo, cuando el número de partículas crece, y que lo ha de hacer pues en eso se basa la potencia del método, esta manera de trabajar es mejorable.

A continuación veremos algunas opciones, pero la idea general es que necesitamos tener una descomposición espacial de nuestro escenario que nos permita encontrar de manera eficiente el vecindario de cada una de nuestras partículas y que tendremos que mantener actualizada a medida que evoluciona el sistema.

Otra opción es la linked-list. Se trata de superponer una malla con un tamaño de celda 2h, de manera que, dada una partícula, solo tendremos que buscar las partículas que se encuentran en las celdas adyacentes. En C++ tenemos las clases list y vector. La primera por nombre y la segunda porque la clase vector puede crecer dinámicamente. Solo podremos utilizarlas si la h es constante para todas las partículas. La complejidad, con un número suficientemente pequeño de partículas por celda, es O(N).

Finalmente, podemos trabajar con árboles. Los árboles permiten complejidades genéricas en búsquedas de O(N \log N) y permiten trabajar con h diferentes para cada partícula. En astrofísica, por ejemplo, podríamos tener un arbol binario de búsqueda, pues tenemos que simular objetos autogravitantes y necesitamos resolver un problema de n-cuerpos, que parece ser que es lo que mas tiempo de CPU consume. Para hacer esto de manera eficiente necesitamos un BST. Se trata de utilizar este mismo para realizar las búsquedas. La clase map de C++ se implementa con un arbol binario balanceado (AVL).

Existen otras estructuras de datos que permiten hacer descomposiones espaciales que podria ser interesante analizar (UB-tree, Octrees, BST, kd-tree, R-tree).

En la siguiente simulación tenemos 1000 partículas. Cada partícula a representa un péndulo de masa m_a y posición (x_a(t),y_a(t),z_a(t)), con z_a(t) constante, sometido a las fuerzas:

  1. F_1 = -m_ag, la fuerza de la gravedad, con g constante.
  2. F_2, la fuerza elática de la suspensión, de sentido opuesta a la posición del péndulo y de magnitud k_1d, con k_1 > 0 y d la variación de la longitud del péndulo respecto de su longitud en reposo l.
  3. F_3, la fuerza de fricción, con sentido opuesta a la velocidad del péndulo y magnitud k_2v, con v = \sqrt{(x_a')^2+(y_a')^2} y k_2 > 0.

La ley de Newton F = ma queda, para cada uno de los péndulos, es decir, para cada partícula, como el sistema de dos ecuaciones diferenciales de segundo orden:

m_a \begin{bmatrix} x_a''(t) \\ y_a''(t) \end{bmatrix} = F_1 + F_2 + F_3

que es:

m_ax_a'' = -k_1 (1-\frac{l}{\sqrt{x_a^2+y_a^2}})x_a - k_2\frac{x'_a}{\sqrt{x_a'^2+y_a'^2}}

m_ay_a'' = -k_1 (1-\frac{l}{\sqrt{x_a^2+y_a^2}})y_a - k_2\frac{y'_a}{\sqrt{x_a'^2+y_a'^2}}

Podemos expresar este sistema como un sistemad de ecuaciones de primer orden introduciendo incongnitas adicionales: u_a = x_a' y v_a = y_a', con lo que nos queda:

\left \{  \begin{array}{rcl}  x_a' &=& u_a \\  y_a' &=& v_a \\  m_au_a' & = & k_1 (\frac{l}{\sqrt{x_a^2+y_a^2}}-1)x_a - k_2\frac{u_a}{\sqrt{u_a^2+v_a^2}} \\  m_av_a' & = & k_1 (\frac{l}{\sqrt{x_a^2+y_a^2}}-1)y_a - k_2\frac{v_a}{\sqrt{u_a^2+v_a^2}} - m_ag \\  \end{array}  \right .

En la siguiente simulación tomamos m_a aleatoria para cada partícula a, por lo que tendremos que resolver el sistema de EDOs numéricamente para cada una de ellas. Además, tendremos l=1, k_1 = 10000 y k_2=1 en todos los péndulos. Cada uno de éstos se distribuye a lo largo de una espilar cilíndrica y las velocidades iniciales son, para todas ellas, u_{a_0}=0, v_{a_0} = -0.7 y v_{z_0} = 0. El intervalo de tiempo es desde t_i = 0 hasta t_f=10 que se divide en 400 subintervalos.

El resultado es el siguiente:

En el artículo A modified SPH approach for fluids with large density differences de F. Ott y E. Schnetter se explica una nueva modificación del método SPH que permite la interacción entre fluidos con densidades muy diferentes.

En la aproximación del SPH estandar tenemos:

\rho_i = \sum_j m_j W_{ij}

con W_{ij}:=W(\boldsymbol{x}_i - \boldsymbol{x}_j), o también:

\frac{d}{dt}\rho_i = \sum_j m_j \boldsymbol{v}_{ij} \cdot \nabla W_{ij}

donde \boldsymbol{v}_{ij}:=(\boldsymbol{v}_i - \boldsymbol{v}_j) y \nabla W_{ij}:=(\nabla W)(\boldsymbol{x}_i - \boldsymbol{x}_j), que tiene la ventaja de que los valores iniciales para \rho_i pueden escogerse libremente permitiendo la modelización de discontinuidades.

Ellos proponen utilizar,

\frac{d}{dt}\rho_i = m_i \sum_j \boldsymbol{v}_{ij} \cdot \nabla W_{ij}

deducida suavizando la densidad de partículas n_i = \frac{1}{V_i} en lugar de la densidad de masa \rho_i, donde V_i = \frac{m_i}{\rho_i} y \frac{d}{dt}m_i = 0.

justificandolo con tests positivos en tres problemas: ecuación de advección, onda sonora y tubo de choque.

En el review que Rosswog hace sobre el método SPH, en especial en sus aplicaciones a la astrofísica, hay un apartado dedicado al tratamiento de los shocks. El tratamiento de los strong shocks es uno de los puntos débiles del método SPH (penetración de partículas parcialmente resuelto mediante XSPH). En este apartado comenta que, básicamente, exiten dos maneras de tratarlos:

  1. Utilizar soluciones (analíticas) de un problema de Riemann entre dos partículas adyacentes.
  2. Ampliar la discontinuidad hasta una longitud resoluble numéricamente, para lo que se añaden terminos disipativos artificiales extra al flujo.

Comenta que la mayoria de las implementaciones de SPH utilizan esta segunda aproximación, consiguiendo esta amplicación o mediante el método de discretización utilizado, utilizando métodos que modifican las ecuaciones originales de una manera que nos puede venir bien, tipo el esquema Lax; o mediante la adición de terminos ad-hoc, algo bastante habitual, de tipo presión.

Pero también existen métodos que implementan la primera, y cita dos artículos:

  1. El artículo Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver de Inutsuka y Shu-Ichiro.
  2. El artículo Implementations and tests of Godunov-type particle hydrodynamics de Cha y Whitworth.

En el primero comenta como consiguen, tras algunos intentos fallidos,  introducir Riemann Solver en las ecuaciones de manera similar a como se consigue en los métodos basados en malla, manteniendo la formulación conservativa de las mismas. Para empezar, comentan que mientras el error en la mayoria de los métodos basados en malla proviene de la truncación de la serie de de Taylor, en SPH, si consideramos:

\langle f \rangle(\boldsymbol{x}) := \int f(\boldsymbol{x'})W(\boldsymbol{x'}-\boldsymbol{x},h)d\boldsymbol{x'}

como la convolución de la función f(x) con el kernel. Mediante Taylor, de f alrededor de x=x_i, podemos ver que:

\langle f\rangle(\boldsymbol{x}_i) - f(\boldsymbol{x}_i) = \frac{h^2_{ef}}{4} \nabla^2 f(\boldsymbol{x}) + O(h^4_{ef}).

donde

h^2_{ef}:=2 \int x^2 W(\boldsymbol{x},h)d\boldsymbol{x},

por lo que el error proviene de la convolución con la función simétrica W.

También se tiene:

\langle \frac{\partial f}{\partial x} \rangle (\boldsymbol{x}) = \int \frac{\partial f(\boldsymbol{x'})}{\partial x'} W(\boldsymbol{x'}-\boldsymbol{x},h)d\boldsymbol{x'}

= - \int f(\boldsymbol{x'}) \frac{\partial}{\partial x'} W(\boldsymbol{x'}-\boldsymbol{x},h)d\boldsymbol{x'}

= \frac{\partial}{\partial x} \int f(x')W(x'-x,h)dx'

De manera que \nabla f = \nabla \langle f \rangle. Si definimos:

\rho(\boldsymbol{x}):=\sum_j m_j W(\boldsymbol{x}-\boldsymbol{x}_j,h),

entonces:

\sum_j \frac{m_j}{\rho(\boldsymbol{x})}W(\boldsymbol{x}-\boldsymbol{x}_j,h) = 1,

\sum_j m_j \nabla \big [ \frac{W(\boldsymbol{x}-\boldsymbol{x}_j,h)}{\rho(\boldsymbol{x})} \big ]= 0

que son la clave del método, pues ahora se puede hacer:

f_i:= \langle f \rangle (\boldsymbol{x}_i) = \int f(\boldsymbol{x'})W(\boldsymbol{x'}-\boldsymbol{x}_i,h)d\boldsymbol{x'}

= \int \sum_j m_j \frac{f(\boldsymbol{x}')}{\rho(\boldsymbol{x}')}W(\boldsymbol{x}'-\boldsymbol{x}_i,h)W(\boldsymbol{x}'-\boldsymbol{x}_j,h)d\boldsymbol{x}'

\sum_j \Delta f_{i,j}

donde:

\Delta f_{i,j}:=\int m_j \frac{f(\boldsymbol{x}')}{\rho(\boldsymbol{x}')}W(\boldsymbol{x}'-\boldsymbol{x}_i,h)W(\boldsymbol{x}'-\boldsymbol{x}_j,h)d\boldsymbol{x}'.

En el método SPH estandar nos encontramos, por una parte, con:

\sum_j \frac{m_j}{\rho(\boldsymbol{x}_j)}W(\boldsymbol{x}-\boldsymbol{x}_j,h) \approx 1

y por otra, con:

\Delta f_{i,j} \approx m_j \frac{f(\boldsymbol{x}_j)}{\rho(\boldsymbol{x}_j)}W(\boldsymbol{x}_i-\boldsymbol{x}_j,h)

y que, evidentemente, son aproximaciones mas pobres a las anteriormente obtenidas.

A continuación derivan las ecuaciones para la posición:

\ddot{\boldsymbol{x}}:=\int \frac{d\boldsymbol{v}(\boldsymbol{x})}{dt}W(\boldsymbol{x}-\boldsymbol{x}_i,h)d\boldsymbol{x}

= - \int \frac{1}{\rho(\boldsymbol{x})} \nabla P(x) W(\boldsymbol{x}-\boldsymbol{x}_i,h)d\boldsymbol{x}

y para la energía

\dot{u}_i := \int \frac{du(\boldsymbol{x})}{dt} W(\boldsymbol{x}-\boldsymbol{x}_i,h)d\boldsymbol{x}

= \sum_j m_j \int \frac{P}{\rho^2}[\boldsymbol{v}-\dot{\boldsymbol{x}}_i] \cdot [\nabla W(\boldsymbol{x}-\boldsymbol{x}_i,h)W(\boldsymbol{x}-\boldsymbol{x}_j,h) - W(\boldsymbol{x}-\boldsymbol{x}_i,h) \nabla W(\boldsymbol{x}-\boldsymbol{x}_j,h) ] d\boldsymbol{x},

que se modifica a

$latex$

para tener en cuenta la disipación

En el segundo artículo define e implementa el método Godunov-type Particle Hydrodynamics (GPH). Básicamente GPH = SPH + Godunov upwind schema.

En el artículo “MAGMA: a three-dimensional, Lagrangian magnetohydrodynamic code for merger applications” de S. Rosswog y D. Price comentan como introducir campos magnéticos en SPH. Las ecuaciones ya discretizadas quedan:

Ecuación de densidad:

\rho = \sum_b m_b W(r-r_b,h)

Ecuación del momento (conh: “grad-h” term, mag: magnetic force term, g: self-gravity and gravitational softening term)

\frac{d}{dt}v_{a,MHD} = \frac{d}{dt} (v_{a,h}+v_{a,h,dis}+v_{a,g}+v_{a,mag}+v_{a,mag,dis})

donde

\frac{d}{dt}v_{a,h} = -\sum_b m_b \big ( \frac{P_a}{\Omega_a \rho_a^2} \nabla_a W_{ab}(h_a)+ \frac{P_b}{\Omega_b \rho_b^2} \nabla_a W_{ab}(h_b) \big )

\frac{d}{dt}v_{a,h,dis} =

\frac{d}{dt}v_{a,g} = -G \sum_b m_b \big [ \frac{\phi'_{ab}(h_a) + \frac{\phi'_{ab}(h_b)}{}}{2} \big ]\hat{e}_{ab}

-\frac{G}{2} \sum_b m_b \big [ \frac{\zeta_a}{\Omega_a} \nabla_a W_{ab}(h_a) + \frac{\zeta_b}{\Omega_b} \nabla_a W_{ab}(h_b) \big ]

con

\phi'_{ab} = \frac{\partial \phi}{\partial|r_a-r_b|}, \zeta_k := \frac{\partial h_k}{\partial \rho_k} \sum_b m_b \frac{\partial \phi_{kb}(h_k)}{\partial h_k}

y

\Omega_a := \big ( 1- \frac{\partial h_a}{\partial \rho_a} \cdot \sum_b m_b \frac{\partial}{\partial h_a} W_{ab}(h_a) \big )

\frac{d}{dt}v_{a,mag} = - \sum_b \frac{m_b}{\mu_0} \big \{ \frac{B_a^2 / 2}{\Omega_a \rho_a^2} \nabla_a W_{ab}(h_a) + \frac{B_b^2 / 2}{\Omega_b \rho_b^2} \nabla_a W_{ab}(h_b) \}

+ \sum_b \frac{m_b}{\mu_0} \big \{ \frac{B_a(B_a \cdot \overline{\nabla_a W_{ab}}) - B_b(B_b \cdot \overline{\nabla_a W_{ab}})}{\rho_a \rho_b} \big \}

con

\overline{\nabla_a W_{ab}} = \frac{1}{2} \big [ \frac{1}{\Omega_a} \nabla_a W_{ab}(h_a) + \frac{1}{\Omega_b} \nabla_a W_{ab}(hb) \big ]

\frac{d}{dt}v_{a,mag,dis} =

Ecuación de la energía (con h: “grad-h” term, AV: Artificial Viscosity term y C: Condutivity):

\frac{d}{dt}u_{a,MDH} = \frac{d}{dt} (du_{a,h} + du_{a,AV} + du_{a,C})

donde

\frac{d}{dt}u_{a,h} = \frac{1}{\Omega_a} \frac{P_a}{\rho_a^2} \sum_b m_b v_{ab} \cdot \nabla_a W_{ab}(h_a)

\frac{d}{dt}u_{a,AV} =

\frac{d}{dt}u_{a,C} =

En [Rosswog 2009] tenemos las ecuaciones de la hidrodinámica en forma Lagrangiana discretizadas y en su forma mas básica. Las partículas avanzaran en el tiempo siguiendo las siguientes ecuaciones:

Para empezar, no hay necesidad de resolver la ecuación de continuidad ya que la masa de las partículas permanece fija. Podemos obtener las densidades mediante:

\rho_a = \sum_b m_b W_{ab}

La ecuación del momento queda:

\frac{d}{dt}\vec{v}_a = - \sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2} + \Pi_{ab} ) \nabla_a W_{ab}

La ecuación de evolución para la energía interna específica puede escribirse como:

\frac{d}{dt} u_a = \sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{1}{2}\Pi_{ab}) \vec{v}_{ab} \nabla_b W_{ab}

Rosswog las llama “vanilla ice” SPH.

En [Monaghan 1992] comenta el caso del método SPH en relatividad especial.

Para empezar asumimos que el fluido está constituido por bariones, por lo que el tensor de energia-momento es:

T^{\mu \nu} = (n m_0 c^2 + n \tau + P) U^\mu U^\nu + P g^{\mu \nu}

donde los indices griegos van de 0 a 3 y los coeficientes de la métrica se definen

g_{00} = -1 y g_{ij} = 1

En estas ecuaciones, n representa la densidad de bariones, P es la presión, \tau la energía térmica, c la velocidad del sonido, U^\nu la 4-velocidad con U_\nu U^\nu = -1 y m_0 la masa.

Las ecuaciones del momento se siguen de:

\frac{\partial}{\partial x^\nu} T^{i \nu} = 0

que es:

\frac{d}{dt} q = - \frac{1}{N} \nabla P

y en forma SPH queda:

\frac{d}{dt}q_a = -\sum_b \nu_b (\frac{P_a}{n_a^2} + \frac{P_b}{N_b^2}) \nabla_a W_{ab}

donde \nu_b es el número de bariones asociados a la partícula b.

Y la de la energía se sigue de:

\frac{\partial}{\partial x^j} T^{0j} = 0

que es:

\frac{d}{dt} \epsilon = - \frac{1}{N} \nabla \cdot (Pv)

y en foma SPH queda:

\frac{d}{dt} \epsilon_a = -\sum_b m_b (\frac{P_a v_a}{N_a^2} + \frac{P_b v_b}{N_b^2}) \nabla W_{ab}

Según Monaghan en [Monaghan, 1992], uno de los padres del método, buscaban un método con el que fuera fácil trabajar y diera buenos resultados y encontraron eso y mucho mas.

Además de no  necesitar malla para calcular las derivadas, puesto que se calculan de manera analítica a partir de una fórmula de interpolación, las ecuaciones de conservación del momento y la energía pasan a ser un conjunto de ecuaciones diferenciales ordinarias fáciles de entender.

En esencia, el método SPH es un método de interpolación que permite expresar una función a partir de los valores en un conjunto desordenado de puntos a los que llamamos particulas.

La idea es que aproximamos cualquier función A(r) de la siguiente manera:

A_I(r) = \int A(r') W(r-r',h) dr'

donde r es el vector posición, W es una función de pesos o kernel (con las propiedades ya comentadas en otros posts) y h es la longitud de suavizado. Este tipo de aproximaciones reciben el nombre de interpolación integral.

Para trabajar de manera numérica, la interpolador integral se aproxima por:

A_S(r) = \sum_b m_b \frac{A_b}{\rho_b} W(r-r_b,h)

donde el índice del sumatorio b identifica a cada partícula y el sumatorio es sobre todas las partículas. La partícula b tiene masa m_b, posición r_b, densidad \rho_b y velocidad v_b. El valor de la función A en r_b es A_b.

En el artículo [Rosswog 2009], Stephan Rosswog hace un repaso detallado del método SPH centrandose especialmente en sus aplicaciones en astrofísica. Repasamos el apartado que hace referencia a las ecuaciones de la hidrodinámica en forma Lagrangiana.

A diferencia de los metodos basados en malla, que son Eulerianos, es decir, métodos donde  describimos el fluido desde un punto fijo del espacio, el Smoothed Particle Hydrodynamics es totalmente Lagrangiano, por lo que describimos el fluido desde un sistema de coordenadas fijado en una particula del fluido en movimiento.

La derivada Lagrangiana o sustancial respecto del tiempo, \frac{d}{dt}, se relaciona con la derivada Euleriana respecto al tiempo, \frac{\partial}{\partial t} de la siguiente manera:

\frac{d}{dt} = \frac{dx^i}{dt} \frac{\partial}{\partial x^i} + \frac{\partial}{\partial t} = \vec{v} \cdot \nabla + \frac{\partial}{\partial t}

Aplicando esta relación a las ecuaciones en forma Euleriana, las ecuaciones de la hidrodinámica en forma Lagrangiana quedan:

  1. Ecuación de continuidad: \frac{d}{dt} \rho = - \rho \nabla \cdot \vec{v}
  2. Ecuacion del momento: \frac{d}{dt} \vec{v} = -\frac{\nabla P}{\rho} + \vec{f}
  3. Ecuación de la energía: \frac{d}{dt}u = \frac{P}{\rho^2} \frac{d}{dt} \rho = - \frac{P}{\rho} \nabla \cdot \vec{v}
  4. Ecuación de estado, que describe la termodinámica del fluido estelar: P = (\gamma -1) \cdot \rho \cdot \epsilon (ecuación del gas ideal)

En el artículo  [Gingold & Monaghan, 1977] se presenta por primera vez el método Smoothed Particle Hydrodynamics. Los autores, originalmente, buscaban un método que permitiera tratar problemas en astrofísica asimétricos (sin simetría esférica, sin simetría axial, etc.) . En estos casos, los métodos de diferencias finitas no se adaptan bien, pues requieren elevar el número de puntos en la malla para seguir con la precisión deseada su evolución, lo cual complica enormemente la evolución de las integrales múltiples.

Lo que pensaron es utilizar la descripción Lagrangiana del flujo del fluido que centra su atención en los elementos del fluido. En la discretización, estos elementos se mueven siguiendo las leyes de Newton con fuerzas debidas a los gradientes de presión y a otras fuerzas de cuerpo como la gravedad, rotación o magnéticas. ¿Qué método utilizar para determinar las fuerzas que actuan en un momento determinado sobre un elemento de fluido?

Para empezar, para elementos de fluido de igual masa, el número de elementos por unidad de volumen debe ser proporcional a la densidad. Además, sin ningún tipo especial de simetría, la posicion de los elementos será aleatoria de acuerdo con la densidad. Para recuperar la densidad de la distribución conocida de los elementos es equivalente a recuperar una distribución de probabilidad a partir de una muestra. Existen dos métodos para conseguir esto que funcionan bien con problemas de fluidos: el smoothing kernel method y  la técnica del spline delta. Ambos métodos se pueden pensar como la aproximación de una integral por el procedimiento de Monte Carlo.

diciembre 2017
L M X J V S D
« Ago    
 123
45678910
11121314151617
18192021222324
25262728293031