SPH

You are currently browsing articles tagged SPH.

Definimos a los kernels como funciones del tipo:

$latex W_{ab}=W(boldsymbol{r}_a – boldsymbol{r}_b,h)$,

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

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

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

Para empezar, $latex 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 $latex q = frac{r}{h}$, siendo $latex r$ la distancia ente las partículas, por lo que:

$latex 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 $latex a$, la función que nos da la distancia de esta a cualquier punto dentro del soporte compacto es:

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

siendo $latex q_a$ lo mismo añadiendo el factor $latex h$.

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

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

$latex = 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:

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

De la misma manera, si tenemos:

$latex 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:

$latex 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:

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

 

Tags: , , , ,

Construimos una malla de $latex 100times 100$ puntos de nuestros kernels en $latex 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:

[youtube http://www.youtube.com/watch?v=toa1Asoo600?rel=0&w=420&h=315]

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

[youtube http://www.youtube.com/watch?v=VMnyetAVyN0?rel=0&w=420&h=315]

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

[youtube http://www.youtube.com/watch?v=YZiuoCBrfOQ?rel=0&w=420&h=315]

Tags: , , , , ,

En $latex 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):

Tags: ,

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:

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

Tags: ,

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 $latex h$.

En la siguiente animación se utiliza una LinkedList:

Tags: , , ,

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

Supongamos que tenemos $latex N$ partículas. Si $latex 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 $latex h$ (all-pair search). La complejidad del algorítmo es $latex 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 $latex 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 $latex h$ es constante para todas las partículas. La complejidad, con un número suficientemente pequeño de partículas por celda, es $latex O(N)$.

Finalmente, podemos trabajar con árboles. Los árboles permiten complejidades genéricas en búsquedas de $latex O(N log N)$ y permiten trabajar con $latex 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).

Tags: , , , , , , ,

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

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

La ley de Newton $latex 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:

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

que es:

$latex m_ax_a” = -k_1 (1-frac{l}{sqrt{x_a^2+y_a^2}})x_a – k_2frac{x’_a}{sqrt{x_a’^2+y_a’^2}}$

$latex m_ay_a” = -k_1 (1-frac{l}{sqrt{x_a^2+y_a^2}})y_a – k_2frac{y’_a}{sqrt{x_a’^2+y_a’^2}}$

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

$latex 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_2frac{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_2frac{v_a}{sqrt{u_a^2+v_a^2}} – m_ag \
end{array}
right .
$

En la siguiente simulación tomamos $latex m_a$ aleatoria para cada partícula $latex a$, por lo que tendremos que resolver el sistema de EDOs numéricamente para cada una de ellas. Además, tendremos $latex l=1$, $latex k_1 = 10000$ y $latex 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, $latex u_{a_0}=0$, $latex v_{a_0} = -0.7$ y $latex v_{z_0} = 0$. El intervalo de tiempo es desde $latex t_i = 0$ hasta $latex t_f=10$ que se divide en $latex 400$ subintervalos.

El resultado es el siguiente:

[youtube http://www.youtube.com/watch?v=emirYBIds6o?rel=0&w=420&h=315]

Tags: , , ,

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:

$latex rho_i = sum_j m_j W_{ij}$

con $latex W_{ij}:=W(boldsymbol{x}_i – boldsymbol{x}_j)$, o también:

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

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

Ellos proponen utilizar,

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

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

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

Tags: , ,

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:

$latex langle f rangle(boldsymbol{x}) := int f(boldsymbol{x’})W(boldsymbol{x’}-boldsymbol{x},h)dboldsymbol{x’}$

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

$latex langle frangle(boldsymbol{x}_i) – f(boldsymbol{x}_i) = frac{h^2_{ef}}{4} nabla^2 f(boldsymbol{x}) + O(h^4_{ef})$.

donde

$latex h^2_{ef}:=2 int x^2 W(boldsymbol{x},h)dboldsymbol{x}$,

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

También se tiene:

$latex langle frac{partial f}{partial x} rangle (boldsymbol{x}) = int frac{partial f(boldsymbol{x’})}{partial x’} W(boldsymbol{x’}-boldsymbol{x},h)dboldsymbol{x’}$

$latex = – int f(boldsymbol{x’}) frac{partial}{partial x’} W(boldsymbol{x’}-boldsymbol{x},h)dboldsymbol{x’}$

$latex = frac{partial}{partial x} int f(x’)W(x’-x,h)dx’$

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

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

entonces:

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

$latex 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:

$latex f_i:= langle f rangle (boldsymbol{x}_i) = int f(boldsymbol{x’})W(boldsymbol{x’}-boldsymbol{x}_i,h)dboldsymbol{x’}$

$latex = 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)dboldsymbol{x}’$

$latex sum_j Delta f_{i,j}$

donde:

$latex 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)dboldsymbol{x}’$.

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

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

y por otra, con:

$latex 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:

$latex ddot{boldsymbol{x}}:=int frac{dboldsymbol{v}(boldsymbol{x})}{dt}W(boldsymbol{x}-boldsymbol{x}_i,h)dboldsymbol{x}$

$latex = – int frac{1}{rho(boldsymbol{x})} nabla P(x) W(boldsymbol{x}-boldsymbol{x}_i,h)dboldsymbol{x}$

y para la energía

$latex dot{u}_i := int frac{du(boldsymbol{x})}{dt} W(boldsymbol{x}-boldsymbol{x}_i,h)dboldsymbol{x}$

$latex = 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) ] dboldsymbol{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.

Tags: , , , ,

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:

$latex 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)

$latex 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

$latex 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 )$

$latex frac{d}{dt}v_{a,h,dis} = $

$latex 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}$

$latex -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

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

y

$latex 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 )$

$latex 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) }$

$latex + 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

$latex 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 ]$

$latex 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):

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

donde

$latex 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) $

$latex frac{d}{dt}u_{a,AV} = $

$latex frac{d}{dt}u_{a,C} = $

Tags: ,

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:

$latex rho_a = sum_b m_b W_{ab}$

La ecuación del momento queda:

$latex 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:

$latex 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.

Tags: , , , , , ,

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:

$latex 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 $latex 0$ a $latex 3$ y los coeficientes de la métrica se definen

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

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

Las ecuaciones del momento se siguen de:

$latex frac{partial}{partial x^nu} T^{i nu} = 0$

que es:

$latex frac{d}{dt} q = – frac{1}{N} nabla P$

y en forma SPH queda:

$latex 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 $latex nu_b$ es el número de bariones asociados a la partícula $latex b$.

Y la de la energía se sigue de:

$latex frac{partial}{partial x^j} T^{0j} = 0$

que es:

$latex frac{d}{dt} epsilon = – frac{1}{N} nabla cdot (Pv)$

y en foma SPH queda:

$latex 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} $

Tags: , , , , , ,

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 $latex A(r)$ de la siguiente manera:

$latex A_I(r) = int A(r’) W(r-r’,h) dr’$

donde $latex r$ es el vector posición, $latex W$ es una función de pesos o kernel (con las propiedades ya comentadas en otros posts) y $latex 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:

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

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

Tags: , ,

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, $latex frac{d}{dt}$, se relaciona con la derivada Euleriana respecto al tiempo, $latex frac{partial}{partial t}$ de la siguiente manera:

$latex 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: $latex frac{d}{dt} rho = – rho nabla cdot vec{v}$
  2. Ecuacion del momento: $latex frac{d}{dt} vec{v} = -frac{nabla P}{rho} + vec{f}$
  3. Ecuación de la energía: $latex 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: $latex P = (gamma -1) cdot rho cdot epsilon$ (ecuación del gas ideal)

Tags: , , , , , , , , , , ,

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.

Tags: , , , , , , , ,

FireStats icon Powered by FireStats