Shocks en SPH

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.

Dejar un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *


¡IMPORTANTE! Responde a la pregunta: ¿Cuál es el valor de 3 10 ?